Diffusion and drift of cosmic rays in highly turbulent magnetic fields

We determine numerically the parallel, perpendicular and antisymmetric diffusion coefficients for charged particles propagating in highly turbulent magnetic fields, by means of extensive Monte Carlo simulations. We propose simple expressions, given in terms of a small set of fitting parameters, to account for the diffusion coefficients as functions of magnetic rigidity and turbulence level, and corresponding to different kinds of turbulence spectra. The results obtained satisfy scaling relations, which make them useful for describing the cosmic ray origin and transport in a variety of different astrophysical environments.


JCAP10(2004)007
Diffusion and drift of cosmic rays in highly turbulent magnetic fields The diffusion and drift of charged particles across highly turbulent magnetic fields are key issues in describing the transport of cosmic rays in different astrophysical environments, e.g. the interplanetary, interstellar and intergalactic media, as well as the efficiency of Fermi acceleration processes at cosmic ray sources. In particular, it has been shown that the inclusion of drift effects in the transport equation leads naturally to an explanation for the knee, for the second knee and for the observed behaviour of the composition and anisotropies between the knee and the ankle [1]- [4]. However, the accuracy of the investigations performed so far is limited by the lack of conclusive results concerning the behaviour of the diffusion tensor under highly turbulent conditions as a function of the particle energy and the relevant magnetic field parameters. In particular, the magnetic fields in the Galaxy are highly turbulent because the mean random field is of the order of the mean regular field. Moreover, since there are reversals in the orientation of the regular field, this implies the existence of regions with negligible regular fields in which the turbulence prevails.
Perturbative studies for low turbulence have been developed for a long time [5]- [7], but these analytic methods cease to be applicable for high turbulence levels, and only recently have the parallel and transverse diffusion coefficients been calculated numerically for regimes with high turbulence [8,9]. The aim of this work is to provide a thorough and more systematic calculation of these coefficients, and to parametrize the results in order to make them useful in a variety of different kinds of applications. Moreover, we present here also a numerical evaluation of the Hall diffusion coefficient that is responsible for the drift effects, which so far has never been evaluated quantitatively under highly turbulent conditions. It should also be remarked that, while in [8,9] only the Kolmogorov spectrum of fluctuations in the random magnetic field was considered, in this work other types of turbulence are studied as well (namely, the Kraichnan and Bykov-Toptygin turbulence spectra, which bracket a wide range of possible turbulence spectra).
A relativistic particle of charge Ze propagating in a uniform regular magnetic field B 0 describes a helical path characterized by a pitch angle θ and a Larmor radius given by The component of the velocity parallel to B 0 is v = c cos θ, while the radius of the helical trajectory is r L sin θ. In the presence of a random magnetic field B r with a maximum scale of turbulence L max , the particles scatter off the magnetic irregularities and change their pitch angle, but not their velocity. The pitch angle scattering proceeds mainly in resonance (i.e., the scattering is dominated by the inhomogeneities with scales of the order of r L ), and hence it is an effective mechanism of isotropization as long as r L < L max . For instance, for the galactic magnetic field, with strength B 0 a few µG and maximum scale of turbulence L max 100 pc, the pitch angle scattering leads to a diffusive regime for protons with energies up to a few 10 17 eV.
In general, the diffusion tensor D ij can be written as where b = B 0 /B 0 is a unit vector along the regular magnetic field, δ ij is the Kronecker delta symbol and ijk is the Levi-Civita fully antisymmetric tensor. The symmetric terms contain the diffusion coefficients parallel and perpendicular to the mean field, D and D ⊥ , which describe diffusion due to small scale turbulence, while the antisymmetric term contains the Hall diffusion coefficient D A . The diffusion along the magnetic field direction is due to the pitch angle scattering and leads to a diffusion coefficient given by where λ is the mean free path in the parallel direction [1]. In this expression, λ depends on the power of the random magnetic field modes at scales of order of r L , i.e.
where dE r /dk is the power spectrum of the random magnetic field energy density. In general, a component of the diffusion tensor is associated with some particular physical processes which may not contribute to the other components. The diffusion transverse to the regular magnetic field is due to both pitch angle scattering, which is the mechanism that prevails in the parallel diffusion, and to the wandering of the magnetic field lines themselves, which drag with them the diffusing particles in the direction perpendicular to B 0 [6,10,11]. When the turbulence level is small, the diffusion in the orthogonal direction is much slower than in the parallel one and is strongly affected by the field line random walk, but in the limit of very high turbulence, the parallel and perpendicular motions become similar.
The Hall coefficient D A describes in turn the macroscopic drift associated with the gradient of the CR density, leading to a current This macroscopic drift is orthogonal to both B 0 and ∇N, and in particular is also present for a constant regular field. The relation between this macroscopic current and the microscopic drift associated with gradients and curvature in B 0 has led to some confusion and controversy in the literature [12]- [15]. A common approach is to identify a guiding centre, associated with the instantaneous radius of curvature of the particle's trajectory, and to define a guiding centre velocity v g . In the so-called first-order orbit theory, an averaging over a particle's gyroperiod is performed assuming that the scale of variations in the field is much smaller than the Larmor radius, and then an ensemble averaging over a given distribution of particles is carried out. Alternatively, v g can be calculated by averaging locally over a distribution of particles in phase space [15]. In the latter case, the results are given in terms of an expansion of the anisotropy of the momenta in spherical harmonics, but in principle no assumptions concerning the scale of variations of the field are required. On the other hand, the average particle velocity v can be calculated from the Vlasov (or collisionless Boltzmann) equation, averaging again over a particle distribution in phase space. As shown in [15], the mean particle and guiding centre velocities are related through their definition, and they differ in a term which is sometimes called diamagnetic drift. A qualitative difference that can be pointed out is that v has a contribution arising from density gradients and can be non-vanishing even in the case of a uniform regular magnetic field, while v g depends only on the field gradients.
Notice that, when taking an average of the particle velocities inside a given volume element, one is including particles whose centres of gyration are outside the volume considered and, inversely, considering the guiding centres inside a given volume may correspond to taking into account particles which are actually outside that volume element. Since the observable quantity is the average particle velocity v and not the abstract concept of the guiding centre motion v g , it is clear that it is the former that will enter into the drift motions.
As commented above, the diffusion coefficients were calculated analytically only in the case of small turbulence levels, while results valid under highly turbulent conditions require instead a numerical approach [8,9]. Following [8], we will consider charged particles propagating in a magnetic field of the form B(r) = B 0ẑ + B r (r), where the first term represents a uniform regular field directed along the z-direction, while the second corresponds to the random component. In order to approximate numerically the isotropic and spatially homogeneous turbulent field, one can sum over a large number (N m ) of plane waves with wavevector direction, polarization and phase chosen randomly [16,8], i.e.
where the two orthogonal polarizationsξ α n (α = 1, 2) are in the plane perpendicular to the wavevector direction (i.e.,ξ α n ⊥ k n , so as to ensure that ∇ · B = 0). The wavenumber distribution is taken according to a constant logarithmic spacing between k min = 2π/L max and k max = 2π/L min , where L min and L max are the minimum and maximum scales of turbulence, respectively. The energy density of the turbulent component is taken as dE r /dk ∝ k −γ , where the spectral index γ is given by the kind of mechanism that builds up the turbulence. Hence, the plane wave amplitudes satisfy A 2 (k n ) = N B 2 r k −γ n (k n − k n−1 ), with N a normalization constant that ensures that n A 2 (k n ) = B 2 r . In this work we will consider in particular three spectra of interest for astrophysical applications, namely a Kolmogorov spectrum with γ = 5/3 (which is the only case studied numerically in the past, and is particularly attractive since, according to observations [17,18], the density fluctuations in the interstellar medium follow this turbulence spectrum), a Kraichnan hydromagnetic spectrum with γ = 3/2 [19] and the Bykov-Toptygin spectrum with γ = 2 [20].
By means of the Kubo formalism [21]- [23], the diffusion coefficients D ij can be computed directly by taking ensemble averages of the decorrelation between different components of the single-particle velocities, i.e.
and where · · · denotes the ensemble average taken over an isotropic distribution of many particles). In [23], it was assumed that the velocity decorrelations were modulated by exponential factors, adopting the following ansätze:

JCAP10(2004)007
Diffusion and drift of cosmic rays in highly turbulent magnetic fields where ω = c/r L is the Larmor angular gyrofrequency and where τ , τ ⊥ and τ A are the decorrelation timescales associated with the different diffusion components. Then, the diffusion coefficients were obtained by integrating the proposed expressions for R ij (t) using equation (7). The diffusion coefficients were found to be and It is interesting to note that expressions of this form, but in which only a single timescale τ appears for the three diffusion coefficients, were obtained also in other analytic approaches that assumed a single scattering process to be responsible for all the decorrelations [24]- [26]. However, it should be pointed out that the expressions proposed in [23] for R ij (t) assume implicitly a small departure from the helical trajectories, and they are no longer adequate for high turbulence levels (see below for further discussions). Moreover, there is no general theory providing the decorrelation timescales, and this requires then further assumptions in order to make equations (11)-(13) useful. These additional assumptions are sometimes obscure and lead to results often at odds with the outcome of numerical simulations, as pointed out in [2]. Alternatively, the parallel and perpendicular diffusion coefficients can be calculated from the asymptotic rate of increase of the mean squared displacements in each direction, namely and In any case, the straightforward method consists basically in generating a sample configuration for the random magnetic field (by choosing randomly the propagation direction, polarization and phase of the N m plane waves in equation (6)), and then following the trajectory of a particle that propagates in the total (regular plus random) field from the origin with a random initial direction. The results should then be averaged over a large number of different field configurations in order to calculate the corresponding diffusion coefficients. In [8], for instance, the results were obtained by integrating firstly the trajectories of 2500 particles for a given field configuration, and then by averaging over 50 different field realizations. We found, however, that the results are more biased by the particular choice of field configuration than by the choice of the initial velocity of the particle. Hence, it turns out to be more convenient to average directly over a large number of field realizations. Indeed, we typically generated

JCAP10(2004)007
Diffusion and drift of cosmic rays in highly turbulent magnetic fields following the trajectory of a single particle in each of them, and used N m = 100 modes in all simulations.
For the numerical integration of the particle trajectories under the influence of the Lorentz force, we adopted a time step ∆t = 0.1 r L /c in the Runge-Kutta routine. In [9], the parallel and perpendicular diffusion coefficients were computed from equations (14) and (15), and this method requires one to follow the particle trajectories for quite long times in order to reach the asymptotic region, typically longer than t = 1000 r L /c. We found it more convenient to instead compute directly the particle decorrelation functions R ij (t) and integrate them through equation (7), since this procedure requires one to follow the particle trajectories for times typically not larger than t = 100 r L /c. The underlying reason for this seems to be that the displacements keep more memory of the initial velocity adopted than the velocities themselves. In addition, this method allows one to compute the antisymmetric diffusion coefficient by means of R yx , while averages such as ∆x∆y vanish for large times, and hence they are not useful for computing D A . Figure 1 shows the time dependence of the decorrelation functions associated with different velocity components, for the Kolmogorov case and different turbulence levels, with the turbulence level defined here 4 by σ 2 ≡ B 2 r /B 2 0 . For low turbulence, it is observed in figure 1(a) that, after a given initial transient period, the correlation functions behave as equations (8)-(10), and hence one can calculate the decorrelation timescales directly. For instance, following the local maxima of the sinusoidal functions one finds τ ⊥ = τ A , as is apparent from figure 2. In [23], these decorrelation timescales were considered to be equal as a simplifying assumption, but it was suggested that new effects could arise in a general case with τ ⊥ = τ A . Here we find that these decorrelation timescales are indeed equal, and actually this could be expected from the fact that the physical origins for the decorrelations are the same in the two cases. For high turbulence levels, however, the amplitude decrease is much more abrupt, as is shown in figure 1(b). The decorrelation functions become vanishingly small already in the transient period, and equations (8)-(10) are hence not valid any longer.
Defining the dimensionless parameter ρ ≡ r L /L max , the results given in terms of D/(cL max ρ) versus ρ are universal and can be scaled in order to calculate the diffusion tensor for different sets of values for the regular field amplitude, random field length scale and particle energies. For instance, the range 0.01 ≤ ρ ≤ 1 investigated in this work can be regarded as corresponding to 10 15 ≤ E/eV ≤ 10 17 for protons propagating in the Galaxy (for B 0 = 1 µG and L max = 100 pc). Alternatively, considering protons propagating in the interplanetary field, with strength B 0 = 50 µG and L max = 0.01 AU, the limit to the diffusion regime at ρ = 1 corresponds to the kinetic energy E k = 1.8 GeV. 5 Since pitch angle scattering proceeds mainly in resonance, the particle propagation is essentially independent of the minimum scale length of turbulence L min , as long as L min r L . For definiteness, we then adopt L min = 0.1 r L . Figure 3 shows the numerical results obtained for the parallel diffusion coefficient corresponding to the Kraichnan spectrum and for different turbulence levels in the range 4 Alternatively, one can define the turbulence level as η ≡ B 2 r /( B 2 r + B 2 0 ) = σ 2 /(1 + σ 2 ) [9]. The definition adopted here follows that in [8]. 5 In [8], however, the power spectrum of the random interplanetary field is considered to continue with constant amplitude from 0.01 AU up to a maximum scale of 1 AU, hence extending considerably the diffusive regime to higher energies.

JCAP10(2004)007
Diffusion and drift of cosmic rays in highly turbulent magnetic fields  (8) and (9)). The curves are guides to the eye.
showing a crossover between the resonant and non-resonant scattering regimes that takes place around ρ 0.2. Indeed, the actual scale length that separates the two regimes is determined by the correlation length of the turbulence spectrum L c , defined as where the point r(L) is displaced with respect to the origin by a distance L along a fixed direction. Then, considering a spectrum of fluctuations of index γ extending between the scale lengths L min and L max , the correlation length is given by [28] Hence, L c /L max 0.2 is a quite representative value for the random field power spectra considered in this work, and this explains the change of regime observed at ρ 0.2. According to these considerations, an appropriate way of fitting the results is to interpolate D between the power laws that characterize the low and high rigidity regimes. A convenient way to achieve this is by means of the expression where the parameters N and ρ are given in table 1. This expression appears to fit our numerical results well up to σ 2 10, which is already a very high turbulence level in most   astrophysical applications of interest. For even higher turbulence levels, the results are best described asymptotically by the parameters N asint and ρ asint , also given in table 1.
In order to check the consistency with the previous numerical results for highly turbulent parallel diffusion (which was studied in [9] assuming a Kolmogorov spectrum of fluctuations), figure 4 shows a comparison between our fit and the results of [9], both corresponding to a Kolmogorov spectrum and for the same turbulence levels. Notice that in [9] the Larmor radius 6 r * L is defined by making the replacement B 0 → B 2 0 + B 2 r in equation (1), hence coupling in r * L the dependence on rigidity and on the turbulence level, while, on the other hand, the dimensionless parameter ρ * is defined as ρ * = 2πr * L /L max . Hence, the results have to be rescaled according to the relations D/(cL max ρ) = D/(r * L c)/ √ 1 + σ 2 and ρ = ρ * √ 1 + σ 2 /2π. As can be seen in figure 4, the two data sets agree reasonably well.
The perpendicular diffusion coefficient is most easily parametrized from the results for the ratio D ⊥ /D , since this ratio exhibits little dependence on the rigidity. In order to avoid the subdiffusive regime (see discussion below), we first restricted ourselves to  Figure 4. Comparison between the fit to D given in this work (formally only valid in the range 0.01 ≤ ρ ≤ 1) and the results of [9], both corresponding to a Kolmogorov spectrum and for the same turbulence levels.  (19)). For comparison, results obtained from fits to the data given in [8] and [9] are also shown.

JCAP10(2004)007
Diffusion and drift of cosmic rays in highly turbulent magnetic fields  ρ ≥ 0.03 and σ 2 ≥ 1. As in previous numerical investigations [8,9], we found the ratio D ⊥ /D to be independent of rigidity in the low rigidity region (namely, for ρ ≤ 0.2), while it scales as ρ −2 in the high rigidity region. Hence, our results can be conveniently parametrized by means of the expression where the parameters N ⊥ and a ⊥ are also given in table 1. In [2]- [4], a fit to the data of [9] (for the low rigidity regime) was given as while the results presented in [8], which correspond to the turbulence range 0.03 ≤ σ 2 ≤ 1, can be accounted for by the expression Figure 5 shows the fit given by equation (19) for the Kolmogorov case and different turbulence levels. For the sake of comparison, the results of [9] (as given by equation (20)) and of [8] (as given by equation (21)) are also shown. For moderate turbulence levels (σ 2 = 1) our results show an excellent agreement with the results of [8], while the agreement with [9] tends to be better for larger turbulence.
For low rigidities and turbulence levels (for instance, ρ 0.01 and σ 2 ≤ 1), we have found evidence for the phenomenon of subdiffusion, as already reported in [9]. Indeed, a low rigidity particle at low turbulence tends to remain attached to a field line, and hence transverse diffusion chiefly proceeds by the transverse random walking of field lines. The mean perpendicular displacement is then expected to evolve more slowly with time than in the case of normal diffusion (i.e., ∆x 2 ⊥ ∝ t m with m < 1). Figure 6 shows (∆x ⊥ ) 2 /2∆t/(cL max ρ) as a function of time, for different values of rigidity and turbulence levels, in the case of a Kolmogorov spectrum of fluctuations. While a plateau is attained at large t for ρ = 0.1 and σ 2 = 1, hence corresponding to the usual diffusion relation ∆x 2 ⊥ ∝ t, for lower values of rigidity and turbulence level the phenomenon of subdiffusion shows up. As an approach to understanding the subdiffusive regime, the so-called compound diffusion has been proposed and investigated, in which particles are assumed to be strictly tied to the field lines, while they scatter back and forth along the lines [29,22,30,31]. For the limiting case of compound diffusion, it turns out that m = 1/2, while in the case of three-dimensional particle transport, m is expected to have a smooth dependence on rigidity and turbulence, such that 1/2 ≤ m(ρ, σ 2 ) < 1. Further investigations aiming at a detailed, quantitative description of subdiffusion are currently in progress.
Concerning the antisymmetric diffusion coefficient, the numerical results can be fitted by the expression where and where the values for the parameter N A are given in table 1. In the limit of very low turbulence, this expression tends to the appropriate value D A cr L /3 [1] (see also equation (13) in the large τ A limit), while it vanishes in the limit of very strong turbulence, as expected. In [2]- [4], the propagation of cosmic rays diffusing in the Galaxy was studied and, due to the lack of an expression like that given by equations (22) and (23) (i.e. valid even under highly turbulent conditions), a simple ρ-independent prescription for D A /cL max ρ was adopted, namely which was inferred for a Kolmogorov spectrum of fluctuations. In figure 7 this ad hoc prescription is compared to the actual fit of equations (22) and (23), calculated for the Kolmogorov case and the rigidities ρ = 0.01, 0.1, 1. One can observe a sound agreement between the prescription of [2]- [4] and the fit corresponding to ρ = 0.01 and, in any case, the dependence of D A /ρ on ρ is apparent only for very high turbulence levels σ 2 1. In summary, in this work we performed extensive Monte Carlo simulations to determine the parallel, transverse and antisymmetric diffusion coefficients that describe the propagation of cosmic rays under highly turbulent conditions. We examined the simple analytical approach proposed in [23], and found that the expressions given there in terms of mean trajectory decorrelations are only meaningful for low turbulence levels. Furthermore, we found that, as long as their approach is valid, the decorrelation timescales associated with transverse and antisymmetric diffusion are indeed equal. We evaluated the diffusion Comparison between the fit to D A given in this work (for the Kolmogorov spectrum and the rigidities ρ = 0.01, 0.1, 1) and the simple ρindependent prescription adopted in [2]- [4] (see equation (24)).
coefficients and parametrized the results by means of simple expressions, which agree with the expected behaviour in the limit of low turbulence levels. Moreover, the results obtained were compared to the previous numerical calculations performed in [8,9]. In this respect, this work extends the previous investigations, since it also takes into consideration other possible turbulence spectra in addition to the Kolmogorov case (namely, the Kraichnan and the Bykov-Toptygin spectra), it provides useful parametrization formulae and it includes the study of the antisymmetric diffusion coefficient. Finally, we compared the new results for the antisymmetric coefficient with the prescription adopted previously in [2]- [4] for explaining the cosmic ray spectrum, composition and anisotropies in the region between the knee and the ankle. We hope that these results will be useful in a variety of different astrophysical scenarios related to the origin and transport of cosmic rays.