Anisotropic cosmic ray diffusion in isotropic Kolmogorov turbulence

ABSTRACT

Understanding the time-scales for diffusive processes and their degree of anisotropy is essential for modelling cosmic ray transport in turbulent magnetic fields. We show that the diffusion time-scales are isotropic over a large range of energy and turbulence levels, notwithstanding the high degree of anisotropy exhibited by the components of the diffusion tensor for cases with an ordered magnetic field component. The predictive power of the classical scattering relation as a description for the relation between the parallel and perpendicular diffusion coefficients is discussed and compared to numerical simulations. Very good agreement for a large parameter space is found, transforming classical scattering relation predictions into a computational prescription for the perpendicular component. We discuss and compare these findings, in particular, the time-scales to become diffusive with the time-scales that particles reside in astronomical environments, the so-called escape time-scales. The results show that, especially at high energies, the escape times obtained from diffusion coefficients may exceed the time-scales required for diffusion. In these cases, the escape time cannot be determined by the diffusion coefficients.

magnetic fields, scattering, turbulence, software: simulations, cosmic rays

1 INTRODUCTION

Scattering by magnetic field inhomogeneities is a fundamental process in cosmic ray transport, whether it be dominated by diffusion (Strong & Moskalenko 1998; Evoli et al. 2008; Kissmann 2014; Merten et al. 2017) or includes a component of self-regulated streaming along the background magnetic field (Kulsrud & Pearce 1969).

Spatial diffusion is described via the tensor

$$\begin{eqnarray*} \boldsymbol {\kappa } = {\begin{pmatrix}\kappa _\perp &\quad \kappa _A &\quad 0 \\ -\kappa _A &\quad \kappa _\perp &\quad 0 \\ 0 &\quad 0 &\quad \kappa _\parallel \end{pmatrix}} \quad \end{eqnarray*}$$

(1)

within a magnetic field B_tot = B + b, whereb is the turbulent component and the z-axis of the coordinate system is aligned with the background magnetic field B. Here, we define the perpendicular and parallel diffusion with respect to the mean magnetic field, hence κ_zz = κ_∥ and κ_xx = κ_yy = κ_⊥ for B = Be_z. Spatial diffusion enters the Parker transport equation

$$\begin{eqnarray*} \frac{\partial n}{\partial t} + {\bf u}\cdot \nabla n &=& \nabla \cdot (\hat{\kappa }\nabla n) + \frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2\kappa _{pp}\frac{\partial n}{\partial p}\right) \nonumber \\ &&+\, \frac{p}{3}(\nabla \cdot {\bf u})\frac{\partial n}{\partial p} + S . \end{eqnarray*}$$

(2)

Here, u denotes advection speed, p is momentum, κ_pp is the scalar momentum diffusion coefficient, and n is the phase-space cosmic ray density averaged over the direction of particle momentum. In most systems, the transport is anisotropic for weak turbulence levels b ≲ B, where the parallel diffusion coefficient κ_∥ is larger than the perpendicular diffusion coefficient κ_⊥. Due to the geometry of many pertinent astrophysical objects, however, the perpendicular component often plays a decisive role in the escape times τ_esc of particles from the system. Examples include those objects where the perpendicular spatial structures have distinctly smaller scales than the structures parallel to the ordered background field. For instance, in flux tubes in the lobes of radio galaxies (Bell et al. 2019) or in elongated jets of active galactic nuclei (AGN), the large-scale field structure is believed to be helical. Thus, the transport perpendicular to the mean-field of the jet is associated with shorter escape times than those associated with particles leaving along the magnetic field lines of the jet. The orientation of AGN jets with respect to the observer determines the escape process relevant for observation, and objects with all jet orientations, even precessing jets, are subject of investigation. The viewing angle toward the Earth, then, determines what signatures can be seen (see e.g. de Bruijn et al. 2020). Blazars are AGN seen along the jet axis, where parallel diffusion is relevant, while the perpendicular escape process can dominate over the parallel one in inclined blazars or Fanaroff–Riley galaxies (see e.g. Becker, Biermann & Rhode 2005; Merten et al. 2021; Tavecchio 2021).

There also exist regions in the Milky Way for which the escape time of charged particles perpendicular to the Galactic plane is more efficient than the time of escape along the radial direction, despite the much smaller perpendicular diffusion, due to the small-scale height of the Galactic disc (Effenberger et al. 2012; Gaggero et al. 2015; Evoli et al. 2017; Reichherzer et al. 2022).

The often complex links between magnetic field geometry and anisotropy of cosmic ray diffusion must be taken into account when simulating cosmic ray transport, as they also affect observable quantities such as the energy spectra of cosmic ray primaries (protons, heavier nuclei) and secondaries (electrons, gamma-rays, neutrinos). As the diffusion tensor changes depending on the propagation regime (set by the energy, turbulence level, turbulence spectrum, intermittency, etc.), these dependencies enter the transport equation, and therefore, the final characteristics of the observed multi-messenger signatures. In particular, the leaky-box model of the Milky Way predicts that the cosmic ray energy spectrum observed at the Earth is steepened by diffusion. In this model, it is assumed that the particle transport is dominated by scalar diffusion and that the system is in steady-state, i.e. ∂n/∂t ≈ 0. Assuming propagation in a fixed scale height, the spectrum is given by the ratio of the source spectrum Q(E) ∝ E^−α and the diffusion coefficient |$\kappa _i(E)\propto E^{\gamma _i}$|⁠, i.e. |$N(E)\propto E^{-\alpha -\gamma _i}$| (Jokipii 1966; Berezinskii et al. 1990; Becker Tjus & Merten 2020), where i represents the perpendicular or the parallel component, depending on the geometry of the system and the dominating component of the diffusion tensor. If we better understand the relationships between the perpendicular and parallel components of the diffusion tensor, observations of one of these components can be used to determine the other one.

This paper is organized as follows. A discussion of the theoretical background of the relation between the diagonal diffusion tensor components is presented in Section 2. Section 3 covers the simulation set-up for diffusion coefficient and time-scale calculations. In Section 4, the parameter space in energy and turbulence levels is examined with respect to the validity of the classical scattering relation (CSR), finding that it provides a good description of numerical results. The paper is concluded by a discussion of the results in the context of perpendicular diffusion and its consequences for the escape time-scales of charged particles.

2 COMPONENTS OF THE DIFFUSION TENSOR

The understanding of purely parallel cosmic ray transport has been advanced significantly over time (Jokipii 1966; Jokipii & Parker 1969; Giacalone & Jokipii 1999; Casse, Lemoine & Pelletier 2002; Shalchi 2009; Buffie, Heesen & Shalchi 2013; Snodin et al. 2016; Reichherzer et al. 2020; Deligny 2021), but, there still are several open questions, such as what turbulence levels allow for the application of quasi-linear theory (QLT) and how the turbulence spectrum and intermittency influence particle transport (Shukurov et al. 2017; Friedrich et al. 2018). Several studies were dedicated to perpendicular diffusion coefficients and their relationship to the parallel diffusion coefficient (Giacalone & Jokipii 1999; Mace, Matthaeus & Bieber 2000; Casse et al. 2002; Matthaeus et al. 2003; Candia & Roulet 2004; Marco, Blasi & Stanev 2007; Minnie et al. 2007; Fatuzzo et al. 2010; Plotnikov, Pelletier & Lemoine 2011; Harari, Mollerach & Roulet 2014, 2015; Snodin et al. 2016; Subedi et al. 2017; Giacinti, Kachelriess & Semikoz 2018; Dundovic et al. 2020; Reichherzer et al. 2020). Knowledge of the relationships between the perpendicular and parallel components of the diffusion tensor are important in determining one component with the knowledge of the other. When measuring diffusion coefficients of cosmic rays in e.g. galaxies (Heesen 2021), the orientation of the galaxy plays an important role, because geometrical arguments often allow the measurement of only one of the components of the diffusion tensor.

Analytical theories must describe the components of the diffusion coefficients over a wide range of turbulence levels and energies. A strong turbulence level (b/B ≫ 1), for example, results in equal perpendicular and parallel components of the diffusion coefficient, as the charged particles propagate through isotropic turbulence (Jokipii & Parker 1969; Bieber & Matthaeus 1997; Giacalone & Jokipii 1999).

Test particle simulations for two-dimensional, slab, or composite (two-dimensional & slab) turbulence can be adequately described using analytical models such as QLT (in some domains), non-linear guiding centre theory, and unified non-linear theory (see e.g. Shalchi 2020 for a review). Despite the successes of these turbulence models, for isotropic three-dimensional turbulence, tension with these theories was found at small reduced rigidities (Dundovic et al. 2020; Reichherzer et al. 2020). So far, no encompassing theory exists capturing the relation between the parallel and perpendicular components that agrees with simulation results over the whole range of turbulence levels b/B for isotropic three-dimensional turbulence.

Instead of considering diffusion models, which, due to their underlying assumptions, only apply to certain parameter ranges or exhibit significant limitations (see e.g. Mertsch 2020), one may employ a more general, non-linear theory, such as the Bieber and Matthaeus (BAM) theory (Bieber & Matthaeus 1997). Therein, fluctuations cause particles to deviate from the ideal helices assumed in QLT, because of continuous change of the pitch angle. Therefore, particle velocities cannot be treated as being correlated over the complete trajectory, as assumed in QLT. Instead, the general, physically motivated assumption of exponentially decaying velocity correlations applies. Consequently, the velocity correlations that are crucial for determining the running diffusion coefficients |$\kappa _{ii}(t) = \langle x_i^2\rangle /(2t)$| read

$$\begin{eqnarray*} \left\langle v_\perp (t)v_\perp (0) \right\rangle &=& \frac{v^2}{3}\, \cos (\Omega t)\, \mathrm{e}^{-t/\tau _\perp }, \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} \left\langle v_\parallel (t)v_\parallel (0) \right\rangle &=& \frac{v^2}{3}\, \mathrm{e}^{-t/\tau _\parallel }, \end{eqnarray*}$$

(4)

with 〈·〉 being the ensemble average, while Ω = v/r_g denotes the angular frequency of the unperturbed orbit, and the effective time-scales for the decorrelation of the trajectories along and perpendicular to the background field are written as τ_∥ and τ_⊥, respectively. Note that the gyroradius r_g of the unperturbed trajectory is calculated only with respect to the ordered background field B. In contrast to the BAM model, we make no assumptions regarding the time-scales for perpendicular diffusion and, in particular, do not use the field-line random walk (FLRW) coefficient as a measure of the perpendicular mean-free path. Avoiding restrictive assumptions is important, since the description of perpendicular transport as FLRW is based on large parallel mean-free paths and requires the absence of resonant scattering, neither of which criterion is fulfilled in all key astrophysical environments for the particle energies considered in the present work (the effects of FLRW (Sonsrettee et al. 2016; Shalchi 2021) are unlikely to be completely negligible, however, a point to which we return in the discussion of Fig. 1). Gyroresonant scattering is defined in QLT as cos Θ = l/(2πr_g), where Θ denotes the pitch angle, l is the wavelength of a turbulent fluctuation. In an approach to generalize this criteria for strong turbulence levels, we employ |$\cos {\Theta } = l / (2\pi \rho \, l_\mathrm{c})$| for defining the different transport regimes, where the unitless reduced rigidity yields ρ = r_g/l_c · B/B_tot = E/(qcB_totl_c). Note that this expression coincides with the QLT expression for b ≪ B, where B ≈ B_tot holds. As −1 ≤ cos Θ ≤ 1, the maximum ρ for which gyroresonance occur is ρ = 5/2π. This is the upper boundary of the resonant-scattering regime (RSR). For larger particle energies, and thus larger gyroradii, the resonant-scattering criterion is only sustained for a fraction of possible pitch angles. This fraction decreases with increasing gyroradius. Beyond the RSR, the system enters the transition regime (TR), followed by the quasi-ballistic regime (QBR) as soon as the ratio 1/ρ becomes negligible. The energy dependence of the diffusion coefficients differs significantly in these regimes (Reichherzer et al. 2022). The description of perpendicular transport in the RSR also has to account for the cross-field diffusion caused by resonant scattering (see Desiati & Zweibel 2014, for a discussion).

Figure 1.

Comparison between simulation results and theoretical predictions from the CSR for κ_⊥/κ_∥ for all investigated simulation parameter combinations of the turbulence level b/B and the energy E. The colour-code indicates the deviation according to equation (11). The theoretical predictions underestimate the ratios κ_⊥/κ_∥, as shown in Figs 2 and 3. Each point comprises the results of 20–50 test particle simulations for the indicated parameter combination as described in Section 3. The magnetic field strength is |$B=1\, \mu$|G and the maximal fluctuation wavelengths are l_max ≈ 82 pc. Using scaling laws, this plot can be applied to other astrophysical environments, e.g. to AGN jet plasmoids, as illustrated with the upper x-axis labels and described in the text.

The diffusion coefficients can be computed according to the Taylor–Green–Kubo (TGK) formalism (Shalchi 2009),

$$\begin{eqnarray*} \kappa _\perp &=& \lim \nolimits _{t \rightarrow \infty } \int \nolimits _0^t \mathrm{d}\tau \, \frac{v^2}{3}\cdot \cos (\Omega \tau)\cdot e^{-\tau /\tau _\perp } = \frac{v\, r_\mathrm{g}}{3}\frac{\Omega \tau _\perp }{1+\Omega ^2\tau _\perp ^2},\nonumber\\ \kappa _\parallel &=& \lim \nolimits _{t \rightarrow \infty } \int \nolimits _0^t \mathrm{d}\tau \, \frac{v^2}{3}\cdot e^{-\tau /\tau _\parallel } = \frac{v^2\tau _\parallel }{3} = \frac{v\, \lambda _\parallel }{3}. \end{eqnarray*}$$

(5)

The relation between the parallel mean-free path and the parallel diffusion coefficient κ_∥ = vλ_∥/3 establishes the relation τ_∥ = λ_∥/v, which defines τ_∥ as the diffusion time-scale in parallel direction. Similarly, we define τ_⊥ as the diffusion time-scale in the perpendicular direction. The relation between the perpendicular and the parallel components of the spatial diffusion coefficient thus reads

$$\begin{eqnarray*} \frac{\kappa _\perp }{\kappa _\parallel }= \frac{\tau _\perp }{\tau _\parallel } \frac{\Omega \, r_\mathrm{g}}{v}\frac{1}{1+ \Omega ^2\, \tau _\perp ^2} = \frac{\tau _\perp }{\tau _\parallel } \frac{1}{1+ v^2 \tau _\perp ^2/r_\mathrm{g}^2}\, . \end{eqnarray*}$$

(6)

Note that there are different assumptions made by several studies about the time-scales τ_⊥ and τ_∥ to further exploit equation (6). In the hard-sphere model introduced by Gleeson (1969), which is also known as the CSR within standard kinetic theory, only one isotropic scattering process is present. The assumption that these two time-scales coincide as the decorrelations are caused by a single scattering process is also present in Isenberg & Jokipii (1979), Balescu, Wang & Misguich (1994), Casse et al. (2002). Without knowledge of the exact contribution of different processes to the decay of the velocity correlation functions, we also employ τ_∥ ≈ τ_⊥. We justify this assumption by the isotropic character of the FLRW in isotropic three-dimensional turbulence. Consequently, the diffusion time-scales for this process are isotropic. We remark that the influence of resonant scattering, which dominates for lower energies, may contradict this assumption. For equal diffusion time-scales, the CSR prediction yields

$$\begin{eqnarray*} \frac{\kappa _\perp }{\kappa _\parallel }= \frac{1}{1+\lambda _\parallel ^2 /r_\mathrm{g}^2}\, . \end{eqnarray*}$$

(7)

Consequently, whenever the time-scales in the perpendicular and parallel directions are identical, the perpendicular diffusion can be expressed as

$$\begin{eqnarray*} \kappa _\perp = \frac{\kappa _\parallel }{1+(3\kappa _\parallel)^2/(c\, r_\mathrm{g})^2}. \end{eqnarray*}$$

(8)

Note that for the highly relativistic particles (Lorentz factor γ ≫ 1) considered throughout the study, we approximate the velocity by the speed of light. Previous investigations (Giacalone & Jokipii 1999) have already shown that equation (8) will not hold everywhere. In this paper, we will systematically investigate the parameter space spanned by particle energy and turbulence level to deduce in which this expression and the underlying assumption τ_⊥ = τ_∥ hold.

3 NUMERICAL SIMULATION SET-UP

Test-particle propagation of highly relativistic protons is simulated with the software crpropa,¹ a publicly available tool for simulations of cosmic ray transport (Alves Batista et al. 2016; Merten et al. 2017; Alves Batista et al. 2021). We performed simulations of highly relativistic protons of at least several PeV, but many conclusions apply to the general case as well, because only the ratio of gyroradius to the correlation length l_c ≈ l_max/5 is relevant for the diffusion-coefficient scaling, which extends over more than two orders of magnitude. Specifically, reducing both B_tot and the particle energy by the same factor results in the same diffusion tensor.

Magnetohydrodynamic turbulence can be simulated numerically, where the spectrum of turbulence is the result of a turbulent energy cascade (Stone, Ostriker & Gammie 1998; Cho, Lazarian & Vishniac 2003; Cohet & Marcowith 2016). In contrast, simplified synthetic turbulence allows for a larger inertial range (Mertsch 2020) and can be constructed in such a way that the turbulence properties exactly match the assumptions on turbulence underlying the theories investigated in this paper (Schlegel et al. 2020).

We follow the synthetic turbulence generation described in Reichherzer et al. (2022). For the isotropic three-dimensional turbulent component, we employ a Kolmogorov-like spectrum G(k) ∝ (k/k_min)^−5/3 for k_min ≤ k ≤ k_max, with k_min = 2π/l_max and k_max = 2π/l_min, and G(k) = 0 otherwise. We choose l_min = 1.7 pc and l_max ≈ 82 pc to span a large inertial range of the turbulence and we store the turbulent field on a grid with 1024³ grid points and a spacing of l_min/2. We note that this adopted inertial-range Kolmogorov scaling for the turbulence is an idealization, as interstellar turbulence is driven on many scales, from as much as 100 pc by superbubbles to kinetic scales by cosmic rays themselves (Beresnyak & Lazarian 2015). Moreover, the Alfvénic turbulence assumed here, unlike hydrodynamic Kolmogorov turbulence, is known to be highly anisotropic, even when considering compressibility effects that can generate an isotropic component (Beresnyak & Lazarian 2015). Note that anisotropy causes inefficient scattering, leading to increased parallel diffusion coefficients (Chandran 2000; Yan & Lazarian 2002).

Our simulations of particle transport in a combined turbulent field b and mean field |$B = 1\, \mu$|G are based on the well-established TGK formalism (Shalchi 2009), (see e.g. Giacalone & Jokipii (1999), Casse et al. (2002), Globus, Allard & Parizot (2008), Marco et al. (2007), Minnie et al. (2007), Snodin et al. (2016), Giacinti et al. (2018), Dundovic et al. (2020), Reichherzer et al. (2022)). We inject 2000 particles at time t = 0 and position x_i = δ(x_{i, 0}), and integrate the particle trajectories in the magnetic field by numerically solving the Lorentz–Newton equation via the energy-conserving Boris-push method (Boris & Shanny 1972). The step-size of the particles in our integrator is s_step = min(r_g/5, l_max/20) in order to sufficiently resolve the gyrations and the turbulent fluctuations.

The running diffusion coefficients are determined according to |$\kappa _{ii}(t) = \langle \Delta x_{i}^{2}\rangle /(2\, t\,)$| by averaging over the particles. After κ_ii(t) converges to a constant value, we calculate the diffusion coefficient κ_ii by averaging the final running diffusion coefficients over time to reduce statistical fluctuations. In order to ensure representative results for isotropic turbulence, we perform 20–50 simulations for each parameter combination, which are only distinguished by randomly changing phases of the turbulent fluctuations with the same statistical and spectral properties.

We perform simulations only for particles in the RSR, TR, and QBR. The latter is characterized by good agreement between numerical simulation results and theoretical predictions. For example, an energy scaling of κ ∝ E² is obtained here, independent of the turbulence level. This is to be expected since diffusion is no longer dominated by resonant scattering. However, the situation is different for the RSR, where resonant scattering becomes dominant according to the resonance criterion.

Further reducing of the gyroradii results into two effects:

• The magnitude of the fluctuations decreases according to the relation |$b \propto (G(k)\, k)^{1/2} \propto k^{-1/3}$|⁠, given by the Kolmogorov-type spectrum for smaller fluctuation wavelengths. As the resonance criterion connects the gyroradii with the fluctuations wavelengths, reducing particle gyroradii results in decreasing magnitude of fluctuations that are relevant for resonant scatterings for these particles. In addition to the decrease in the magnitude of the fluctuations, the scales of the changes of the turbulent magnetic field are considerably larger for small gyration radii, so that an effective directed magnetic field is established on the relevant scales. Thus, an effective b/B ≪ 1 exists on small scales.

• According to the resonant-scattering criterion and the limited resolution of the turbulence (l_min > spacing/2), the range of pitch angles capable of resonant interactions decreases. As the needed RAM to store the turbulence on the grid is very large, the resolution is severely constrained and limits the range of fluctuations to a few orders of magnitude. The missing resonant-interaction possibilities due to the constrained fluctuation wavelength range lead to a numerically increased diffusion coefficient.

Whereas effect (i) establishes the prerequisites for weak turbulence, missing resonant-scattering possibilities must be taken into account due to effect (ii). Since we are both interested in the influence of the turbulence level and want to keep the numerical effort tractable, we omit the low energies. This also justifies neglecting the effect of cosmic ray self-confinement, which is deemed ineffective at energies above ∼100–200 GeV (Zweibel 2017).

4 RELATION BETWEEN PERPENDICULAR AND PARALLEL COMPONENTS

In this section, we carry out a comprehensive parameter study to substantially expand on the results obtained in the different investigations (Jokipii & Parker 1969; Giacalone & Jokipii 1999) and extract relevant insights on the validity of the application of the equations (7) and (8). For this purpose, a parameter space over three orders of magnitude in particle energy and two orders of magnitude in turbulence level is examined; thus, covering different energy regimes, such as the RSR, QBR, and TR.

The calculation of the diffusion coefficients via the TGK mechanism described in Section 3 (see equation (5)), using the temporal convergence of the running diffusion coefficient, facilitates the determination of the parallel mean-free paths. Therefore, instead of determining λ_∥ directly, we determine it by means of equation (5) via the diffusion coefficient κ_∥ as

$$\begin{eqnarray*} \lambda _\parallel = \frac{3\kappa _\parallel }{c}. \end{eqnarray*}$$

(9)

This approach simplifies the analytical prediction for the ratio of the diffusion coefficient components in the CSR, see equation (7), to

$$\begin{eqnarray*} \left(\frac{\kappa _\perp }{\kappa _\parallel }\right)_\mathrm{CSR} = \frac{1}{1+\lambda _{\parallel , \mathrm{sim}}^2 /r_\mathrm{g}^2} = \frac{1}{1+(3\kappa _{\parallel , \mathrm{sim}})^2/(c\, r_\mathrm{g})^2}. \end{eqnarray*}$$

(10)

The deviation between simulation results (labelled sim) and the CSR prediction (labelled CSR) for the ratio κ_⊥/κ_∥ yields

$$\begin{eqnarray*} \varepsilon = \left| \left(\frac{\kappa _\perp }{\kappa _\parallel } \right)_\mathrm{sim} - \left(\frac{\kappa _\perp }{\kappa _\parallel } \right)_\mathrm{CSR} \right| \cdot \left(\frac{\kappa _\perp }{\kappa _\parallel } \right)_\mathrm{CSR}^{-1} \end{eqnarray*}$$

(11)

$$\begin{eqnarray*} \hphantom{\varepsilon} = \left|\left(\frac{\kappa _\perp }{\kappa _\parallel } \right)_\mathrm{sim} - \frac{1}{1+(3\kappa _{\parallel ,\mathrm{sim}})^2/(c\, r_\mathrm{g})^2}\right| \cdot \left(1 + \frac{(3\kappa _{\parallel ,\mathrm{sim}})^2}{(c r_\mathrm{g})^2}\right). \end{eqnarray*}$$

(12)

Fig. 1 presents this deviation ε with different colours between simulated and CSR-predicted ratios of the perpendicular and parallel diffusion coefficients. The displayed data points represent all parameter combinations considered in the present study.

Since only the ratio r_g/l_c is relevant for the diffusion coefficients, this plot is universally applicable as indicated by the axis labels below and above the figure. The labels of the lower x-axis of Fig. 1 correspond precisely to the simulation parameters. With an arbitrary factor, the range of particle energy can be scaled if, at the same time, the product of magnetic field strength and the correlation length is scaled in the same way. For illustration, the upper x-axis labels are obtained with a magnetic field 10⁶ times larger than that on the lower x-axis labels. At the same time, the correlation length is scaled down by a factor of 5 · 10⁷. Effectively, therefore, the energy range has decreased by a factor of 50 and represents the particle energies in the set-up of an exemplary AGN jet plasmoid (Hoerbe et al. 2020; Reichherzer et al. 2021). Accordingly, the CSR relation describes the transport of PeV protons, which are a relevant contributor to multi-messenger processes in hadronic blazar models (Böttcher 2019).

As no encompassing theory currently exists capturing the relation between the parallel and perpendicular components that agrees with simulation results over the whole range of turbulence levels b/B and energies for isotropic three-dimensional turbulence (see Dundovic et al. 2020), it is not surprising that the CSR fails for some parameter combinations. Fig. 1 demonstrates only poor agreement (deviations between |$10^2{{\ \rm per\ cent}} \lesssim \varepsilon \lesssim 10^3{{\ \rm per\ cent}}$|⁠) between CSR and numerical simulations for small energies and weak turbulence levels. Particles in the zone of poor agreement satisfy r_g ≲ l_max ≲ λ_∥, which suggests that although these particles are scattered by resonant waves, a significant part of the perpendicular displacement is due to FLRW, a point to which we return in the discussion of Figs 2 and 4. A precise characterization of FLRW for our turbulence model is currently in progress and will be described in future publication. Fig. 1 shows that a large parameter space is well-described by this simple theoretical prediction from the CSR. As we have shown here, an all-encompassing theory must also describe the anisotropic diffusion time-scales in the RSR for weak turbulence levels, with otherwise isotropic diffusion time-scales even for strongly anisotropic diffusion coefficients at high energies. Note that in many astrophysical environments (e.g. galaxies; Jansson & Farrar 2012; Kleimann et al. 2019; Shukurov et al. 2019), similar orders of magnitude of b and B are present at the injection scales of the turbulence. However, scaling of the turbulent spectra according to Kolmogorov or Kraichnan leads to b ≪ B on the smaller fluctuation scales relevant to lower-energy particles. In the following, we discuss the reasons for the good agreement between simulations and CSR predictions toward high energies (see Section 4.1) and toward strong turbulence levels (see Section 2) that were identified in the overview plot.

Figure 2.

Ratios of the perpendicular to the parallel diffusion coefficient from test particle simulations as functions of the ratio of the parallel mean-free paths and the gyroradii of the particles (top), defined in equation (13). Each point comprises the results of 20–50 test particle simulations and indicates, via its colour-code, the fixed particle energy used in these simulations. The simulation set-up and parameters are described in Section 3. Theoretical predictions obtained from equation (7) are shown as the grey line that is labelled theory. In the lower panel, the deviation ε between simulation and theory is shown in per cent and is determined via the definition presented in equation (11).

Figure 3.

Same data as in Fig. 2, except that the colour-coding indicates the turbulence level in each simulation. In contrast to Fig. 2, here the deviations for large λ_∥/r_g are visible in the log–log scaling. In the lower panel, the deviation ε between simulation and theory is shown in per cent and is determined via the definition presented in equation (11). Good agreement between theory and simulations for strong turbulence levels is found.

Figure 4.

Parallel mean-free path versus unperturbed gyroradii normalized to l_max, respectively. The same colour-coding from Fig. 1 indicates the deviation between numerical simulations and CSR predictions. The parameter space for which r_g ≲ l_max ≲ λ_∥ holds, is highlighted in grey. The upper x-axis provides the corresponding energies of particles in AGN jet plasmoids with the given parameters. The dotted, straight, and dashed lines represent the predicted scalings for fixed turbulence levels in the limits for QLT (RSR), Bohm (RSR), and high energies (QBR), respectively.

4.1 Range of validity in particle energy

Fig. 2 shows simulation results for the ratios of the diffusion coefficient components κ_⊥/κ_∥ as functions of the ratio between the parallel mean-free paths and the gyroradii λ_∥/r_g. Note that, from equation (9),

$$\begin{eqnarray*} \frac{\lambda _\parallel }{ r_\mathrm{g}} = \frac{\kappa _\parallel }{c\, r_\mathrm{g}}. \end{eqnarray*}$$

(13)

The particle energy used for each simulation is colour-coded. Each data point in Fig. 2 comprises the results of 20–50 test particle simulations and indicates, via its colour-code, the fixed particle energy used in these simulations. Theoretical predictions obtained by equation (7) are shown as the grey line labelled theory. The deviation ε between theory and simulation results, defined in equation (11), is presented in the lower panel.

It should be noted that there is no straightforward relationship between λ_∥/r_g and energy, since the dependence of λ_∥ on energy depends on the transport regime. For this reason, the colour of the points, according to the particle energy, is essential to visualize the energy dependence. This figure, especially in the lower panel, as well as Fig. 1, shows better agreement between theory and simulation results for higher particle energies. Particles that show large deviations (⁠|${\gtrsim}100{{\ \rm per\ cent}}$|⁠) between simulation results and CSR predictions satisfy r_g ≲ l_max ≲ λ_∥, as demonstrated in the lower panel of Fig. 2. These low-energy particles (purple) correspond to small gyroradii, while the parallel mean-free paths are large due to the weak turbulence level. Fig. 4 visualizes this claim by showing the simulation results for the parallel mean-free paths versus the unperturbed gyroradii normalized to l_max, respectively. The colour-coding provides the deviation between QBR and simulation results according to equation (11). The grey area shows the parameter space spanned by r_g ≲ l_max ≲ λ_∥. The large deviations lie primarily in this parameter space. The scaling with energy and r_g follows the predictions from the high-energy theory proposed in Subedi et al. (2017) as discussed in more detail in Reichherzer et al. (2022).

Other definitive general statements based on these results alone are not possible because the turbulence level varies in addition to the particle energy. This will be discussed in more detail in the following subsection, where we investigate the agreement between theory and simulations for fixed turbulence levels, and thus, separate the influence of the turbulence level and the energy.

4.2 Range of applicable turbulence levels

In order to simplify the investigation at λ_∥/r_g ≫ 1, where deviations between theory and simulations are apparent in the lower panel of Fig. 2, the same data are shown in log–log representation in Fig. 3. To enable investigation of the turbulence-level dependence of the deviation ε between theory and simulation, the simulation points are colour-coded according to b/B.

The differences between the diagonal elements of the diffusion tensor diminish with increasing b/B, so that κ_⊥/κ_∥ ≈ 1 in the limit b ≫ B. Here, the ratios λ_∥/r_g are expected to be small, as λ_∥ decreases with increasing turbulence level. Since the CSR predicts exactly this, the agreement with the simulation results for this parameter range, which is characterized by strong turbulence levels, is good.

To further investigate the turbulence-level dependence of the agreement between theory and simulations, we parametrize the CSR as

$$\begin{eqnarray*} \frac{\kappa _\perp }{\kappa _\parallel } = \frac{a_1}{1+a_2 (\lambda _\parallel /r_\mathrm{g})^{a_3}} \stackrel{\lambda _\parallel \gg r_\mathrm{g}}{\approx } \frac{a_1}{a_2}\left(\frac{\lambda _\parallel }{r_\mathrm{g}}\right)^{-a_3}, \end{eqnarray*}$$

(14)

where the latter expression yields as a good approximation in the limit of large λ_∥, and fit the simulation data for turbulence levels. The best-fitting results for the parameters a₁, a₂, and a₃ for all turbulence levels presented in this study are shown as functions of b/B in Fig. 5. The CSR predictions are a₁ = 1, a₂ = 1, and a₃ = 2 as indicated by the horizontal solid lines in the plot. As it is apparent in Fig. 3, the approximation presented in equation (14) for large λ_∥ is used for the fits for weak turbulence levels, where all simulation points approximately obey |$\kappa _\perp /\kappa _\parallel \propto (\lambda _\parallel / r_\mathrm{g})^{a_3}$|⁠. Theoretically, this criterion is justified by the fact that the parallel diffusion coefficients scale with (b/B)^α (with α = −2 for QLT). This scaling is directly transferable to the parallel mean-free paths, since λ_∥ = 3κ_∥/c. Weaker turbulence levels, thus, lead to larger λ_∥. Tests have shown the value of 2 for the turbulence level to be practical as a switch between both expressions presented in equation (14).

Figure 5.

Fit parameters a₃ and the ratio a₁/a₂ as functions of the turbulence level b/B. For strong turbulence levels, where λ_∥ ≫ r_g holds, we fit the approximation in equation (14). However, for unified representation, we also show only the ratios of a₁ and a₂ for weaker turbulence levels, although individual values of the parameters are available. The dash–dotted horizontal line indicates the prediction of the CSR for the ratio a₁/a₂, whereas the solid grey horizontal line illustrates the theoretical prediction of the CSR for a₃ according to equation (5). The fits were performed to the subset of simulations shown in Fig. 3 that match the turbulence level b/B to which the points correspond in the present plot.

For strong turbulence levels, the fits agree well with these CSR predictions. A difficulty arises at weak turbulence levels, where, due to the dependence of the parallel diffusion on the turbulence level, simulation data only exist for λ_∥ ≫ r_g. This leads to large uncertainties in the fits. As we commented below equation (11), in these cases, there may also be a contribution from FLRW.

Furthermore, we examine the predictive power of the CSR for the perpendicular component computed from κ_∥ according to the expression from equation (8). The perpendicular diffusion coefficients (circles) obtained from the CSR description and the numerical simulation results for κ_∥, are directly compared with κ_⊥ from the numerical simulations. For this purpose, Fig. 6 shows the perpendicular diffusion coefficients as functions of the reduced rigidity. Note that here the gyration radius is calculated with respect to the total magnetic field and differs from the definition of the gyration radius of the unperturbed orbit r_g, since the latter depends only on the ordered magnetic field component B.

Figure 6.

Perpendicular diffusion coefficient as functions of reduced rigidity ρ. The markers show the association of the data points of the numerical simulations in the respective transport regimes. The circles indicate the values expected from the CSR for the perpendicular diffusion coefficients. These are calculated using the numerical simulation results for κ_∥ according to equation (8). The lower panel shows the relative deviation according to equation (15) between simulations and RSR predictions.

The energy-independent perpendicular diffusion coefficient characteristic of FLRW diffusion in QBR is apparent in Fig. 6. The lower panel provides the deviations between our numerical simulations and the CSR predictions

$$\begin{eqnarray*} \xi = \left| \frac{\kappa _{\perp ,\mathrm{sim}} - \kappa _{\perp ,\mathrm{CSR}}}{\kappa _{\perp ,\mathrm{CSR}}}\right| = \frac{\kappa _{\perp ,\mathrm{sim}}}{\kappa _{\parallel ,\mathrm{sim}}} \left(1 + \frac{9 \kappa _{\parallel ,\mathrm{sim}}^2}{c^2\, r_\mathrm{g}^2} \right) . \end{eqnarray*}$$

(15)

As follows from the discussion above, there is poor predictive power for the perpendicular diffusion coefficients in the RSR for weak turbulence levels. In contrast, there is good agreement in QBR that improves with energy, where the parallel mean-free path increases as λ_∥ ∝ E² and suppresses perpendicular diffusion according to

$$\begin{eqnarray*} \kappa _\perp = \frac{\kappa _\parallel }{1+\lambda _\parallel ^2 /r_\mathrm{g}^2} \stackrel{\lambda _\parallel \gg r_\mathrm{g}}{\approx } \frac{c\, r_g^2}{3 \lambda _\parallel } ~ \stackrel{\mathrm{QBR}}{\propto }~ \rho ^0. \end{eqnarray*}$$

(16)

The λ_∥ ≫ r_g approximation holds for a large parameter space according to Fig. 4.

Note that a related set-up involves the spatial diffusion of beam-injected ions and alpha particles in fusion plasmas, where energy-independent diffusion due to magnetic fluctuations is predicted by decorrelation theory (Hauff et al. 2009; Pueschel 2012; Pueschel et al. 2012). Future work, accounting for pitch-angle-dependence, will assess to which degree this theory can be seen as a (non-relativistic) extension of the present results.

4.3 Implications for time-scales of diffusion

In the previous subsections, we found that the CSR describes the diffusion coefficients accurately for a wide range of parameters. We emphasize that this covers important astrophysical environments involving relativistic plasmoids in blazar jets. In Figs 1 and 4, we have shown that the transport of PeV protons can be described well across turbulence levels using the CSR theory. Furthermore, this finding is complemented by our results on the isotropy of the diffusion time-scales. We have shown that the diffusion time-scales coincide in both parallel and perpendicular directions over a very large parameter space, even if the diffusion coefficients are anisotropic (κ_⊥ ≪ κ_∥) due to strong coherent background fields. Two effects cause the isotropy in the diffusion time-scales:

• Increasing the turbulence level isotropizes the diffusion process and, therefore also the time-scales needed to correctly describe transport diffusively in the parallel and perpendicular directions. This manifests itself in the fact that the CSR is able to describe the data more accurately with increasing turbulence level, independent of energy, as this theory applies τ_⊥ ≈ τ_∥.

• Isotropization of the diffusion time-scales also arises for increased particle energy, especially for r_g ≳ l_c, as the contribution from resonant scattering decreases. It is important to emphasize that equal parallel and perpendicular time-scales for diffusion do not lead to equal diffusion coefficient components for this case. The finding of isotropic diffusion time-scales in combination with energy-independent perpendicular diffusion indicates that high-energy particles follow field lines, a behaviour commonly seen for low-energy particles with small gyroradii. The identification of field-line random walk as the cause for diffusion of high-energy particles is in agreement with Minnie et al. (2009).

Furthermore, we want to emphasize that our comparisons reveal poor agreement of the numerical simulation results with the predictions of the CSR for the parameter range r_g ≲ l_max ≲ λ_∥. The reasons for this can be manifold:

• The isotropization of diffusion time-scales assumed to derive equation (7) may fail in RSR. In the RSR, resonant scattering dominates the diffusion behaviour. It is known that resonant scattering depends on the pitch angle and consequently on the orientation of the ordered background magnetic field. This may lead to different diffusion time-scales, along and perpendicular to the ordered magnetic field.

• The simplistic BAM model may fail to model physical processes that become relevant in the aforementioned parameter range that also describes the Galactic CR transport. Therefore, more refined theories, such as the non-linear guiding centre theory and the unified non-linear theory (see e.g. Shalchi 2020 for a review), have been developed recently. As noted before, these theories apply to numerous turbulence models, except for isotropic Kolmogorov turbulence, where deviations between predictions and numerical simulations are apparent.

When diffusion times are shorter than acceleration times and good agreement between CSR and numerical simulation results is apparent, important consequences for the escape times perpendicular to the ordered field arise. Here, it should be noted that the presented findings are only valid for acceleration times smaller than the escape times of the particles; otherwise, particles would already become diffusive during acceleration. Assuming transport to be in the diffusive limit, the escape time of charged particles to reach distance d yields

$$\begin{eqnarray*} \tau _{\mathrm{esc},\perp } \approx \frac{d^2}{2\kappa _\perp } \stackrel{\lambda _\parallel \gg r_\mathrm{g}}{\approx } \frac{3d^2\lambda _\parallel }{2cr_\mathrm{g}^2}. \end{eqnarray*}$$

(17)

Here, we use the λ_∥ ≫ r_g approximation, which, according to Fig. 4, applies to a large parameter space spanning all transport regimes and encompassing many turbulence levels in the RSR as well. In QBR, with |$\lambda _\parallel \propto r_\mathrm{g}^2$| for high-energy particles, the scaling of the perpendicular escape time is again energy-independent. Assuming equal diffusion time-scales yields

$$\begin{eqnarray*} \tau _{\perp } \approx \tau _{\parallel } \approx \frac{\lambda _\parallel }{c}. \end{eqnarray*}$$

(18)

For l_c ≲ r_g, the diffusive assumption used to derive the perpendicular escape time can be seen to break down when computing the ratio of the time-scales for diffusion and escape

$$\begin{eqnarray*} \frac{\tau _{\perp }}{\tau _{\mathrm{esc},\perp }} \approx \frac{2r_\mathrm{g}^2}{3d^2}. \end{eqnarray*}$$

(19)

For the diffusive approach to be applicable, this ratio must be less than 1, from which the condition |$r_\mathrm{g} \lesssim \sqrt{3/2} d$| follows. This condition for r_g should only be understood as an upper limit, since in perpendicular diffusion, there is the peculiarity that the maximum of the perpendicular running diffusion coefficient in weak turbulence is already reached in the first gyration (max(⁠|$\kappa _{\perp }(t)) = c\, r_\mathrm{g} / \pi)$| and does not coincide with the final converged κ_⊥. Using the temporal maximum of the perpendicular diffusion coefficient as a more realistic approach in equation (17) for the goal of deriving a realistic condition for the validity of the perpendicular escape time results in

$$\begin{eqnarray*} \tau _{\mathrm{esc},\perp } \approx \frac{d^2}{2\kappa _\perp } \approx \frac{d^2\, \pi }{2cr_\mathrm{g}}, \end{eqnarray*}$$

(20)

and therefore,

$$\begin{eqnarray*} \frac{\tau _{\perp }}{\tau _{\mathrm{esc},\perp }} \approx \frac{2r_\mathrm{g} \lambda _\parallel }{3\pi \, d^2}. \end{eqnarray*}$$

(21)

Interpreting requirements more strictly, the diffusive limit should only be applied under the condition

$$\begin{eqnarray*} r_\mathrm{g} \lesssim \frac{3\pi d^2}{2\lambda _\parallel }. \end{eqnarray*}$$

(22)

The practicality of this criterion is attributed to the fact that it depends only on well-studied parallel diffusion, making it easy to evaluate. Especially as high-energy particles in weak turbulence have large parallel mean-free paths λ_∥, equation (22) poses a restrictive criterion for the usage of the perpendicular escape times based on the perpendicular diffusion coefficient.

This criterion can be used to check if cosmic rays have enough time to become diffusive in astrophysical environments. The condition becomes particularly useful and, at the same time, restrictive for high energy particles (in the QBR), which are in compact structures with small values for d. In this context, examples include perpendicular diffusion in elongated jets with relatively small cross-sectional scales with respect to the extent along the jet axis.

5 CONCLUSION

Diffusive propagation and its degree of anisotropy is essential for modelling cosmic ray transport in turbulence. In this study, we examine the ratios of perpendicular to parallel diffusion coefficients in isotropic three-dimensional Kolmogorov turbulence superimposed on a uniform background field over a parameter space spanning from the resonant scattering transport regime to the quasi-ballistic regime, for weak (QLT limit) and strong turbulence levels (b/B ≳ 1). Test-particle simulations in synthetic magnetic fields and analysis, as defined by the Taylor–Green–Kubo formalism, are the basis of this study. We show that a simple analytical approach, the CSR, yields a very good description of the simulation results throughout a large portion of parameter space, which covers several transport regimes and turbulence levels. This includes high-energy particles in the quasi-ballistic regime with large mean-free paths that follow field lines and exhibit energy-independent perpendicular diffusion coefficients that are characteristic of FLRW diffusion. Where resonant scattering is dominant over diffusive transport and the influence of FLRW is subordinate, however, significant discrepancies between CSR predictions and our simulation results arise. These discrepancies are only apparent in the resonant scattering regime for weak turbulence levels. For all other parameter regimes, good agreement between CSR and simulations is found, from which the following findings emerge:

• The perpendicular diffusion coefficients follow directly from the parallel component via

  $$\begin{eqnarray*} \kappa _\perp = \frac{\kappa _\parallel }{1+(3\kappa _\parallel)^2/(c\, r_\mathrm{g})^2}. \end{eqnarray*}$$

  (23)

  The description of perpendicular diffusion via the CSR relation is beneficial, in that one may recover the difficult-to-evaluate perpendicular diffusion from the better-understood and more accessible parallel diffusion using only a few physically well-motivated assumptions.

• Furthermore, the diffusion time-scales for parallel and perpendicular diffusion are identical and can be obtained directly from the parallel diffusion coefficients via τ_⊥ ≈ τ_∥ ≈ 3κ_∥/c².

Since the denominator in equation (23) is greater than or equal to 1 for all parameter combinations (and much greater for weak turbulence levels), the parallel diffusion coefficients are greater than the perpendicular counterparts. Viewed in the context of finding (ii), equation (23) provides crucial insights into when particles may be considered diffusive. Whereas the time-scales of parallel diffusion are proportional to the parallel diffusion coefficient, the time-scales of perpendicular diffusion increase as the perpendicular diffusion coefficient decreases. In exploring this behaviour, we present a condition (see equation (22)) on the time-scale of perpendicular diffusion, permitting estimation of whether diffusion perpendicular to the directed magnetic field can occur in given astrophysical environments. We anticipate that this criterion will be especially useful in describing the propagation of high energy cosmic rays, which interact primarily with large-scale isotropic turbulence in the Galactic disc and beyond.

ACKNOWLEDGEMENTS

We acknowledge funding from the DFG, SFB1491. This work is supported by the ‘ADI 2019’ project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02 (PR). PR also acknowledges support from the German Academic Exchange Service and by the RUB Research School via the Project International funding. PR acknowledge support by Institut Pascal at Université Paris-Saclay during the Paris-Saclay Astroparticle Symposium 2021, with the support of the P2IO Laboratory of Excellence (programme ‘Investissements d’avenir’ ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), the P2I research departments of the Paris-Saclay university, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master projet UCMN and EuCAPT. LM acknowledges financial support from the Austrian Science Fund (FWF) under grant agreement number I 4144-N27. EZ gratefully acknowledges support by the US National Science Foundation through grant AST2007323. We thank Leander Schlegel for usefull discussions.

DATA AVAILABILITY

Simulations were performed with the publicly available tool crpropa (Merten et al. 2017) (the specific version used for the simulations is CRPropa 3.1-f6f818d36a64), supported by various analysis tools: numpy (van der Walt, Colbert & Varoquaux 2011), matplotlib (Hunter 2007), pandas (McKinney 2010), and jupyter notebooks (Kluyver et al. 2016).

The simulation data presented in this paper are available to interested researchers upon reasonable request.

Footnotes

1

The specific version used for the simulations is CRPropa 3.1-f6f818d36a64.