EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 555 (July 2013) A&A, 555 (2013) A101 Full HTML
Free Access
Issue
A&A
Volume 555, July 2013
Article Number A101
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201220966
Published online 09 July 2013

© ESO, 2013

1. Introduction

Understanding the acceleration processes of cosmic rays has been a long-standing problem in astrophysics (e.g., Drury 1983; Blandford & Eichler 1987). It is currently accepted that high-energy particles with energies up to 1015 eV are accelerated in supernova remnants (e.g., Abdo et al. 2010; Hartquist & Morfill 1983) via diffusive shock acceleration, thus producing a power-law 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; Letessier-Selvon & Stanev 2011) and with laser experiments (e.g., Kuramitsu et al. 2011). Particles with low energies up to 1010 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). Energy-dependent 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 outward-propagating Alfvénic turbulence, which has been named the slab component. However, Alfvénic waves have also been assumed for the investigation of cosmic-ray 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 cosmic-ray 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, large-scale 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 & Siemieniec-Ozi¸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 quasi-linear theory (Shalchi & Schlickeiser 2004).

In the presence of shock waves (Schlickeiser 1989b), the first-order 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 test-particle 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 energy-dependence 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 cosmic-ray spectrum (Fisk & Gloeckler 2006). Therefore, we intent to demonstrate that Monte-Carlo test-particle simulations are a suitable tool for tracking the time evolution of a particle velocity distribution or an energy-spectrum.

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 test-particle code PADIAN is introduced and the simulations results are presented, respectively. In Sect. 5, an analytical solution of the momentum-diffusion 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 large-scale 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 bend-over 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 three-dimensional 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 B0 as δB = δB(z) for ; (ii) the composite model, which superposes the slab field by a two-dimensional complement as δB = δBslab(z) + δB2D(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 three-dimensional turbulence, something that is beyond the scope of the present investigation.

3. Monte-Carlo simulations

For the test-particle simulations, the test-particle Monte-Carlo 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 kn 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 Newton-Lorentz 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 RA ≡ vA/(ℓ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 time-dependent 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 so-called 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 RA/R, the sooner Dp attains an asymptotic value1.

3.3. Rigidity distribution

In most test-particle simulations, a number of discrete initial particle rigidities are assumed which are then investigated separately, which is acceptable when one attempts to investigate the rigidity-dependence 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 cosmic-ray observations (see Letessier-Selvon & 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 high-energy 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, Rmin and Rmin are chosen so that particles with rigidities comparable to and much larger than RA are present but, at the same time, the rigidity interval is kept moderately narrow.

Table 1

Parameters used for the Monte-Carlo simulation.

thumbnail 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 RA = 10-4 (upper panel) and RA = 10-3 (lower panel) have been used. (This figure is available in color in electronic form.)

thumbnail Fig. 2

Energy spectral index for the cases of RA = 10-4 (upper panel) and RA = 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 power-law fits (dashed lines) to the simulation data (thin lines). (This figure is available in color in electronic form.)

thumbnail Fig. 3

Maximum of the distribution function as the time increases. Again the cases of RA = 10-4 (upper panel) and RA = 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 power-law 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 power-law distribution (upper panel) and, for larger Alfvén velocities, a distribution that shows a power-law 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 life-time as distinguished by continuous versus single injection.

thumbnail Fig. 4

Numerically calculated running momentum diffusion coefficient for the cases of RA = 10-4 (solid line) and RA = 10-3 (dashed line). We note that for the last part of the lower curve, only 103 particles could be used, while the rest was obtained from a simulation with 104 particles. The two diamonds illustrate the quasi-linear 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 vAt/ℓ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 vAt/ℓ0 = 102 corresponds to Ωt = 105 and Ωt = 106, respectively, which illustrates that for the parameters chosen here, scattering in momentum takes an extremely long time until the diffusive limit is reached2.

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 Fokker-Planck approach (Parker 1965; Schlickeiser 2002) or at least a diffusion-convection equation. Usually, the terms that describe spatial diffusion are of higher magnitude than momentum-diffusion 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 & Siemieniec-Ozi¸ebło 1997), it is sometimes convenient to neglect spatial diffusion altogether. Under certain conditions3 the diffusion-convection equation can then be expressed as (10)which is a linear parabolic partial differential equation. In principle, Dp 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 RL the Larmor radius, so that Dp ∝ ps. 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 Dp = D0pq and the initial distribution written in the form f(p,a) = f0p−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 self-similar 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. Self-similar solution

Consider again Eq. (12). For a self-similar 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 C1, C2, and C3 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)3H0, 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 pmin and pmax.

thumbnail Fig. 5

Evolution of the momentum distribution function as obtained from the self-similar 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 C3 = −3C1. Then from Eq. (21)with Λ = −a one has C1 = D0a.

The solution for the distribution function is then in . From the values at t = 0 one therefore has ζmin,max = pmin,max(t = 0).

From the normalization condition (Eq. (23c)) one can choose E0 = 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 so-called hard-sphere 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

thumbnail 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., ⟨p0⟩ = ⟨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 time-dependent turbulent component in the solar wind. Throughout, radiative losses (Schlickeiser & Lerche 2007; Stawarz & Petrosian 2008; Schaefer-Rolffs et al. 2009) have been neglected as only wave-particle interactions were included via the Lorentz force (cf. Tautz 2010a).

Two independent methods were used:

  • A Monte-Carlo simulation, which, in the presence of artificialturbulence obtained from the superposition of slab Alfvénwaves, integrates the Newton-Lorentz 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 momentum-diffusion equation has been analyzed using a self-similar 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 momentum-diffusion 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 pick-up ion distributions (Fahr & Fichtner 2011), termination-shock 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 self-similar solution based on a momentum diffusion coefficient of the form Dp ∝ pq with q = 2 for reasons given in the text. This spectral index change requires a departure from a scale-independent 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 momentum-diffusion coefficient Dp for which a simplified form is assumed. As detailed in Appendix B, a more general, non-linear approach shows that Dp 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 self-similar 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 quasi-equilibrium 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.

1

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 diffusion-like behavior the corresponding equivalent diffusion coefficient is time-dependent and reaches an asymptotic constant value only later in time.

2

Of course, “long” is a relative term as 106 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 104 fewer gyro periods – that time scale is indeed extremely long.

3

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

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Science, 327, 1103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Achterberg, A. 1981, A&A, 97, 259 [NASA ADS] [Google Scholar]
  3. Antecki, T., Schlickeiser, R., & Zhang, M. 2013, ApJ, 764, 89 [NASA ADS] [CrossRef] [Google Scholar]
  4. Belcher, J. W., & Davis, L. 1971, J. Geophys. Res., 76, 3534 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bieber, J. W., Matthaeus, W. H., Smith, C. W., et al. 1994, ApJ, 420, 294 [Google Scholar]
  6. Blandford, R. D., & Eichler, D. 1987, Phys. Rep., 154, 1 [Google Scholar]
  7. Bruno, R., & Carbone, V. 2005, Liv. Rev. Sol. Phys., 2, 1 [Google Scholar]
  8. Büsching, I., Kopp, A., Pohl, M., et al. 2005, ApJ, 619, 314 [NASA ADS] [CrossRef] [Google Scholar]
  9. Buts, V. A., Manuilenko, O. V., Tolstoluzhsky, A. P., & Turkin, Y. A. 2001, Prob. Atomic Sci. Tech., 2001, 92 [Google Scholar]
  10. Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 [NASA ADS] [CrossRef] [Google Scholar]
  11. Dosch, A., & Shalchi, A. 2010, Adv. Space Res., 46, 1208 [NASA ADS] [CrossRef] [Google Scholar]
  12. Dröge, W. 2005, Adv. Space Res., 35, 532 [NASA ADS] [CrossRef] [Google Scholar]
  13. Drury, L. O. 1983, Rep. Prog. Phys., 46, 973 [NASA ADS] [CrossRef] [Google Scholar]
  14. Duffy, P., & Blundell, K. M. 2005, Plasma Phys. Contr. Fusion, 47, B667 [CrossRef] [Google Scholar]
  15. Fahr, H., & Fichtner, H. 2011, A&A, 533, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Fatuzzo, M., & Melia, F. 2011, MNRAS, 410, L23 [NASA ADS] [CrossRef] [Google Scholar]
  17. Fatuzzo, M., & Melia, F. 2012, ApJ, 750, 21 [NASA ADS] [CrossRef] [Google Scholar]
  18. 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]
  19. Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  20. Fisk, L. A. 1996, J. Geophys. Res., 101, 15547 [Google Scholar]
  21. Fisk, L. A., & Gloeckler, G. 2006, ApJ, 640, L79 [NASA ADS] [CrossRef] [Google Scholar]
  22. Fisk, L. A., & Gloeckler, G. 2007a, Space Sci. Rev., 130, 153 [Google Scholar]
  23. Fisk, L. A., & Gloeckler, G. 2007b, Proc. Nat. Acad. Sci., 104, 5749 [Google Scholar]
  24. Fisk, L. A., & Gloeckler, G. 2008, ApJ, 686, 1466 [NASA ADS] [CrossRef] [Google Scholar]
  25. Fisk, L. A., Gloeckler, G., & Schwadron, N. A. 2010, ApJ, 720, 533 [NASA ADS] [CrossRef] [Google Scholar]
  26. Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204 [Google Scholar]
  27. Guyer, J., Wheeler, D., & Warren, J. A. 2009, Comp. Sci. Eng., 11, 6 [Google Scholar]
  28. Hamilton, R. J., & Petrosian, V. 1992, ApJ, 398, 350 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hartquist, T. W., & Morfill, G. E. 1983, ApJ, 266, 271 [NASA ADS] [CrossRef] [Google Scholar]
  30. Heinbach, U., & Simon, M. 1995, ApJ, 441, 209 [NASA ADS] [CrossRef] [Google Scholar]
  31. Jokipii, J. R. 1966, ApJ, 146, 480 [Google Scholar]
  32. Jokipii, J. R., & Lee, M. A. 2010, ApJ, 713, 475 [NASA ADS] [CrossRef] [Google Scholar]
  33. Klein, M. J. 1955, Phys. Rev. Lett., 97, 1466 [Google Scholar]
  34. Kolmogorov, A. N. 1991, Proc. Royal Soc. London, Ser. A: Math. Phys. Sci., 434, 9 [Google Scholar]
  35. Kuramitsu, Y., Nakanii, N., Kondo, K., et al. 2011, Phys. Plasmas, 18, 010701 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lee, M. A., & Völk, H. J. 1975, ApJ, 198, 485 [NASA ADS] [CrossRef] [Google Scholar]
  37. Letessier-Selvon, A., & Stanev, T. 2011, Rev. Mod. Phys., 83, 907 [NASA ADS] [CrossRef] [Google Scholar]
  38. Litvinenko, Y. E., & Schlickeiser, R. 2011, ApJ, 732, L31 [NASA ADS] [CrossRef] [Google Scholar]
  39. Longair, M. S. 2011, High Energy Astrophysics (Cambridge: University Press) [Google Scholar]
  40. Malkov, M. A., Diamond, P., Drury, L. O., & Sagdeev, R. Z. 2010, ApJ, 721, 750 [NASA ADS] [CrossRef] [Google Scholar]
  41. Matthaeus, W. H., Goldstein, M. L., & Roberts, D. A. 1990, J. Geophys. Res., 95, 20673 [Google Scholar]
  42. Michałek, G., Ostrowski, M., & Schlickeiser, R. 1999, Sol. Phys., 184, 339 [NASA ADS] [CrossRef] [Google Scholar]
  43. Moraal, H., Caballero-Lopez, R. A., & Ptuskin, V. S. 2008, in Proc. 3th International Cosmic Ray conference, eds. R. Caballero, J. C. D’Olivio, G. Medina-Tanco, L. Nellen, F. A. Sánchez, & J. F. Valdés-Galicia (Mexico City, Mexico: Universidad Nacional Autónoma de México), 1, 865 [Google Scholar]
  44. Nagano, M., & Watson, A. A. 2000, Rev. Mod. Phys., 72, 689 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  45. 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]
  46. Ostrowski, M., & Siemieniec-Ozi¸ebło, G. 1997, Astropart. Phys., 6, 271 [NASA ADS] [CrossRef] [Google Scholar]
  47. O’Sullivan, S., Reville, B., & Taylor, A. M. 2009, MNRAS, 400, 248 [NASA ADS] [CrossRef] [Google Scholar]
  48. Park, J., & Petrosian, V. 1995, ApJ, 446, 699 [NASA ADS] [CrossRef] [Google Scholar]
  49. Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
  50. Parker, E. N. 1965, Planet. Space Sci., 13, 9 [Google Scholar]
  51. Parker, E. N., & Tidman, D. A. 1958, Phys. Rev., 111, 1206 [NASA ADS] [CrossRef] [Google Scholar]
  52. Pohl, M. 2009, Phys. Rev. D, 79, 041301 [NASA ADS] [CrossRef] [Google Scholar]
  53. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes (Cambridge: University Press) [Google Scholar]
  54. Rausch, M., & Tautz, R. C. 2012, MNRAS, 428, 2333 [Google Scholar]
  55. Schaefer-Rolffs, U., Lerche, I., & Tautz, R. C. 2009, J. Phys. A: Math. Theor., 42, 105501 [NASA ADS] [CrossRef] [Google Scholar]
  56. Schlickeiser, R. 1985, A&A, 143, 431 [NASA ADS] [Google Scholar]
  57. Schlickeiser, R. 1989a, ApJ, 336, 243 [NASA ADS] [CrossRef] [Google Scholar]
  58. Schlickeiser, R. 1989b, ApJ, 336, 246 [NASA ADS] [Google Scholar]
  59. Schlickeiser, R. 1994, ApJS, 20, 929 [NASA ADS] [CrossRef] [Google Scholar]
  60. Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin: Springer) [Google Scholar]
  61. Schlickeiser, R., & Lerche, I. 2007, A&A, 476, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Schlickeiser, R., Artmann, S., & Dröge, W. 2009, Open Plasma Phys. J., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
  63. Shalchi, A. 2005, Phys. Plasmas, 12, 052905 [Google Scholar]
  64. Shalchi, A. 2006, MNRAS, 371, 1898 [NASA ADS] [CrossRef] [Google Scholar]
  65. Shalchi, A., & Schlickeiser, R. 2004, A&A, 420, 799 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Shalchi, A., Bieber, J. W., Matthaeus, W. H., & Qin, G. 2004, ApJ, 616, 617 [NASA ADS] [CrossRef] [Google Scholar]
  67. Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
  68. Stawarz, Ł., & Petrosian, V. 2008, ApJ, 681, 1725 [NASA ADS] [CrossRef] [Google Scholar]
  69. Tautz, R. C. 2010a, Comput. Phys. Commun., 81, 71 [Google Scholar]
  70. Tautz, R. C. 2010b, Plasma Phys. Contr. Fusion, 52, 045016 [Google Scholar]
  71. Tautz, R. C. 2012, in Turbulence: Theory, Types and Simulation, ed. R. J. Marcuso (New York: Nova Publishers), 365 [Google Scholar]
  72. Tautz, R. C., & Lerche, I. 2010, Phys. Lett. A, 374, 4573 [NASA ADS] [CrossRef] [Google Scholar]
  73. Tautz, R. C., & Shalchi, A. 2010, J. Geophys. Res., 115, A03104 [NASA ADS] [CrossRef] [Google Scholar]
  74. Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2006a, J. Phys. G: Nuclear Part. Phys., 32, 809 [NASA ADS] [CrossRef] [Google Scholar]
  75. Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2006b, J. Phys. G: Nuclear Part. Phys., 32, 1045 [NASA ADS] [CrossRef] [Google Scholar]
  76. Tautz, R. C., Shalchi, A., & Schlickeiser, R. 2008, ApJ, 685, L165 [NASA ADS] [CrossRef] [Google Scholar]
  77. Tautz, R. C., Dosch, A., & Lerche, I. 2012, A&A, 545, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Tverskoi, B. A. 1967, J. Exp. Theor. Phys., 25, 317 [NASA ADS] [Google Scholar]
  79. Vainio, R., & Schlickeiser, R. 1999, A&A, 343, 303 [NASA ADS] [Google Scholar]
  80. Vocks, C., Salem, C., & Lin, R. P. 2005, ApJ, 627, 540 [NASA ADS] [CrossRef] [Google Scholar]
  81. Vocks, C., Mann, G., & Rausche, G. 2008, A&A, 480, 527 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Weinhorst, B., & Shalchi, A. 2010, MNRAS, 403, 287 [NASA ADS] [CrossRef] [Google Scholar]
  83. Zank, G., Li, G., Florinski, V., et al. 2006, J. Geophys. Res., 111, A06108 [NASA ADS] [CrossRef] [Google Scholar]
  84. Zhang, M., & Lee, M. A. 2011, Space Sci. Rev., DOI: 10.1007/s11214-011-9754-3 [Google Scholar]

Appendix A: Fourier solution of the momentum-diffusion 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 self-similar 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ω/D0 [−(q + 1)] −(q + 4)/(q + 1) and m = (q + 4)/(q + 1) one has (A.3)which has the solution (A.4)where Jn(z) is the Bessel function of the first kind of order n and where F0(ω) is constant in ξ. Then (A.5)We note that 2 − m = (q − 2)/(q + 1).

Suppose that F0(ω) is a simple power in ω, i.e., F0(ω) = f0(−iω)γ. Now set α = α0(−iω) so that (A.6)Set −iωξ2−m ≡ σ2 so that 2σ dσ = −idω ξ2−m and σ = ± ξ(2−m)/2e3iπ/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 power-law F0(ω) ∝ ωγ. 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 F0(ω) = ∑ fnω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 Gram-Schmidt orthonormalization process one has the solution by superposition.

Appendix B: Non-linear momentum diffusion

thumbnail Fig. B.1

Quasi-linear result for the momentum diffusion coefficient (solid line). The dotted line is the power-law fit, which yields Dp ∝ p0.4s.

With ε ≡ RA/R, the general expression for the Fokker-Planck 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 quasi-linear 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 wavenumber-dependent 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 Dp ∝ p0.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 non-linear theory (Shalchi et al. 2004) showed that, compared to previous, quasi-linear results (Shalchi & Schlickeiser 2004), the velocity-dependence of the momentum diffusion coefficients is flattened. Here, the second-order quasi-linear theory (SOQLT, see Shalchi 2005; Tautz et al. 2008; Tautz & Lerche 2010) is used to derive the non-linear 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  ⟨ z2(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 pitch-angles close to 90°, simple analytical expressions can be obtained for  ⟨ z2(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 momentum-diffusion coefficient, the result reads (B.8a)where (B.8b)is the second-order 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 Dp is shown in Fig. B.2 in comparison to results obtained from a test-particle simulation and the quasi-linear 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 Fokker-Planck coefficient to be considerably higher in this range. In the calculation, the Alfvén rigidity is set to RA = 10-3 in the calculation, hence the approximation is no longer valid in this regime.

thumbnail Fig. B.2

Comparison of the momentum diffusion coefficient Dp as obtained from numerical test-particle simulations (triangles), quasi-linear theory (dashed line), and second-order quasi-linear 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 ε = RA/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

Table 1

Parameters used for the Monte-Carlo simulation.

In the text

All Figures

thumbnail 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 RA = 10-4 (upper panel) and RA = 10-3 (lower panel) have been used. (This figure is available in color in electronic form.)
In the text
thumbnail Fig. 2

Energy spectral index for the cases of RA = 10-4 (upper panel) and RA = 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 power-law fits (dashed lines) to the simulation data (thin lines). (This figure is available in color in electronic form.)
In the text
thumbnail Fig. 3

Maximum of the distribution function as the time increases. Again the cases of RA = 10-4 (upper panel) and RA = 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
thumbnail Fig. 4

Numerically calculated running momentum diffusion coefficient for the cases of RA = 10-4 (solid line) and RA = 10-3 (dashed line). We note that for the last part of the lower curve, only 103 particles could be used, while the rest was obtained from a simulation with 104 particles. The two diamonds illustrate the quasi-linear result from Eq. (B.2)in Appendix B.
In the text
thumbnail Fig. 5

Evolution of the momentum distribution function as obtained from the self-similar 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
thumbnail 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., ⟨p0⟩ = ⟨p(t = 0)⟩.
In the text
thumbnail Fig. B.1

Quasi-linear result for the momentum diffusion coefficient (solid line). The dotted line is the power-law fit, which yields Dp ∝ p0.4s.
In the text
thumbnail Fig. B.2

Comparison of the momentum diffusion coefficient Dp as obtained from numerical test-particle simulations (triangles), quasi-linear theory (dashed line), and second-order quasi-linear theory (crosses).
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.