Direct Simulations of Particle Acceleration in Fluctuating Electromagnetic Field across a Shock
Abstract
We simulate the acceleration processes of collisionless particles in a shock structure with magnetohydrodynamical (MHD) fluctuations. The electromagnetic field is represented as a sum of MHD shock solution () and torsional Alfven modes spectra (). We represent fluctuation modes in logarithmic wavenumber space. Since the electromagnetic fields are represented analytically, our simulations can easily cover as large as eight orders of magnitude in resonant frequency, and do not suffer from spatial limitations of box size or grid spacing. We deterministically calculate the particle trajectories under the Lorenz force for time interval of up to ten years, with a time step of . This is sufficient to resolve Larmor frequencies without a stochastic treatment. Simulations show that the efficiency of the first order Fermi acceleration can be parametrized by the fluctuation amplitude . Convergence of the numerical results is shown by increasing the number of wave modes in Fourier space while fixing .
Efficiency of the first order Fermi acceleration has a maximum at . The acceleration rate depends on the angle between the shock normal and , and is highest when the angle is zero. Our method will provide a convenient tool for comparing collisionless turbulence theories with, for example, observations of bipolar structure of super nova remnants (SNRs) and shelllike synchrotronradiating structure.
1 Introduction
Cosmic rays have the spectrum of up to the socalled kneeenergy of . Cosmic ray propagation theories suggest energy spectra at the cosmic ray acceleration sites (e.g. Strong et al., 2007).
The current description of cosmic ray acceleration up to knee energy is the well known firstorder Fermi acceleration (Axford et al., 1977; Bell, 1978). In the firstorder Fermi acceleration model, magnetic turbulence is an important agent for particle acceleration. Turbulence makes the particle momenta isotropic, thus portion of the particles cross the shock front many times. Expectation value of the kinetic energy after times of shock crossing is . On the other hand, the probability for a particle to survive shock crossing can be roughly estimated as . This gives us the powerlaw spectrum of .
Ellison et al. (1996), Lucek & Bell (2000), and Bell & Lucek (2001) have done simulations to describe the selfconsistent generation of turbulence, with approximations such as gyrocenter approximation, random walk approximation, or lowering the dimension. On the other hand, the recent development of the particlein cell simulation has made it possible to describe the particle acceleration in electronpositron plasma selfconsistently (e.g. Spitkovsky, 2008).
In this letter, we propose an alternative approach to the simulation of cosmic ray acceleration. We have calculated the motion of particles deterministically, solving the particles’ cyclotron motion from Larmor radii of thermal particles to that of knee energy particles . According to the theories, we assume turbulence spectrum in space. This allows us to cover a large dynamicrange of space and energy, which enables us a direct comparison of the accelerated cosmic ray spectra with the observations.
2 Numerical Scheme
2.1 Representation of Turbulence
Upstream and Downstream Regions
In our method, the electromagnetic field and velocity field of a continuous region are given by
(1)  
(2)  
(3) 
where the amplitude and the wavenumber of the th mode are
(4)  
(5)  
(6)  
(7) 
and the initial phase of the th mode is .
Here is the spectral index that reflects the nature of the turbulence and is the total number of the modes (), is the amplitude for each mode, is its wavenumber, and are two mutually perpendicular unit vectors that are perpendicular to . We choose to be either parallel or antiparallel to .Equation (6) means that are logarithmically discrete.
We use as the measure of the strength of the fluctuation, independent of . Because increasing while keeping (1) keeps the magnetic energy in fluctuation mode, and (2) keeps the expectation value of the fluctuation field the same, if are independent. We will confirm these properties in section 3.
The argument to derive in space is summarized below: Variables in space are marked by tilde. The power law energy spectrum is in Kolmogorov turbulence case. This energy spectrum is in linear bin. In log energy bin the spectral power is ; and since , . Thus, we get for our discretization of the turbulent magnetic field.
Junction Conditions
We assumed strong shock junction condition with low plasma limit at the shock front: , , and , where is the shock compression ratio, and are components of the normal and tangential to the shock, respectively.
2.2 Initial Condition and Equation of Motion
Initial Condition
For each set of initial condition we introduce electromagnetic fields described in section 2.1. We choose a set of initial turbulence phase , sign of and from uniform distribution. Then we put protons in Boltzmann distribution of temperature at the upstream side of the shock.
fluid speed in shock frame  
Alfven speed in fluid frame  
unperturbed magnetic field strength  
torsional Alfven mode energy  
angle of to shock normal  
maximum wavelength of turbulence  
minimum wavelength of turbulence  
temperature of the particles  
particle mass  
particle charge  

Evolution
We make each turbulence mode propagate at Alfven velocity of the uniform field as in Equation 1  3, and updated the particle with Lorenz force, with 4th order RungeKutta method. We choose time discretization for each particle at every timestep, so that and always hold. Typical time step is whereas the Larmor period of thermal particle for is .
(8)  
(9) 
3 Result
We have made following examinations for the results of our method. First, we have traced the particles’ energy as the function of shock crossing number (Figure 3). The inclination of the curves match the inclinationof the theoretical prediction, . Secondly, we have traced the spatial location where the particles gained their kinetic energy. We have found that 94% of the final kinetic energy have been earned within final Larmor radius away from the shock. This is consistent with the firstorder Fermi acceleration picture. Thirdly, we have studied the validity of our Fourier representation in space. We have kept the physical parameters and increased the number of modes per decade ; we see that the spectra converge, and converge to the theoretical powerlaw spectrum (Figure 3). This justifies our use of space discretization.
We have done a large number of simulations while varying the background conditions, from to , , the ratio of magnetic energy in fluctuation mode to that in background field from to , , the angle between the background field and the shock normal from to . Figure 3 shows the time evolution of the energy spectrum for years. The highenergy end of the spectrum gradually grows, and reaches by years.
In our simulation all the particles start its motion in the given time. Since we don’t include the backreaction from the particles to the electromagnetic field in our simulations, timeintegral of energy spectra at each timeslice gives the steady state energy spectra. This steady state spectrum is also shown in Figure 3 with thick curve. The nonthermal component has powerlaw spectrum that meets the observational requirement mentioned in section 1. We can also estimate the “injection rate” to be the proportion of particles that have more than of energy after 1 years. For , , , and , the injection rate was , , , and , respectively.The other parameters are , and .
In Figure 4 we show for all the parameter range the ratio of the particle numbers that were accelerated to have energy greater than . We see that the acceleration is most efficient at polar region () when . We can understand this dependence of the spectra with background fluid parameters as follows; particles are trapped in Larmor motion and tend to move in direction of . Thus particles more easily cross the shockfront when is parallel to shock normal. If the turbulence amplitude is much weaker, fewer particles get reflected by pitch angle scattering, and Fermi acceleration is suppressed. The injection is more efficient for smaller , because more energy is distributed to modes with smallest wavelengths which are resonant with the thermal particles.
If a spherical shock emerges in a uniform mean magnetic field, there are two polar region where the mean magnetic field is parallel to the shock normal, and the equatorial region has the mean magnetic field perpendicular to the shock normal. Thus, we expect the Fermi acceleration process to be only active in the pair of polar region. This might explain the bipolar structure we see at SN 1006.
We have also checked the acceleration rate in threedimensional(isotropic), rather than onedimensional(anisotropic) distribution of . We have found that less significant dependence of injection rate on with larger . We can interpret this as follows; if the turbulence spectrum is isotropic and the maximum turbulence wavelength is large, turbulence modes with largest wavelength and strongest amplitude play the role of local . Thus we observe almost isotropic Fermi acceleration. This might explain the many SNRs with no typical orientation.
4 Discussion
Some might question the validity of value much greater than unity. However, Uchiyama et al. (2007) observed extremely fast varying Xray images at SNR RX J1713.73946. Their observation may indicate that magnetic field is locally enhanced up to in year in SNR, which corresponds to case in our model. Our simulations suggest that such fastvariating spots in SNRs might be the sites of galactic cosmic ray acceleration.
Although we have ignored many of the Fourier modes by adopting space discretization of the turbulence spectrum, the validity of the approximation can be argued in many ways. Most importantly we have confirmed that our measure in space lead to convergence. Figure 3 shows energy spectra for , with increasing . Turbulent cascade, from which the very turbulence arises, is by nature a logarithmic process: a mode of a certain wavelength couples with the mode of half the wavelength by nonlinear term of Euler equation. First order Fermi acceleration also is a logarithmic process: particles gain energy as an exponential function of shock crossing number . All these reason combined, waves in logarithmically discretized wavenumber space act as a sufficient ladder to carry up cosmic ray particles.
References
 Axford et al. (1977) Axford, W., Leer, E., & Skadron, G. 1977, in International Cosmic Ray Conference, Vol. 11, International Cosmic Ray Conference, 132
 Bamba et al. (2003) Bamba, A., Yamazaki, R., Ueno, M., & Koyama, K. 2003, ApJ, 589, 827
 Bell (1978) Bell, A. 1978, MNRAS, 182, 147
 Bell & Lucek (2001) Bell, A., & Lucek, S. 2001, MNRAS, 321, 433
 Ellison et al. (1996) Ellison, D., Baring, M., & Jones, F. 1996, ApJ, 473, 1029
 Lucek & Bell (2000) Lucek, S., & Bell, A. 2000, MNRAS, 314, 65
 Spitkovsky (2008) Spitkovsky, A. 2008, ApJ, 682, L5
 Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle Science, 57, 285
 Uchiyama et al. (2007) Uchiyama, Y., Aharonian, F., Tanaka, T., Takahashi, T., & Maeda, Y. 2007, Nature, 449, 576