Issue 
A&A
Volume 555, July 2013



Article Number  A101  
Number of page(s)  10  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201220966  
Published online  09 July 2013 
Modification of cosmicray energy spectra by stochastic acceleration
^{1}
Zentrum für Astronomie und Astrophysik, Technische Universität
Berlin,
Hardenbergstraße 36,
10623
Berlin,
Germany
email: rct@gmx.eu
^{2}
Institut für Geowissenschaften, Naturwissenschaftliche Fakultät
III, MartinLutherUniversität Halle, 06099
Halle,
Germany
email: lercheian@yahoo.com
Received:
20
December
2012
Accepted:
22
May
2013
Context. Typical space plasmas contain spatially and temporally variable turbulent electromagnetic fields. Understanding the transport of energetic particles and the acceleration mechanisms for charged particles is an important goal of today’s astroparticle physics.
Aims. To understand the acceleration mechanisms at the particle source, subsequent effects have to be known. Therefore, the modification of a particle energy distribution, due to stochastic acceleration, needs to be investigated.
Methods. The diffusion in momentum space was investigated by using both a MonteCarlo simulation code and by analytically solving the momentumdiffusion equation. For simplicity, the turbulence was assumed to consist of onedimensional Alfvén waves.
Results. Using both methods, it is shown that, on average, all particles with velocities comparable to the Alfvén speeds are accelerated. This influences the energy distribution by significantly increasing the energy spectral index.
Conclusions. Because of electromagnetic turbulence, a particle energy spectrum measured at Earth can drastically deviate from its initial spectrum. However, for particles with velocities significantly above the Alfvén speed, the effect becomes negligible.
Key words: plasmas / magnetic fields / turbulence / acceleration of particles / solar wind / cosmic rays
© ESO, 2013
1. Introduction
Understanding the acceleration processes of cosmic rays has been a longstanding problem in astrophysics (e.g., Drury 1983; Blandford & Eichler 1987). It is currently accepted that highenergy particles with energies up to 10^{15} eV are accelerated in supernova remnants (e.g., Abdo et al. 2010; Hartquist & Morfill 1983) via diffusive shock acceleration, thus producing a powerlaw energy spectrum. The underlying processes date back to Fermi (1949) and have since been improved (Duffy & Blundell 2005) by incorporating many additional effects (see Longair 2011, for an overview). The resulting energy spectrum agrees both with measurements (e.g., Nagano & Watson 2000; LetessierSelvon & Stanev 2011) and with laser experiments (e.g., Kuramitsu et al. 2011). Particles with low energies up to 10^{10} eV, instead, originate from the Sun via coronal mass ejections.
In either case, however, there are many processes that may change the energy spectrum of the particles on their way to Earth (Malkov et al. 2010). Energydependent loss and cooling processes, for example, modify the energy distribution function (Buts et al. 2001), generally because the interstellar (and the interplanetary) medium contains turbulent electric fields. Based on Mariner 2 measurements, it was found (Belcher & Davis 1971) that, to some extent, the heliospheric plasma is dominated by outwardpropagating Alfvénic turbulence, which has been named the slab component. However, Alfvénic waves have also been assumed for the investigation of cosmicray acceleration at the galactic center (Fatuzzo & Melia 2011, 2012), for example.
Now the presence of electromagnetic turbulence gives rise to reacceleration processes (Heinbach & Simon 1995) so that particles gain additional energy while propagating through an otherwise homogeneous medium. The anomalous cosmicray component is also influenced by stochastic acceleration processes (Moraal et al. 2008), thereby giving rise (Fisk & Gloeckler 2007a; Fisk et al. 2010) to a suprathermal tail in the distribution function with −5 as the energy spectral index. We note, however, that thermodynamic arguments and the assumption of compressional turbulence may be sufficient to derive a similar spectral index (Fisk & Gloeckler 2007b). Likewise, adiabatic cooling is a more efficient process than the interaction with turbulent electric fields (Zhang & Lee 2011; Litvinenko & Schlickeiser 2011; Tautz et al. 2012). In anisotropic helical turbulence, largescale electric fields may also give rise to particle acceleration (Fedorov & Stehlik 2008).
For these reasons, much effort has been focused on the general properties of momentum diffusion (e.g., Michałek et al. 1999; Lee & Völk 1975; Achterberg 1981) in slab and isotropic turbulence with Alfvén and oblique propagating fast magnetosonic waves, respectively. Furthermore, electrons near the Sun may also be accelerated thanks to the interaction with Whistler waves (Hamilton & Petrosian 1992; Vocks et al. 2005, 2008). Momentum diffusion (Schlickeiser 2002, 1994) can be considered as the underlying process of stochastic acceleration, where particles mostly gain significant amounts of kinetic energy. If spatial diffusion is to be neglected, stochastic acceleration can be reduced to pure momentum diffusion (Ostrowski & SiemieniecOzi¸ebło 1997). The connection of stochastic acceleration and momentum diffusion has been investigated analytically using magnetohydrodynamic (MHD) waves (e.g., Skilling 1975; Schlickeiser 1989a) by applying quasilinear theory (Shalchi & Schlickeiser 2004).
In the presence of shock waves (Schlickeiser 1989b), the firstorder Fermi acceleration is usually dominant compared to stochastic acceleration in the downstream region (Vainio & Schlickeiser 1999). In contrast, in this article the stochastic acceleration outside a shock wave is investigated by means of testparticle simulations. Such codes neglect the influence of the particles on the turbulence (Giacalone & Jokipii 1999; Michałek et al. 1999; Tautz 2010a), which is generally justified if concentration is focused on a tenuous component such as cosmic rays. Because the turbulence consists of MHD plasma waves that carry an electric field, momentum diffusion occurs so that the average particle energy is changed. In most cases, particles with discrete energies are assumed so that the energydependence of transport parameters such as the mean free path can be investigated. Here, in contrast, an initial distribution function is assumed, the time evolution of which is then traced so that the influence of the stochastic acceleration can be evaluated. Furthermore, the two cases will be distinguished by (i) an ensemble that evolves without supplying new particles; and (ii) continuous particle injection according to the prescribed initial distribution. Case (i) has been applied to solar energetic particle events (Dröge 2005; Schlickeiser et al. 2009). Case (ii) corresponds to the generation of the anomalous cosmicray spectrum (Fisk & Gloeckler 2006). Therefore, we intent to demonstrate that MonteCarlo testparticle simulations are a suitable tool for tracking the time evolution of a particle velocity distribution or an energyspectrum.
The paper is organized as follows: in Sect. 2, the basic properties of the electromagnetic turbulence model are presented. In Sects. 3 and 4, the numerical testparticle code PADIAN is introduced and the simulations results are presented, respectively. In Sect. 5, an analytical solution of the momentumdiffusion equation is derived, and is compared to the numerical results. A brief summary and a discussion of the results with regard to future applications are given in Sect. 6.
2. Astrophysical turbulence properties
In astrophysical plasmas such as the heliosphere, turbulent electromagnetic fields are generated by instabilities or by largescale motion cascaded towards smaller scales. Measurements in the solar wind (e.g., Bruno & Carbone 2005) confirm the Kolmogorov (1991) power law, which gave rise to a spectral turbulent energy distribution of the form (e.g., Bieber et al. 1994; Tautz et al. 2006b) (1)where ℓ_{0} = 0.03 a.u. denotes the turbulence bendover scale.
In addition to the turbulence spectrum, the particle behavior is also influenced by the turbulence geometry, i.e., the angular distribution with respect to a preferred direction, which is given by an ambient magnetic field such as the solar magnetic field (Parker 1958; Fisk 1996). Apart from fully threedimensional anisotropic models (Matthaeus et al. 1990; Weinhorst & Shalchi 2010; Rausch & Tautz 2012), there are three models that have been used over the years: (i) the slab model, which assumes that the turbulent field depends solely on the coordinate along the mean magnetic field B_{0} as δB = δB(z) for ; (ii) the composite model, which superposes the slab field by a twodimensional complement as δB = δB_{slab}(z) + δB_{2D}(x,y); and (iii) the isotropic model, which does not assume a preferred direction for the turbulent field. For reasons of simplicity, the slab model is used here.
For the propagating plasma waves, which enter Eq. (3)via the wave frequency ω, specifically linearly polarized, undamped slab Alfvén waves are chosen. In that case, the dispersion relation reads (2)where the Alfvén speed is with ρ the mass density. Other possible wave types such as slow and fast magnetosonic waves or Whistler waves require threedimensional turbulence, something that is beyond the scope of the present investigation.
3. MonteCarlo simulations
For the testparticle simulations, the testparticle MonteCarlo code PADIAN (Tautz 2010a) is used, which solves the trajectories of a large number of particles as they are being scattered by electromagnetic turbulence. Dimensionless variables are introduced as τ = Ωt for the time with Ω = qB/(mc) the Larmor frequency, and R = γv/(Ωℓ_{0}) as a dimensionless rigidity per magnetic field with ℓ_{0} a characteristic length scale (see below).
3.1. Numerical turbulence generation
Numerically, the turbulent electromagnetic fields are obtained from the superposition of a number of propagating MHD plasma waves as (3)where the wavenumbers k_{n} are distributed logarithmically in the interval and where β is a random phase angle. For the amplitude and the polarization vector, one has and , respectively, with the primed coordinates determined by a rotation matrix with random angles.
With these assumptions and for the variables introduced above, the NewtonLorentz equation can be expressed as (4)where and are unit vectors in the directions of the turbulent and the background magnetic fields, respectively. The unit vector for the electric field is obtained using Faraday’s induction law (Schlickeiser 2002). Furthermore, the parameter R_{A} ≡ v_{A}/(ℓ_{0}Ω) is called the Alfvén rigidity.
3.2. Diffusion coefficients
In the normalized coordinates, spatial diffusion coefficients are defined through the mean square displacements in the directions parallel and perpendicular to the ambient magnetic field as (Tautz 2010a) (5)Here, κ is timedependent and hence is frequently denoted the running diffusion coefficient.
Using the same general derivation, a dimensionless running momentum diffusion coefficient can be defined through (cf. Ostrowski & Michałek 1993; Tautz 2010b) (6)From the inverse momentum diffusion coefficient, the socalled acceleration time scale (e.g., Zank et al. 2006; O’Sullivan et al. 2009; Dosch & Shalchi 2010) can be constructed as (7)Whether or not the momentum diffusion coefficient is truly diffusive, i.e., whether it has a finite constant value for large times, is not entirely clear at present (cf. Tautz 2010b; Tautz & Shalchi 2010). Generally, the larger the ratio R_{A}/R, the sooner D_{p} attains an asymptotic value^{1}.
3.3. Rigidity distribution
In most testparticle simulations, a number of discrete initial particle rigidities are assumed which are then investigated separately, which is acceptable when one attempts to investigate the rigiditydependence of the diffusion coefficients. Here, however, the goal is to trace the evolution of a given rigidity distribution function, while the particles are scattered in momentum space.
In accordance with cosmicray observations (see LetessierSelvon & Stanev 2011, for recent observations), the initial rigidity distribution is modeled in the form of a power law. The requirement that the distribution function be normalized to unity leads to (8a)in and f(R,a) = 0 otherwise. We note that Eq. (8a)is valid for a ≠ 1, whereas for a = 1 one has the special form (8b)again in . For the numerical simulations, the spectral index has to be chosen so that, on the one hand, a relatively steep initial energy spectrum is obtained. On the other hand, the steeper the spectrum, the fewer particles will be found at the highenergy end of the spectrum, thus requiring a larger total number of particles. As a compromise, here the value a = 1.5 will be used. The same holds true for the minimum and maximum rigidity values. Whereas the interpretation of the numerical simulations requires a large number of particles per rigidity interval, R_{min} and R_{min} are chosen so that particles with rigidities comparable to and much larger than R_{A} are present but, at the same time, the rigidity interval is kept moderately narrow.
Parameters used for the MonteCarlo simulation.
Fig. 1 Evolution of the particle distribution as a function of the rigidity R as the normalized time Ωt increases. For the Alfvén rigidity, values of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel) have been used. (This figure is available in color in electronic form.) 
Fig. 2 Energy spectral index for the cases of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel). Shown is the energy distribution function at the beginning (black lines) and at the end of the simulation for the cases of single injection (red lines) and continuous particle injection (blue lines). Also shown are the spectral indices obtained from powerlaw fits (dashed lines) to the simulation data (thin lines). (This figure is available in color in electronic form.) 
Fig. 3 Maximum of the distribution function as the time increases. Again the cases of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel) are shown. The black and red lines denote the cases of single injection and continuous particle injection, respectively. The thick lines simply show cubic fits of the actual simulation data. Additionally, the average particle rigidity is shown in the lower panel (dotted line). (This figure is available in color in electronic form.) 
In the simulation code, particle rigidities are randomly drawn from the distribution function (Eq. (8a)) using the transformation method (e.g., Press et al. 2007, p. 362). The method requires uniform random numbers x ∈ [0,1] to be fed to the inverse cumulant F of the distribution function, which reads (9a)for a ≠ 1 and (9b)if a = 1, thereby returning rigidity values R with the probability prescribed by Eq. (8a).
4. Simulation results
Two basic sets of simulations need to be distinguished: (i) single injection, where all particles are injected at time t = 0 and their trajectories are monitored; and (ii) continuous injection, where, throughout the simulation, new particles are injected with rigidities determined according to the distribution function from Eq. (8a). Technically; case (ii) is realized by choosing the maximum simulation time as a random number so that, throughout the simulation, particles with random life times are present. All simulation parameters are summarized in Table 1.
In Fig. 1, the evolution of the distribution function is traced as time increases. Initially, a powerlaw distribution according to Eq. (8a)is assumed. As the simulation time increases, a maximum in the distribution function is quickly established. Because of the stochastic acceleration of the particles, the maximum is then constantly shifted towards higher energies. In the case of a low Alfvén velocity, the process proceeds relatively slowly. We note that the results shown in Fig. 1 correspond to case (i).
The quantitative variation of the energy spectral index is illustrated in Fig. 2, which illustrates the time evolution of the velocity distribution by comparing the initial and the final particle spectra. Whereas initially all particles are distributed according to the power law in Eq. (8a), this is later changed towards a broken powerlaw distribution (upper panel) and, for larger Alfvén velocities, a distribution that shows a powerlaw feature only in a very narrow energy range. In all simulation runs, the spectral indices have been fitted to the steepest range. The notable result is that in the case of continuous injection the spectrum is only moderately flattened even though newly injected particles are present at all times.
The maximum particle rigidity is depicted in Fig. 3, showing a steady increase due to stochastic acceleration. The only exception is found for low Alfvén rigidity in the case of continuous injection, where the average rigidity remains constant for some time. Again, it is the Alfvén rigidity that has the prominent influence on the resulting average particle energy, at least qualitatively, and not the average particle lifetime as distinguished by continuous versus single injection.
Fig. 4 Numerically calculated running momentum diffusion coefficient for the cases of R_{A} = 10^{4} (solid line) and R_{A} = 10^{3} (dashed line). We note that for the last part of the lower curve, only 10^{3} particles could be used, while the rest was obtained from a simulation with 10^{4} particles. The two diamonds illustrate the quasilinear result from Eq. (B.2)in Appendix B. 
The associated running momentum diffusion coefficient determined through Eq. (6)is shown in Fig. 4. We note that the time is normalized as v_{A}t/ℓ_{0}, which illustrates that for smaller Alfvén velocities larger times are needed until the momentum diffusion coefficient reaches its asymptotic limit. It is precisely the fact that almost a 1/t time dependence can be seen before finally a transition to a constant (diffusive) value occurs that makes it so numerically difficult to obtain a reliable value for the momentum diffusion coefficient. The maximum simulation time v_{A}t/ℓ_{0} = 10^{2} corresponds to Ωt = 10^{5} and Ωt = 10^{6}, respectively, which illustrates that for the parameters chosen here, scattering in momentum takes an extremely long time until the diffusive limit is reached^{2}.
5. Analytical stochastic acceleration
In principle, the stochastic variations in the particle energy may be positive or negative. However, as shown in Sect. 4, the net effect was a significant energy gain of the particle ensemble. To verify this finding, we present an approach to this problem by directly solving the equation of momentum diffusion. The behavior of charged particles in turbulent electromagnetic fields has to be described using a FokkerPlanck approach (Parker 1965; Schlickeiser 2002) or at least a diffusionconvection equation. Usually, the terms that describe spatial diffusion are of higher magnitude than momentumdiffusion terms. For particular parameter values analytical solutions are feasible (e.g., Büsching et al. 2005; Pohl 2009), whereas a full solution requires numerical tools (Guyer et al. 2009).
However, if concentration is to be focused on stochastic energy changes only (e.g., Sect. 3.3 in Tverskoi 1967; Blandford & Eichler 1987; Ostrowski & SiemieniecOzi¸ebło 1997), it is sometimes convenient to neglect spatial diffusion altogether. Under certain conditions^{3} the diffusionconvection equation can then be expressed as (10)which is a linear parabolic partial differential equation. In principle, D_{p} depends also on time, which however will be neglected here for reasons of simplicity. We note that, by defining a dimensionless rigidity diffusion coefficient , Eq. (10)can be written in dimensionless form by using the normalized variables introduced in Sect. 3.
For Alfvén waves, the diffusion coefficient (see, e.g., Schlickeiser 1985, 1989a; Stawarz & Petrosian 2008) reads (11)with R_{L} the Larmor radius, so that D_{p} ∝ p^{s}. It should be noted however that Eq. (11)was derived only for a simple turbulence spectrum of the form G(k) ∝ k^{−s} with s ∈ [1,2] (for a more detailed discussion see, e.g., Schlickeiser 2002, Sect. 13.1.3).
Consider the diffusion coefficient written as D_{p} = D_{0}p^{q} and the initial distribution written in the form f(p,a) = f_{0}p^{−a}. Then Eq. (10)is given as (12)In general, there are two options. First, when f is initially prescribed as the power law in Eq. (8a), then a Fourier approach can be used as shown in Appendix A. This allows the reduction of the momentum diffusion equation to an ordinary differential equation but has the disadvantage that the Fourier transform of the initial distribution function must be known. Second, when f is initially unknown, a selfsimilar solution can be sought, which is based on the separation f(t,R) = S(t)H(p/E(t)) and will be discussed in the next subsection.
5.1. Selfsimilar solution
Consider again Eq. (12). For a selfsimilar solution one seeks a solution in the form (13)where S, H, and E need to be determined. With ζ = p/E(t), Eq. (12)can then be written (14)which, by using Eq. (13), yields
5.1.1. Solution for general q ≠ 2
Collecting terms, one has (15)which is separable when (the coefficients C_{1}, C_{2}, and C_{3} are constants) so that (16c)Then additionally one must have (16d) which integrates trivially so that one obtains S(t) in the form (17)and one has (18)For further simplification use the fact that the distribution function is normalized, i.e., (19)so that S(t)^{1} = E(t)^{3}H_{0}, yielding Ṡ/S = −3Ė/E and Finally, the resulting differential equation for H(t) reads (20)which needs a numerical solution except for special values of q such as q = 2.
5.1.2. Solution for q = 2
For instance when q = 2, set H = ζ^{Λ} so that (21)yielding (22)Now at t = 0, the distribution function f is discontinuous; therefore, let H = ζ^{−a} in and H = 0 otherwise, where ζ_{min} and ζ_{max} are fixed values to be related eventually to p_{min} and p_{max}.
Fig. 5 Evolution of the momentum distribution function as obtained from the selfsimilar solution of the momentum diffusion equation for the case q = 2. The three lines illustrate the distribution at the times t = 0 (black solid line), t = 1 (blue dashed line), and t = 2 (red dotdashed line). (This figure is available in color in electronic form.) 
The normalization requires Hence, S(t)E(t)^{3} must be constant, which again relates the coefficients C_{3} = −3C_{1}. Then from Eq. (21)with Λ = −a one has C_{1} = D_{0}a.
The solution for the distribution function is then in . From the values at t = 0 one therefore has ζ_{min,max} = p_{min,max}(t = 0).
From the normalization condition (Eq. (23c)) one can choose E_{0} = 1 so that (25)Then the distribution function f(t,p) has the same normalization for all times (as must be the case); however, the minimum and maximum momenta increase as time increases. In Fig. 5, the distribution function is shown for different times, illustrating the shift of the particles toward higher momenta.
We note that Eq. (12)shows that for q = 2 there is scale independence of the momentum (as can easily be seen by replacing p by Lp where L is an arbitrary constant). Equation (12)is unaltered by such a replacement. Thus if the initial distribution has no scale, as with a power law, then the distribution never develops a scale and retains the power law. If, however, the initial distribution is terminated at specific values of momenta then these end values must adjust with time to maintain the scale independence and also to maintain a constant particle population. This anomalous result happens only for q = 2; for any other value of q there is an overall scale for the momentum in Eq. (12)so that the particle evolution with time will end up not only with adjusted end point values but also with a change in the spectral index, as indeed is shown by the numerical illustrations. Furthermore, we note that the case q = 2 agrees with the socalled hardsphere scattering (e.g., Parker & Tidman 1958; Park & Petrosian 1995) and thus possesses a tangible physical interpretation.
Finally, the average momentum gain per time can be obtained using (26)where The relative momentum change is shown in Fig. 6, thereby illustrating that particles are efficiently accelerated in agreement with the simulation results shown in Sect. 4. Therefore, it has been confirmed that, on average, particles indeed gain energy, thus confirming the numerical ansatz.
6. Discussion and conclusion
Fig. 6 Relative average momentum gain as described through Eq. (27)as a function of time. The reference value is the initial average momentum, i.e., ⟨p_{0}⟩ = ⟨p(t = 0)⟩. 
In this paper, the stochastic acceleration of energetic charged particles, due to momentum diffusion in Alfvénic electromagnetic turbulence, has been investigated. The initial particle distribution function was prescribed as a power law starting with the Alfvén velocity and covering three orders of magnitude. Even though it was found (Schlickeiser 1989a) that the dominant contribution to the stochastic acceleration of charged particles is provided by fast magnetosonic waves and not by Alfvén waves, here the case of Alfvén waves was chosen because these waves constitute the dominant timedependent turbulent component in the solar wind. Throughout, radiative losses (Schlickeiser & Lerche 2007; Stawarz & Petrosian 2008; SchaeferRolffs et al. 2009) have been neglected as only waveparticle interactions were included via the Lorentz force (cf. Tautz 2010a).
Two independent methods were used:

A MonteCarlo simulation, which, in the presence of artificialturbulence obtained from the superposition of slab Alfvénwaves, integrates the NewtonLorentz equation for a largenumber of test particles. Additionally, the cases weredistinguished of all particles being injected (i) at the initial timet = 0; or (ii) continuously throughout the simulation with probabilities given by the initial distribution function.

By neglecting spatial diffusion, the momentumdiffusion equation has been analyzed using a selfsimilar approach. The solution allows one to predict the time evolution of the distribution function. A general method to obtain a solution by Fourier transforming the momentumdiffusion equation is presented in Appendix A.
By monitoring the average momentum value, it was confirmed in both cases that particles are efficiently accelerated. Furthermore, the simulation results show that, for a distribution function that is initially of the form f ∝ p^{−a}, the spectral index a is significantly increased for later times. The steepened spectrum is in remarkable agreement with p^{5} spectra observed in various scenarios such as pickup ion distributions (Fahr & Fichtner 2011), terminationshock particles (Jokipii & Lee 2010), and the solar corona (quiet solar wind particles; Fisk & Gloeckler 2008; Antecki et al. 2013). We note that, because steep spectra take considerably more computing time, a relatively flat initial spectrum has been chosen; however, the agreement with observations is of itself no guarantee that this agreement cannot also be achieved with a steep injected spectrum. Demonstrating the validity of this point is, however, best delayed to a future communication. In addition, here the acceleration is due to homogeneous turbulent electromagnetic fields generated by (unstable) plasma waves. This steepening of the spectrum could not be reproduced in the analytical selfsimilar solution based on a momentum diffusion coefficient of the form D_{p} ∝ p^{q} with q = 2 for reasons given in the text. This spectral index change requires a departure from a scaleindependent momentum behavior for Eq. (12), something that occurs for all values of q except q = 2. Therefore, the latter case has been chosen for reasons of simplicity.
Furthermore, the results depend on the details of the momentumdiffusion coefficient D_{p} for which a simplified form is assumed. As detailed in Appendix B, a more general, nonlinear approach shows that D_{p} exhibits a more complicated dependence on the Alfvén velocity and on the turbulence power spectrum. In addition, the cases of single and continuous particle injection have been considered in the simulation. Whereas the second case might be more realistic, the first case is useful in that, by superposition of different time injection solutions, a continuous injection model can be constructed where the source term varies with time. This, however, is beyond the scope of the present investigation.
Future work should attempt to modify the analytical selfsimilar solution so that the effect of a variable energy spectral index is incorporated as predicted by the numerical simulations. Generally, a much larger number of test particles than is currently available is required if the evolution of a distribution function is to be investigated for a broader momentum range and/or for a considerably steeper initial spectrum. This is especially true when the particles are continuously injected so that the evolution of a quasiequilibrium can be monitored. Furthermore, real physical scenarios such as particle acceleration at the outer regions of the solar system require the inclusion of spatial structures, thereby giving rise to spatially variable diffusion coefficients. The investigation of this problem will be deferred to a future paper.
At first glance is might seem puzzling that diffusion coefficients are not constant even though there is turbulence from the start. A long time ago, when comparing random particle motions with diffusive processes, Chandrasekhar (1943) found that truly diffusionlike behavior the corresponding equivalent diffusion coefficient is timedependent and reaches an asymptotic constant value only later in time.
Of course, “long” is a relative term as 10^{6} gyro periods is a short time compared to the dynamical time scales in most astronomical sources. However, both computationally and with respect to spatial diffusion – which takes about a factor of 10^{4} fewer gyro periods – that time scale is indeed extremely long.
For example, the principle of detailed balance has to be fulfilled, which requires an equilibrium state (e.g., Klein 1955). Specifically, in the present case it states that the probability Ψ(p,p′) for changing the particle momentum from p to p′ satisfies the condition Ψ(p,p′) = Ψ(p′,p).
Acknowledgments
F.K. thanks D. Breitschwerdt for support during his M. Sc. thesis. We are also appreciative of comments by a referee that helped us sharpen some otherwise fuzzy points.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Science, 327, 1103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Achterberg, A. 1981, A&A, 97, 259 [NASA ADS] [Google Scholar]
 Antecki, T., Schlickeiser, R., & Zhang, M. 2013, ApJ, 764, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Belcher, J. W., & Davis, L. 1971, J. Geophys. Res., 76, 3534 [NASA ADS] [CrossRef] [Google Scholar]
 Bieber, J. W., Matthaeus, W. H., Smith, C. W., et al. 1994, ApJ, 420, 294 [Google Scholar]
 Blandford, R. D., & Eichler, D. 1987, Phys. Rep., 154, 1 [Google Scholar]
 Bruno, R., & Carbone, V. 2005, Liv. Rev. Sol. Phys., 2, 1 [Google Scholar]
 Büsching, I., Kopp, A., Pohl, M., et al. 2005, ApJ, 619, 314 [NASA ADS] [CrossRef] [Google Scholar]
 Buts, V. A., Manuilenko, O. V., Tolstoluzhsky, A. P., & Turkin, Y. A. 2001, Prob. Atomic Sci. Tech., 2001, 92 [Google Scholar]
 Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Dosch, A., & Shalchi, A. 2010, Adv. Space Res., 46, 1208 [NASA ADS] [CrossRef] [Google Scholar]
 Dröge, W. 2005, Adv. Space Res., 35, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. O. 1983, Rep. Prog. Phys., 46, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Duffy, P., & Blundell, K. M. 2005, Plasma Phys. Contr. Fusion, 47, B667 [CrossRef] [Google Scholar]
 Fahr, H., & Fichtner, H. 2011, A&A, 533, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fatuzzo, M., & Melia, F. 2011, MNRAS, 410, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Fatuzzo, M., & Melia, F. 2012, ApJ, 750, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Fedorov, Y. I., & Stehlik, M. 2008, in Proc. 21st European Cosmic Ray Symposium, eds. P. Király, K. Kudela, M. Stehlík, & A. W. Wolfendale (Košice, Slovakia: Inst. Exp. Phys., Slovak Academy of Sciences), 241 [Google Scholar]
 Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Fisk, L. A. 1996, J. Geophys. Res., 101, 15547 [Google Scholar]
 Fisk, L. A., & Gloeckler, G. 2006, ApJ, 640, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Fisk, L. A., & Gloeckler, G. 2007a, Space Sci. Rev., 130, 153 [Google Scholar]
 Fisk, L. A., & Gloeckler, G. 2007b, Proc. Nat. Acad. Sci., 104, 5749 [Google Scholar]
 Fisk, L. A., & Gloeckler, G. 2008, ApJ, 686, 1466 [NASA ADS] [CrossRef] [Google Scholar]
 Fisk, L. A., Gloeckler, G., & Schwadron, N. A. 2010, ApJ, 720, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204 [Google Scholar]
 Guyer, J., Wheeler, D., & Warren, J. A. 2009, Comp. Sci. Eng., 11, 6 [Google Scholar]
 Hamilton, R. J., & Petrosian, V. 1992, ApJ, 398, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Hartquist, T. W., & Morfill, G. E. 1983, ApJ, 266, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Heinbach, U., & Simon, M. 1995, ApJ, 441, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Jokipii, J. R. 1966, ApJ, 146, 480 [Google Scholar]
 Jokipii, J. R., & Lee, M. A. 2010, ApJ, 713, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Klein, M. J. 1955, Phys. Rev. Lett., 97, 1466 [Google Scholar]
 Kolmogorov, A. N. 1991, Proc. Royal Soc. London, Ser. A: Math. Phys. Sci., 434, 9 [Google Scholar]
 Kuramitsu, Y., Nakanii, N., Kondo, K., et al. 2011, Phys. Plasmas, 18, 010701 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, M. A., & Völk, H. J. 1975, ApJ, 198, 485 [NASA ADS] [CrossRef] [Google Scholar]
 LetessierSelvon, A., & Stanev, T. 2011, Rev. Mod. Phys., 83, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Litvinenko, Y. E., & Schlickeiser, R. 2011, ApJ, 732, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Longair, M. S. 2011, High Energy Astrophysics (Cambridge: University Press) [Google Scholar]
 Malkov, M. A., Diamond, P., Drury, L. O., & Sagdeev, R. Z. 2010, ApJ, 721, 750 [NASA ADS] [CrossRef] [Google Scholar]
 Matthaeus, W. H., Goldstein, M. L., & Roberts, D. A. 1990, J. Geophys. Res., 95, 20673 [Google Scholar]
 Michałek, G., Ostrowski, M., & Schlickeiser, R. 1999, Sol. Phys., 184, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Moraal, H., CaballeroLopez, R. A., & Ptuskin, V. S. 2008, in Proc. 3th International Cosmic Ray conference, eds. R. Caballero, J. C. D’Olivio, G. MedinaTanco, L. Nellen, F. A. Sánchez, & J. F. ValdésGalicia (Mexico City, Mexico: Universidad Nacional Autónoma de México), 1, 865 [Google Scholar]
 Nagano, M., & Watson, A. A. 2000, Rev. Mod. Phys., 72, 689 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ostrowski, M., & Michałek, G. 1993, in 23rd International Cosmic Ray Conference, held 19–30 July, at University of Calgary, Alberta, Canada, eds. D. A. Leahy, R. B. Hickws, & D. Venkatesan (Singapore: World Scientific), 2, 309 [Google Scholar]
 Ostrowski, M., & SiemieniecOzi¸ebło, G. 1997, Astropart. Phys., 6, 271 [NASA ADS] [CrossRef] [Google Scholar]
 O’Sullivan, S., Reville, B., & Taylor, A. M. 2009, MNRAS, 400, 248 [NASA ADS] [CrossRef] [Google Scholar]
 Park, J., & Petrosian, V. 1995, ApJ, 446, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
 Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
 Parker, E. N., & Tidman, D. A. 1958, Phys. Rev., 111, 1206 [NASA ADS] [CrossRef] [Google Scholar]
 Pohl, M. 2009, Phys. Rev. D, 79, 041301 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes (Cambridge: University Press) [Google Scholar]
 Rausch, M., & Tautz, R. C. 2012, MNRAS, 428, 2333 [Google Scholar]
 SchaeferRolffs, U., Lerche, I., & Tautz, R. C. 2009, J. Phys. A: Math. Theor., 42, 105501 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1985, A&A, 143, 431 [NASA ADS] [Google Scholar]
 Schlickeiser, R. 1989a, ApJ, 336, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1989b, ApJ, 336, 246 [NASA ADS] [Google Scholar]
 Schlickeiser, R. 1994, ApJS, 20, 929 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin: Springer) [Google Scholar]
 Schlickeiser, R., & Lerche, I. 2007, A&A, 476, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R., Artmann, S., & Dröge, W. 2009, Open Plasma Phys. J., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A. 2005, Phys. Plasmas, 12, 052905 [Google Scholar]
 Shalchi, A. 2006, MNRAS, 371, 1898 [NASA ADS] [CrossRef] [Google Scholar]
 Shalchi, A., & Schlickeiser, R. 2004, A&A, 420, 799 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shalchi, A., Bieber, J. W., Matthaeus, W. H., & Qin, G. 2004, ApJ, 616, 617 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Stawarz, Ł., & Petrosian, V. 2008, ApJ, 681, 1725 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C. 2010a, Comput. Phys. Commun., 81, 71 [Google Scholar]
 Tautz, R. C. 2010b, Plasma Phys. Contr. Fusion, 52, 045016 [Google Scholar]
 Tautz, R. C. 2012, in Turbulence: Theory, Types and Simulation, ed. R. J. Marcuso (New York: Nova Publishers), 365 [Google Scholar]
 Tautz, R. C., & Lerche, I. 2010, Phys. Lett. A, 374, 4573 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., & Shalchi, A. 2010, J. Geophys. Res., 115, A03104 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2006a, J. Phys. G: Nuclear Part. Phys., 32, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2006b, J. Phys. G: Nuclear Part. Phys., 32, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2008, ApJ, 685, L165 [NASA ADS] [CrossRef] [Google Scholar]
 Tautz, R. C., Dosch, A., & Lerche, I. 2012, A&A, 545, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tverskoi, B. A. 1967, J. Exp. Theor. Phys., 25, 317 [NASA ADS] [Google Scholar]
 Vainio, R., & Schlickeiser, R. 1999, A&A, 343, 303 [NASA ADS] [Google Scholar]
 Vocks, C., Salem, C., & Lin, R. P. 2005, ApJ, 627, 540 [NASA ADS] [CrossRef] [Google Scholar]
 Vocks, C., Mann, G., & Rausche, G. 2008, A&A, 480, 527 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weinhorst, B., & Shalchi, A. 2010, MNRAS, 403, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Zank, G., Li, G., Florinski, V., et al. 2006, J. Geophys. Res., 111, A06108 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, M., & Lee, M. A. 2011, Space Sci. Rev., DOI: 10.1007/s1121401197543 [Google Scholar]
Appendix A: Fourier solution of the momentumdiffusion equation
The analytical solution of Eq. (10)for the initial distribution from Eq. (8a)proceeds as follows. There are two alternatives: (i) a direct solution via Fourier transformation; and (ii) a selfsimilar solution approach. While the second approach is discussed in the text, here some basic considerations are shown regarding the first approach.
Set dξ = p^{−(q + 2)}dp so that ξ = p^{−(q + 1)}/(q + 1) when Eq. (12)takes the form (A.1)With the Fourier transform dependence in time proportional to e^{−iωt}, Eq. (A.1)yields (A.2)where F(ω,ξ) is the Fourier transform of f. Now so that Eq. (A.2)can be written With α = −iω/D_{0} [−(q + 1)] ^{−(q + 4)/(q + 1)} and m = (q + 4)/(q + 1) one has (A.3)which has the solution (A.4)where J_{n}(z) is the Bessel function of the first kind of order n and where F_{0}(ω) is constant in ξ. Then (A.5)We note that 2 − m = (q − 2)/(q + 1).
Suppose that F_{0}(ω) is a simple power in ω, i.e., F_{0}(ω) = f_{0}(−iω)^{γ}. Now set α = α_{0}(−iω) so that (A.6)Set −iωξ^{2−m} ≡ σ^{2} so that 2σ dσ = −idω ξ^{2−m} and σ = ± ξ^{(2−m)/2}e^{3iπ/4}ω^{1/2}. Careful inspection of the integration limits shows that The resulting inverse Fourier integral of the distribution function then reads (A.8)which, at t = 0, has a ξ dependence as f(0,ξ) ∝ ξ^{(2γ + 1)(m − 2)}. Now originally, f(0,p) is given in the form f ∝ p^{−a} and ξ ∝ p^{−(q + 1)}. From Eq. (A.8)it then follows that (A.9)which connects the expression for f with the frequency powerlaw F_{0}(ω) ∝ ω^{γ}. And so one has the solution, provided no limits are set on the p^{−a} behavior. Furthermore, for finite times, t > 0, one then has a relatively simple integral for f(t,ξ) from Eq. (A.8). We note also that, if F_{0}(ω) = ∑ f_{n}ω^{n} then superposition of inverse Fourier integrals of the above type gives the general solution. Thus for any f(p) that is expandible in orthonormal functions expressed as power combinations through a GramSchmidt orthonormalization process one has the solution by superposition.
Appendix B: Nonlinear momentum diffusion
Fig. B.1 Quasilinear result for the momentum diffusion coefficient (solid line). The dotted line is the powerlaw fit, which yields D_{p} ∝ p^{0.4s}. 
With ε ≡ R_{A}/R, the general expression for the FokkerPlanck coefficient of momentum diffusion is of the form (B.1)which was derived under the assumption of slab linearly polarized, undamped Alfvén waves (Schlickeiser 2002).
Using quasilinear theory (Jokipii 1966), the momentum diffusion coefficient for Alfvén turbulence with vanishing cross helicity is given through (Schlickeiser 2002) (B.2)where (B.3a)ensures the normalization of the spectrum from Eq. (1), from which only the wavenumberdependent part (B.3b)is used here. For the spectral index, s = 5/3 is chosen in accordance with Eq. (1). In the limit R ≪ 1, the resulting momentum diffusion coefficient can be fitted by a power law D_{p} ∝ p^{0.4s} (see Fig. B.1).
To obtain an average momentum diffusion coefficient for a distribution function from Eq. (8a), an expectation value has to be calculated as (B.4)For the parameters given in Table 1, one has (B.5)which is in approximate agreement with the simulation results shown in Fig. 4.
A recent investigation of Shalchi (2006) using the weakly nonlinear theory (Shalchi et al. 2004) showed that, compared to previous, quasilinear results (Shalchi & Schlickeiser 2004), the velocitydependence of the momentum diffusion coefficients is flattened. Here, the secondorder quasilinear theory (SOQLT, see Shalchi 2005; Tautz et al. 2008; Tautz & Lerche 2010) is used to derive the nonlinear equivalent to judge the resulting differences.
While in QLT, particles are assumed to follow the unperturbed spiral trajectory, SOQLT takes into account the distribution of particles around their unperturbed orbits through resonance broadening. For the Fourier transform, the new approach involves deviations from the unperturbed orbit, , as (B.6)where a Gaussian shape is assumed for the distribution function so that (B.7)with the (vanishing) mean of the unperturbed parallel coordinate.
The mean square parallel coordinate ⟨ z^{2}(t) ⟩ that enters the width of the Gaussian as is obtained using QLT, thereby justifying the designation of the resulting theory as being of the second order. In the limit of large times t ≫ Ω^{1} and for pitchangles close to 90°, simple analytical expressions can be obtained for ⟨ z^{2}(t) ⟩ , thereby modifying the resonance function. This has a deep impact on some transport parameters. For instance, the mean free path in isotropic magnetostatic turbulence (Tautz et al. 2006a, 2008) is infinitely large in QLT, but in SOQLT agrees well with simulations (Tautz 2012).
For the momentumdiffusion coefficient, the result reads (B.8a)where (B.8b)is the secondorder resonance function. The integrals in Eq. (B.8a)cannot be solved analytically so that a numerical solution, which in this case is intricate, has to be evoked.
The result for D_{p} is shown in Fig. B.2 in comparison to results obtained from a testparticle simulation and the quasilinear result from Eq. (B.2). At first view, the SOQLT result seems to be an improvement compared to QLT for rigidities larger than R = 0.0015 but one has to be careful before coming to any conclusions. A conspicuous feature of the values obtained with SOQLT can be seen at the lower rigidity range. The values exhibit a discrepancy in comparison to the simulation, because the slope of the curve is positive for rigidities smaller than R = 0.0015. As was pointed out in the preceding, the simulation in this parameter regime provides the most reliable results, since for small rigidities it reached almost diffusive behavior. Therefore, one would expect a decreasing slope everywhere and the values for the FokkerPlanck coefficient to be considerably higher in this range. In the calculation, the Alfvén rigidity is set to R_{A} = 10^{3} in the calculation, hence the approximation is no longer valid in this regime.
Fig. B.2 Comparison of the momentum diffusion coefficient D_{p} as obtained from numerical testparticle simulations (triangles), quasilinear theory (dashed line), and secondorder quasilinear theory (crosses). 
There is an additional feature of the simulations that has to be taken into account. As shown in Fig. 4, the momentum diffusion coefficient decreases for a longer time if the ratio ε = R_{A}/R is decreased, i.e., for larger initial particle rigidities. Therefore, longer simulation times would be required in principle to ensure that the diffusive regime is reached for all initial rigidities. Therefore, the values of the momentum diffusion coefficient are overestimated for larger initial rigidities. However, quantitative statements about the expected reduction are not possible; instead, simulations with longer running times would be required, which are beyond what is currently available.
All Tables
All Figures
Fig. 1 Evolution of the particle distribution as a function of the rigidity R as the normalized time Ωt increases. For the Alfvén rigidity, values of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel) have been used. (This figure is available in color in electronic form.) 

In the text 
Fig. 2 Energy spectral index for the cases of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel). Shown is the energy distribution function at the beginning (black lines) and at the end of the simulation for the cases of single injection (red lines) and continuous particle injection (blue lines). Also shown are the spectral indices obtained from powerlaw fits (dashed lines) to the simulation data (thin lines). (This figure is available in color in electronic form.) 

In the text 
Fig. 3 Maximum of the distribution function as the time increases. Again the cases of R_{A} = 10^{4} (upper panel) and R_{A} = 10^{3} (lower panel) are shown. The black and red lines denote the cases of single injection and continuous particle injection, respectively. The thick lines simply show cubic fits of the actual simulation data. Additionally, the average particle rigidity is shown in the lower panel (dotted line). (This figure is available in color in electronic form.) 

In the text 
Fig. 4 Numerically calculated running momentum diffusion coefficient for the cases of R_{A} = 10^{4} (solid line) and R_{A} = 10^{3} (dashed line). We note that for the last part of the lower curve, only 10^{3} particles could be used, while the rest was obtained from a simulation with 10^{4} particles. The two diamonds illustrate the quasilinear result from Eq. (B.2)in Appendix B. 

In the text 
Fig. 5 Evolution of the momentum distribution function as obtained from the selfsimilar solution of the momentum diffusion equation for the case q = 2. The three lines illustrate the distribution at the times t = 0 (black solid line), t = 1 (blue dashed line), and t = 2 (red dotdashed line). (This figure is available in color in electronic form.) 

In the text 
Fig. 6 Relative average momentum gain as described through Eq. (27)as a function of time. The reference value is the initial average momentum, i.e., ⟨p_{0}⟩ = ⟨p(t = 0)⟩. 

In the text 
Fig. B.1 Quasilinear result for the momentum diffusion coefficient (solid line). The dotted line is the powerlaw fit, which yields D_{p} ∝ p^{0.4s}. 

In the text 
Fig. B.2 Comparison of the momentum diffusion coefficient D_{p} as obtained from numerical testparticle simulations (triangles), quasilinear theory (dashed line), and secondorder quasilinear theory (crosses). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.