# Particle transport in intense small-scale magnetic turbulence with a mean field

I. Plotnikov, G. Pelletier, M. Lemoine

► To cite this version:

I. Plotnikov, G. Pelletier, M. Lemoine. Particle transport in intense small-scale magnetic turbulence with a mean field. Astronomy and Astrophysics - A&A, 2011, 532, 10.1051/0004-6361/201117182. insu-03624335

HAL Id: insu-03624335

https://insu.hal.science/insu-03624335

Submitted on 30 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Particle transport in intense small-scale magnetic turbulence with a mean field

I. Plotnikov<fnref n="1" />, G. Pelletier<fnref n="1" />, and M. Lemoine<fnref n="2" />

Received 3 May 2011 / Accepted 9 June 2011

## ABSTRACT

Various astrophysical studies have motivated the investigation of the transport of high energy particles in magnetic turbulence, either in the source or en route to the observation sites. For strong turbulence and large rigidity, the pitch-angle scattering rate is governed by a simple law involving a mean free path that increases proportionally to the square of the particle energy. In this paper, we show that perpendicular diffusion deviates from this behavior in the presence of a mean field. We propose an exact theoretical derivation of the diffusion coefficients and show that a mean field significantly changes the transverse diffusion even in the presence of a stronger turbulent field. In particular, the transverse diffusion coefficient is shown to reach a finite value at large rigidity instead of increasing proportionally to the square of the particle energy. Our theoretical derivation is corroborated by a dedicated Monte Carlo simulation. We briefly discuss several possible applications in astrophysics.

**Key words.** magnetic fields – diffusion – scattering – turbulence – methods: analytical – methods: numerical

## 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 $\rho$, the level of turbulence $\eta$, and the power spectrum of magnetic fluctuations in three dimensions (3D) $S_{3d}(\mathbf{k})$. These quantities are defined as

$$ \rho \equiv \frac{\bar{r}_L}{\ell_c} = \frac{\epsilon}{e\bar{B}\ell_c}, \qquad (1) $$

<footnote n="1">¹ UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 38041 Grenoble, France
e-mail: illya plotnikov@obs.ujf-grenoble.fr</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, CNRS – UPMC, 98bis boulevard Arago, 75014 Paris, France</footnote>
where $\bar{r}_L$ denotes the Larmor radius of the particle in the total (mean $\mathbf{B}_0$ and turbulent $\delta\mathbf{B}$) field $\bar{B}$ where $\bar{B}^2 = \mathbf{B}_0^2 + \delta\mathbf{B}^2$, $\epsilon$ the energy of the particle, and $\ell_c$ the coherence length of the fluctuations.

The turbulence level $\eta$ is defined as

$$ \eta \equiv \frac{\langle \delta \mathbf{B}^2 \rangle}{\langle \delta \mathbf{B}^2 \rangle + \mathbf{B}_0^2}, \qquad (2) $$

where $\eta \to 0$ corresponds to weak turbulence and $\eta \to 1$ corresponds to pure turbulence with no mean field.

The correlation function $C(r)$ of the random field

$$ C(r) \equiv \frac{\langle \delta \mathbf{B}(x+r) \delta \mathbf{B}(x) \rangle}{\langle \delta \mathbf{B}^2 \rangle}, \qquad (3) $$

can be written in terms of the one-dimensional power spectrum $S(k) \propto k^2 S_{3d}(\mathbf{k})$

$$ C(r) = \frac{\int dk S(k) \sin(kr)/(kr)}{\int dk S(k)}. \qquad (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) \propto (k/k_{\min})^{-\beta}$ for $k_{\min} \le k \le k_{\max}$ and zero otherwise, one finds for the coherence length $\ell_c \approx 0.77 k_{\min}^{-1}$. Alternatively, one can define the coherence length as we do here, to be

$$ \ell_c \equiv \int_0^{+\infty} dr C(r), \qquad (5) $$

where one then derives in a straightforward way

$$ \ell_c = \frac{\pi}{2} \frac{1}{\eta} \int_0^{+\infty} dk k^{-1} S(k), \qquad (6) $$

and the presence of $1/\eta$ results from our choice of normalization for the power spectrum

$$ \int_0^{+\infty} dk S(k) \equiv \eta, \qquad (7) $$

where in practice the spectrum is bounded between $k_{\min}$ and $k_{\max}$. Both definitions for $\ell_c$ coincide to within a factor close to unity. As a function of the spectrum index $\beta$, the coherence length is close to either $k_{\min}^{-1}$ on large scale for $\beta > 1$, or to $k_{\max}^{-1}$ on small scales for $\beta < 1$.

The scattering frequency $v_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

$$ v_s \approx \frac{\pi c}{3 \bar{r}_L^2} \frac{\int_{k \bar{r}_L > 1} k^{-1} S(k) dk}{\int S(k) dk}, \qquad (8) $$

an expression that extends to the strong turbulence regime the results of the quasi-linear theory. This leads to the scalings

$$ \begin{aligned} v_s &\approx \frac{2}{3} \eta \frac{c}{\ell_c} \rho^{\beta-2} \quad (\rho \ll 1) \\ v_s &\approx \frac{2}{3} \eta \frac{c}{\ell_c \rho^2} \quad (\rho \gg 1). \end{aligned} \qquad (9) $$

The Bohm scaling holds only in the very special case where $\beta = 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

$$ \tau_c \equiv \int_0^{+\infty} C(|\Delta x(\tau)|) d\tau, \qquad (10) $$

where $\Delta x(\tau)$ represents the displacement after a time $\tau$ 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 $\rho \gg 1$, then $\tau_c \sim \ell_c/c$. This correlation time is much shorter than the scattering time $v_c^{-1} \sim \eta^{-1} \rho^2 \ell_c/c$ in this regime. The correlation time $\tau_c$ can be recovered from Eq. (10) by using the ballistic approximation $\Delta x(\tau) \approx c\tau$, which is appropriate in this regime $\rho \gg 1$, in which case Eq. (4) leads to $\tau_c = (\pi/2)(\beta - 1)\beta^{-1}k_{\min}^{-1}/c \sim \ell_c/c$. We note that in the special case where the power-law index of turbulence $\beta = 1$ (Bohm regime) Eqs. (4) and (5) lead to $\ell_c = (\lambda_{\min}/4) \log(\lambda_{\max}/\lambda_{\min})$, where $\lambda_{\min}$ and $\lambda_{\max}$ are the shortest and the longest wavelengths of turbulence.

If a particle experiences a chaotic motion on a length-scale smaller than $\ell_c$, corresponding to the regime $\rho \ll 1$, then the estimate is more complicated to obtain but one finds that $\tau_c \sim \rho^\beta \ell_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} \rho^{-1}$, which give rise to gyroresonant interactions with the particle of rigidity $\rho$. The modes with wavelengths longer than the Larmor radius (i.e., $k < k_{\min} \rho^{-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 $\rho$. The heuristic estimate for $\rho < 1$ is consistent with quasi-linear theory when $\eta \ll 1$ and with numerical results in the strong turbulence regime (Casse et al. 2002) can then be written as

$$ \tau_c \approx \frac{1}{\eta c} \int_{k > k_{\min} \rho^{-1}} dk k^{-1} S(k), \qquad (11) $$

which bears some resemblance to the case discussed before for $\rho \gg 1$, except that $\rho$ 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} \rho^{-1}$ for the reasons given above. The calculation then implies that $\tau_c \sim \rho^\beta \ell_c/c$ as announced. The particle trajectory undergoes decoherence before traveling $\ell_c$ because of the large number of wiggles in the random field.

Thus, except for $\eta \sim 1$ and $\rho \sim 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_\parallel = c^2/(3v_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_\parallel$ (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_\perp = \eta^{2.3} D_\parallel$ 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 $\rho \gg 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 $\rho \gg 1$, the velocity changes by $1/\rho$ 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{d}{dt}v = [\hat{\Omega}_0 + \delta\hat{\Omega}(t)] \cdot v. \qquad (12) $$

The quantities $\hat{\Omega}_0$ and $\delta\hat{\Omega}(t)$ are rotation operators developed as linear combinations of the generators of the Lie algebra of the rotation group, $\hat{L}_1, \hat{L}_2, \hat{L}_3$

$$ \hat{L}_1 := \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix}, \quad \hat{L}_2 := \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix}, \quad \hat{L}_3 := \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (13) $$

In detail, $\hat{\Omega}_0 = \Omega_0 B_0^i \hat{L}_i / B_0$, where $B_0^i$ denotes the ith component of $B_0$ and $\Omega_0 \equiv c/r_{L,0}$ the Larmor pulsation defined with respect to the mean field. With this notation, $\hat{\Omega}_0 \cdot v = \Omega_0 v \times B_0 / B_0$. The operator $\delta\hat{\Omega}(t)$ is decomposed in a similar way as the generators of the rotation group, and $\delta\Omega \equiv c/r_L$, where $r_L$ is now measured relatively to $\delta B$.

To solve the equation of motion, one uses an auxiliary variable $u$ that is defined as

$$ v(t) \equiv \hat{R}_0(t) \cdot u(t), \qquad (14) $$

where

$$ \hat{R}_0(t) \equiv \exp(t\hat{\Omega}_0). \qquad (15) $$

We then define

$$ \hat{\Omega}(t) \equiv \hat{R}_0(t)^{-1} \cdot \delta\hat{\Omega}(t) \cdot \hat{R}_0(t), \qquad (16) $$

one finds that $u(t)$ obeys

$$ \frac{d}{dt}u = \hat{\Omega}(t) \cdot u. \qquad (17) $$

This equation is solved as

$$ u(t) = \mathcal{T} \exp \left[ \int_0^t \hat{\Omega}(t') dt' \right] \cdot u(0). \qquad (18) $$

Because the operator in the exponent is time dependent, to preserve the exponential character of the solution, a time-ordering operator $\mathcal{T}$ has to be introduced, as we now explain.

We note that $u(0) = v(0)$, thus the solution for $v$ is given by

$$ v(t) = \hat{R}_0(t) \cdot \mathcal{T} \exp \left[ \int_0^t \hat{\Omega}(t') dt' \right] \cdot v(0). \qquad (19) $$

The regular part of the field generates the regular rotation matrix $\hat{R}_0(t)$, while the exponential accounts for the effect of the turbulent part. The time-ordering operator $\mathcal{T}$ maintains the chronological order of the products in the non-commuting $\hat{\Omega}(t_k)$ in the expansion of the exponential operator, i.e.

$$ \begin{aligned} \mathcal{T} \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) &= \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \quad \text{if } t_1 > t_2, \\ &= \hat{\Omega}(t_2) \cdot \hat{\Omega}(t_1) \quad \text{if } t_2 > t_1, \end{aligned} \qquad (20) $$

and so on for higher order products. Alternatively, the time-ordered expansion can be written as a Dyson series

$$ \mathcal{T} \exp \left[ \int_0^t \hat{\Omega}(t') dt' \right] \equiv 1 + \sum_{n=1}^{n=+\infty} \int_0^t dt_1 \dots \int_0^{t_{n-1}} dt_n \hat{\Omega}(t_1) \dots \hat{\Omega}(t_n). \qquad (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

$$ \left\langle \mathcal{T} \exp \left[ \int_0^t \hat{\Omega}(t') dt' \right] \right\rangle = \mathcal{T} \exp \left[ \frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \left\langle \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \right\rangle \right]. \qquad (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_0$, with $\delta B \cdot B_0 = 0$.

We define the projection operators $\hat{\pi}_\perp$ on the plane transverse to $B_0$ and the projection operator $\hat{\pi}_\parallel$ along $B_0$. We now define the correlation function of the random rotation matrices $\delta\hat{\Omega}$

$$ \langle \delta\hat{\Omega}(t_1) \delta\hat{\Omega}(t_2) \rangle = \langle \delta\Omega^i(t_1) \delta\Omega^j(t_2) \rangle \hat{L}_i \hat{L}_j, \qquad (23) $$

where

$$ \langle \delta\Omega^i(t_1) \delta\Omega^j(t_2) \rangle = 2\tau_c \delta(t_1 - t_2) \left[ \frac{1}{2} \langle \delta\Omega_\perp^2 \rangle \hat{\pi}_\perp^{ij} + \langle \delta\Omega_\parallel^2 \rangle \hat{\pi}_\parallel^{ij} \right]. \qquad (24) $$

The scalars $\langle \delta\Omega_\perp^2 \rangle$ and $\langle \delta\Omega_\parallel^2 \rangle$ characterize the relative strengths of the turbulence in the parallel (to $B_0$) and perpendicular directions. In particular, for 3D isotropic turbulence, $\langle \delta\Omega_\perp^2 \rangle = 2\langle \delta\Omega_\parallel^2 \rangle$, in which case the above correlator becomes proportional to the identity. Then, using the properties of $\hat{\pi}_\perp, \hat{\pi}_\parallel$ and the $\hat{L}_i$, one finds

$$ \langle \delta\hat{\Omega}(t_1) \delta\hat{\Omega}(t_2) \rangle = -2\tau_c \delta(t_1 - t_2) \times \left[ \langle \delta\Omega^2 \rangle \hat{\mathbf{1}} - \langle \delta\Omega_\parallel \rangle^2 \hat{\pi}_\parallel - \frac{1}{2} \langle \delta\Omega_\perp \rangle^2 \hat{\pi}_\perp \right], \qquad (25) $$

where $\langle \delta\Omega^2 \rangle = \langle \delta\Omega_\parallel \rangle^2 + \langle \delta\Omega_\perp \rangle^2$. We note that the above correlation function holds for $\delta\hat{\Omega}$, which should not be confused with
$\hat{\Omega}$, the latter being the quantity of relevance for calculating the transport properties, as expressed in Eq. (22). However,

$$ \langle \hat{\Omega}(t_1) \hat{\Omega}(t_2) \rangle = e^{-t_1 \hat{\Omega}_0} \langle \delta\hat{\Omega}(t_1) e^{-(t_2-t_1)\hat{\Omega}_0} \delta\hat{\Omega}(t_2) \rangle e^{t_2 \hat{\Omega}_0}, \quad (26) $$

and, because $[\hat{\pi}_{\parallel}, e^{t\hat{\Omega}_0}] = [\hat{\pi}_{\perp}, e^{t\hat{\Omega}_0}] = 0$, the correlation function for $\hat{\Omega}$ is the same as that for $\delta\hat{\Omega}$.

Using Eq. (22), one then finds the solution for $v$:

$$ \langle v(t) \rangle = \hat{R}_0(t) \cdot \exp \left[ \frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \langle \hat{\Omega}(t_1) \hat{\Omega}(t_2) \rangle \right] v(0), \quad (27) $$

where the average is taken over the possible realizations of the turbulent field. This leads to

$$ \langle v(t) \rangle = \hat{R}_0(t) \cdot \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \langle \delta\Omega_{\parallel}^2 \rangle \hat{\pi}_{\parallel} - \frac{1}{2} \langle \delta\Omega_{\perp}^2 \rangle \hat{\pi}_{\perp} \right) \right] v(0). \quad (28) $$

Using the properties of $\hat{\pi}_{\parallel}$ and $\hat{\pi}_{\perp}$, this can be rewritten as

$$ \langle v(t) \rangle = \hat{R}_0(t) \cdot \left\{ \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \langle \delta\Omega_{\parallel}^2 \rangle \right) \right] \hat{\pi}_{\parallel} \right. \\ \left. + \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \frac{1}{2} \langle \delta\Omega_{\perp}^2 \rangle \right) \right] \hat{\pi}_{\perp} \right\} v(0). \quad (29) $$

Therefore, one derives the general results

$$ \langle v_{\parallel}(0) v_{\parallel}(t) \rangle = \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \langle \delta\Omega_{\parallel}^2 \rangle \right) \right] \langle v_{\parallel}(0)^2 \rangle. \quad (30) $$

In the transverse direction,

$$ \begin{aligned} \langle v_{\perp}(0) \cdot v_{\perp}(t) \rangle &= \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \frac{1}{2} \langle \delta\Omega_{\perp}^2 \rangle \right) \right] \\ & \quad \times {}^{\mathrm{T}}v_{\perp}(0) \cdot \hat{R}_0 \cdot v_{\perp}(0) \\ &= \exp \left[ -t \tau_c \left( \langle \delta\Omega^2 \rangle - \frac{1}{2} \langle \delta\Omega_{\perp}^2 \rangle \right) \right] \\ & \quad \times \cos(\Omega_0 t) \langle v_{\perp}(0)^2 \rangle. \end{aligned} \quad (31) $$

The last equality follows from developing the exponential $\hat{R}_0 = \exp(t\hat{\Omega}_0)$, noting that $\hat{\Omega}_0 = \Omega_0 \hat{L}_3$ for $B_0$ oriented along $z$, $\hat{L}_3^{2n} = (-1)^n \hat{\pi}_{\perp}, \hat{L}_3^{2n+1} = (-1)^n \hat{L}_3, {}^{\mathrm{T}}v_{\perp}(0) \cdot \hat{\pi}_{\perp} \cdot v_{\perp}(0) = \langle v_{\perp}(0)^2 \rangle$, and ${}^{\mathrm{T}}v_{\perp}(0) \cdot \hat{L}_3 \cdot v_{\perp}(0) = 0$.

The parallel $D_{\parallel}$ and perpendicular $D_{\perp}$ diffusion coefficients are directly obtained from the correlation functions of the velocity components after averaging over the initial velocities

$$ \begin{aligned} D_{\parallel} &= \int_{0}^{+\infty} dt \langle v_{\parallel}(0) v_{\parallel}(t) \rangle, \\ D_{\perp} &= \frac{1}{2} \int_{0}^{+\infty} dt \langle v_{\perp}(0) \cdot v_{\perp}(t) \rangle. \end{aligned} \quad (32) $$

Using Eqs. (30) and (31), this leads to

$$ \begin{aligned} D_{\parallel} &= \frac{1}{3} \frac{c^2}{v_s}, \\ D_{\perp} &= \frac{1}{3} c^2 \frac{v_{\perp}}{v_{\perp}^2 + \Omega_0^2}, \end{aligned} \quad (33) $$

where

$$ \begin{aligned} v_s &= \tau_c \langle (\delta\Omega^2) - (\delta\Omega_{\parallel}^2) \rangle, \\ v_{\perp} &= \tau_c \langle (\delta\Omega^2) - \frac{1}{2} \langle \delta\Omega_{\perp}^2 \rangle \rangle. \end{aligned} \quad (34) $$

These expressions for $D_{\perp}$ 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, $\langle \delta\Omega_{\perp}^2 \rangle = 2\langle \delta\Omega_{\parallel}^2 \rangle = \frac{2}{3}\langle \delta\Omega^2 \rangle$, so that

$$ v_{\perp} = v_s = \frac{2}{3} \langle \delta\Omega^2 \rangle \tau_c = \frac{2}{3} \frac{\eta c}{\rho^2 \ell_c}. \quad (35) $$

One may note that the expression for $v_s$ matches that derived from a random walk argument for pitch angle diffusion. We also note that the above calculation for $v_s$ may be applied to the regime $\rho \ll 1$, as long as the correlation time is shorter than the scattering time. This is true in the case of $v_s = (2/3)\eta\rho^{\delta-2}$, which is the standard quasi-linear theory result. The result for the perpendicular coefficient cannot, of course, be extended to the regime $\rho \ll 1$, as the above calculation does not account for field line wandering.

In case (B), for 2D transverse isotropic turbulence, $\langle \delta\Omega_{\parallel}^2 \rangle = 0$, $\langle \delta\Omega_{\perp}^2 \rangle = \langle \delta\Omega^2 \rangle$, hence

$$ v_s = 2v_{\perp} = \langle \delta\Omega^2 \rangle \tau_c. \quad (36) $$

This demonstrates that the transverse diffusion coefficients follows the scalings, which we express here for case (A), i.e., isotropic turbulence

$$ \begin{aligned} D_{\perp} &\approx D_{\parallel} \approx \frac{1}{2} c \ell_c \rho^2 / \eta \quad (1 \ll \rho \ll \bar{B} / B_0), \\ D_{\perp} &\approx \frac{2}{9} c \ell_c \frac{\bar{B}^2}{B_0^2} \quad (\bar{B} / B_0 \ll \rho). \end{aligned} \quad (37) $$

The transition between these two regimes takes place at $\rho \sim \bar{B}/B_0 \approx \eta^{1/2}(1-\eta)^{-1/2}$, corresponding to $v_s \sim \Omega_0$. At larger rigidities, the perpendicular diffusion coefficient remains constant, while the parallel diffusion coefficient continues to increase as $\rho^2$.

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_0 + \delta B$ as before. The regular field is oriented along $z$, and $\delta 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 $\delta B$ by summing over plane wave modes ($N_{\text{mod}}$) with turbulent wavelengths extending from $L_{\text{min}} = 1 \equiv 2\pi/k_{\text{max}}$ to $L_{\text{max}} \equiv 2\pi/k_{\text{min}}$, the power spectrum following a truncated power law between $k_{\text{min}}$ and $k_{\text{max}}$. In detail

$$ \delta B(\mathbf{x}) = \sum_n G_n(k_n) \xi_n \cos(k_n \cdot \mathbf{x} + \beta_n). \quad (38) $$
With Fourier modes of amplitude $G_n$, and wave vectors $k_n = \frac{2\pi}{L_n}e_k$ isotropically distributed, the unitary vector $\xi_n$ is perpendicular to $k_n$ in order to ensure that $\nabla \cdot \delta\mathbf{B}=0$, and $\beta_n \in [0, 2\pi]$ represents the random phase. The power spectrum is normalized by the turbulence parameter $\eta$ introduced earlier such that $\langle \delta\mathbf{B}^2 \rangle = B_0^2\eta/(1-\eta)$. For definiteness, the mode amplitudes are constructed according to a Kolmogorov cascade with logarithmic spacing between wavenumbers: $G_n \propto k_n^{-5/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 $\tilde{r}_L$ 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 $100\rho\tilde{r}_L/c$, 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 $\gamma$

$$ \frac{dv}{dt} = \frac{q}{\gamma mc} v \times (\mathbf{B}_0 + \delta\mathbf{B}). \quad (39) $$

At this point, we define the numerical rigidity $\rho' \equiv 2\pi\tilde{r}_L/L_{max}$, which differs from the previous physical definition by a numerical factor of order unity, as discussed earlier. The exact relation between $\ell_c$ and $L_{max}/2\pi$ depends on the dynamic range and the power-law index of turbulence. In the following, the conversion factor between both rigidities is derived using $\ell_c \approx 0.1 L_{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{aligned} D_{ij}(t) &= \frac{1}{2N_p N_{field}} \sum_{n=1}^{N_{field}} \sum_{k=1}^{N_p} \frac{(x_i(t) - x_i(t_0))(x_j(t) - x_j(t_0))_{n,k}}{t - t_0} \\ &= \frac{\langle \Delta x_i \Delta x_j \rangle}{2\Delta t}. \end{aligned} \quad (40) $$

The average is performed over different particle trajectories and different field realizations. For each value of $\rho$, we take $N_{field} \times N_p = 10^3$ different trajectories with random initial velocity directions. The asymptotic value for $t \to \infty$ (plateau)

<figure><img src="image_2.png" /><figcaption>Fig. 1. The diffusion coefficient variation is plotted in units of $cL_{max}/(2\pi)$ as a function of rigidity $\rho'$ in pure turbulence ($B_0 = 0$ or $\eta = 1$). The dashed line is drawn as a reference for a scaling $D_{iso} \propto \rho'^2$. For $\rho' > 1$, $D_{iso}$ is indeed proportional to $\rho^2$.</figcaption></figure>

is roughly constant and defines the actual diffusion regime. It gives the diffusion coefficient as $D_{ij}$ as $t \to \infty$, precisely when $t \gg v_s^{-1}$. This method of coefficient estimation appears precise enough for an integration involving $10^3$ 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 = 0$) and a strong turbulent case with $\delta B \gg B_0$. Results are presented in the following sub-sections.

## 3.2. Pure turbulence $B_0 = 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 $\rho' \gg 1$.

Here we set $\eta = 1$ and $\delta 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

$$ D_{iso} = \frac{\langle \Delta x^2 \rangle + \langle \Delta y^2 \rangle + \langle \Delta z^2 \rangle}{6\Delta t}. \quad (41) $$

Figure 1 shows numerical values calculated for $\rho'$ going from 1 to 100. The diffusion coefficient is plotted in units of $cL_{max}/(2\pi)$ as a function of rigidity $\rho'$. A power law is observed for $1 < \rho' < 100$, as predicted by the theory, namely $D_{iso} \propto \rho'^2 \propto \epsilon^2$. One may be able to discern a slight deviation at $\rho' \approx 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} \sim 1$. We retain however the current field configuration, taking this effect into account when interpreting the results.

## 3.3. Weak mean field $B_0 < \delta B$

We now consider the case where a constant weak mean field $B_0$ along the $z$ direction is present. In this case, two different
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The parallel diffusion coefficient $D_\parallel$ plotted in units of $(cL_{\max}/(2\pi))$ as a function of $\rho'$ for different degrees of turbulence $\delta B^2/B_0^2 \in [1, 9999]$. Here $D_\parallel \propto \rho'^2$, as in the case of purely isotropic turbulence without mean field. As long as $\delta B^2/B_0^2 \gg 1$, the strength of the turbulence does not influence the normalization of $D_\parallel$.</figcaption></figure>

diffusion coefficients are defined $D_\parallel = D_{zz}$ and $D_\perp = (D_{xx} + D_{yy})/2$. Overall, we explored five different levels of turbulence $\eta = \{0.5, 0.9, 0.99, 0.999, 0.9999\}$, spanning five orders of magnitude in $\delta B^2/B_0^2$. The rigidity $\rho'$ ranges from 1 to 100 for each value of $\eta$. At each calculation point $\{\eta, \rho'\}$, the coefficients are evaluated by averaging over $10^3$ particles (10 particles $\times$ 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_\parallel \propto \rho'^2$. For $\eta > 0.5$, the turbulence level has almost no influence on the value of $D_\parallel \approx 0.9c\bar{\tau}_L\rho'$. Therefore as expected, the mean field, as long as it remains weak enough, has no influence on the diffusion of particles along its direction.

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 $\rho'$ for different degrees of turbulence. In each case, its value saturates to a constant value when $\rho' \sim \delta B/B_0$. This value behaves proportionally to the turbulence degree; in detail, $D_\perp \approx 0.13c(L_{\max}/2\pi)\delta B^2/B_0^2$, in excellent agreement with our theoretical prediction from Eq. (37). Individual particle trajectories reveal a weakly perturbed helical path when $\rho' \gg 1$. Therefore, a strong small-scale turbulence acts as a collection of small-scale scattering centers, each producing a small deflection.

According to the theory, $D_\parallel$ is the limit of a function $c^2g_\parallel(t)/3$ as $t \to \infty$, precisely as $t > t_s$, the function being

$$g_\parallel(t) = \frac{1 - e^{-v_s t}}{v_s}. \qquad (42)$$

In a similar way $D_\perp$ is the limit of a function $c^2g_\perp(t)/3$ as $t \to \infty$, in addition to when $t > t_s$, the function being

$$g_\perp(t) = \frac{v_s}{\Omega_0^2 + v_s^2} \left\{ 1 - e^{-v_s t} \left[ \cos(\Omega_0 t) - \frac{\Omega_0}{v_s} \sin(\Omega_0 t) \right] \right\}. \quad (43)$$

The numerical simulation reproduces these types of behavior, although the transverse evolution departs slightly from the above formula before reaching the scattering time $\tau_s$. Nevertheless, the agreement between the theory and the numerical simulation holds during the linear growth at the beginning of the evolution

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** The transverse diffusion coefficient $D_\perp$ plotted in units of $cL_{\max}/(2\pi)$ as a function of $\rho'$ for different degrees of turbulence $\delta B^2/B_0^2 \in [1, 9999]$. The diffusion coefficient saturates at $\rho' \sim \delta B/B_0$. Below this value, its behavior is similar to the parallel diffusion coefficient. Beyond $\rho'$, its value becomes independent of particle rigidity.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Transition toward parallel and perpendicular diffusion. Before reaching its asymptotic value for $t > \tau_s$, the transverse diffusion rate decreases as in a sub-diffusive regime.</figcaption></figure>

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_\perp/D_\parallel$ can be seen in Fig. 5 as a function of $\rho'$. When the turbulence level $\eta$ is close to 1 and $\rho'$ is not too large, the transport appears isotropic $D_\perp/D_\parallel \approx 1$. At higher rigidities, its behavior follows the law $\propto \rho'^{-2}$ for all turbulence levels, illustrating the saturation of the transverse coefficient and in agreement with the theoretical prediction.

### 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Anisotropy ratio $D_{\perp}/D_{\parallel}$ as function of $\rho'$ for different levels of turbulence $\delta B^2/B_0^2 \in [1, 9999]$, as indicated by the various symbols. The dashed line provides a guide for a $\rho'^{-2}$ scaling.</figcaption></figure>

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 ($\delta B^2 \sim B_0^2$). Their simulations provided results for $\rho \le 1$ and $\eta \le 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}/\ell_c = 10$ for moderate values of $\delta B^2/B_0^2$ (Fig. 6 of their paper) in which $D_{\perp}/D_{\parallel}$ is plotted as a function of $\lambda_{\parallel}/r_{L,0}$ ($\lambda_{\parallel}$ denoting the mean free path in the parallel direction). By inspecting their figure, one can see that they varied $\delta B^2/B_0^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 ($\delta B^2 \ll B_0^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_0$ 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 $\delta B^2/B_0^2 \in [0.1, 99]$. An FFT algorithm was used to construct the magnetic field. For $\rho' > 1$, these authors found evidence of anisotropic scattering $D_{\perp}/D_{\parallel} < 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_{\parallel}$ and $D_{\perp}$ with the present results.

Parizot (2004) presented simulations of particle propagation in pure isotropic turbulence. The results in the regime $r_L \gg \ell_c$ leads to a diffusion coefficient with a quadratic scaling, $D \propto E^2$, in agreement with our results from Sect. 3.2.

Shalchi & Dosch (2009) derived an analytical expression for the diffusion anisotropy ratio $D_{\perp}/D_{\parallel}$ in the framework of a non-linear guiding centre theory. They assume an isotropic turbulence $\delta B$ with a mean field $B_0$. 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_{\perp}/D_{\parallel}$ (Eq. (15) in

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Ratio $D_{\perp}/D_{\parallel}$ as a function of $\rho'$ for $\eta = 0.99$, compared to theoretical 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_{\parallel}$ from simulations (see Fig. 2). Dashed curve: present theoretical prediction with analytical $D_{\parallel} = c^2/(3v_s)$. Dot-dashed curve: analytical prediction from Shalchi & Dosch (2009), their Eq. (15).</figcaption></figure>

their work) that depends on two parameters was obtained. The first parameter corresponds to the ratio of the mean free path ($\lambda_{\parallel}$) along the mean field direction to the coherence length $\ell_c$ of the turbulent field. The second parameter is the turbulence level $\delta B^2/B_0^2 = \eta/(1-\eta)$. Shalchi & Dosch (2009) thus find that the transport becomes highly anisotropic, meaning $D_{\perp}/D_{\parallel} \ll 1$ when $\lambda_{\parallel}/\ell_c \gg 1$ and/or $\delta B^2/B_0^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 $\lambda_{\parallel}$ as a function of $\rho'$, which could be done by using our results of $D_{\parallel}$ for which $\rho' = (4\pi\eta/30)^{1/2}(\lambda_{\parallel}/\ell_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 $\rho'$ 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 $\delta B^2/B^2$, 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 $\delta B^2/B_0^2 = 1$ and $-2.4$ for $\delta B^2/B_0^2 = 10^4$. At this point, it could be argued that our definition of $\lambda_{\parallel}$ as a function of $\rho'$ is inaccurate. However, on physical grounds, the scaling $\lambda_{\parallel} \propto \epsilon^2$ when $\bar{r}_L \gg \ell_c$ 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 $\epsilon^{-2}$, whereas the Larmor frequency decreases like $\epsilon^{-1}$. Instead of finding a single diffusion coefficient that increases like $\epsilon^2$, 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 $\langle\epsilon\rangle \sim \gamma_{sh}m_p c^2$ as revealed clearly by particle-in-cell simulations (e.g., Sironi & Spitkovsky 2011). There is a single plasma frequency $\omega_{p*} \sim \omega_{pi}$, where $\omega_{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 $\sigma$ is here defined as the flux of magnetic energy crossing the shock over the flux of matter energy, $\sigma \equiv B_0^2 \sin^2 \theta_B / (4\pi\rho_u c^2)$ (where $\theta_B$ is the angle of the background magnetic field with the shock normal, and $\rho_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 $\delta_* = c/\omega_{p*}$ and that a fraction $\epsilon_B$ of shock dissipated energy is converted into electromagnetic turbulence, i.e.

$$ \frac{\langle\delta B^2\rangle}{8\pi} = 2\epsilon_B\gamma_{sh}^2\rho_u c^2, \qquad (44) $$

where the rigidity of shock accelerated particles of energy $\epsilon$ is given by

$$ \rho \approx \epsilon_B^{-1/2} \frac{\delta_*}{\ell_c} \frac{\epsilon}{\langle\epsilon\rangle}. \qquad (45) $$

Current simulations indicate values of $\epsilon_B \sim 0.01-0.1$, hence $\rho > 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 $2\sqrt{2}\gamma_{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} \approx 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{c}{3} t_{\text{ret}} = \sqrt{2D_{\perp} t_{\text{ret}}}. \qquad (46) $$

Therefore $t_{\text{ret}} = 18D_{\perp}/c^2$ and Fermi cycles are possible until $t_{\text{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_{\text{ret}} \ge \tau_s$. Using the second limit to constrain diffusive returns, one obtain $D_{\perp}/D_{\parallel} \ge 1/6$, equivalent to $v_s \ge \sqrt{5}\Omega_0$ when $D_{\perp}$ is replaced by its expression from Eq. (33). Fiducial values for a relativistic shock in the interstellar medium provide an energy limit $E_{\text{lim}} \sim 10^{19}$ eV. This limit is somewhat irrelevant because $t_{\text{ret}} \gg R_{\text{acc}}/c$ at this energy, where $R_{\text{acc}}/c$ is the shock dynamics timescale. Hence, the returns appear to be efficient when the condition $\tau_s < t_{\text{ret}} \ll R_{\text{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 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 $\ell_c \sim 10-100$ pc, a mean field intensity of $3\mu G$ approximately and a turbulent field of the same order, the rigidity of particles of energy $E$ is given by $\rho \approx 2(E/10^{17} \text{eV})(\ell_c/10 \text{pc})^{-1}(\bar{B}/5 \mu\text{G})^{-1}$, while the Larmor radius $r_L \approx 20 \text{ pc}(E/10^{17} \text{eV})(\bar{B}/5 \mu\text{G})^{-1}$. Assuming $\eta \approx 0.5$ in the Galaxy and using Eq. (33), the perpendicular mean free path is then of order $\lambda_{\perp} \sim 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<sup>17</sup> 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.

## Appendix A: Average of the time ordered exponential

We solve the differential equation $\dot{u} = \hat{\Omega} \cdot u$ in successive iterations that leads to a Dyson series, the average of which is composed of products of the form

$$ \hat{A}_{2p}(t) = \int_0^t dt_1 \int_0^{t_1} dt_2 \dots \int_0^{t_{2p-1}} dt_{2p} \langle \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \dots \hat{\Omega}(t_{2p}) \rangle, \quad (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 $\langle X(t_i)X(t_l) \rangle \langle X(t_j)X(t_k) \rangle$ with $t_i \ge t_j \ge t_k \ge t_l$, while crossed terms are of the form $\langle X(t_i)X(t_k) \rangle \langle X(t_j)X(t_l) \rangle$ with $t_i \ge t_j \ge t_k \ge 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

$$ \hat{A}_{2p}(t) = \int_0^t dt_1 \int_0^{t_1} dt_2 \dots \int_0^{t_{2p-1}} dt_{2p} \langle \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \dots \hat{\Omega}(t_{2p-1}) \cdot \hat{\Omega}(t_{2p}) \rangle. \quad (A.2) $$

We introduce the short-hand notation $\langle \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \rangle = \hat{C}(t_1 - t_2) = 2\tau_c\delta(t_1 - t_2)\hat{C}_0$. Then one can calculate $\hat{A}_{2p}(t)$ by recursion, starting from the last double integral in the product

$$ \hat{A}_{2p}(t) = \frac{(\tau_c t)^p}{p!} \hat{C}_0^p. \quad (A.3) $$

We consider now the integral of the second order moment

$$ \hat{K}(t) = \int_0^t dt_1 \int_0^t dt_2 \langle \hat{\Omega}(t_1) \cdot \hat{\Omega}(t_2) \rangle = 2\tau_c t \hat{C}_0. \quad (A.4) $$

Therefore

$$ \hat{A}_{2p}(t) = \frac{1}{2^p p!} \hat{K}^p(t), \quad (A.5) $$

hence summing all the terms of the series,

$$ \sum_{p=0}^{p=+\infty} \hat{A}_{2p}(t) = \exp\left[\frac{1}{2}\hat{K}(t)\right]. \quad (A.6) $$

Further details can be found in Frisch (1966) and Pelletier (1977).

## References

* Aloisio, R., & Berezinsky, V. 2004, ApJ, 612, 900

* Bell, A. R. 2004, MNRAS, 353, 550

* Bell, A. R., & Lucek, S. G. 2001, MNRAS, 321, 433

* Candia, J., & Roulet, E. 2004, JCAP, 10, 007

* Casse, F., Lemoine, M., & Pelletier, G. 2002, Phys. Rev. D, 65, 023002

* Fatuzzo, M., Melia, F., Todd, E., Adams, F. C., et al. 2010, ApJ, 725, 515

* Frisch, U. 1966, Annales d'Asphysique, 29, 645

* Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204

* Jokipii, J. R. 1966, ApJ, 146, 480

* Jokipii, J. R. 1973, ApJ, 183, 1029

* Lemoine, M., & Pelletier, G. 2010, MNRAS, 402, 321

* Lemoine, M., & Pelletier, G. 2011, MNRAS [arXiv:1102.1308]

* Lemoine, M., Pelletier, G., & Revenu, B. 2006, ApJL, 645, 129

* Lyutikov, M., & Uzdensky, D. 2003, ApJ, 589, 893

* Niemiec, A., Ostrowsky, M., & Pohl, M. 2006, ApJ, 650, 1020

* Parizot, E. 2004, Nuc. Phys. B, 136, 169

* Pelletier, G. 1977, J. Plasma Phys., 18, 49

* Pelletier, G., Lemoine, M., & Marcowith, A. 2009, MNRAS, 393, 587

* Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1986, Numerical Recipes (Cambridge: Cambridge Univ. Press)

* Shalchi, & Dosch 2009, Phys. Rev. D, 79, 083001

* Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75