© ESO, 2011

1. Introduction

The scattering and the spatial diffusion of high energy particles off magnetic turbulence play a crucial role in many fields of astrophysics. For instance, they are key ingredients of Fermi acceleration processes because they directly control the efficiency and the rate of particle acceleration. They determine the properties of the confinement of astrophysical objects from jets to galaxies and clusters of galaxies, and governs the transport of the particles through interplanetary, interstellar, or intergalactic space. Diffusion has long been described by a quasi-linear theory approach (Jokipii 1966, 1973), which allows us to calculate the diffusion coefficients when the turbulent field is significantly weaker than the background field. However, in many circumstances the level of turbulence turns out to be large so that this standard picture requires extension. Several studies have examined the transport properties in strong turbulence by means of numerical simulations, e.g. Giacalone & Jokipii (1999), Casse et al. (2002), Candia & Roulet (2004), and Fatuzzo et al. (2010). Most of these investigations have focused on the situation in which a large-scale turbulence cascades toward small dissipative scales – as in the Kolmogorov scheme – and in which particles interact with gyroresonant modes of the turbulence spectrum. From the point of view of the particle, the turbulence therefore occurs on large scales, as the coherence length of the magnetic field corresponds roughly to the maximal scale of the turbulent spectrum.

However, in a variety of physical situations, the Larmor radius of the particle can exceed the coherence scale of the turbulence. The transport of particles downstream of a relativistic shock wave provides a clear example of this situation. The mean field is there mostly transverse to the flow because of Lorentz tranform effects and shock compression, and the turbulence that is excited in the shock precursor is generated on microscopic plasma skin depth scales. In this case, perpendicular diffusion at high (possibly very high) rigidity controls the transport of the particles back and forth from the shock. More generally, the high rigidity regime likely plays an important role in the deconfinement process of particles of high energy, when their Larmor radius exceeds the size of the astrophysical system. However, this high rigidity regime has received little attention so far, except for the pioneering study of Shalchi & Dosch (2009). The pitch angle scattering rate is known to increase in proportion to the square of the particle energy in this limit, but the behavior of the transverse diffusion coefficient, which is crucial in the above contexts deserves a careful analysis. This analysis is the objective of the present paper. It will be found in particular that even a weak mean field, as measured relatively to the turbulent component, can affect the scaling of the perpendicular diffusion coefficient.

The present paper describes both a theoretical and a numerical study of diffusion at high rigidity. The theoretical aspects are discussed in Sect. 2, while the numerical simulations are presented in Sect. 3. Finally in Sect. 4 we summarize our results and discuss some applications.

2. Transport of high rigidity particles with a mean field

2.1. Notations and summary of previous results

The transport of particles in magnetostatic turbulence is characterized by the reduced rigidity ρ, the level of turbulence η, and the power spectrum of magnetic fluctuations in three dimensions (3D) S_3d(k). These quantities are defined as $\mathit{\rho }\mathrm{\equiv }\frac{\mathit{r̅}\mathrm{L}}{{\mathit{\ell }}_{\mathrm{c}}}\mathrm{=}\frac{\mathit{ϵ}}{\mathit{e}\mathit{B̅}{\mathit{\ell }}_{\mathrm{c}}}\mathit{,}$(1)where denotes the Larmor radius of the particle in the total (mean B₀ and turbulent δB) field where , ϵ the energy of the particle, and ℓ_c the coherence length of the fluctuations.

The turbulence level η is defined as $\mathit{\eta }\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\frac{\mathrm{⟨}{{\delta B}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}{{\delta B}}^{\mathrm{2}}\mathrm{⟩}\mathrm{+}{{{B}}_{\mathbf{0}}}^{\mathrm{2}}}\mathit{,}$(2)where η → 0 corresponds to weak turbulence and η → 1 corresponds to pure turbulence with no mean field.

The correlation function C(r) of the random field $\mathit{C}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\frac{\mathrm{⟨}{\delta B}\mathrm{(}{x}\mathrm{+}{r}\mathrm{)}{\delta B}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{⟩}}{\mathrm{⟨}{{\delta B}}^{\mathrm{2}}\mathrm{⟩}}\mathit{,}$(3)can be written in terms of the one-dimensional power spectrum S(k) ∝ k²S_3d(k) $\mathit{C}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{\int }\mathrm{d}\mathit{k}\hspace{0.17em}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\sin \mathrm{\left(}\mathit{kr}\mathrm{\right)}\mathit{/}\mathrm{\left(}\mathit{kr}\mathrm{\right)}}{\mathrm{\int }\mathrm{d}\mathit{k}\hspace{0.17em}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{·}$(4)Casse et al. (2002) defined the coherence length as the scale at which C(r) is maximum; if the power spectrum takes the form of a broad-band truncated power-law S(k) ∝ (k/k_min) ^− β for k_min ≤ k ≤ k_max and zero otherwise, one finds for the coherence length ${\mathit{\ell }}_{\mathrm{c}}\mathrm{\simeq }\mathrm{0.77}{\mathit{k}}_{\min }^{-1}$. Alternatively, one can define the coherence length as we do here, to be ${\mathit{\ell }}_{\mathrm{c}}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{r}\hspace{0.17em}\mathit{C}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(5)where one then derives in a straightforward way ${\mathit{\ell }}_{\mathrm{c}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\pi }}{\mathrm{2}}\frac{\mathrm{1}}{\mathit{\eta }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{k}\hspace{0.17em}{\mathit{k}}^{-1}\hspace{0.17em}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(6)and the presence of 1/η results from our choice of normalization for the power spectrum ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{k}\hspace{0.17em}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\mathit{\eta ,}$(7)where in practice the spectrum is bounded between k_min and k_max. Both definitions for ℓ_c coincide to within a factor close to unity. As a function of the spectrum index β, the coherence length is close to either ${\mathit{k}}_{\min }^{-1}$ on larges scale for β > 1, or to ${\mathit{k}}_{\max }^{-1}$ on small scales for β < 1.

The scattering frequency ν_s is defined as the reciprocal of the decorrelation time of the pitch angle of the particle, the latter being defined relative to the direction of the mean field. As discussed in Casse et al. (2002), the scattering frequency can be written ${\mathit{\nu }}_{\mathrm{s}}\hspace{0.17em}\mathrm{\approx }\hspace{0.17em}\frac{\mathit{\pi }}{\mathrm{3}}\frac{\mathit{c}}{\mathit{r̅}\begin{array}{}\mathrm{2}\\ \mathrm{L}\end{array}}\frac{{\mathrm{\int }}_{\mathit{k}\mathit{r̅}\mathrm{L}\mathit{>}\mathrm{1}}\hspace{0.17em}{\mathit{k}}^{-1}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{k}}{\mathrm{\int }\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{k}}\mathit{,}$(8)an expression that extends to the strong turbulence regime the results of the quasi-linear theory. This leads to the scalings $\begin{array}{ccc}{\mathit{\nu }}_{\mathit{s}}& \mathrm{\simeq }& \frac{\mathrm{2}}{\mathrm{3}}\mathit{\eta }\frac{\mathit{c}}{{\mathit{\ell }}_{\mathrm{c}}}{\mathit{\rho }}^{\mathit{\beta }\mathrm{-}\mathrm{2}} \mathrm{(}\mathit{\rho }\mathrm{\ll }\mathrm{1}\mathrm{)}\\ {\mathit{\nu }}_{\mathit{s}}& \mathrm{\simeq }& \frac{\mathrm{2}}{\mathrm{3}}\mathit{\eta }\frac{\mathit{c}}{{\mathit{\ell }}_{\mathrm{c}}{\mathit{\rho }}^{\mathrm{2}}} \mathrm{(}\mathit{\rho }\mathrm{\gg }\mathrm{1}\mathrm{)}\mathit{.}\end{array}$(9)The Bohm scaling holds only in the very special case where β = 1. In addition to these quantities, the notion of correlation time also plays an important role because it measures the time beyond which a particle experiences a force that is decorrelated from the initial one, along the particle trajectory. It is then defined as ${\mathit{\tau }}_{\mathrm{c}}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{C}\mathrm{\left(}\mathrm{\right|}\mathrm{\Delta }\mathit{x}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{\left|}\mathrm{\right)}\mathrm{d}\mathit{\tau ,}$(10)where Δx(τ) represents the displacement after a time τ in the turbulence. In quasi-linear theory, only the unperturbed trajectory is inserted into this definition, although one can extend that definition with a diffusive trajectory as we later indicate.

If a relativistic particle travels over a coherence length of the turbulent field without having displayed any wiggle, corresponding to the regime ρ ≫ 1, then τ_c ~ ℓ_c/c. This correlation time is much shorter than the scattering time ${\mathit{\nu }}_{\mathrm{s}}^{-1}\mathrm{~}{\mathit{\eta }}^{-1}{\mathit{\rho }}^{\mathrm{2}}{\mathit{\ell }}_{\mathrm{c}}\mathit{/}\mathit{c}$ in this regime. The correlation time τ_c can be recovered from Eq. (10) by using the ballistic approximation Δx(τ) ≃ cτ, which is appropriate in this regime ρ ≫ 1, in which case Eq. (4) leads to ${\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathit{\beta }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{\beta }}^{-1}{\mathit{k}}_{\min }^{-1}\mathit{/}\mathit{c}\mathrm{~}{\mathit{\ell }}_{\mathrm{c}}\mathit{/}\mathit{c}$. We note that in the special case where the power-law index of turbulence β = 1 (Bohm regime) Eqs. (4) and (5) lead to ℓ_c = (λ_min/4)log (λ_max/λ_min), where λ_min and λ_max are the shortest and the longest wavelengths of turbulence.

If a particle experiences a chaotic motion on a length-scale smaller than ℓ_c, corresponding to the regime ρ ≪ 1, then the estimate is more complicated to obtain but one finds that τ_c ~ ρ^βℓ_c/c as follows. Since the correlation time remains shorter than the scattering time, Casse et al. (2002) proposed a heuristic estimate in which decorrelation arises out of the small-scale modes with wavenumber k > k_minρ^-1, which give rise to gyroresonant interactions with the particle of rigidity ρ. The modes with wavelengths longer than the Larmor radius (i.e., k < k_minρ^-1) construct the field line to which the particle is attached, hence do not cause decorrelation on timescales shorter than the scattering time. The above correlation time is indeed shorter than the scattering time and increases with ρ. The heuristic estimate for ρ < 1 is consistent with quasi-linear theory when η ≪ 1 and with numerical results in the strong turbulence regime (Casse et al. 2002) can then be written as ${\mathit{\tau }}_{\mathrm{c}}\mathrm{\simeq }\frac{\mathrm{1}}{\mathit{\eta c}}{\mathrm{\int }}_{\mathit{k}\mathit{>}{\mathit{k}}_{\min }{\mathit{\rho }}^{-1}}\mathrm{d}\mathit{k}\hspace{0.17em}{\mathit{k}}^{-1}\mathit{S}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(11)which bears some resemblance to the case discussed before for ρ ≫ 1, except that ρ explicitly enters the sinc function, since one must now follow the orbit of the particle around the field line, and the integral is limited to k > k_minρ^-1 for the reasons given above. The calculation then implies that τ_c ~ ρ^βℓ_c/c as announced. The particle trajectory undergoes decoherence before traveling ℓ_c because of the large number of wiggles in the random field.

Thus, except for η ~ 1 and ρ ~ 1 for which the correlation time becomes comparable to the scattering time, a Markovian theory of the scattering process is appropriate, even if the turbulence is strong, stronger even than the mean field. This is an essential key for the present discussion.

Independently of the rigidity, the parallel diffusion coefficient is always given by D_∥ = c²/(3ν_s), even in the strong regime of turbulence. As for the transverse diffusion coefficient, in the strong regime at low rigidities, it does not follow a law similar to the quasi-linear result but is proportional to D_∥ (Casse et al. 2002) because of the magnetic field line wandering that transmits parallel diffusion in the transverse direction. Casse et al. (2002) found in particular that D_⊥ = η^2.3D_∥ at small rigidities, which rules out the conjecture of Bohm’s diffusion. In the next section, we discuss the transverse diffusion in the large rigidity regime.

2.2. Transverse diffusion at large rigidity

As mentioned previously, in the large rigidity regime ρ ≫ 1, the correlation time is (much) shorter than the scattering time, hence we expect to derive the parallel and transverse diffusion coefficients using a Markovian description of the trajectory. In particular when ρ ≫ 1, the velocity changes by 1/ρ only over a correlation time. This implies that significant changes in the velocity occur on timescales that are much longer than the correlation time. Therefore we can assimilate the effect of small-scale fluctuations to a fully decorrelated white noise on the relevant timescales.

To calculate the particle transport in a random field, one has to use the solution of the differential equation that governs the evolution of the particle velocity v$\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{v}\mathrm{=}\left[\mathbf{\Omega ̂}\mathbf{0}\mathrm{+}{\delta }\mathbf{\Omega ̂}\mathbf{\left(}{t}\mathbf{\right)}\right]\mathrm{·}{v}\mathit{.}$(12)The quantities $\mathbf{\Omega ̂}{\mathbf{0}}_{}$ and are rotation operators developed as linear combinations of the generators of the Lie algebra of the rotation group, ${L̂}{\mathbf{1}}_{}$, ${L̂}{\mathbf{2}}_{}$, ${L̂}{\mathbf{3}}_{}$$\mathit{L̂}\mathrm{1}\hspace{0.17em}\mathrm{:}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -1\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right)\mathit{,}\mathit{L̂}\mathrm{2}\hspace{0.17em}\mathrm{:}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -1& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,}\mathit{L̂}\mathrm{3}\hspace{0.17em}\mathrm{:}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{0}& -1& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{.}$(13)In detail, $\mathbf{\Omega ̂}{}_{\mathbf{0}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{B}}_{\mathrm{0}}^{\mathit{i}}{L̂}{}_{{i}}\mathit{/}{\mathit{B}}_{\mathrm{0}}$, where ${\mathit{B}}_{\mathrm{0}}^{\mathit{i}}$ denotes the ith component of B₀ and Ω₀ ≡ c/r_L,0 the Larmor pulsation defined with respect to the mean field. With this notation, . The operator is decomposed in a similar way as the generators of the rotation group, and δΩ ≡ c/r_L, where r_L is now measured relatively to δB.

To solve the equation of motion, one uses an auxiliary variable u that is defined as ${v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}{R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{u}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(14)where ${R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\exp \mathrm{(}\mathit{t}{\mathbf{\Omega ̂}\mathbf{0}}^{\mathrm{)}}\mathit{.}$(15)We then define $\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\equiv }{R̂}\mathbf{0}\mathrm{(}\mathit{t}{\mathrm{)}}^{-1}\mathrm{·}{\delta }\mathbf{\Omega ̂}\mathrm{(}\mathit{t}\mathrm{)}\mathrm{·}{R̂}\mathbf{0}\mathrm{(}\mathit{t}\mathrm{)}\mathit{,}$(16)one finds that u(t) obeys $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{u}\mathrm{=}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{u}\mathit{.}$(17)This equation is solved as ${u}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{𝒯}\exp \left[{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\right]\mathrm{·}{u}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(18)Because the operator in the exponent is time dependent, to preserve the exponential character of the solution, a time-ordering operator has to be introduced, as we now explain.

We note that u(0) = v(0), thus the solution for v is given by ${v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}\mathrm{𝒯}\exp \left[{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\right]\mathrm{·}{v}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(19)The regular part of the field generates the regular rotation matrix , while the exponential accounts for the effect of the turbulent part. The time-ordering operator maintains the chronological order of the products in the non-commuting $\widehat{\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}$ in the expansion of the exponential operator, i.e. $\begin{array}{ccc}\mathrm{𝒯}\mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\mathrm{if}\hspace{0.17em}\hspace{0.17em}{\mathit{t}}_{\mathrm{1}}\mathit{>}{\mathit{t}}_{\mathrm{2}}\mathit{,}\\ & \mathrm{=}& \mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{·}\mathbf{\Omega ̂}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\mathrm{if}\hspace{0.17em}\hspace{0.17em}{\mathit{t}}_{\mathrm{2}}\mathit{>}{\mathit{t}}_{\mathrm{1}}\mathit{,}\end{array}$(20)and so on for higher order products. Alternatively, the time-ordered expansion can be written as a Dyson series $\begin{array}{ccc}\mathrm{𝒯}\exp \left[{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\right]& \mathrm{\equiv }& \mathrm{1}\\ & \mathrm{+}& \sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{n}\mathrm{=}\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{1}}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathit{n}\mathrm{-}\mathrm{1}}}\mathrm{d}{\mathit{t}}_{\mathit{n}}\hspace{0.17em}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{...}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\mathit{.}\\ & & \end{array}$(21)We now use the following theorem that holds for a Gaussian stationary random process in the white noise limit. As discussed in detail in the Appendix, this is a direct generalization to any Lie algebra of a well-known result for a scalar random process, with no other restriction than the white noise assumption $\begin{array}{ccc}& & ⟨\mathrm{𝒯}\exp \left[{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\right]⟩\mathrm{=}\\ & & \end{array}$(22)Various properties of the turbulent field can be considered i.e. that is either isotropic with no helicity, isotropic with helicity, or anisotropic with rotation invariance in the transverse direction. All these cases can be easily treated, although we focus on two relevant cases: (A) 3D isotropic turbulence and (B) two-dimensional (2D) isotropic turbulence in the plane transverse to B₀, with δB·B₀ = 0.

We define the projection operators ${\pi ̂}{\mathbf{\perp }}_{}$ on the plane transverse to B₀ and the projection operator ${\pi ̂}{\mathbf{\parallel }}_{}$ along B₀. We now define the correlation function of the random rotation matrices $\widehat{{\delta }\mathbf{\Omega }}$$\mathrm{⟨}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{⟩}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathit{i}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{\delta }{\mathrm{\Omega }}^{\mathit{j}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}{L̂}{i}{L̂}{j}\mathit{,}$(23)where $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathit{i}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathit{\delta }{\mathrm{\Omega }}^{\mathit{j}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{2}{\mathit{\tau }}_{\mathrm{c}}\mathit{\delta }\mathrm{(}{\mathit{t}}_{\mathrm{1}}\mathrm{-}{\mathit{t}}_{\mathrm{2}}\mathrm{)}\hspace{0.17em}\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}{{\pi ̂}\mathbf{\perp }}^{\mathit{ij}}\mathrm{+}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}{{\pi ̂}\mathbf{\parallel }}^{\mathit{ij}}\right]\mathit{.}$(24)The scalars $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}$ and $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}$ characterize the relative strengths of the turbulence in the parallel (to B₀) and perpendicular directions. In particular, for 3D isotropic turbulence, $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{2}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}$, in which case the above correlator becomes proportional to the identity. Then, using the properties of ${\pi ̂}{\mathbf{\perp }}_{}$, ${\pi ̂}{\mathbf{\parallel }}_{}$ and the ${L̂}{{i}}_{}$, one finds $\begin{array}{ccc}\mathrm{⟨}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{⟩}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \mathrm{-}\mathrm{2}{\mathit{\tau }}_{\mathrm{c}}\mathit{\delta }\mathrm{(}{\mathit{t}}_{\mathrm{1}}\mathrm{-}{\mathit{t}}_{\mathrm{2}}\mathrm{)}\\ & & \hspace{0.17em}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\left[\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathbf{1̂}\mathrm{-}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}{\mathrm{⟩}}^{\mathrm{2}}{\pi ̂}\mathbf{\parallel }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}{\pi ̂}\mathbf{\perp }\right]\mathit{,}\end{array}$(25)where $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\equiv }\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}{\mathrm{⟩}}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}$. We note that the above correlation function holds for $\widehat{{\delta }\mathbf{\Omega }}$, which should not be confused with $\widehat{\begin{array}{}{\mathbf{˜}}_{}\\ \mathbf{\Omega }\end{array}}$, the latter being the quantity of relevance for calculating the transport properties, as expressed in Eq. (22). However, $\mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{e}}^{\mathrm{-}{\mathit{t}}_{\mathrm{1}}\mathbf{\Omega ̂}\mathbf{0}}\mathrm{⟨}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}{\mathit{e}}^{\mathrm{-}\mathrm{(}{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}\mathrm{)}\mathbf{\Omega ̂}\mathbf{0}}\widehat{{\delta }\mathbf{\Omega }}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}{\mathit{e}}^{{\mathit{t}}_{\mathrm{2}}\mathbf{\Omega ̂}\mathbf{0}}\mathit{,}$(26)and, because , the correlation function for $\widehat{\begin{array}{}{\mathbf{˜}}_{}\\ \mathbf{\Omega }\end{array}}$ is the same as that for $\widehat{{\delta }\mathbf{\Omega }}$.

Using Eq. (22), one then finds the solution for v: $\mathrm{⟨}{v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}\exp \left[\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{2}}\hspace{0.17em}\mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\right]{v}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{,}$(27)where the average is taken over the possible realizations of the turbulent field. This leads to $\mathrm{⟨}{v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}\exp \left[\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}\left(\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}{\pi ̂}\mathbf{\parallel }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}{\pi ̂}\mathbf{\perp }\right)\right]{v}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(28)Using the properties of ${\pi ̂}{\mathbf{\parallel }}_{}$ and ${\pi ̂}{\mathbf{\perp }}_{}$, this can be rewritten as $\begin{array}{ccc}\mathrm{⟨}{v}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& {R̂}\mathbf{0}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}\{\exp \mathrm{[}\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}{\mathrm{(}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{)}}}^{\mathrm{]}}{\pi ̂}\mathbf{\parallel }\\ & & \mathrm{+}\exp \left[\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}\left(\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\right)\right]{\pi ̂}\mathbf{\perp }\}{v}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}\end{array}$(29)Therefore, one derives the general results $\mathrm{⟨}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\exp \mathrm{[}\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}{\mathrm{(}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{)}}}^{\mathrm{]}}\mathrm{⟨}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{(}\mathrm{0}{\mathrm{)}}^{\mathrm{2}}\mathrm{⟩}\mathit{.}$(30)In the transverse direction, $\begin{array}{ccc}\mathrm{⟨}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{·}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& \exp \left[\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}\left(\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\right)\right]\hspace{0.17em}\\ & & \mathrm{\times }{\hspace{0.17em}}^{\mathrm{T}}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{·}{R̂}\mathbf{0}\mathrm{·}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\\ & \mathrm{=}& \exp \left[\mathrm{-}\mathit{t}{\mathit{\tau }}_{\mathrm{c}}\left(\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\right)\right]\\ & & \end{array}$(31)The last equality follows from developing the exponential , noting that for B₀ oriented along z, , , , and .

The parallel D_∥ and perpendicular D _⊥  diffusion coefficients are directly obtained from the correlation functions of the velocity components after averaging over the initial velocities $\begin{array}{ccc}{\mathit{D}}_{\mathrm{\parallel }}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \int \begin{array}{}\mathrm{+}\mathrm{\infty }\\ \mathrm{0}\end{array}\mathrm{d}\mathit{t}\hspace{0.17em}\mathrm{⟨}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{v}}_{\mathrm{\parallel }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\mathit{,}\\ {\mathit{D}}_{\mathrm{\perp }}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{d}\mathit{t}\hspace{0.17em}\mathrm{⟨}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{·}{{v}}_{\mathbf{\perp }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\mathit{.}\end{array}$(32)Using Eqs. (30) and (31), this leads to $\begin{array}{ccc}{\mathit{D}}_{\mathrm{\parallel }}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \frac{\mathrm{1}}{\mathrm{3}}\frac{{\mathit{c}}^{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{s}}}\mathit{,}\\ {\mathit{D}}_{\mathrm{\perp }}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& \end{array}$(33)where $\begin{array}{ccc}{\mathit{\nu }}_{\mathrm{s}}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& {\mathit{\tau }}_{\mathrm{c}}\mathrm{(}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{)}}\mathit{,}\\ {\mathit{\nu }}_{\mathrm{\perp }}& \hspace{0.17em}\mathrm{=}\hspace{0.17em}& {\mathit{\tau }}_{\mathrm{c}}\left(\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\right)\mathit{.}\end{array}$(34)These expressions for D _⊥  are formally similar to the results of the so-called classical diffusion theory, although they are obtained here under different physical assumptions; in particular, a strong turbulence situation is assumed.

In case (A), for 3D isotropic turbulence, $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{2}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}$, so that ${\mathit{\nu }}_{\mathrm{\perp }}\mathrm{=}{\mathit{\nu }}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}{\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathit{\eta }}{{\mathit{\rho }}^{\mathrm{2}}}\frac{\mathit{c}}{{\mathit{\ell }}_{\mathrm{c}}}\mathrm{·}$(35)One may note that the expression for ν_s matches that derived from a random walk argument for pitch angle diffusion. We also note that the above calculation for ν_s may be applied to the regime ρ ≪ 1, as long as the correlation time is shorter than the scattering time. This is true in the case of ν_s = (2/3)ηρ^β − 2, which is the standard quasi-linear theory result. The result for the perpendicular coefficient cannot, of course, be extended to the regime ρ ≪ 1, as the above calculation does not account for field line wandering.

In case (B), for 2D transverse isotropic turbulence, $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{0}$, $\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}$, hence ${\mathit{\nu }}_{\mathrm{s}}\mathrm{=}\mathrm{2}{\mathit{\nu }}_{\mathrm{\perp }}\mathrm{=}\mathrm{⟨}\mathit{\delta }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{⟩}{\mathit{\tau }}_{\mathrm{c}}\mathit{.}$(36)This demonstrates that the transverse diffusion coefficients follows the scalings, which we express here for case (A), i.e., isotropic turbulence $\begin{array}{ccc}{\mathit{D}}_{\mathrm{\perp }}& \mathrm{\simeq }& {\mathit{D}}_{\mathrm{\parallel }}\hspace{0.17em}\mathrm{\simeq }\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}}\mathit{c}{\mathit{\ell }}_{\mathrm{c}}{\mathit{\rho }}^{\mathrm{2}}\mathit{/}\mathit{\eta } \mathrm{(}\mathrm{1}\mathrm{\ll }\mathit{\rho }\mathrm{\ll }\mathit{B̅}\mathit{/}{\mathit{B}}_{\mathrm{0}}\mathrm{)}\mathit{,}\\ {\mathit{D}}_{\mathrm{\perp }}& \mathrm{\simeq }& \end{array}$(37)The transition between these two regimes takes place at , corresponding to ν_s ~ Ω₀. At larger rigidities, the perpendicular diffusion coefficient remains constant, while the parallel diffusion coefficient continues to increase as ρ².

This result is supported by the numerical simulation that we now discuss.

3. Numerical simulation of the transport with a mean field for high rigidities

3.1. Numerical set up

A Monte Carlo strategy is adopted to measure the diffusion coefficients by integrating a large number of particle trajectories in given turbulent magnetic field configurations. Averages are then performed and statistical values of the diffusion coefficients deduced. The numerical set up is presented hereafter: we first discuss the construction of the magnetic field, then the integration of particle motion from Lorentz-Newton equation, and finally the estimates of the diffusion coefficients.

The total magnetic field is expressed as B = B₀ + δB as before. The regular field is oriented along z, and δB is assumed to be isotropic in the three dimensions. An algorithm similar to Giacalone & Jokipii (1999) is used to construct the turbulent component of magnetic field δB by summing over plane wave modes (N_mod) with turbulent wavelengths extending from L_min = 1 ≡ 2π/k_max to L_max ≡ 2π/k_min, the power spectrum following a truncated power law between k_minand k_max. In detail

$\begin{array}{ccc}{\delta B}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}}{\mathit{G}}_{\mathit{n}}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}{{\xi }}_{{n}}\cos \left({{k}}_{{n}}\mathrm{·}{x}\mathrm{+}{\mathit{\beta }}_{\mathit{n}}\right)\mathit{.}& & \end{array}$(38)With Fourier modes of amplitude G_n, and wave vectors ${{k}}_{{n}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{L}}_{\mathit{n}}}{{e}}_{{k}}$ isotropically distributed, the unitary vector ξ_n is perpendicular to k_n in order to ensure that ∇·δB=0, and β_n ∈  [0,2π] represents the random phase. The power spectrum is normalized by the turbulence parameter η introduced earlier such that $\mathrm{⟨}{{\delta B}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{\eta }\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}$. For definiteness, the mode amplitudes are constructed according to a Kolmogorov cascade with logarithmic spacing between wavenumbers: ${\mathit{G}}_{\mathit{n}}\mathrm{\propto }{\mathit{k}}_{\mathit{n}}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}}$. We note that the details of the inertial range of the turbulence are not important because we are interested in the scattering properties at large rigidities, when the particle Larmor radius is larger than all turbulent length-scales. For a detailed presentation of the numerical turbulent magnetic field construction, the reader is referred to Sect. 2.B of Casse et al. (2002) and Sect. 3 of Giacalone & Jokipii (1999).

Several tests of the dynamic range of turbulence L_max/L_min and the magnetic wave-modes N_mod were performed. The main difficulty is that the scattering timescales increase as a square of particle rigidity. For large rigidities, it is thus difficult to preserve the accuracy with time when achieving particle diffusion together with a realistic magnetic field model. To develop a simulation that operates over a few scattering times, one needs to achieve an integration time of at least , as in our simulations. One must also strike a compromise with the number of plane wave modes to save computational time; values of order 200–300 have emerged as a satisfactory compromise between accuracy and calculation time. To save computational time, and because the small scales of the turbulent cascade are of little influence, the dynamic range has been shortened to L_min/L_max = 0.1. Tests performed with a larger dynamic range have provided similar results; the highest accuracy is obtained when modes are concentrated on the largest scale. It is explained physically by the high energy particles interacting only with the largest magnetic structures.

Particle motion is solved using the Lorentz-Newton equation of motion that preserves its energy, hence its Lorentz factor γ$\begin{array}{ccc}\frac{\mathrm{d}{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{q}}{\mathit{\gamma mc}}{v}\mathrm{\times }\left({{B}}_{\mathbf{0}}\mathrm{+}{\delta B}\right)\mathit{.}& & \end{array}$(39)At this point, we define the numerical rigidity , which differs from the previous physical definition by a numerical factor of order unity, as discussed earlier. The exact relation between ℓ_c and L_max/2π depends on the dynamic range and the power-law index of turbulence. In the following, the conversion factor between both rigidities is derived using ℓ_c ≃ 0.1L_max, a good approximation for a Kolmogorov-type spectrum.

The numerical integration of Eq. (39) is performed using a Bulirsch-Stoer schema (Press et al. 1986). Once a large number of particle trajectories were calculated and stored, statistical averages instead being performed. Given the number of particles N_p for each field realization and the number of field realizations N_field, the diffusion tensor coefficient (i,j) is evaluated as $\begin{array}{ccc}{\mathit{D}}_{\mathit{ij}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}{\mathit{N}}_{\mathrm{p}}{\mathit{N}}_{\mathrm{field}}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{field}}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{p}}}\frac{\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{x}}_{\mathit{i}}\mathrm{(}{\mathit{t}}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right)}\mathrm{\left(}{\mathit{x}}_{\mathit{j}}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{x}}_{\mathit{j}}\mathrm{(}{\mathit{t}}_{\mathrm{0}}\mathrm{)}{\mathrm{)}}_{\mathit{n,k}}}{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}\\ & \mathrm{=}& \frac{\mathrm{⟨}\mathrm{\Delta }{\mathit{x}}_{\mathit{i}}\mathrm{\Delta }{\mathit{x}}_{\mathit{j}}\mathrm{⟩}}{\mathrm{2}\mathrm{\Delta }\mathit{t}}\mathrm{·}\end{array}$(40)The average is performed over different particle trajectories and different field realizations. For each value of ρ, we take N_field × N_p = 10³ different trajectories with random initial velocity directions. The asymptotic value for t → ∞ (plateau) is roughly constant and defines the actual diffusion regime. It gives the diffusion coefficient as D_ij as t → ∞, precisely when $\mathit{t}\mathrm{\gg }{\mathit{\nu }}_{\mathrm{s}}^{-1}$. This method of coefficient estimation appears precise enough for an integration involving 10³ particles. A complementary technique consists of evaluating time correlations between velocities over particle trajectories. With 1000 particles in the transport regime studied here, this method is affected by numerical noise for the velocity correlation function, hence is not presented.

Two different cases were investigated numerically: a pure turbulence situation (B₀ = 0) and a strong turbulent case with δB ≫ B₀. Results are presented in the following sub-sections.

3.2. Pure turbulence B₀ = 0

These simulations were performed to test the correctness and accuracy of the code. On theoretical grounds (see the appendix of Casse et al. 2002; Aloisio et al. 2004; Pelletier et al. 2009) and previous numerical works (Parizot 2004), we expect the diffusion coefficient to evolve as the square of energy (e.g. rigidity) when ρ′ ≫ 1.

Here we set η = 1 and δB isotropically distributed by construction, so that the three space directions are equivalent. The equivalence of the three directions was numerically verified in our simulations. The diffusion coefficient is evaluated as ${\mathit{D}}_{\mathrm{iso}}\mathrm{=}\frac{\mathrm{⟨}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\mathrm{\Delta }{\mathit{y}}^{\mathrm{2}}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\mathrm{\Delta }{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{6}\mathrm{\Delta }\mathit{t}}\mathrm{·}$(41)Figure 1 shows numerical values calculated for ρ′ going from 1 to 100. The diffusion coefficient is plotted in units of cL_max/(2π) as a function of rigidity ρ′. A power law is observed for 1 < ρ′ < 100, as predicted by the theory, namely D_iso ∝ ρ^′2 ∝ ϵ². One may be able to discern a slight deviation at ρ′ close to 100. This purely numerical effect disappears when taking a larger number of magnetic wave-modes on scales close to L_max by defining L_min/L_max ~ 1. We retain however the current field configuration, taking this effect into account when interpreting the results.

Fig. 1 The diffusion coefficient variation is plotted in units of cLmax/(2π) as a function of rigidity ρ′ in pure turbulence (B0 = 0 or η = 1). The dashed line is drawn as a reference for a scaling Diso ∝ ρ′2. For ρ′ > 1, Diso is indeed proportional to ρ2.

3.3. Weak mean field B₀ < δB

We now consider the case where a constant weak mean field B₀ along the z direction is present. In this case, two different diffusion coefficients are defined D_∥ = D_zz and D _⊥  = (D_xx + D_yy)/2. Overall, we explored five different levels of turbulence η =  { 0.5,0.9,0.99,0.999,0.9999 } , spanning five orders of magnitude in $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$. The rigidity ρ′ ranges from 1 to 100 for each value of η. At each calculation point  { η,ρ } , the coefficients are evaluated by averaging over 10³ particles (10 particles  ×  100 field realizations), as before.

As shown in Fig. 2, the parallel diffusion coefficient retains the same dependence on rigidity as in pure turbulence, D_∥ ∝ ρ^′2. For η > 0.5, the turbulence level has almost no influence on the value of ${\mathit{D}}_{\mathrm{\parallel }}\mathrm{\simeq }\mathrm{0.9}\mathit{c}\mathit{r̅}{}_{\mathrm{L}}\mathit{\rho }^{\mathrm{\prime }}$. Therefore as expected, the mean field, as long as it remains weak enough, has no influence on the diffusion of particles along its direction.

Fig. 2 The parallel diffusion coefficient D∥ plotted in units of (cLmax/(2π)) as a function of ρ′ for different degrees of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$. Here D∥ ∝ ρ2, as in the case of purely isotropic turbulence without mean field. As long as $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\gg }\mathrm{1}$, the strength of the turbulence does not influence the normalization of D∥.

The picture is different for the transverse coefficient when the particle rigidity becomes large. In Fig. 3, the simulated transverse diffusion coefficient is plotted as a function of rigidity ρ′ for different degrees of turbulence. In each case, its value saturates to a constant value when ρ′ ~ δB/B₀. This value behaves proportionally to the turbulence degree; in detail, ${\mathit{D}}_{\mathrm{\perp }}\mathrm{\simeq }\mathrm{0.13}\mathit{c}\mathrm{(}{\mathit{L}}_{\max }\mathit{/}\mathrm{2}\mathit{\pi }\mathrm{)}\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$, in excellent agreement with our theoretical prediction from Eq. (37). Individual particle trajectories reveal a weakly perturbed helical path when ρ′ ≫ 1. Therefore, a strong small-scale turbulence acts as a collection of small-scale scattering centers, each producing a small deflection.

Fig. 3 The transverse diffusion coefficient D ⊥  plotted in units of cLmax/(2π) as a function of ρ′ for different degrees of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$. The diffusion coefficient saturates at ρ′ ~ δB/B0. Below this value, its behavior is similar to the parallel diffusion coefficient. Beyond ρ′, its value becomes independent of particle rigidity.

According to the theory, D_∥ is the limit of a function c²g_∥(t)/3 as t → ∞, precisely as t > t_s, the function being ${\mathit{g}}_{\mathrm{\parallel }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{-}{\mathit{\nu }}_{\mathrm{s}}\mathit{t}}}{{\mathit{\nu }}_{\mathrm{s}}}\mathrm{·}$(42)In a similar way D _⊥  is the limit of a function c²g _⊥ (t)/3 as t → ∞, in addition to when t > t_s, the function being ${\mathit{g}}_{\mathrm{\perp }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\nu }}_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{\nu }}_{\mathrm{s}}^{\mathrm{2}}}\left\{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{-}{\mathit{\nu }}_{\mathrm{s}}\mathit{t}}\left[\cos \left({\mathrm{\Omega }}_{\mathrm{0}}\mathit{t}\right)\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{0}}}{{\mathit{\nu }}_{\mathrm{s}}}\sin \left({\mathrm{\Omega }}_{\mathrm{0}}\mathit{t}\right)\right]\right\}\mathrm{·}$(43)

Fig. 4 Transition toward parallel and perpendicular diffusion. Before reaching its asymptotic value for t > τs, the transverse diffusion rate decreases as in a sub-diffusive regime.

The numerical simulation reproduces these types of behavior, although the transverse evolution departs slightly from the above formula before reaching the scattering time τ_s. Nevertheless, the agreement between the theory and the numerical simulation holds during the linear growth at the beginning of the evolution and when the evolution approaches the asymptotic behavior. The numerical results confirm the theory we proposed in the previous section for the asymptotic regime. The scattering time is clearly the time beyond which spatial diffusion takes place. We can also note that there is a sub-diffusion regime before the settlement of the transverse diffusion regime.

The anisotropy ratio D _⊥ /D_∥ can be seen in Fig. 5 as a function of ρ′. When the turbulence level η is close to 1 and ρ′ is not too large, the transport appears isotropic D _⊥ /D_∥ ≃ 1. At higher rigidities, its behavior follows the law  ∝ ρ^′ − 2 for all turbulence levels, illustrating the saturation of the transverse coefficient and in agreement with the theoretical prediction.

Fig. 5 Anisotropy ratio D ⊥ /D∥ as function of ρ′ for different levels of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$, as indicated by the various symbols. The dashed line provides a guide for a ρ′-2 scaling.

3.4. Comparisons with previous results

Transverse diffusion at high rigidity, as far as we know, has been poorly studied in the literature. However, we can compare our results with several previous numerical and theoretical studies with different limits.

The seminal study of Giacalone& Jokipii (1999) focused on the propagation of mildly relativistic particles (E = 1 MeV to 1 GeV) in the interplanetary magnetic field ($\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathrm{~}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$). Their simulations provided results for ρ ≤ 1 and η ≤ 0.5. However, they also performed several simulations in which the particle energy and the coherence length remained fixed, while the turbulence level was varied. In particular, they examined the case r_L,0/ℓ_c = 10 for moderate values of $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ (Fig. 6 of their paper) in which D _⊥ /D_∥ is plotted as a function of λ_∥/r_L,0 (λ_∥ denoting the mean free path in the parallel direction). By inspecting their figure, one can see that they varied $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ from 0.05 to 30. As a result, they found a classical scattering theory scaling but no physical explanation was proposed. Strictly speaking, the classical theory is valid only for weak turbulence ($\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathrm{\ll }{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$), which clearly does not apply to those simulations. The present theoretical framework provides a clear explanation of this result, which we confirmed with additional detailed numerical simulations. It is found, for instance, that particles with large rigidities do not interact directly with the magnetic field lines but experience an overall magnetic topology dominated by the mean field with “infinite” coherence length. As a result, the particles execute regular orbits around B₀ and undergo random deflections on the coherence length-scale.

The simulations of Casse et al. (2002) investigated weak as well as strong turbulence regimes where $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{0.1}\mathit{,}\mathrm{99}\mathrm{\right]}$. An FFT algorithm was used to construct the magnetic field. For ρ′ > 1, these authors found evidence of anisotropic scattering D _⊥ /D_∥ < 1 for all turbulence levels. However, only three simulations points were computed in the high rigidity range and the estimate of the power law slope was inaccurate. Nevertheless, a reasonable agreement is obtained when comparing values of D_∥ and D _⊥  with the present results.

Parizot (2004) presented simulations of particle propagation in pure isotropic turbulence. The results in the regime r_L ≫ ℓ_c leads to a diffusion coefficient with a quadratic scaling, D ∝ E², in agreement with our results from Sect. 3.2.

Fig. 6 Ratio D⊥/D∥ as a function of ρ′ for η = 0.99, compared to theorectical predictions and other numerical simulations. Filled diamonds: our simulation results. Star symbols: results from Casse et al. 2002. Solid curve: present theoretical prediction with best-fit D∥ from simulations (see Fig. 2). Dashed curve: present theoretical prediction with analytical D∥ = c2/(3νs). Dot-dashed curve: analytical prediction from Shalchi & Dosch (2009), their Eq. (15).

Shalchi & Dosch (2009) derived an analytical expression for the diffusion anisotropy ratio D _⊥ /D_∥ in the framework of a non-linear guiding centretheory. They assume an isotropic turbulence δB with a mean field B₀. No assumption was made about either the level of turbulence or about particle energy, hence their result should be valid for any particle rigidity and turbulent field strength. An expression of D _⊥ /D_∥ (Eq. (15) in their work) that depends on two parameters was obtained. The first parameter corresponds to the ratio of the mean free path (λ_∥) along the mean field direction to the coherence length ℓ_c of the turbulent field. The second parameter is the turbulence level $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathit{\eta }\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}$. Shalchi & Dosch (2009) thus find that the transport becomes highly anisotropic, meaning D _⊥ /D_∥ ≪ 1 when λ_∥/ℓ_c ≫ 1 and/or $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ is not too large (see Figs. 1 and 2 of their work). Therefore, our present conclusions agree with theirs, at least at a qualitative level. A detailed comparison would require us to define λ_∥ as a function of ρ′, which could be done by using our results of D_∥ for which ρ′ = (4πη/30)^1/2(λ_∥/ℓ_c)^1/2. With this substitution, we can directly compare their predictions to our results. In Fig. 6, we plot the ratio of diffusion coefficients as a function of ρ′ from our numerical simulations and compare these results to both the predictions of Shalchi & Dosch (2009) and the theoretical model developed in Sect. 2. Good agreement is found between the simulation results (diamond symbols) and our theory (solid curve and dashed curve); however, the predictions of Shalchi & Dosch (2009) disagree with the numerical results, increasingly so as the rigidity increases. In particular, their analysis predicts a scaling with a slope  − 2.4 instead of the value of  − 2 observed here. Repeating the same comparisons for each simulated value of δB²/B², we were unable to find agreement between the predictions of Shalchi & Dosch (2009) and our simulations; the predicted values always lie below the numerical results, with a different power-law scaling, comprised between –2.5 for $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ and –2.4 for $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}{\mathrm{10}}^{\mathrm{4}}$. At this point, it could be argued that our definition of λ_∥ as a function of ρ′ is inaccurate. However, on physical grounds, the scaling λ_∥ ∝ ϵ² when remains robust. Therefore, the discrepancy between the power law scalings should not be affected by uncertainties in the numerical prefactors. We think that the “guiding center” assumption in their work is questionable.

4. Summary and some astrophysical applications

4.1. Summary

Our investigation of the diffusion process in small-scale turbulence with a mean field revealed that, despite its smallness, the mean field plays a role in transverse diffusion because the scattering frequency decreases like ϵ^-2, whereas the Larmor frequency decreases like ϵ^-1. Instead of finding a single diffusion coefficient that increases like ϵ², we found an anisotropic diffusion with a transverse coefficient that reaches a limit value at large rigidities. The theory we proposed is based on a single assumption, namely that the correlation time is much smaller than the scattering time, which is valid for both small and large rigidities. The only regime where the theory fails is for a rigidity close to 1 and a high turbulence level; however, the interpolation is obvious. The theory allows us to derive a correct pitch-angle scattering rate and a correct parallel diffusion coefficient for every rigidity. It provides a transverse diffusion coefficient similar to the classical scattering theory formula, despite the arbitrary level of turbulence, which is a correct result for large rigidity. At low rigidity, the present theory is incorrect because it does not take into account the effect of field line wandering described in Casse et al. (2002).

4.2. Particle transport in relativistic shock environments

One major application of the diffusion theory in small-scale turbulence is the transport of supra-thermal particles in the vicinity of a relativistic shock. By crossing the shock transition, electrons and protons reach more or less the same characteristic energy  ⟨ ϵ ⟩  ~ γ_shm_pc² as revealed clearly by particle-in-cell simulations (e.g., Sironi & Spitkovsky 2011). There is a single plasma frequency ω_p∗ ~ ω_pi, where ω_pi is the ion plasma frequency in the upstream or unshocked plasma. This length-scale characterizes the typical length scale of the microturbulence excited in the shock precursor, as transmitted downstream of the shock transition and viewed in the downstream rest frame. The generation of short scale intense micro-turbulence is possible only at low magnetizations of the upstream plasma (Sironi & Spitkovsky 2011), where the magnetization parameter σ is here defined as the flux of magnetic energy crossing the shock over the flux of matter energy, $\mathit{\sigma }\mathrm{\equiv }{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{B}}\mathit{/}\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{u}}{\mathit{c}}^{\mathrm{2}}$ (where θ_B is the angle of the background magnetic field with the shock normal, and ρ_u the unshocked plasma mass density). However, this same level of magnetization also permits the efficient acceleration of particles through a first-order Fermi process at the shock front (Lemoine & Pelletier 2010, 2011). For larger magnetizations – the exact level depending on the shock Lorentz factor, see the above references – the Fermi process cannot develop because of a lack of efficient scattering in the microturbulence (Lemoine et al. 2006; Niemiec et al. 2006; Pelletier et al. 2009). In brief, the development of the Fermi process hinges on the development of micro-turbulence, which itself requires (in the absence of external sources of turbulence) a sufficiently low magnetization level. The situation in which particles are accelerated is by far the most interesting as it should produce directly observable signatures, in the form of radiation and possibly neutrinos.

The transport properties of these accelerated particles is then directly governed by the parallel and perpendicular diffusion coefficients in the limit of large rigidity, as discussed above. We assume that the microturbulence has a typical length-scale close to δ_∗ = c/ω_p ∗  and that a fraction ϵ_B of shock dissipated energy is converted into electromagnetic turbulence, i.e. $\frac{\mathrm{⟨}\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{8}\mathit{\pi }}\mathrm{=}\mathrm{2}{\mathit{ϵ}}_{\mathit{B}}{\mathit{\gamma }}_{\mathrm{sh}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{u}}{\mathit{c}}^{\mathrm{2}}\mathit{,}$(44)where the rigidity of shock accelerated particles of energy ϵ is given by $\mathit{\rho }\mathrm{\approx }{\mathit{ϵ}}_{\mathit{B}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{\delta }}_{\mathrm{\ast }}}{{\mathit{\ell }}_{\mathrm{c}}}\frac{\mathit{ϵ}}{\mathrm{⟨}\mathit{ϵ}\mathrm{⟩}}\mathrm{·}$(45)Current simulations indicate values of ϵ_B ~ 0.01 − 0.1, hence ρ > 1 and all the more so at high energy.

In this regime, the perpendicular diffusion coefficient that we discussed in the previous section becomes particularly relevant, as the mean magnetic field is mostly perpendicular to the shock normal in the downstream frame, since the transverse components (relatively to the shock normal) are increased by $\mathrm{2}\sqrt{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{sh}}$, while the parallel component remains the same as in the upstream frame. Therefore, perpendicular diffusion at high rigidity plays an essential role in the transport of particles in the downstream flow of a relativistic shock.

We consider the diffusive behavior of particles in the downstream rest frame. In this frame the shock front appears to move away with velocity V_shock ≃ c/3. Achieving Fermi cycles requires the particle to return to the shock front. The return time is then measured by identifying shock front with the particle mean displacements $\frac{\mathit{c}}{\mathrm{3}}{\mathit{t}}_{\mathrm{ret}}\mathrm{=}\sqrt{\mathrm{2}{\mathit{D}}_{\mathrm{\perp }}{\mathit{t}}_{\mathrm{ret}}}\mathit{.}$(46)Therefore t_ret = 18D _⊥ /c² and Fermi cycles are possible until t_ret is neither large nor too short. While the first case is constrained by confinement in the acceleration site, the second one is related to the diffusive approximation that is valid only when t_ret ≥ τ_s. Using the second limit to constrain diffusive returns, one obtain D _⊥ /D_∥ ≥ 1/6, equivalent to ${\mathit{\nu }}_{\mathrm{s}}\mathrm{\ge }\sqrt{\mathrm{5}}{\mathrm{\Omega }}_{\mathrm{0}}$ when D _⊥  is replaced by its expression from Eq. (33). Fiducial values for a relativistic shock in the interstellar medium provide an energy limit E_lim ~ 10¹⁹ eV. This limit is somewhat irrelevant because t_ret ≫ R_acc/c at this energy, where R_acc/c is the shock dynamics timescale. Hence, the returns appear to be efficient when the condition τ_s < t_ret ≪ R_acc/c is satisfied. Further investigation would require us to solve a kinetic equation taking into account acceleration, scattering and energy losses processes. Diffusion coefficients obtained in this work may be relevant to providing more realistic results. Previous works assumed Bohm diffusion or isotropic pitch-angle scattering.

Detailed discussion of the performance of the relativistic Fermi process is beyond of the scope of the present paper and is left to future work.

In certain astrophysical settings, the transverse diffusion may play a key role in the transport of particles upstream of a relativistic shock, most particularly if the shock propagates in a wind with a dominant toroidal field at large distances. These circumstances can be encountered in particular when a gamma-ray burst explodes in the wind of the progenitor, or at the termination shock of a pulsar wind.

4.3. High-energy cosmic rays

The above result about transverse diffusion has a broader application than Fermi acceleration at shocks, as it governs the confinement properties of any relativistic flow containing a small-scale turbulence, where “small” is measured relatively to the Larmor radius of the test particles propagating in this flow. This concerns in particular the propagation of very high-energy cosmic rays in our Galaxy. Assuming a coherence length of interstellar turbulence ℓ_c ~ 10 − 100   pc, a mean field intensity of 3   μG approximately and a turbulent field of the same order, the rigidity of particles of energy E is given by , while the Larmor radius . Assuming η ≃ 0.5 in the Galaxy and using Eq. (33), the perpendicular mean free path is then of order λ _⊥  ~ 6   pc with these values of energy and magnetic field. This implies that the escape, or transport across the disk magnetic field of particles of energy  ≥ 10¹⁷   eV is governed by the perpendicular diffusion in the high rigidity regime discussed above. Quite interestingly, this energy range presumably corresponds to the transition between the Galactic and extragalactic cosmic-ray components in the all-particle spectrum.

Finally, one could mention another application of the present discussion, to the field of magnetic reconnection. There, transverse diffusion likely plays a role in the control of particle diffusion across the field lines with small-scale turbulence being associated with the dissipation of magnetic energy. The reconnection rate depends on two fundamental parameters (Lyutikov & Uzdensky 2003): magnetization and the Lundquist number that involves diffusion across field lines. In general, one assumes Bohm diffusion for simplicity but the present work provides the grounds for a more accurate estimate.

References

Appendix A: Average of the time ordered exponential

We solve the differential equation $\dot{{u}}\mathrm{=}\widehat{\begin{array}{}{}_{\mathbf{˜}}\\ \mathbf{\Omega }\end{array}}\mathrm{·}{u}$ in successive iterations that leads to a Dyson series, the average of which is composed of products of the form ${{Â}}_{\mathrm{2}\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathrm{1}}}\mathrm{d}{\mathit{t}}_{\mathrm{2}}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}}\mathrm{d}{\mathit{t}}_{\mathrm{2}\mathit{p}}\hspace{0.17em}\mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathit{...}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}\mathit{p}}\mathrm{\right)}\mathrm{⟩}\mathit{,}$(A.1)which can be compared to Eq. (21).

For a Gaussian process, each average of order 2p products can be divided into a sum of p products of second-order moments, the sum containing (2p − 1) !  !  terms. We assume a stationary random process such that the second order moment is an even function of the time difference. In the white noise limit, the “nested” and “crossed” averages vanish, only the “unconnected” averages remaining in the expansion. Nested terms contain products of the form  ⟨ X(t_i)X(t_l) ⟩  ⟨ X(t_j)X(t_k) ⟩  with t_i ≥ t_j ≥ t_k ≥ t_l, while crossed terms are of the form  ⟨ X(t_i)X(t_k) ⟩  ⟨ X(t_j)X(t_l) ⟩  with t_i ≥ t_j ≥ t_k ≥ t_l. These terms vanish as the various delta functions associated with the second order moments cancel each other as a result of the time ordering in the upper bounds of the integrals. Thus, only the unconnected average remains at each order $\begin{array}{ccc}{{Â}}_{\mathrm{2}\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}\mathit{t}\\ \mathrm{0}\end{array}\mathrm{d}{\mathit{t}}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathrm{1}}}\mathrm{d}{\mathit{t}}_{\mathrm{2}}\mathit{...}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}}\mathrm{d}{\mathit{t}}_{\mathrm{2}\mathit{p}}\hspace{0.17em}\\ & & \mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathit{...}\mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}\mathit{p}\mathrm{-}\mathrm{1}}\mathrm{\right)}\mathrm{·}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}\mathit{p}}\mathrm{\right)}\mathrm{⟩}\mathit{.}\end{array}$(A.2)

We introduce the short-hand notation $\mathrm{⟨}\widehat{\begin{array}{}{}_{\mathbf{˜}}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\widehat{\begin{array}{}{}_{\mathbf{˜}}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{\equiv }{Ĉ}\mathrm{(}{\mathit{t}}_{\mathrm{1}}\mathrm{-}{\mathit{t}}_{\mathrm{2}}\mathrm{)}\mathrm{=}\mathrm{2}{\mathit{\tau }}_{\mathrm{c}}\mathit{\delta }\mathrm{(}{\mathit{t}}_{\mathrm{1}}\mathrm{-}{\mathit{t}}_{\mathrm{2}}\mathrm{)}{Ĉ}{}_{\mathbf{0}}$. Then one can calculate Â_2p(t) by recursion, starting from the last double integral in the product ${{Â}}_{\mathrm{2}\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}{\mathit{\tau }}_{\mathrm{c}}\mathit{t}{\mathrm{)}}^{\mathit{p}}}{\mathit{p}\mathrm{!}}{{Ĉ}}_{\mathbf{0}}^{\mathit{p}}\mathit{.}$(A.3)We consider now the integral of the second order moment ${K̂}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{d}{\mathit{t}}_{\mathrm{2}}\hspace{0.17em}\mathrm{⟨}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\widehat{\begin{array}{}\mathbf{˜}\\ \mathbf{\Omega }\end{array}}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{2}{\mathit{\tau }}_{\mathrm{c}}\mathit{t}{Ĉ}\mathbf{0}\mathit{.}$(A.4)Therefore ${{Â}}_{\mathrm{2}\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{2}}^{\mathit{p}}\mathit{p}\mathrm{!}}{{K̂}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(A.5)hence summing all the terms of the series, $\sum _{\mathit{p}\mathrm{=}\mathrm{0}}^{\mathit{p}\mathrm{=}\mathrm{+}\mathrm{\infty }}{{Â}}_{\mathrm{2}\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\exp \left[\frac{\mathrm{1}}{\mathrm{2}}{K̂}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right]\mathit{.}$(A.6)Further details can be found in Frisch (1966) and Pelletier (1977).

All Figures

	Fig. 1 The diffusion coefficient variation is plotted in units of cLmax/(2π) as a function of rigidity ρ′ in pure turbulence (B0 = 0 or η = 1). The dashed line is drawn as a reference for a scaling Diso ∝ ρ′2. For ρ′ > 1, Diso is indeed proportional to ρ2.
In the text

Fig. 2 The parallel diffusion coefficient D∥ plotted in units of (cLmax/(2π)) as a function of ρ′ for different degrees of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$. Here D∥ ∝ ρ2, as in the case of purely isotropic turbulence without mean field. As long as $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\gg }\mathrm{1}$, the strength of the turbulence does not influence the normalization of D∥.

In the text

Fig. 3 The transverse diffusion coefficient D ⊥  plotted in units of cLmax/(2π) as a function of ρ′ for different degrees of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$. The diffusion coefficient saturates at ρ′ ~ δB/B0. Below this value, its behavior is similar to the parallel diffusion coefficient. Beyond ρ′, its value becomes independent of particle rigidity.

In the text

	Fig. 4 Transition toward parallel and perpendicular diffusion. Before reaching its asymptotic value for t > τs, the transverse diffusion rate decreases as in a sub-diffusive regime.
In the text

	Fig. 5 Anisotropy ratio D ⊥ /D∥ as function of ρ′ for different levels of turbulence $\mathit{\delta }{\mathit{B}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\in }\mathrm{\left[}\mathrm{1}\mathit{,}\mathrm{9999}\mathrm{\right]}$, as indicated by the various symbols. The dashed line provides a guide for a ρ′-2 scaling.
In the text

Fig. 6 Ratio D⊥/D∥ as a function of ρ′ for η = 0.99, compared to theorectical predictions and other numerical simulations. Filled diamonds: our simulation results. Star symbols: results from Casse et al. 2002. Solid curve: present theoretical prediction with best-fit D∥ from simulations (see Fig. 2). Dashed curve: present theoretical prediction with analytical D∥ = c2/(3νs). Dot-dashed curve: analytical prediction from Shalchi & Dosch (2009), their Eq. (15).

In the text