# Diffusion of relativistic charged particles and field lines in isotropic
turbulence

Marco Kuhlen<EMAIL_ADDRESS>Philipp Mertsch
<EMAIL_ADDRESS>Vo Hong Minh Phan<EMAIL_ADDRESS>aachen.de Institute for Theoretical Particle Physics and Cosmology (TTK),
RWTH Aachen University, 52056 Aachen, Germany

###### Abstract

The transport of non-thermal particles across a large-scale magnetic field in
the presence of magnetised turbulence has been a long-standing issue in high-
energy astrophysics. Of particular interest is the dependence of the parallel
and perpendicular mean free paths $\lambda_{\parallel}$ and $\lambda_{\perp}$
on rigidity $\mathcal{R}$. We have revisited this important issue with a view
to applications from the transport of Galactic cosmic rays to acceleration at
astrophysical shocks. We have run test particle simulations of cosmic ray
transport in synthetic, isotropic Kolmogorov turbulence at unprecedentedly low
reduced rigidities $r_{\text{g}}/L_{\text{c}}\simeq 10^{-4}$, corresponding to
$\mathcal{R}\simeq 10\,\text{TV}$ for a turbulent magnetic field of
$B_{\text{rms}}=4\,\mu\text{G}$ and correlation length
$L_{\text{c}}=30\,\text{pc}$. Extracting the (asymptotic) parallel and
perpendicular mean free paths $\lambda_{\parallel}$ and $\lambda_{\perp}$, we
have found $\lambda_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ as
expected for a Kolmogorov turbulence spectrum. In contrast, $\lambda_{\perp}$
has a faster dependence on $r_{\text{g}}/L_{\text{c}}$ for $10^{-2}\lesssim
r_{\text{g}}/L_{\text{c}}\lesssim 1$, but for $r_{\text{g}}/L_{\text{c}}\ll
10^{-2}$, also $\lambda_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$. Our
results have important implications for the transport of Galactic cosmic rays.

††preprint: TTK-22-35

##  I Introduction

There is overwhelming observational evidence that the transport of high-
energy, non-thermal particles in the Galaxy is predominantly diffusive: For
Galactic cosmic rays (CRs), the grammage, that is the inferred amount of
matter that CRs experience on average, implies that CRs cross the gaseous disk
of the Galaxy many times between acceleration and observation. The small
anisotropies in the observed arrival directions of CRs also point to an
efficient randomisation of CR directions. (See, e.g. Ref. Gabici  _et al._
(2019) for a review on Galactic CRs.) For collisionless plasmas like CRs, the
interactions with turbulent magnetic fields act as the sole agent of such
randomisation.

The first attempt at clarifying the interaction of waves and particles was the
so-called Quasi-Linear Theory (QLT), a perturbative approach that considers a
small turbulent magnetic field $\bm{\delta B}$ on top of the regular
background field $\bm{B}_{0}$ Jokipii (1966); Kennel and Engelmann (1966);
Hall and Sturrock (1967); Hasselmann and Wibberenz (1970). The scattering rate
is estimated by approximating the force experienced by particles on
unperturbed trajectories, that is solutions of the equation of motion in the
background field $\bm{B}_{0}$ only. Magnetic turbulence is oftentimes modelled
as a spectrum of different waves. It can be shown that the interactions
between the turbulent field and the charged particles are gyro-resonant,
meaning that particles with a gyroradius $r_{\text{g}}$ are scattered by modes
with a wavenumber $k\sim 1/r_{\text{g}}$. Under this condition, a coherent
force accumulates over many gyrations resulting in scattering of the
particle’s pitch-angle, i.e. the angle between the particle momentum and
$\bm{B}_{0}$. Scattering by an ensemble of wave-like perturbations, a
particle’s pitch-angle becomes diffuse and once particles have reversed
direction, the spatial transport along the background field also becomes
diffusive Hasselmann and Wibberenz (1970). Due to the resonant nature of the
interactions, the mean free path $\lambda_{\parallel}$ will depend on particle
rigidity 111Rigidity $\mathcal{R}$ is defined as the ratio of particle
momentum $p$ and charge $q$, that is $\mathcal{R}\equiv pc/q$ with the speed
of light $c$. and the turbulence power spectrum. Here, we introduce the
reduced rigidity $r_{\text{g}}{}/L_{\text{c}}$ where
$r_{\text{g}}{}=pc/(eB_{\text{rms}})=\mathcal{R}/B_{\text{rms}}$ is the
relativistic gyroradius, $L_{\text{c}}$ denotes the correlation length, $p$ is
the momentum and $e$ the charge of the particle. Generically, for a $k^{-5/3}$
Kolmogorov Kolmogorov (1931) power spectrum for instance,
$\lambda_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$.

Applying QLT to slab turbulence, where the wavevector
$\bm{k}\parallel\bm{B}_{0}$, the pitch-angle diffusion coefficient vanishes
for pitch-angles close to $90^{\circ}$. This is known as the $90^{\circ}$
problem Voelk (1975). The parallel mean free path, $\lambda_{\parallel}$,
which is proportional to the inverse of the pitch-angle diffusion coefficient,
however, remains finite if the turbulence spectrum is not too steep Shalchi
(2009). For isotropic turbulence, instead, $\lambda_{\parallel}$ diverges
Tautz  _et al._ (2006a). Various modifications of QLT have been suggested to
cure this problem. (See Ref. Shalchi (2020) for a recent review.) The most
promising suggestions are dynamical turbulence Tautz  _et al._ (2006b) and
non-linear theories Shalchi (2009), both of which lead to a broadening of the
resonance condition and can hence cure the vanishing of the pitch-angle
diffusion coefficient. Transport in the perpendicular direction, however, is a
different matter altogether. It is believed that it depends both on the
transport of particles along field lines as well as the transport of the field
lines itself. For the composite turbulence models popular in the space science
communities, various theories have been developped Bieber and Matthaeus
(1997); Matthaeus  _et al._ (2003); Shalchi  _et al._ (2004). See again Ref.
Shalchi (2020) for an account of successes and failures of such theories.

In simulations, the equations of motions are solved for a large number of
particles in individual realisations of the turbulent magnetic field
$\bm{\delta B}$. This field is either the result of MHD simulations Beresnyak
_et al._ (2011); Xu and Yan (2013); Cohet and Marcowith (2016) or has been
created as “synthetic turbulence”, that is according to a heuristic turbulence
model (e.g. isotropic turbulence with a power law power spectrum) Giacalone
and Jokipii (1994); Michalek and Ostrowski (1997, 1998); Giacalone and Jokipii
(1999); Casse  _et al._ (2001); Qin  _et al._ (2002); DeMarco  _et al._
(2007); Hussein  _et al._ (2015); Giacinti  _et al._ (2018); Reichherzer  _et
al._ (2020). While MHD turbulence is more realistic, its dynamical range is
limited; conversely, synthetic turbulence offers a larger dynamical range, but
it does not capture the full dynamics of MHD turbulence. However, even in
synthetic turbulence, the simulations have not been performed at GV or TV
rigidities for which the most precise direct observations of CRs are
available. Extrapolations from larger rigidities are of limited value due to
the possibility that the transport regime changes with rigidity. Various
features observed in the spectra of TV CRs could be taken as hints for such
transitions. See Ref. Mertsch (2020) for a review of test particle
simulations.

In the absence of a reliable theory and of simulations at the relevant
rigidities, the simplest expectation is for the parallel and perpendicular
mean free path to scale with rigidity in the same way while the dependence on
$\delta B^{2}/B_{0}^{2}$ should be the opposite Melrose (1980),

$\displaystyle\lambda_{\parallel}$
$\displaystyle\propto\left(\frac{r_{\text{g}}}{L_{\text{c}}}\right)^{1/3}\left(\frac{\delta
B^{2}}{B_{0}^{2}}\right)^{-1}\,,$ (1) $\displaystyle\lambda_{\perp}$
$\displaystyle\propto\left(\frac{r_{\text{g}}}{L_{\text{c}}}\right)^{1/3}\left(\frac{\delta
B^{2}}{B_{0}^{2}}\right)\,.$ (2)

Numerical simulations of isotropic turbulence, however, show a different
scaling of $\lambda_{\parallel}$ and $\lambda_{\perp}$. While
$\lambda_{\parallel}$ is in agreement with the
$(r_{\text{g}}/L_{\text{c}})^{1/3}$ behaviour generically expected for gyro-
resonant interactions with Kolmogorov turbulence, the scaling of
$\lambda_{\perp}$ with $r_{\text{g}}/L_{\text{c}}$ is faster. This behaviour
had been indicated already in the simulations by De Marco et al. DeMarco  _et
al._ (2007). The simulations by Snodin el al. Snodin  _et al._ (2016) and
Giacinti et al. Giacinti  _et al._ (2018) show a similar trend, even though
the respective authors refrain from a verdict on the reliability of this
result. More recently, Dundovic et al. Dundovic  _et al._ (2020) investigated
this behaviour in more detail, finding a ratio
$\lambda_{\perp}/\lambda_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{\Delta
s}$ with $\Delta s\simeq 0.2$. We also note that at the reduced rigidities
where simulation results are available, the scaling of $\lambda_{\perp}$ with
$(\delta B/B_{0})$ indicated in eq. (1) is not observed either, see for
instance Fig. 6 of Ref. Mertsch (2020). So far, neither the scaling with
rigidity nor with $(\delta B/B_{0})$ has been understood.

Note that for a correlation length of $30\,\text{pc}$ and
$B_{\text{rms}}=4\,\mu\text{G}$, reduced rigidities
$r_{\text{g}}{}/L_{\text{c}}\simeq 10^{-2}\mathellipsis 10^{-1}$ correspond to
rigidities $\mathcal{R}\simeq(1\mathellipsis 10)\,\text{PV}$. If extrapolated
to GV rigidities, this difference could have important phenomenological
consequences. Diffuse emission at GeV energies is sensitive to the spectra and
spatial distribution of CR nuclei and leptons of GV and TV rigidities. Some of
the deviations between isotropic diffusion models and observations have been
suggested to be due to a different rigidity-dependence of
$\lambda_{\parallel}$ and $\lambda_{\perp}$ Evoli  _et al._ (2012); Cerri  _et
al._ (2017); Reichherzer  _et al._ (2020). Anisotropic diffusion might also
play an important role in the confinement of CRs in different parts of the
Galaxy and in the overall escape at large rigidities Giacinti  _et al._
(2018).

We have revisited the important question on the nature of perpendicular
transport with a particular focus on the rigidity-dependence of the
perpendicular mean free path $\lambda_{\perp}$ in isotropic turbulence. To
this end, we have run a large suite of numerical test-particle simulations
reaching unprecedentedly low rigidities, $r_{\text{g}}{}/L_{\text{c}}\simeq
10^{-4}$ which for $L_{\text{c}}=30\,\text{pc}$ and
$B_{\text{rms}}=4\,\mu\text{G}$ corresponds to rigidities of $\mathcal{R}\sim
10\,\text{TV}$. These small rigidities have been made possible by the use of
graphics processing units (GPUs) for solving the equations of motions. (See
Ref. Tautz (2016) for an earlier use of GPUs in test particle simulations.) We
have found that the perpendicular mean free path $\lambda_{\perp}$ scales
differently than the parallel one for $10^{-2}\lesssim
r_{\text{g}}/L_{\text{c}}\lesssim 1$, but for $r_{\text{g}}/L_{\text{c}}\ll
10^{-2}$ the same rigidity-dependence is recovered. If the transition from the
standard $(r_{\text{g}}/L_{\text{c}})^{1/3}$-dependence at low values of
$r_{\text{g}}/L_{\text{c}}$ to the faster dependence on
$r_{\text{g}}/L_{\text{c}}$ at intermediate values of
$r_{\text{g}}/L_{\text{c}}$ depended on $(\delta B^{2}/B_{0}^{2})$, this could
explain why at intermediate values of $r_{\text{g}}/L_{\text{c}}$, a $(\delta
B^{2}/B_{0}^{2})$-dependence different from eq. (2) was observed.

The outline of this paper is as follows. In Sec. II, we explain the
preparation of an isotropic, turbulent magnetic field on a computer, describe
the test particle simulations and introduce the various correlation functions
and diffusion coefficients we determined from the simulation results. Our
results are presented and discussed in Sec. III. Sec. IV contains a discussion
of various sources of systematic uncertainties. We put our results into
context in Sec. V, offering some thoughts on the phenomenological consequences
of our findings, and conclude in Sec. VI.

##  II Methods

###  II.1 Test particle simulations

In order to study the perpendicular transport of cosmic rays in the presence
of a turbulent magnetic field, we will be making use of test particle
simulations, meaning that the back-reaction of cosmic rays onto the magnetic
field is neglected (see Ref. Mertsch (2020) for a review). In this limit the
magnetic field is defined once, before a set of test particle trajectories is
computed by solving the Newton-Lorentz equation

$\frac{\mathrm{d}}{\mathrm{d}t}\bm{p}=q\left[\bm{E}(\bm{r},t)+\frac{\bm{v}}{c}\times\bm{B}(\bm{r},t)\right]\,,$
(3)

where $q$ is the particle charge, $\bm{v}$ and $\bm{p}$ are the velocity and
momentum vectors, and $\bm{E}$ and $\bm{B}$ represent the electric and
magnetic vector fields. We will be assuming magnetostatic turbulence
throughout, meaning that $\partial\bm{B}/\partial t=0$. The electric field is
neglected since the high mobility of charges in typical astrophysical plasmas
is effectively preventing the presence of any large scale electric field
$\bm{E}_{0}$. On scales larger than the Debye length the electric charges
effectively shield any charge overdensity producing an electric potential
Debye and Hückel (1923). While a small scale electric field $\bm{\delta E}$
could in principle be present, it follows from Faraday’s induction law that
for Alfvén-wave type magnetic field perturbations, $|\bm{\delta
E}|\sim(v_{A}/c)|\bm{\delta B}|$ with $v_{A}$ denoting the Alfvén velocity.
Since in most astrophysical environments the Alfvén velocity is much smaller
than the speed of light, i.e. $v_{A}/c\simeq 10^{-5}-10^{-3}$ for typical
interstellar medium conditions Ferrière (2001), the small scale electric field
can also be neglected here. Note that in the absence of an electric field, the
particle energy is conserved. This has to hold for the test particle
trajectories in the simulations also.

The computed test particle trajectories can then be used to infer parameters
of the transport model, such as the diffusion coefficient Giacalone and
Jokipii (1994); Michalek and Ostrowski (1997, 1998); Giacalone and Jokipii
(1999); Casse  _et al._ (2001); Qin  _et al._ (2002); DeMarco  _et al._
(2007); Giacinti  _et al._ (2018); Reichherzer  _et al._ (2020), predict
spectra or large-scale anisotropies Giacinti  _et al._ (2014); Savchenko  _et
al._ (2015); Génolini and Ahlers (2021), or investigate effects beyond the
standard picture of cosmic ray transport such as the generation of small-scale
anisotropies in the arrival directions Giacinti and Sigl (2012); Ahlers and
Mertsch (2015); López-Barquero  _et al._ (2016); Ahlers and Mertsch (2017);
Mertsch and Ahlers (2019); Kuhlen  _et al._ (2022a).

####  II.1.1 Synthetic turbulence

We model the total magnetic field as the superposition of a large-scale,
coherent field $\bm{B}_{0}$ and a small-scale, turbulent field $\bm{\delta
B}$. Realistic magnetic field turbulence can be generated by solving the
equations of magnetohydrodynamics (see Ref. Schekochihin (2020) for a review).
There has been some effort towards using the magnetic fields generated by
these magnetohydrodynamic simulations and propagating cosmic rays through them
Beresnyak  _et al._ (2011); Xu and Yan (2013); Cohet and Marcowith (2016), but
the dynamical range that can be simulated limits these simulations to
normalized rigidities $r_{\text{g}}/L_{\text{c}}>10^{-2}$ Cohet and Marcowith
(2016). For a turbulent magnetic field of correlation length
$L_{\text{c}}=30\,\text{pc}$ and $B_{\text{rms}}=4\,\mu\text{G}$, this
corresponds to rigidities larger than $1\,\text{PV}$. Simulating particles
with magnetic rigidities small enough to compare to observations at tens of TV
instead is very challenging as these particles will interact with resonant
waves that have very small wavelengths. In order to simulate particles at
smaller rigidities one has to compromise between retaining as much of the
relevant physics as possible and making the model simple enough such that the
required dynamical range can be covered. To this end, we will use synthetic
magnetic field turbulence.

The synthetic, turbulent magnetic field considered in the following will be a
Gaussian random field that has the properties of a magnetic field and as such
fulfils Maxwell’s equations. In addition, the form of the power spectrum, that
is the Fourier transform of the two-point correlation function,

$\displaystyle
P_{ij}(\bm{k})\equiv\int\mathrm{d}^{3}r\,e^{\imath\bm{k}\cdot\bm{r}}\langle\delta
B_{i}(\bm{r}_{0})\delta B_{j}(\bm{r}_{0}+\bm{r})\rangle\,,$ (4)

can be motivated from analytical models of magnetic field turbulence. (We use
angled brackets to denote averaging over an ensemble of turbulent magnetic
fields.) In order to satisfy the divergencelessness condition, the most
general isotropic correlation tensor Batchelor and Press (1953) must take the
form

$P_{ij}(\bm{k})=g^{3D}(k)\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}+\imath\sigma(k)\epsilon_{ijm}\frac{k_{m}}{k}\right)$
(5)

with the $k$-dependent 3D power spectrum $g^{3D}(k)$. It can be shown that
this is also the most general form of the correlation tensor for axisymmetric
turbulence Matthaeus and Smith (1981) if $\sigma(k)$ is replaced with
$\sigma(k_{\parallel},k_{\perp})$. Here, we will only consider linearly
polarized waves, therefore the helicity function is set to be zero, i.e.
$\sigma(k)=0$. The normalization of the 3D power spectrum has to be such that

$\displaystyle\delta B^{2}$
$\displaystyle\equiv\int\mathrm{d}^{3}k~{}\left[P_{xx}(k)+P_{yy}(k)+P_{zz}(k)\right]$
(6)
$\displaystyle=8\pi\int_{0}^{\infty}\mathrm{d}k\,k^{2}g^{3D}(k)=8\pi\int_{0}^{\infty}\mathrm{d}k~{}g(k),$

where the equivalent one-dimensional power spectrum $g(k)=k^{2}g^{3D}(k)$ has
been introduced. In the following, a simple power law for the magnetic power
spectrum $g(k)\propto k^{-\gamma}$ with a spectral index $\gamma=5/3$ for
Kolmogorov turbulence Kolmogorov (1941) is assumed for
$2\pi/L_{\text{max}}\leq k\leq 2\pi/L_{\text{min}}$ and zero otherwise. A
power law magnetic power spectrum with the correct normalization is then given
by

$g(k)=\delta
B^{2}k^{-\gamma}\frac{(\gamma-1)(2\pi/L_{\text{max}})^{\gamma-1}}{1-(L_{\text{min}}/L_{\text{max}})}.$
(7)

We characterise the overall strength of the turbulent component with respect
to the total rms B-field through the turbulence level $\eta$,

$\eta\equiv\frac{\delta B^{2}}{B_{0}^{2}+\delta B^{2}}\,.$ (8)

This can be inverted to give

$\frac{\delta B^{2}}{B_{0}^{2}}=\frac{\eta}{1-\eta}\,.$ (9)

The correlation length of the magnetic field $L_{\text{c}}$ can be defined as
Harari  _et al._ (2002)

$\int_{-\infty}^{\infty}\mathrm{d}L~{}\langle\bm{\delta
B}(\bm{r}_{0})\cdot\bm{\delta B}(\bm{r}_{0}+\bm{\Delta r}(L))\rangle\equiv
L_{\text{c}}\delta B^{2},$ (10)

where $(\bm{r}_{0}+\bm{\Delta r}(L))$ is displaced with respect to
$\bm{r}_{0}$ by a distance $L$ along a fixed direction. For a power law
magnetic power spectrum, the correlation length evaluates to

$L_{\text{c}}=\frac{1}{2}L_{\text{max}}\frac{\gamma-1}{\gamma}\frac{1-(L_{\text{min}}/L_{\text{max}})^{\gamma}}{1-(L_{\text{min}}/L_{\text{max}})^{\gamma-1}}\,.$
(11)

For a Kolmogorov turbulence spectrum, assuming $L_{\rm max}\gg L_{\rm min}$,
one finds $L_{\text{c}}\simeq L_{\text{max}}/5$.

In our simulations, we are not sensitive to the individual parameter values of
particle rigidity $\mathcal{R}$, correlation length $L_{\text{c}}$ (or
equivalently outer scale $L_{\text{max}}$) and rms B-field, but only to the
combination $r_{\text{g}}/L_{\text{c}}$. Therefore, our simulations can be
applied to different combinations of these parameters. However, when
considering the use for a particular physics case, e.g. transport of Galactic
CRs, we oftentimes adopt particular fiducial parameter values. Specifically,
we consider an outer scale $L_{\text{max}}=150\,\text{pc}$ Chepurnov  _et al._
(2010) which for Kolmogorov turbulence corresponds to a correlation length of
$L_{\text{c}}\simeq 30\,\text{pc}$. Smaller values have been inferred for
interarm regions Haverkorn  _et al._ (2008). With an rms value of
$4\,\mu\text{G}$, the gyroradius evaluates to
$r_{\text{g}}=0.270\,\text{pc}(\mathcal{R}/\text{PV})(B_{\text{rms}}/4\,\mu\text{G})^{-1}$.

Finally, it is important to recognise that in individual realisations of the
turbulent magnetic field, modes with wavelengths much longer than the
gyroradius under consideration also create a large-scale, coherent
contribution to the total magnetic field. However, in contrast to background
field $\bm{B}_{0}$, the strength and direction of this contribution is not
fixed and varies between different realisations of the turbulent magnetic
field. We will refer to the sum of the background field $\bm{B}_{0}$ and these
additional contribution as the effective background field
$\bm{B}_{\text{eff}}$. In Sec. III.1 and III.2, we will highlight some effects
due the effective background field.

###  II.2 Generating synthetic turbulence

In the following, we will review the generation of synthetic, isotropic
magnetic field turbulence on a computer, using the so-called harmonic method,
grid method and nested grid method.

####  II.2.1 Harmonic Method

In the harmonic method, the turbulent magnetic field is computed for each
position $\bm{r}$ as the superposition of $N$ plane waves with (random)
amplitudes $A_{n}$, phases $\beta_{n}$ and polarisations $\hat{\bm{\xi}}_{n}$.
The wave vectors are split into a normalization and a random direction
$\bm{k}_{n}=k_{n}\hat{\bm{k}}_{n}$. The superposition of plane waves can then
be written as

$\bm{\delta
B}(\bm{r})=\sqrt{2}\sum_{n=0}^{N-1}A_{n}\hat{\bm{\xi}}_{n}\cos{(k_{n}\hat{\bm{k}}_{n}\cdot\bm{r}+\beta_{n})}.$
(12)

For isotropic turbulence, the $\hat{\bm{k}}_{n}$ are isotropically
distributed,

$\hat{\bm{k}}_{n}=\begin{pmatrix}\sqrt{1-\eta_{n}^{2}}\cos{\phi_{n}}\\\
\sqrt{1-\eta_{n}^{2}}\sin{\phi_{n}}\\\ \eta_{n}\end{pmatrix},$ (13)

with $\phi_{n}$ sampled uniformly in $[0,2\pi[$ and $\eta_{n}$ sampled
uniformly in $[-1,1]$. Note that this is not equivalent to sampling
$\theta_{n}$ uniformly in $[0,\pi]$ and setting $\eta_{n}=\cos{\theta_{n}}$,
as this does not sample the three-dimensional unit sphere uniformly. To fulfil
Gauss’ law Giacalone and Jokipii Giacalone and Jokipii (1994) originally
proposed a polarization vector of the form

$\hat{\bm{\xi}}_{n}=\begin{pmatrix}-\imath\sin{\alpha_{n}}\sin{\phi_{n}}+\cos{\alpha_{n}}\cos{\phi_{n}}\eta_{n}\\\
\imath\sin{\alpha_{n}}\cos{\phi_{n}}+\cos{\alpha_{n}}\sin{\phi_{n}}\eta_{n}\\\
-\sqrt{1-\eta_{n}^{2}}\cos{\alpha_{n}}\end{pmatrix},$ (14)

together with using a complex exponential of the form

$\bm{\delta
B}(\bm{r})=\sum_{n=0}^{N-1}A_{n}\hat{\bm{\xi}}_{n}\exp{\imath(k_{n}\hat{\bm{k}}_{n}\cdot\bm{r}+\beta_{n})}.$
(15)

instead of Eq. (12). Tautz and Dosch Tautz and Dosch (2013) made a slight
modification to this vector to avoid using complex numbers altogether. In
practice, this yields a small runtime improvement. The polarization vector
they suggested is given by

$\hat{\bm{\xi}}_{n}=\begin{pmatrix}-\sin{\alpha_{n}}\sin{\phi_{n}}+\cos{\alpha_{n}}\cos{\phi_{n}}\eta_{n}\\\
\sin{\alpha_{n}}\cos{\phi_{n}}+\cos{\alpha_{n}}\sin{\phi_{n}}\eta_{n}\\\
-\sqrt{1-\eta_{n}^{2}}\cos{\alpha_{n}}\end{pmatrix},$ (16)

where $\alpha_{n}$ is drawn from a uniform distribution between $[0,2\pi[$. In
total, this method therefore requires saving six real numbers per mode $n$.
The number of wave modes typically adopted for a dynamical range of
$k_{\text{min}}/k_{\text{max}}\sim 10^{4}$ is
$N=\mathcal{O}(100)-\mathcal{O}(1000)$ Mertsch (2020).

####  II.2.2 Grid method

The major drawback of the harmonic method is that while in principle it is
very accurate, at each position, for every particle the sum has to be
evaluated up to large $N$ which is computationally very expensive. To
circumvent this problem grid-based methods can be used Qin  _et al._ (2002);
DeMarco  _et al._ (2007); Snodin  _et al._ (2016).

In the grid method, the magnetic field is generated only once on a large grid
that is then repeated periodically to span the simulation volume. For each
particle and time step, the magnetic field is then interpolated between the
closest grid points. This method benefits greatly from the well developed Fast
Fourier Transform (FFT) routine Dongarra and Sullivan (2000). To make optimal
use of these routines, however, the number of grid points should be a power of
two Johnson and Frigo (2008a).

We consider a cubic spatial domain of side length $L_{\text{max}}$ and choose
for the grid to have $N^{3}$ grid points. In the harmonic domain, we use a
multi-index $(i,j,k)$ corresponding to a wave vector

$\bm{k}_{(i,j,k)}=\frac{2\pi}{L_{\text{max}}}\begin{pmatrix}i\\\ j\\\
k\end{pmatrix}.$ (17)

We also need the pair of discrete three-dimensional Fourier transforms for the
field component $\delta B_{l}$ and $\delta\tilde{B}_{l}$

$\begin{array}[]{r l}&\delta\tilde{B}_{l}(\bm{k}_{(i,j,k)})\\\
&=\displaystyle\sum_{m=0}^{N-1}\sum_{n=0}^{N-1}\sum_{p=0}^{N-1}e^{-2\pi\imath\left(\frac{im}{N}+\frac{jn}{N}+\frac{kp}{N}\right)}\delta
B_{l}(\bm{r}_{(m,n,p)})\end{array}$ (18)

and

$\displaystyle\delta B_{l}(\bm{r}_{(m,n,p)})$ (19)
$\displaystyle=\frac{1}{N^{3}}\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}e^{2\pi\imath\left(\frac{im}{N}+\frac{jn}{N}+\frac{kp}{N}\right)}\delta
B_{l}(\bm{k}_{(i,j,k)}).$ (20)

The hermiticity condition $\bm{\delta B}({\bm{k}_{(i,j,k)}})=\bm{\delta
B}^{*}(\bm{k}_{(N-i,N-j,N-k)})$ ensures that the magnetic field in position
space is real. This condition reduces the number of independent components by
half. It is the convention used by complex-to-real Fourier transform methods
to reduce the input size in _the last_ dimension to $(N/2+1)$ modes. To
account for negative frequencies any index, $i$ for example, larger than $N/2$
has to be replaced by its negative equivalent in Fourier space, $-(N-i)$.

For each point on the three-dimensional $k$-grid, $\bm{k}_{(i,j,k)}$, a random
direction $\delta\hat{\bm{B}}(\bm{k}_{(i,j,k)})$ of the magnetic field
amplitude orthogonal to $\bm{k}_{(i,j,k)}$ is chosen in order to ensure that
$\nabla\cdot\bm{B}=0$. The amplitude is scaled according to the magnetic power
spectrum $|\bm{\delta
B}(\bm{k}_{(i,j,k)})|\propto|\bm{k}_{(i,j,k)}|^{-\gamma/2-1}$; the additional
$-1$ in the exponent is needed to account for the difference between the
three-dimensional power spectrum $g^{3D}(k)$ and the effective one-dimensional
power spectrum $g(k)=k^{2}g^{3D}(k)$. The condition of a vanishing mean
magnetic field $\langle\bm{\delta B}\rangle=0$ is fulfilled by requiring
$\bm{\delta B}(\bm{k}_{(0,0,0)})=\bm{0}$. The magnetic field can then be
transformed from Fourier space to position space by using the discrete inverse
Fourier transform in Eq. (19). Finally, the variance has to be scaled to the
correct value $\langle\delta B^{2}\rangle$.

####  II.2.3 Nested Grid Method in Isotropic Turbulence

Since the magnetic field in three-dimensional isotropic turbulence depends on
three spatial components the memory requirements of the grid-based method grow
as the number of grid points cubed. Simulating a turbulent field with large
enough dynamical range for the particle rigidities considered can thus quickly
exceed the computer memory. It has therefore been suggested Giacinti  _et al._
(2012) to use a number of nested grids, each grid corresponding to a certain
range in wave number. The magnetic field can be calculated as the
superposition of the different magnetic field grids in position space where
the grids are repeated periodically over the entire simulation volume.

While evaluating the magnetic field in the grid-based method is
computationally much less expensive, it requires a prohibitively large amount
of memory to cover a large dynamical range. This memory requirement depends
strongly on the dynamical range that is needed for the problem. We will assume
a particle rigidity of $10\,\text{TV}$ and an outer scale of turbulence of
$L_{\text{max}}=150\,\text{pc}$. The smallest wavelength on the grid has to be
at least as small as the gyroradius at this energy. Assuming a field strength
of $B_{\text{rms}}=4\,\mu\text{G}$ this means $r_{g}/L_{\text{max}}=1.78\times
10^{-5}$. To avoid artefacts due to periodicity, the grid should be larger
than $L_{\text{max}}$ by a factor of a few, let’s say $4$. Also, it has proven
necessary in practice to have a smallest wavelength that is at least a factor
10 smaller than the gyroradius of the particle.

This means that the number of grid points required to span this dynamical
range is at least $N=10\times 4\times 2/(1.78\times 10^{-5})\simeq 4,494,382$.
Accounting for the requirements of the FFT algorithms Cooley and Tukey (1965);
Johnson and Frigo (2008b) that the number of grid points should be a power of
$2$ this increases $N$ to $N=2^{23}=8,388,608$ per dimension. For isotropic
turbulence this results in a total $N^{3}=590,295,810,358,705,651,712$ grid
points for each field direction. The total memory requirement, assuming
$8\,\text{B}$ per double precision floating point number, then evaluates to
$14,167,100\,\text{PB}$ which makes such simulations currently unfeasible.

To get around this constraint Ref. Giacinti  _et al._ (2012) proposed to
superimpose multiple grids with different sizes and therefore different
resolutions. The small wavelengths are resolved on the smallest grid with the
finest resolution while the modes with large wavelengths are resolved on a
larger grid with larger grid spacing.

To implement this method of generating the turbulent magnetic field the
dynamical range $[k_{\text{min}},k_{\text{max}}]$ is split into $G$ different
grids that are set up as described previously with each grid having a
dynamical range $[k_{g},k_{g+1}]$, where
$k_{0}=k_{\text{min}}=2\pi/L_{\text{max}}$ and
$k_{G}=k_{\text{max}}=2\pi/L_{\text{min}}$. In addition to the modes that
contribute to the power spectrum, each grid also includes some modes without
power as padding. This means that the grids are both larger and have a finer
resolution than the theoretical minimum requirement to resolve the modes. This
is necessary to avoid artefacts due to discreteness effects and grid
periodicity.

It is well known that the largest wavenumber that can theoretically be
resolved on a grid with finite resolution is given by the Nyquist wavenumber
$k_{\text{Nyquist}}=(1/2)2\pi/\Delta x$ Nyquist (1928). Waves with wavenumbers
larger than this will, due to aliasing effects, look like waves with smaller
wavenumbers traveling in the opposite direction Condon and Ransom (2016).
Additionally, a factor of a few is included to ensure that the modes with the
smallest wavelengths are always resolved by enough grid points. To avoid grid
periodicity effects one should also make sure that the largest wavelength fits
on the grid at least a few times.

For example, grids at each level $g$ have sizes $L_{g}$ and $N=128$ grid
points, implying $\Delta x_{g}=L_{g}/N$. However, we choose for only the modes
with wavelengths between $L_{g,\text{min}}=L_{g}/20$ and
$L_{g,\text{max}}=L_{g}/8$ to have non-vanishing amplitude. The next grid
level ($g+1$) is defined such that $L_{g+1,\text{max}}=L_{g,\text{min}}$.
Specifically, for the use case with $r_{\text{g}}/L_{\text{c}}=8.91\times
10^{-5}$, we have chosen 16 grid levels. At level $g=0$, the grid has a size
$L_{0}=40\,L_{\text{c}}$, $L_{0,\text{max}}=5\,L_{\text{c}}$ and
$L_{0,\text{min}}=2\,L_{\text{c}}$. At level $g=1$, $L_{1}=16\,L_{\text{c}}$,
$L_{1,\text{max}}=2\,L_{\text{c}}$ and $L_{1,\text{min}}=0.8\,L_{\text{c}}$
and so on. Finally, at level $g=15$, $L_{15}=4.3\times 10^{-5}\,L_{\text{c}}$,
$L_{15,\text{max}}=5.4\times 10^{-6}\,L_{\text{c}}$ and
$L_{15,\text{min}}=2.1\times 10^{-6}\,L_{\text{c}}$. The dynamic range
$L_{0,\text{max}}/L_{15,\text{min}}$ covered is thus more than six orders of
magnitude.

To ensure a continuous power spectrum, the amplitudes on each grid $g$ have to
be rescaled such that they contribute the correct fraction of the total power
$\delta B^{2}$ to the power spectrum. It follows that for a pure power law
power spectrum the variance $\delta B_{g}^{2}$ of an individual grid
containing modes in $[L_{g,\text{min}},L_{g,\text{max}}]$ needs to be
rescaled,

$\delta B_{g}^{2}\propto\delta
B^{2}\frac{L_{g,\text{min}}^{\gamma-1}-L_{g,\text{max}}^{\gamma-1}}{L_{\text{max}}^{\gamma-1}-L_{\text{min}}^{\gamma-1}}\,.$
(21)

Note that in particular for small reduced rigidities, at the end of our
simulations, only a very small fraction of particles will have travelled
further than a few times $L_{0}$. We therefore consider artefacts due to the
grid periodicity to be negligible.

####  II.2.4 ODE solver

To calculate the test particle trajectories the Newton-Lorentz equation (3)
with a vanishing electric field has to be solved. Formally, this is a second
order differential equation in time. To solve this numerically it can be
rewritten as a set of two first order differential equations,

$\displaystyle\frac{\mathrm{d}\bm{r}}{\mathrm{d}t}$ $\displaystyle=\bm{v}\,,$
(22) $\displaystyle\frac{1}{\Omega}\frac{\mathrm{d}\hat{\bm{v}}}{\mathrm{d}t}$
$\displaystyle=\hat{\bm{v}}\times\hat{\bm{B}}\,.$ (23)

Here, we have used that $\gamma=\text{const.}$ and introduced the relativistic
gyrofrequency $\Omega\equiv qB/(\gamma mc)$.

The generic problem of solving a first order ordinary differential equation is
a well studied problem Press  _et al._ (1992), with a multitude of
differential equation solvers available. The most widely used solver for
ordinary differential equations is probably the Runge-Kutta solver Runge
(1895); Kutta (1901) that is already implemented for example in the odeint
package of the boost library Ahnert and Mulansky (2011). It should be noted,
however, that in this method errors on the particle velocity can accumulate
and over a large number of steps the particle energy can change significantly.

The Hamiltonian of a charged particle in a stationary magnetic field is given
by Griffiths (1995)

$H(\bm{r},\bm{p})=\frac{1}{2}\left(\bm{p}-q\bm{A}(\bm{r})\right)^{2}\,,$ (24)

with the vector potential $\bm{A}(\bm{r})$ defined by
$\bm{B}(\bm{r})=\bm{\nabla}\times\bm{A}(\bm{r})$. It is easy to see that since
this Hamiltonian does not explicitly depend on time, the total energy of the
particles should be conserved. While there are symplectic differential
equation solvers that inherently conserve particle energy, they typically
require that the Hamiltonian is separable
$H(\bm{r},\bm{p})=T(\bm{p})+V(\bm{r})$. Due to its ease of implementation and
inherent ability to conserve energy, the differential equation solver that has
become the standard in charged particle propagation is the Boris solver Boris
(1970); Birdsall and Langdon (2005). It has been pointed out that, even though
this algorithm is technically not symplectic, it conserves phase-space volume,
which is a necessary but not a sufficient condition for symplectic algorithms
Qin  _et al._ (2013). It can be shown that the error on the energy has a
global upper bound for all time steps Higuera and Cary (2018). For the Runge-
Kutta method, we have used the Runge-Kutta-Dopri5 solver of the odeint package
in the boost library Ahnert and Mulansky (2011). A comparison of different
differential equation solvers for this problem can be found in Ref. Ripperda
_et al._ (2018). We also compare results for the Boris-Push and for the Runge-
Kutta solver in Sec. IV.2.

In order to study the transport of magnetic field lines, we have also solved
the equations for a field line $\bm{r}_{\text{FL}}(s)=(x(s),y(s),z(s))^{T}$,

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}$ $\displaystyle=\frac{\delta
B_{x}}{\sqrt{B_{0}^{2}+\delta B^{2}}}\,,$ (25)
$\displaystyle\frac{\mathrm{d}y}{\mathrm{d}s}$ $\displaystyle=\frac{\delta
B_{y}}{\sqrt{B_{0}^{2}+\delta B^{2}}}\,,$ (26)
$\displaystyle\frac{\mathrm{d}z}{\mathrm{d}s}$
$\displaystyle=\frac{B_{0}+\delta B_{z}}{\sqrt{B_{0}^{2}+\delta B^{2}}}\,.$
(27)

For the field-line integration, we have employed the Runge-Kutta-Dopri5
solver.

###  II.3 Observables

####  II.3.1 Pitch-angle auto-correlation function

Particle motion becomes diffusive as soon as individual particle velocities
become uncorrelated with their initial velocities. For diffusion along the
background magnetic field, the relevant velocity component is the parallel
one, $v_{\parallel}=v\mu$, $\mu$ being the cosine of the pitch-angle. The
decorrelation of the pitch-angle (cosine) due to interactions with the
turbulent magnetic field is known as pitch-angle scattering.

The decorrelation of the parallel velocity component can be investigated by
considering the pitch-angle autocorrelation function $C(t)$ Casse  _et al._
(2001),

$C(t)\equiv\frac{\langle\mu_{0}\mu(t)\rangle}{\langle\mu_{0}^{2}\rangle}\,,$
(28)

where $\mu_{0}\equiv\mu(0)$. The angled brackets $\langle\mathellipsis\rangle$
here denote averages over initial particle directions and over the ensemble of
turbulent magnetic fields. From this the scattering time can be defined as the
coherence time of the pitch-angle cosine $\mu$,

$\tau_{\text{s}}\equiv\int_{0}^{\infty}\mathrm{d}\tau\,C(\tau)\,.$ (29)

It has previously been argued Casse  _et al._ (2001) that the pitch-angle
autocorrelation function $C(t)$ takes an exponential form,
$C(t)=\exp[-t/\tau]$. For this form, the time constant $\tau$ is immediately
identified as the scattering time. In the following, we will generalise this
to $C(t)=A\exp[-t/\tau]$ with $A\leq 1$. For $A<1$, the scattering time
$\tau_{\text{s}}$ will be suppressed with respect to $\tau$, that is
$\tau_{\text{s}}=A\tau<\tau$.

####  II.3.2 Particle diffusion coefficients

The nature of particle transport, in particular the nature of transport of
charged CRs in turbulent magnetic fields, can be characterised by the time-
dependence of the mean-square displacement in different directions $i$,
$\langle(\Delta r_{i})^{2}\rangle$. If this is of power law form, that is

$\langle(\Delta
r_{i})^{2}\rangle=\left\langle\left(r_{i}(t)-r_{i}(0)\right)^{2}\right\rangle\propto
t^{\alpha}\,,$ (30)

transport is called sub-diffusive if $0<\alpha<1$, diffusive for $\alpha=1$,
super-diffusive for $1<\alpha<2$ and ballistic for $\alpha=2$ Shalchi (2009);
Mertsch (2020). The diffusion coefficient $\kappa$, that is the constant of
proportionality in Eq. (30) for the case of diffusive transport, plays a
central role in the transport theory of charged CRs. In this section, we will
describe various methods for how it can be derived, for instance from test
particle simulations.

We define the running diffusion coefficient with respect to direction $i$ as

$d_{ii}(t)\equiv\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\langle(r_{i}(t)-r_{i}(0))^{2}\rangle\,.$
(31)

Specifically, since the isotropy of the space is broken by the presence of a
background magnetic field, here assumed to point in the $z$-direction, we
distinguish the parallel and the perpendicular running diffusion coefficients,

$\displaystyle d_{\parallel}(t)$
$\displaystyle\equiv\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\langle(z(t)-z(0))^{2}\rangle\,,$
(32) $\displaystyle d_{\perp}(t)$
$\displaystyle\equiv\frac{1}{4}\frac{\mathrm{d}}{\mathrm{d}t}\left(\langle(x(t)-x(0))^{2}\rangle+\langle(y(t)-y(0))^{2}\rangle\right)\,.$
(33)

If transport is diffusive at late times, those converge towards the asymptotic
diffusion coefficients

$\displaystyle\kappa_{\parallel}\equiv\lim_{t\to\infty}d_{\parallel}(t)\,,$
(34) $\displaystyle\kappa_{\perp}\equiv\lim_{t\to\infty}d_{\perp}(t)\,.$ (35)

The TGK formalism Taylor (1922); Green (1954); Kubo (1957) allows to relate
$\kappa_{\parallel}$ and $\kappa_{\perp}$ to the velocity autocorrelation
functions. In particular,

$\langle(\Delta
r_{i})^{2}\rangle=2\int_{0}^{t}\mathrm{d}t^{\prime}\,(t-t^{\prime})\langle
v_{i}(0)v_{i}(t^{\prime})\rangle\,,$ (36)

and for the asymptotic diffusion coefficients

$\displaystyle\kappa_{\parallel}$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}t^{\prime}\,\langle
v_{z}(0)v_{z}(t^{\prime})\rangle\,,$ (37) $\displaystyle\kappa_{\perp}$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}t^{\prime}\,\langle
v_{x}(0)v_{x}(t^{\prime})\rangle=\int_{0}^{\infty}\mathrm{d}t^{\prime}\,\langle
v_{y}(0)v_{y}(t^{\prime})\rangle.$ (38)

For an isotropic particle distribution $\langle\mu^{2}\rangle=1/3$ and so we
can relate the scattering time $\tau_{\text{s}}$, eq. (29), and the parallel
diffusion coefficient $\kappa_{\parallel}$ via the pitch-angle autocorrelation
function $C(t)=3\langle\mu(t)\mu_{0}\rangle=3/v^{2}\langle
v_{z}(0)v_{z}(t)\rangle$,

$\tau_{\text{s}}=\\!\\!\int_{0}^{\infty}\mathrm{d}\tau\,C(\tau)=\frac{3}{v^{2}}\\!\int_{0}^{\infty}\mathrm{d}t^{\prime}\langle
v_{z}(0)v_{z}(t^{\prime})\rangle=\frac{3\kappa_{\parallel}}{v^{2}}.$ (39)

This coincides with the time at which a ballistic running diffusion
coefficient, $d_{\parallel}(t)=(v^{2}/3)t$ intersects with the asymptotic
diffusion coefficient $\kappa_{\parallel}$.

Before concluding this section, we note that an alternative definition of the
running diffusion coefficients is sometimes employed in the literature,

$\displaystyle d_{\parallel,\text{frac}}(t)\\!$
$\displaystyle\equiv\\!\frac{1}{2t}\langle(z(t)-z(0))^{2}\rangle\,,$ (40)
$\displaystyle d_{\perp,\text{frac}}(t)\\!$
$\displaystyle\equiv\\!\frac{1}{4t}\left(\langle(x(t)-x(0))^{2}\rangle\\!+\\!\langle(y(t)-y(0))^{2}\rangle\right)\\!.$
(41)

While the derivative definition of eqs. (32) and (33) agrees with the
fractional definition of eqs. (40) and (41) for the diffusive limit, in
general the running diffusion coefficients differ. In particular, the
fractional definition tends to smoothen features. In the following, we will
adopt the derivative definition although this will be more susceptible to
noise, but we will investigate possible differences in Sec. IV.3.

####  II.3.3 Field line diffusion coefficients

In the same way that we have characterised particle transport through the
power law dependence on time $t$ of the mean-square displacement, we can also
characterise the transport of magnetic field lines through the power law
dependence of the mean-square displacement of the field line $\bm{r}(s)$ on
arc-length $s$,

$\langle(\Delta
r_{i})^{2}\rangle=\left\langle\left(r_{i}(s)-r_{i}(0)\right)^{2}\right\rangle\propto
s^{\alpha}\,.$ (42)

Specifically, we say that field line transport is sub-diffusive if
$0<\alpha<1$, diffusive for $\alpha=1$, super-diffusive for $1<\alpha<2$ and
ballistic for $\alpha=2$.

The running (perpendicular) field line diffusion coefficient with respect to
$s$ is defined as

$d_{\text{FL}}(s)\\!\equiv\\!\frac{1}{4}\frac{\mathrm{d}}{\mathrm{d}s}\left(\langle(x(s)-x(0))^{2}\rangle\\!+\\!\langle(y(s)-y(0))^{2}\rangle\right)\\!.$
(43)

with the asymptotic value

$\kappa_{\text{FL}}\equiv\lim_{s\to\infty}d_{\text{FL}}(s)\,.$ (44)

Analogously, a parallel field line diffusion coefficient could be defined, but
we will not be needing this.

Note that oftentimes low turbulence levels are assumed such that the field
lines are almost parallel to the $z$-direction and $s\approx z$. In that case,
the running (perpendicular) field line diffusion coefficient with respect to
$z$ is defined as

$d^{\prime}_{\text{FL}}(z)\\!\equiv\\!\frac{1}{4}\frac{\mathrm{d}}{\mathrm{d}z}\\!\left(\langle(x(z)-x(0))^{2}\rangle\\!+\\!\langle(y(z)-y(0))^{2}\rangle\right)\\!.$
(45)

with the asymptotic value

$\kappa^{\prime}_{\text{FL}}\equiv\lim_{z\to\infty}d^{\prime}_{\text{FL}}(z)\,.$
(46)

Given that we extract field line positions
$\bm{r}_{\text{FL}}(s)=(x(s),y(s),z(s))^{T}$ for fixed values of $s$ from the
simulations, we have computed the vertical field line displacements as
functions of $z$, that is $x(z)$ and $y(z)$, by computing the mean $x$ and $y$
for given intervals in $z^{2}$. Specifically, we consider $51$ bins,
logarithmically spaced between $z^{2}=1.6\times 10^{-6}L_{\text{c}}^{2}$ and
$1.6\times 10^{6}L_{\text{c}}^{2}$. While this introduces a certain degree of
smoothing into the data, we judge the loss of precision to be acceptable.

##  III Results

###  III.1 Pitch-angle auto-correlation

We start the discussion of our results by looking at the pitch-angle
autocorrelation $C(t)$, see eq. (28). In Fig. 2 we show $C(t)$ as a function
of $t$ for different turbulence levels $\eta$, see eq. (8) at the fixed
reduced rigidity $r_{\text{g}}{}/L_{\text{c}}=8.9\times 10^{-3}$ which for our
fiducial parameters ($L_{\text{max}}=150\,\text{pc}$,
$B_{\text{rms}}=4\,\mu\text{G}$, see Sec. II.1.1) corresponds to a rigidity of
$1\,\text{PeV}$. We have chosen this rigidity since for
$r_{\text{g}}{}/L_{\text{c}}=8.9\times 10^{-3}$, most particles can resonantly
interact with the turbulent magnetic field.

Figure 1: Pitch-angle autocorrelation function $C(t)$ for
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$ and different turbulence levels
$\eta$. The dashed vertical line indicates $\Omega t=2\pi$. The dotted lines
are fits of the form $A\exp[-t/\tau]$ to the different $C(t)$ for $\Omega
t>2\pi$.

Figure 2: Pitch-angle autocorrelation function $C(t)$ for $\eta=0.5$ and
different reduced rigidities $r_{\text{g}}/L_{\text{c}}$. The dashed vertical
line indicates $\Omega t=2\pi$.

Three separate phases can be distinguished in the time-dependence of the
pitch-angle autocorrelation.

####  III.1.1 Early times ($\Omega t\ll 1$)

For $\Omega t\ll 1$, the particles are moving ballistically into their initial
direction and are not yet affected by the magnetic field. The correlation of
all velocity components, in particular of the $z$-component
$\mu=v_{z}/|\bm{v}|$ is therefore one,

$C(t)=\frac{\langle\mu_{0}\mu(t)\rangle}{\langle(\mu_{0})^{2}\rangle}\simeq\frac{\langle\mu_{0}\mu_{0}\rangle}{\langle(\mu_{0})^{2}\rangle}=1\,.$
(47)

####  III.1.2 Intermediate times ($1\ll\Omega t\ll\Omega\tau_{\text{s}}$)

For $1\ll\Omega t\ll\Omega\tau_{\text{s}}$, particles have started to gyrate
in the effective background magnetic field. This results in a plateau in the
pitch-angle autocorrelation function: For low turbulence levels, the
suppression is moderate, for large turbulence levels it is much stronger.
Depending on the initial pitch-angle and the angle $\theta_{\mathrm{eff}}$
between $\bm{B}_{\text{eff}}$ and the $z$-direction, the correlation between
$\mu_{0}$ and $\mu(t)$ oscillates between a suppressed value and $1$. Imagine,
for instance, the time-dependence of the $z$-component of the particle
velocity, $v_{z}(t)$, for the extreme case that the effective background
magnetic field is pointing along the $x$-direction. (This could be the case
for a large turbulence level where $B_{0}$ is subdominant.) In that case, the
$z$-component would oscillate between its initial value $v_{z}(0)$ and the
negative opposite, $-v_{z}(0)$. In the average over all initial particle
directions and over the ensemble of turbulent magnetic fields, this would
result in a suppression of the pitch-angle autocorrelation $C(t)$.

We can formalise this reasoning as follows. Consider the initial (normalised)
velocity
$\hat{\bm{v}}_{0}=(\sqrt{1-\mu_{0}^{2}}\cos\phi_{0},\sqrt{1-\mu_{0}^{2}}\sin\phi_{0},\mu_{0})^{T}$,
parametrised by the initial pitch-angle $\mu_{0}$ and azimuth $\phi_{0}$.
Without loss of generality, we take the effective background field to point in
the direction $(\sin\theta_{\mathrm{eff}},0,\cos\theta_{\mathrm{eff}})^{T}$.
The gyration around this direction by an angle $\psi=\Omega t$ results in the
pitch-angle cosine

$\displaystyle\mu(t)=$
$\displaystyle\sqrt{1-\mu_{0}^{2}}\sin\theta_{\mathrm{eff}}\sin\phi_{0}\sin\Omega
t$ (48)
$\displaystyle-\sqrt{1-\mu_{0}^{2}}\sin\theta_{\mathrm{eff}}\cos\theta_{\mathrm{eff}}\cos\phi_{0}(\cos\Omega
t-1)$ $\displaystyle+\mu_{0}(\cos^{2}\theta_{\mathrm{eff}}+\cos\Omega
t\sin^{2}\theta_{\mathrm{eff}})\,.$

The pitch-angle correlation then computes as

$\displaystyle C(t)$
$\displaystyle=3\left\langle\mu_{0}\mu(t)\right\rangle_{0}$ (49)
$\displaystyle=3\frac{1}{4\pi}\int_{0}^{2\pi}\mathrm{d}\phi_{0}\,\int_{-1}^{1}\mathrm{d}\mu_{0}\,\mu_{0}\mu(t)$
(50)
$\displaystyle=\cos^{2}\theta_{\mathrm{eff}}+\sin^{2}\theta_{\mathrm{eff}}\cos\Omega
t\,,$ (51)

where the angled brackets $\langle\cdot\rangle_{0}$ denote the average over
the initial directions.

For the distribution of effective background field directions, we adopt a von
Mises-Fisher distribution Fisher (1953),

$F(\theta_{\mathrm{eff}})\equiv\sigma_{\mathrm{eff}}^{-2}\operatorname{csch}\left(\sigma_{\mathrm{eff}}^{-2}\right)\exp\left[\sigma_{\mathrm{eff}}^{-2}\cos\theta_{\mathrm{eff}}\right]/2\,,$
(52)

an extension of a 2D Gaussian to the sphere. We estimate the dispersion
$\sigma_{\mathrm{eff}}$ of the $\theta_{\mathrm{eff}}$-distribution to be

$\tan^{2}\sigma_{\mathrm{eff}}=\frac{\delta B^{2}}{B_{0}^{2}}\,.$ (53)

The ensemble-averaged pitch-angle autocorrelation is then given by the Fisher-
distribution-weighted average of eq. (51),

$\displaystyle\langle C(t)\rangle_{\theta_{\mathrm{eff}}}$
$\displaystyle\equiv\int_{0}^{\pi}\mathrm{d}\theta_{\mathrm{eff}}\sin\theta_{\mathrm{eff}}F(\theta_{\mathrm{eff}})C(t)$
(54) $\displaystyle=1-2\sigma_{\mathrm{eff}}^{4}\cos(\Omega
t)+2\sigma_{\mathrm{eff}}^{4}$ (55)
$\displaystyle+2\sigma_{\mathrm{eff}}^{2}\coth\left(\frac{1}{\sigma_{\mathrm{eff}}^{2}}\right)(\cos(\Omega
t)-1)\,.$ (56)

Finally, averaging over gyroperiods,

$\langle
C(t)\rangle_{\theta_{\mathrm{eff}},t}=1+2\sigma_{\mathrm{eff}}^{4}-2\sigma_{\mathrm{eff}}^{2}\coth\left(\frac{1}{\sigma_{\mathrm{eff}}^{2}}\right)\,.$
(57)

For the turbulence levels adopted in Fig. 2, that is $\eta=0.1,0.2,0.5$ and
$0.8$, we find $A=\langle
C(t)\rangle_{\theta_{\mathrm{eff}},t}=0.81,0.66,0.43$ and $0.36$,
respectively.

####  III.1.3 Late times ($\Omega t\gg\Omega\tau_{\text{s}}$)

Finally, for $\Omega t\gg\Omega\tau_{\text{s}}$, the pitch-angle
autocorrelation is exponentially suppressed due to pitch-angle scattering. For
low turbulence levels, the time constant, that is the scattering time
$\tau_{\text{s}}$ is large, for large turbulence levels, it becomes smaller,
as expected.

####  III.1.4 Comparison with simulation results

We have tested our understanding of the plateau for $\Omega t\gg 1$ by fitting
an exponential function $A\exp[-t/\tau_{\text{s}}]$ to the pitch-angle
autocorrelation function $C(t)$. These fits are indicated by the dashed lines
in Fig. 2. The parameter values $A$ for $\eta=0.1,0.2,0.5$ and $0.8$ are
$0.84,0.70,0.40$ and $0.20$ which compares favourably with the above
predictions.

In Fig. 2, we also compare the pitch-angle auto correlation functions for the
fixed turbulence level $\eta=0.5$ and different reduced rigidities. We note
that the suppression for $1\ll\Omega t\ll\Omega\tau_{\text{s}}$ does not
depend on the reduced rigidity, as predicted by our reasoning above, as long
as $r_{\text{g}}/L_{\text{c}}\ll 1$. The scattering time, $\tau_{\text{s}}$,
however grows with reduced rigidity, as expected.

###  III.2 Running diffusion coefficients

####  III.2.1 Running parallel diffusion coefficient

Figure 3:  Running parallel diffusion coefficient $d_{\parallel}(t)$ at
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$ and for various turbulence
levels $\eta$. The solid lines show the mean over the ensemble of turbulent
magnetic fields, the shaded band indicate the standard mean error. The dashed
lines indicate the suppressed ballistic growth. The dotted horizontal lines
show the asymptotic diffusion coefficients $\kappa_{\parallel}$ and the
vertical dot-dashed lines show the scattering time $\tau_{\text{s}}^{\prime}$,
increased due to the suppression of the ballistic growth.

Figure 4:  Same as Fig. 4, but for fixed $\eta=0.5$ and various reduced
rigidities $r_{\text{g}}/L_{\text{c}}$. The grey dashed line indicates the
initial, unsuppressed ballistic phase.

For the same fixed reduced rigidity as before,
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$, we illustrate the dependence of
the running parallel diffusion coefficient $d_{\parallel}(t)$ on the
turbulence level $\eta$ in Fig. 4. One can identify the same three stages as
before and the $\eta$-dependence of the suppression of the pitch-angle auto-
correlation function at intermediate times ($1\ll\Omega
t\ll\Omega\tau_{\text{s}}$) discussed in Sec. III.1.

  1. 1.

For $\Omega t\ll 1$, the particles are moving ballistically into their initial
direction and are not yet affected by the magnetic field. The mean square
displacement, e.g. in the $z$-direction is
$\langle(z-z(0))^{2}\rangle=(1/3)v^{2}t^{2}$. Therefore,
$d_{\parallel}(t)=v^{2}t/3$ or $d_{\parallel}/(vL_{\text{c}})=(1/3)\Omega
t\,r_{\text{g}}/L_{\text{c}}$.

  2. 2.

For $1\ll\Omega t\ll\Omega\tau_{\text{s}}$, particles have started to gyrate
in the effective background magnetic field. While the mean-square displacement
still grows $\propto t^{2}$, it does so with the rate reduced by the same
factor as for the pitch-angle auto-correlation function; for instance, for
$\eta=0.5$, this factor is $\approx 0.4$.

  3. 3.

Finally, for $\Omega t\gg\Omega\tau_{\text{s}}$, pitch-angle scattering is
reducing the mean square displacement to diffusive behaviour,
$\langle(z-z(0))^{2}\rangle\propto t$. Consequently, the running diffusion
coefficient, $d_{\parallel}(t)$ flattens out. The asymptotic value
$\kappa_{\parallel}(t)$ and the transition time, of course, depend on the
reduced rigidity $r_{\text{g}}/L_{\text{c}}$. We will investigate this
rigidity-dependence in more detail in Sec. III.3.

With decreasing turbulence level $\eta$, the scattering time increases and so
does the asymptotic diffusion coefficient.

In Fig. 4 we show the running parallel diffusion coefficient
$d_{\parallel}(t)$, see eq. (32), for different reduced rigidities
$r_{\text{g}}/L_{\text{c}}$ and fixed turbulence level $\eta=0.5$. Note that
we have rescaled the time $\Omega t$ with the reduced rigidities, in order to
have all running diffusion coefficients overlap during the ballistic phase
despite the different $\Omega$.

Figure 5: Running perpendicular diffusion coefficient $d_{\perp}(t)$ at
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$ and for various turbulence
levels $\eta$. The solid lines show the mean over the ensemble of turbulent
magnetic fields, the shaded band indicate the standard mean error. The dashed
lines indicate the suppressed ballistic growth. The dotted horizontal lines
show the asymptotic diffusion coefficients $\kappa_{\perp}$ and the vertical
dot-dashed lines show the scattering time $\tau_{\text{s}}^{\prime}$,
increased due to the suppression of the ballistic growth.

Figure 6:  Same as Fig. 6, but for fixed $\eta=0.5$ and various reduced
rigidities $r_{\text{g}}/L_{\text{c}}$. The grey dashed line indicates the
initial, unsuppressed ballistic phase.  Figure 7: Running parallel and
perpendicular diffusion coefficients $d_{\parallel}(t)$ and $d_{\perp}(t)$,
respectively, for $\eta=0.5$ and for various reduced rigidities
$r_{\text{g}}/L_{\text{c}}$. The dotted and dash-dotted lines indicate the
asymptotic diffusion coefficients $\kappa_{\parallel}$ and $\kappa_{\perp}$,
respectively.

####  III.2.2 Running perpendicular diffusion coefficient

In Fig. 6 we show the running perpendicular diffusion coefficient
$d_{\perp}(t)$, see eq. (33), at fixed reduced rigidities
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$ and for different turbulence
level $\eta$. The time-dependence is more complicated than that of the pitch-
angle auto-correlation function or of the running parallel diffusion
coefficient in that we can distinguish four separate phases which are
separated by the times $\Omega t=1,\,\Omega\tau_{\text{s}}$ and
$\Omega\tau_{\text{c}}$:

  1. 1.

For $\Omega t\ll 1$, the particles are not affected by the magnetic field and
move ballistically. The mean-square displacement in the $x$\- and
$y$-directions does not differ from the one in the $z$-direction and so the
running diffusion coefficients are both
$d_{\parallel}(t)=d_{\perp}(t)=v^{2}t/3$ or
$d_{\parallel}/(vL_{\text{c}})=d_{\perp}/(vL_{\text{c}})=(1/3)\Omega
t\,r_{\text{g}}/L_{\text{c}}$.

  2. 2.

For $1\ll\Omega t\ll\Omega\tau_{\text{s}}$, particle motion does get affected
by the effective background magnetic field and particles start gyrating about
the effective background field direction. This is best seen at low turbulence
levels, e.g. $\eta=0.1$ or $0.2$. There is also a suppression of the
perpendicular running diffusion coefficient, but in the opposite direction:
Whereas for $d_{\parallel}$, the suppression was stronger for larger
turbulence level $\eta$, for $d_{\perp}$ the suppression is larger for small
$\eta$. This can be understood in the following way: In the limit $\delta B\to
0$, the motion in the $x$\- and $y$-directions would be merely a gyration. A
finite $\delta B$ misaligns the effective $B$-field direction from the
$z$-direction such that part of the gyration in the perpendicular plane now
contributes to the parallel motion; in turn, some of the parallel ballistic
motion, contributes to the transport in the $x$\- and $y$-directions. The
larger $\eta$, the stronger is this effect.

  3. 3.

For $\Omega\tau_{\text{s}}\ll\Omega t\ll\Omega\tau_{\text{c}}$, perpendicular
transport is sub-diffusive, $d_{\perp}\propto t^{\alpha-1}$ with $\alpha<1$.
We note that this is reminiscent of so-called compound subdiffusion Webb  _et
al._ (2006), an effect due to diffusive particle transport along diffusive
field lines. However, while compound subdiffusion predicts $d_{\perp}\propto
t^{1/2}$, we do not find this scaling for all $\eta$, see the dotted line in
Fig. 6. We note that the transition to this regime sets in later for smaller
$\eta$ as $\tau_{\text{s}}$ is larger.

  4. 4.

For $\Omega\tau_{\text{c}}\ll\Omega t$, finally, perpendicular transport also
is diffusive and $d_{\perp}$ attains its asymptotic value, $\kappa_{\perp}$.

Figure 8:  Asymptotic parallel and perpendicular mean free paths
$\lambda_{\parallel}$ and $\lambda_{\perp}$ as a function of reduced rigidity
$r_{\text{g}}/L_{\text{c}}$ for different turbulence levels
$\eta=0.2,\,0.5,\,0.8$ and $1$. We also compare the results of our test
particle simulations with the data of Giacinti et al. Giacinti and Kirk (2017)
(empty circles, squares, and pentagons) as well as Lopez-Coto & Giacinti
López-Coto and Giacinti (2018) (empty hexagons). We have indicated power laws
$\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ for
$r_{\text{g}}/L_{\text{c}}<10^{-1}$ and
$\propto(r_{\text{g}}/L_{\text{c}})^{2}$ for $r_{\text{g}}/L_{\text{c}}>1$.
Figure 9:  Ratio $\lambda_{\perp}/\lambda_{\parallel}$ of the asymptotic
perpendicular and parallel mean free paths as a function of reduced rigidity
$r_{\text{g}}/L_{\text{c}}$ for different turbulence levels $\eta$. We have
indicated a constant ratio at low rigidities and a
$(r_{\text{g}}/L_{\text{c}})^{0.2}$ dependence at intermediate rigidities, as
suggested by Dundovic et al. Dundovic  _et al._ (2020).

In Fig. 6 we show the running perpendicular diffusion coefficient
$d_{\perp}(t)$, see eq. (33), for different reduced rigidities
$r_{\text{g}}/L_{\text{c}}$ and fixed turbulence level $\eta=0.5$. For
$r_{\text{g}}/L_{\text{c}}\lesssim 1$, the behaviour is qualitatively the same
as described above for the case of $r_{\text{g}}/L_{\text{c}}=8.9\times
10^{-3}$. For larger values of $r_{\text{g}}/L_{\text{c}}$, the behaviour
becomes markedly different. Here, we can at most make out three phases:

  1. 1.

For $\Omega t\ll 1$, the particles move ballistically in the same way as for
lower rigidities.

  2. 2.

For $1\ll\Omega t\ll\Omega\tau_{\text{s}}$, particles gyrate around the
background field direction. The direction of the effective background field
does not play any role anymore since the particles at these rigidities cannot
be in resonance with any of the turbulence modes anymore. Therefore, particles
in different realisations gyrate around the original $z$-direction. These
gyrations are clearly visible in Fig. 6.

  3. 3.

For $\Omega t\gg\Omega\tau_{\text{s}}$, the gyrations get damped as particles
perform small-angle scattering and the asymptotic perpendicular diffusion
coefficient $\kappa_{\perp}$ is attained. Note that the value becomes
independent of the reduced rigidity $r_{\text{g}}/L_{\text{c}}$.

Having established the subdiffusive phase in the perpendicular diffusion
coefficient in Fig. 6, we can revisit Fig. 4 and also identify a subdiffusive
phase in the parallel diffusion coefficient, see $d_{\parallel}$ for
$\eta=0.8$ and $100\lesssim\Omega t\lesssim 1000$. While this might sound
surprising, this is unavoidable at high enough turbulence levels; after all in
the limit $\eta\to 1$, both running diffusion coefficients must agree.

A first hint of the non-standard dependence on $r_{\text{g}}/L_{\text{c}}$ of
the perpendicular diffusion coefficient $\kappa_{\perp}$ can already be
gleaned from the spacing of the dotted horizontal lines in Figs. 4 and 6. Note
that those lines refer to reduced rigidities $r_{\text{g}}/L_{\text{c}}$ which
are equally spaced in $\log(r_{\text{g}}/L_{\text{c}})$. Constant distances
between lines belonging to neighbouring $r_{\text{g}}/L_{\text{c}}$ values
hence indicate power law behaviour. This is for example the case for the
lowest four lines in Fig. 4, meaning that the parallel diffusion coefficient
follows a power law between $r_{\text{g}}/L_{\text{c}}\sim 10^{-4}$ and
$10^{-1}$. The spacing of the other lines gets wider, meaning that the
dependence on $r_{\text{g}}/L_{\text{c}}$ gets faster. For the perpendicular
diffusion coefficient, however, only the first two lines have the same spacing
as the first lines for the parallel diffusion coefficient. Subsequent lines
have wider spacing, meaning that the diffusion coefficient will have a faster
dependence on $r_{\text{g}}/L_{\text{c}}$. Finally for
$r_{\text{g}}/L_{\text{c}}\gtrsim 1$, all the perpendicular diffusion
coefficients coincide.

Fig. 7 compares the running parallel and perpendicular diffusion coefficients
$d_{\parallel}$ and $d_{\perp}$ of individual rigidities in the same panels.
While the data presented is the same as in Figs. 4 and 6, Fig. 7 allows
comparing $d_{\parallel}$ and $d_{\perp}$ more directly. In particular this
allows directly identifying the time where $d_{\perp}$ starts to be sub-
diffusive (for $r_{\text{g}}/L_{\text{c}}\ll 1$) and where the gyrations in
$d_{\perp}$ get damped (for $r_{\text{g}}/L_{\text{c}}\gtrsim 1$) with the
scattering time $\tau_{\text{s}}^{\prime}$.

###  III.3 Diffusion coefficients as function of rigidity

Figure 10:  Asymptotic parallel and perpendicular mean free paths
$\lambda_{\parallel}$ and $\lambda_{\perp}$ as a function of turbulence levels
$\eta$ for a subset of reduced rigidities $r_{\text{g}}/L_{\text{c}}$. The
solid and dotted lines show the fitted $\lambda_{\parallel}$ and
$\lambda_{\perp}$, respectively, see eqs. (58) and 59.

In Fig. 8 we show the mean free paths $\lambda_{\parallel}$ and
$\lambda_{\perp}$ as functions of rigidity, for the different turbulence
levels $\eta=0.2,\,0.5,\,0.8$ and $1$. Note that we have derived the
asymptotic values by averaging the values $d_{\parallel}$ and $d_{\perp}$
within $\Omega t_{\rm max}/2\lesssim\Omega t\leq\Omega t_{\rm max}$ where
$t_{\rm max}$ is the maximum run time for the simulations. The errors of these
estimates are also obtained by averaging the standard errors of the mean
within this time period. We have also plotted some lines with constant power
law indices. For reduced rigidities $r_{\text{g}}/L_{\text{c}}\ll 1$, the
asymptotic parallel diffusion coefficient $\kappa_{\parallel}$ exhibits the
$(r_{\text{g}}/L_{\text{c}})^{1/3}$-dependence expected for gyro-resonant
interactions due to turbulence with a power spectrum $g(k)\propto k^{-5/3}$.
At rigidities $r_{\text{g}}/L_{\text{c}}\gg 1$, the
$(r_{\text{g}}/L_{\text{c}})^{2}$-dependence of small-angle scattering is
visible (e.g. Gruzinov (2018)). We note that the normalisation of the
diffusion coefficient decreases with increasing turbulence level $\eta$. We
have also plotted the results from Giacinti et al. (2018) Giacinti and Kirk
(2017) which uses the same turbulence model. The agreement is excellent in the
range where the simulation data overlap. Note that Ref. Sonsrettee  _et al._
(2016) has considered rigidities as low as $\sim 2\times 10^{-2}$; however,
the turbulence spectrum considered differed, so we cannot compare directly. By
and large, they are in agreement though.

As far as the perpendicular diffusion coefficient is concerned, its rigidity-
dependence is more complicated. Again, we show a power law
$\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ that matches the simulation results
at low rigidities. However, it appears that the data follow this dependence
only over a limited rigidity range. At intermediate rigidities,
$\kappa_{\perp}$ grows faster than $(r_{\text{g}}/L_{\text{c}})^{1/3}$; for
$r_{\text{g}}/L_{\text{c}}\gtrsim 1$, $\kappa_{\perp}$ becomes constant. The
behaviour at intermediate and large rigidities is in line with the behaviour
seen in Ref. Dundovic  _et al._ (2020). The normalisation at low rigidities
shows the expected ordering in that the diffusion coefficient is smaller for
smaller $\eta$. The deviation from the
$(r_{\text{g}}/L_{\text{c}})^{1/3}$-behaviour sets in at rigidities
$r_{\text{g}}/L_{\text{c}}\sim 10^{-2}$ for $\eta=0.8$ and at
$r_{\text{g}}/L_{\text{c}}\sim 3\times 10^{-3}$ for $\eta=0.5$. We have not
actually observed the $(r_{\text{g}}/L_{\text{c}})^{1/3}$-scaling for
$\eta=0.2$ for the rigidities for which we were able to run simulations. (For
low turbulence levels the scattering time $\tau_{\text{s}}^{\prime}$ for which
$d_{\parallel}$ becomes constant and the time $\tau_{\text{c}}$ for which
$d_{\perp}$ become constant become prohibitively large.) Overall, the
agreement with Giacinti et al. (2018) Giacinti and Kirk (2017) is good.

$r_{\text{g}}/L_{\text{c}}$ | $\eta=0.2$ | $\eta=0.5$ | $\eta=0.8$
---|---|---|---
$\lambda_{\parallel}\left[L_{c}\right]$ | $\lambda_{\perp}\left[L_{c}\right]$ | $\lambda_{\parallel}\left[L_{c}\right]$ | $\lambda_{\perp}\left[L_{c}\right]$ | $\lambda_{\parallel}\left[L_{c}\right]$ | $\lambda_{\perp}\left[L_{c}\right]$
$8.90\times 10^{-5}$ | $-$ | $-$ | $\left(1.614\pm 0.001\right)\times 10^{-1}$ | $\left(3.096\pm 0.006\right)\times 10^{-3}$ | $-$ | $-$
$4.45\times 10^{-4}$ | $-$ | $-$ | $\left(2.621\pm 0.009\right)\times 10^{-1}$ | $\left(5.331\pm 0.045\right)\times 10^{-3}$ | $\left(7.527\pm 0.018\right)\times 10^{-2}$ | $\left(1.040\pm 0.008\right)\times 10^{-2}$
$8.90\times 10^{-4}$ | $-$ | $-$ | $\left(3.030\pm 0.001\right)\times 10^{-1}$ | $\left(6.057\pm 0.048\right)\times 10^{-3}$ | $\left(9.070\pm 0.011\right)\times 10^{-2}$ | $\left(1.202\pm 0.001\right)\times 10^{-2}$
$2.67\times 10^{-3}$ | $-$ | $-$ | $\left(4.203\pm 0.002\right)\times 10^{-1}$ | $\left(9.156\pm 0.039\right)\times 10^{-3}$ | $\left(1.285\pm 0.001\right)\times 10^{-1}$ | $\left(1.904\pm 0.001\right)\times 10^{-2}$
$4.45\times 10^{-3}$ | $\left(1.947\pm 0.050\right)\times 10^{0}$ | $\left(3.945\pm 0.060\right)\times 10^{-3}$ | $\left(5.148\pm 0.118\right)\times 10^{-1}$ | $\left(1.234\pm 0.009\right)\times 10^{-2}$ | $\left(1.584\pm 0.020\right)\times 10^{-1}$ | $\left(2.486\pm 0.033\right)\times 10^{-2}$
$7.12\times 10^{-3}$ | $\left(2.181\pm 0.014\right)\times 10^{0}$ | $\left(4.890\pm 0.048\right)\times 10^{-3}$ | $\left(5.827\pm 0.052\right)\times 10^{-1}$ | $\left(1.584\pm 0.010\right)\times 10^{-2}$ | $\left(1.817\pm 0.022\right)\times 10^{-1}$ | $\left(3.191\pm 0.013\right)\times 10^{-2}$
$8.90\times 10^{-3}$ | $\left(2.310\pm 0.019\right)\times 10^{0}$ | $\left(5.367\pm 0.054\right)\times 10^{-3}$ | $\left(6.265\pm 0.046\right)\times 10^{-1}$ | $\left(1.767\pm 0.012\right)\times 10^{-2}$ | $\left(1.980\pm 0.016\right)\times 10^{-1}$ | $\left(3.499\pm 0.024\right)\times 10^{-2}$
$1.78\times 10^{-2}$ | $\left(2.890\pm 0.026\right)\times 10^{0}$ | $\left(7.503\pm 0.057\right)\times 10^{-3}$ | $\left(7.872\pm 0.061\right)\times 10^{-1}$ | $\left(2.489\pm 0.016\right)\times 10^{-2}$ | $\left(2.613\pm 0.020\right)\times 10^{-1}$ | $\left(5.092\pm 0.018\right)\times 10^{-2}$
$2.67\times 10^{-2}$ | $\left(3.270\pm 0.025\right)\times 10^{0}$ | $\left(9.303\pm 0.048\right)\times 10^{-3}$ | $\left(9.172\pm 0.096\right)\times 10^{-1}$ | $\left(3.111\pm 0.017\right)\times 10^{-2}$ | $\left(3.116\pm 0.022\right)\times 10^{-1}$ | $\left(6.378\pm 0.043\right)\times 10^{-2}$
$4.45\times 10^{-2}$ | $\left(3.968\pm 0.036\right)\times 10^{0}$ | $\left(1.275\pm 0.005\right)\times 10^{-2}$ | $\left(1.133\pm 0.008\right)\times 10^{0}$ | $\left(4.158\pm 0.024\right)\times 10^{-2}$ | $\left(3.879\pm 0.020\right)\times 10^{-1}$ | $\left(8.822\pm 0.065\right)\times 10^{-2}$
$5.34\times 10^{-2}$ | $\left(4.208\pm 0.028\right)\times 10^{0}$ | $\left(1.453\pm 0.007\right)\times 10^{-2}$ | $\left(1.215\pm 0.007\right)\times 10^{0}$ | $\left(4.703\pm 0.021\right)\times 10^{-2}$ | $\left(4.282\pm 0.024\right)\times 10^{-1}$ | $\left(9.951\pm 0.055\right)\times 10^{-2}$
$7.12\times 10^{-2}$ | $\left(4.718\pm 0.025\right)\times 10^{0}$ | $\left(1.721\pm 0.008\right)\times 10^{-2}$ | $\left(1.373\pm 0.015\right)\times 10^{0}$ | $\left(5.596\pm 0.047\right)\times 10^{-2}$ | $\left(4.925\pm 0.030\right)\times 10^{-1}$ | $\left(1.196\pm 0.006\right)\times 10^{-1}$
$8.90\times 10^{-2}$ | $\left(5.245\pm 0.030\right)\times 10^{0}$ | $\left(1.997\pm 0.013\right)\times 10^{-2}$ | $\left(1.499\pm 0.011\right)\times 10^{0}$ | $\left(6.467\pm 0.041\right)\times 10^{-2}$ | $\left(5.723\pm 0.028\right)\times 10^{-1}$ | $\left(1.410\pm 0.006\right)\times 10^{-1}$
$1.78\times 10^{-1}$ | $\left(7.097\pm 0.061\right)\times 10^{0}$ | $\left(3.125\pm 0.021\right)\times 10^{-2}$ | $\left(2.116\pm 0.011\right)\times 10^{0}$ | $\left(1.024\pm 0.006\right)\times 10^{-1}$ | $\left(8.572\pm 0.058\right)\times 10^{-1}$ | $\left(2.329\pm 0.012\right)\times 10^{-1}$
$2.67\times 10^{-1}$ | $\left(8.781\pm 0.067\right)\times 10^{0}$ | $\left(4.131\pm 0.015\right)\times 10^{-2}$ | $\left(2.676\pm 0.025\right)\times 10^{0}$ | $\left(1.331\pm 0.004\right)\times 10^{-1}$ | $\left(1.175\pm 0.005\right)\times 10^{0}$ | $\left(3.225\pm 0.016\right)\times 10^{-1}$
$4.45\times 10^{-1}$ | $\left(1.126\pm 0.009\right)\times 10^{1}$ | $\left(5.784\pm 0.025\right)\times 10^{-2}$ | $\left(3.551\pm 0.038\right)\times 10^{0}$ | $\left(1.893\pm 0.010\right)\times 10^{-1}$ | $\left(1.715\pm 0.010\right)\times 10^{0}$ | $\left(4.813\pm 0.020\right)\times 10^{-1}$
$7.12\times 10^{-1}$ | $\left(1.668\pm 0.013\right)\times 10^{1}$ | $\left(7.138\pm 0.037\right)\times 10^{-2}$ | $\left(5.778\pm 0.043\right)\times 10^{0}$ | $\left(2.405\pm 0.014\right)\times 10^{-1}$ | $\left(2.992\pm 0.015\right)\times 10^{0}$ | $\left(6.917\pm 0.022\right)\times 10^{-1}$
$8.90\times 10^{-1}$ | $\left(2.159\pm 0.016\right)\times 10^{1}$ | $\left(7.651\pm 0.048\right)\times 10^{-2}$ | $\left(7.562\pm 0.055\right)\times 10^{0}$ | $\left(2.582\pm 0.016\right)\times 10^{-1}$ | $\left(4.129\pm 0.028\right)\times 10^{0}$ | $\left(7.969\pm 0.038\right)\times 10^{-1}$
$1.78\times 10^{0}$ | $\left(6.249\pm 0.036\right)\times 10^{1}$ | $\left(7.071\pm 0.037\right)\times 10^{-2}$ | $\left(2.410\pm 0.017\right)\times 10^{1}$ | $\left(2.501\pm 0.016\right)\times 10^{-1}$ | $\left(1.543\pm 0.012\right)\times 10^{1}$ | $\left(9.461\pm 0.055\right)\times 10^{-1}$
$2.67\times 10^{0}$ | $\left(1.413\pm 0.006\right)\times 10^{2}$ | $\left(6.203\pm 0.032\right)\times 10^{-2}$ | $\left(5.371\pm 0.043\right)\times 10^{1}$ | $\left(2.535\pm 0.014\right)\times 10^{-1}$ | $\left(3.402\pm 0.017\right)\times 10^{1}$ | $\left(1.002\pm 0.004\right)\times 10^{0}$
$4.45\times 10^{0}$ | $\left(3.768\pm 0.020\right)\times 10^{2}$ | $\left(6.321\pm 0.039\right)\times 10^{-2}$ | $\left(1.489\pm 0.009\right)\times 10^{2}$ | $\left(2.545\pm 0.014\right)\times 10^{-1}$ | $\left(9.386\pm 0.056\right)\times 10^{1}$ | $\left(1.013\pm 0.004\right)\times 10^{0}$
Table 1: Asymptotic parallel and perpendicular mean free paths $\lambda_{\parallel}$ and $\lambda_{\perp}$ for different reduced rigidities $r_{\text{g}}/L_{\text{c}}$ and at the three turbulence levels $\eta=0.2,\,0.5,\,{\rm and}\,0.8$. $r_{\text{g}}/L_{\text{c}}$ | $\eta=1$ | $\lambda_{0,\parallel}[L_{\text{c}}]$ | $\lambda_{0,\perp}[L_{\text{c}}]$
---|---|---|---
$\lambda_{1}[L_{\text{c}}]$
$4.45\times 10^{-4}$ | $\left(1.449\pm 0.022\right)\times 10^{-2}$ | $\left(2.457\pm 0.012\right)\times 10^{-1}$ | $\left(8.571\pm 0.276\right)\times 10^{-3}$
$2.67\times 10^{-3}$ | $\left(2.950\pm 0.031\right)\times 10^{-2}$ | $\left(3.929\pm 0.009\right)\times 10^{-1}$ | $\left(1.330\pm 0.019\right)\times 10^{-2}$
$2.67\times 10^{-2}$ | $\left(1.035\pm 0.007\right)\times 10^{-1}$ | $\left(8.077\pm 0.089\right)\times 10^{-1}$ | $\left(4.216\pm 0.040\right)\times 10^{-2}$
$8.90\times 10^{-2}$ | $\left(2.418\pm 0.016\right)\times 10^{-1}$ | $\left(1.265\pm 0.011\right)\times 10^{0}$ | $\left(8.721\pm 0.089\right)\times 10^{-2}$
$2.67\times 10^{-1}$ | $\left(6.535\pm 0.038\right)\times 10^{-1}$ | $\left(2.035\pm 0.024\right)\times 10^{0}$ | $\left(1.706\pm 0.011\right)\times 10^{-1}$
$8.90\times 10^{-1}$ | $\left(3.415\pm 0.021\right)\times 10^{0}$ | $\left(4.329\pm 0.060\right)\times 10^{0}$ | $\left(2.875\pm 0.021\right)\times 10^{-1}$
$2.67\times 10^{0}$ | $\left(2.896\pm 0.011\right)\times 10^{1}$ | $\left(2.687\pm 0.031\right)\times 10^{1}$ | $\left(2.545\pm 0.012\right)\times 10^{-1}$
$4.45\times 10^{0}$ | $\left(7.944\pm 0.020\right)\times 10^{1}$ | $\left(7.240\pm 0.077\right)\times 10^{1}$ | $\left(2.550\pm 0.013\right)\times 10^{-1}$
Table 2: Asymptotic mean free paths $\lambda_{1}$ for different reduced
rigidities $r_{\text{g}}/L_{\text{c}}$ at $\eta=1$ and the fit parameters
$\lambda_{0,\parallel}$ and $\lambda_{0,\perp}$ for eqs. 58 and 59.

In Fig. 9 we show the ratio of the perpendicular and parallel mean free paths,
$\lambda_{\perp}/\lambda_{\parallel}$ as functions of rigidity, for the
different turbulence levels $\eta=0.2,0.5$ and $0.8$. In the range between
$r_{\text{g}}/L_{\text{c}}\simeq 10^{-2}$ and $1$, the ratio is increasing
with a behaviour close to the $(r_{\text{g}}/L_{\text{c}})^{0.2}$ suggested by
Dundovic et al. Dundovic  _et al._ (2020). This is due to the faster than
$(r_{\text{g}}/L_{\text{c}})^{1/3}$-dependence of $\lambda_{\perp}$ at these
rigidities. At lower rigidities, however, the ratio becomes constant, meaning
that both $\lambda_{\parallel}$ and $\lambda_{\perp}$ scale in the same way.
At rigidities $\gtrsim 1$, the ratio drops steeply due to fact that
$\lambda_{\perp}$ becomes constant while
$\lambda_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{2}$.

In Fig. 10, we have plotted the asymptotic parallel and perpendicular mean
free paths $\lambda_{\parallel}$ and $\lambda_{\perp}$ as a function of
turbulence level $\eta$ for a subset of the reduced rigidities. We have also
investigated the scaling with $\delta B^{2}/B_{0}^{2}$ that had been
anticipated at small turbulence levels, see eqs. (1) and (2). As $\eta\to 1$,
of course, both $\lambda_{\parallel}$ and $\lambda_{\perp}$ must coincide,
$\lambda_{\parallel},\lambda_{\perp}\to\lambda_{1}$. A simple way to achieve
this scaling is to extend the ($\delta B^{2}/B_{0}^{2}$)-dependence of eqs.
(1) and (2) to

$\displaystyle\lambda_{\parallel}$
$\displaystyle=\left(\lambda_{1}+\lambda_{0,\parallel}\left(\frac{\delta
B^{2}}{B_{0}^{2}}\right)^{-1}\right)\,,$ (58) $\displaystyle\lambda_{\perp}$
$\displaystyle=\left(\lambda_{1}^{-1}+\lambda_{0,\perp}^{-1}\left(\frac{\delta
B^{2}}{B_{0}^{2}}\right)^{-1}\right)^{-1}\,.$ (59)

For those reduced rigidities, for which we have simulation data at the
turbulence level $\eta=1$, we know $\lambda_{1}$ and hence, we have only the
free parameters $\lambda_{0,\parallel}$ and $\lambda_{0,\perp}$ for each
reduced rigidity. We have fitted eqs. (58) and (59) to the simulation data and
show the resulting fits in Fig. 10. We find that eqs. (58) and (59) give an
excellent representation of the data.

Our asymptotic mean free paths are provided in Tbl. 1 (for $\eta=0.2,0.5$ and
$0.8$) and in Tbl. 2 (for $\eta=1$).

###  III.4 Field line running diffusion coefficients

Figure 11: Running field line diffusion coefficients $d_{\text{FL}}(s)$ with
respect to arc length $s$ (purple) and $d_{\text{FL}}(z)$ with respect to the
$z$-coordinate (gray) for a turbulence level $\eta=0.091$. The lines show the
mean over the ensemble of turbulent magnetic fields and the shaded bands
indicate the standard mean error. Figure 12: Same as Fig. 12, but for a
turbulence level $\eta=0.5$.

In Fig. 11, we show the running field line diffusion coefficients
$d_{\text{FL},s}(s)$ with respect to arc length $s$ and $d_{\text{FL},z}(z)$
with respect to the $z$-coordinate for a turbulence level $\eta=0.091$. At
this low turbulence level, the field lines are approximately aligned with the
$z$-direction and so the difference between the definition with respect to arc
length $s$ or $z$-coordinate are very small. The two curves seem to differ
mostly in the level of fluctuations, but this is caused by the fact that we
needed to smoothen the $d_{\text{FL}}$ defined with respect to $z$ for
technical reasons, see Sec. II.3.3. As a function of $s$ and $z$,
respectively, both running diffusion coefficients grow close to $\propto s$
and $\propto z$, initially, that is for $s\ll L_{\text{c}}$ and $z\ll
L_{\text{c}}$. Once the correlation length has been exceeded, the running
diffusion coefficients start to grow more slowly. However, the diffusive
regime is not attained immediately, but rather via some sub-diffusive phase
between $\sim L_{\text{c}}$ and $\sim 10L_{\text{c}}$. The diffusive field-
line wandering regime is finally attained for very large $s$ or $z$, that is
beyond $\sim 100L_{\text{c}}$. We note that a similar subdiffusive phase can
be observed in the field-line simulations of Sonsrettee et al. (2016)
Sonsrettee  _et al._ (2016) for the case of a small turbulence level, see
their Fig. 14. We caution however, that a somewhat different definition of the
running field line diffusion coefficient was adopted there, so their values
cannot be compared directly with our values.

At larger turbulence levels, the differences between the definitions of the
field line diffusion coefficients with respect to arc length $s$ and
$z$-coordinate become more pronounced. This is visible in Fig. 12, where we
show the running field line diffusion coefficients $d_{\text{FL},s}(s)$ with
respect to arc length $s$ and $d_{\text{FL},z}(z)$ with respect to the
$z$-coordinate for a turbulence level $\eta=0.5$. We remind the reader that
$\eta=0.5$ means that the turbulent field $\bm{\delta B}$ has the same energy
density as the background field $\bm{B}_{0}$. Qualitatively, both diffusion
coefficients show a similar behaviour as at low turbulence level, in
particular the field line transport is first ballistic, then sub-diffusive and
eventually diffusive. However, $d_{\text{FL,z}}$ is almost always larger in
value than $d_{\text{FL,s}}$. In particular, in the ballistic phase, it is
larger by almost an order of magnitude and in the late diffusive phase,
$d_{\text{FL,z}}\to\kappa_{\text{FL,z}}\simeq 0.1L_{\text{c}}$ while
$d_{\text{FL,s}}\to\kappa_{\text{FL,s}}\simeq 0.076L_{\text{c}}$. At least
qualitatively, this difference can be understood as follows: For finite
turbulence level $\eta$, the field lines are not generally aligned with the
$z$-direction. Part of the arc length $s$ is therefore accumulated also in the
$x$\- and $y$-directions and so $\langle(\Delta z)^{2}\rangle<(\Delta
s)^{2}=\langle(\Delta x)^{2}+(\Delta y)^{2}+(\Delta z)^{2}\rangle$. Therefore,
$\mathrm{d}\langle(\Delta x)^{2}\rangle/\mathrm{d}z>\mathrm{d}\langle(\Delta
x)^{2}\rangle/\mathrm{d}s$.

##  IV Systematic uncertainties

In the following, we present some results obtained by varying the methodology
used in the test particle simulations and in the computation of the diffusion
coefficients.

###  IV.1 Harmonic vs. Grid Method

Figure 13:  Comparison of the running diffusion coefficients using the
harmonic and the grid-based methods for setting up the turbulent magnetic
fields for $r_{\text{g}}/L_{\text{c}}=0.089$ and $\eta=0.5$. The parallel
(perpendicular) running diffusion coefficient are indicated by dashed (solid)
lines, their asymptotic values by dotted (dot-dashed) lines; shaded bands
indicate the respective standard mean errors. The asymptotic values have been
computed by averages beyond the time indicated by the dotted vertical lines.

Figure 14: Same as Fig. 14, but for $r_{\text{g}}/L_{\text{c}}=0.89$.

In Sec. II.2, we reviewed three methods for setting up the turbulent magnetic
field in test particle simulations: the harmonic method, the grid method and a
variant of the grid method, the so-called nested grid method which we employed
in our work. As discussed, only the nested grid method allows covering the
dynamical range needed to simulate particles at TV rigidities, but it requires
periodically replicating the various grids in order to cover the distances
that particles travel during the simulations. We have attempted to prevent the
introduction of possible artefacts due to this periodicity effect by allowing
for generous amounts of padding. Yet, it is desirable to compare the results
from the nested grid method with those from the harmonic method. Possible
deviations could be interpreted as a sign of systematic issues in either of
the methods.

However, given the computational expense of the harmonic method, such a
comparison is only possible at relatively large rigidities. Here, we perform
said check at reduced rigidities of $r_{\text{g}}/L_{\text{c}}=0.089$ and
$0.89$. While this is far from the low rigidities that are accessible in the
nested grid method, $r_{\text{g}}/L_{\text{c}}=0.089$ is in a range of
particular interest where the rigidity-dependence of the perpendicular
diffusion coefficient deviates from that of the parallel one, see, e.g. Fig.
9. Furthermore, at $r_{\text{g}}/L_{\text{c}}=0.089$, particles interact
resonantly with turbulent waves as witnessed by the rigidity-scaling of
$\lambda_{\parallel}$ close to $\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$, see
Fig. 8.

For the nested grid method, the same grid configurations as discussed in Sec.
II.2.3 were used. For the harmonic method, we have logarithmically spaced
$10^{5}$ modes between $k_{\text{min}}=2\pi/L_{\text{max}}$ and
$k_{\text{max}}=10^{5}k_{\text{min}}$. In Figs. 14 and 14 we compare the
running parallel and perpendicular diffusion coefficients.

The running diffusion coefficients show the same time-dependences and overall
they agree very well: At early times, $\Omega t\ll 1$, the running diffusion
coefficients grow ballistically. The parallel running diffusion coefficient
transitions to a diffusive behaviour beyond the diffusion time
$\tau_{\text{s}}^{\prime}$. The perpendicular running diffusion coefficient
instead experiences some subdiffusive behaviour before also becoming diffusive
for $t\gg\tau_{\text{c}}$. Note that for $r_{\text{g}}/L_{\text{c}}=0.89$
there are some features related to gyrations of particles around the effective
background field visible. The differences between the running diffusion
coefficients obtained with both methods almost always remain within the
uncertainties introduced by the sample variance which in Figs. 14 and 14 are
indicated by the shaded bands. The differences in the asymptotic diffusion
coefficients are $6\,\%$ and $7\,\%$ for $r_{\text{g}}/L_{\text{c}}=0.089$ and
$r_{\text{g}}/L_{\text{c}}=0.89$, respectively, which is larger than the
uncertainties since the asymptotic diffusion coefficients are based on an
average of the running diffusion coefficients over a large number of time
steps. Unfortunately, we have not been able to reach better convergence of
both methods and so we consider this as an estimate of the systematic
uncertainty of our test particle simulations. However, for the purpose of this
analysis, this level of accuracy is certainly sufficient.

Note that we expect periodicity effects, if present, to be even weaker at
smaller rigidities since for the maximum runtimes considered in our
simulations, particles generally travel over smaller distances. We thus
conclude that it is very unlikely that our results are affected by periodicity
effects due to the nested grid method.

###  IV.2 Differential Equation Solvers

Figure 15:  Comparison of the running diffusion coefficients computed with the
Boris-Push solver and an adaptive fourth-order Runge-Kutta solver. We have
adopted a reduced rigidity of $r_{\text{g}}/L_{\text{c}}=0.089$ and the
turbulence level $\eta=0.1$. In the top panel, the running parallel diffusion
coefficient $d_{\parallel}$ is shown, in the bottom panel the running parallel
diffusion coefficient $d_{\perp}$. The colours blue, orange and green refer to
three turbulent field realisations, with the results for the Boris-Push method
shown with solid lines and dark colours and those for the Runge-Kutta solve
with dashed lines and light colours. Note that the same random realisations
have been used for the Boris-Push and the Runge-Kutta methods. We have used a
linear ordinate axis in order to be able to appreciate the similarity between
the results for both methods.

Another potential source of systematic errors in test particle simulations is
the numerical solver for the equations of motion. In the literature, a number
of solvers is being employed, as discussed in Sec. II.2.4 and we have used the
Boris-Push method for all our results so far. Here, we compare the results for
running diffusion coefficients for one combination of rigidity and turbulence
level obtained with the Boris-Push method and with a fourth-order Runge-Kutta
solver.

Generally speaking, the adaptive Runge-Kutta method is computationally
significantly more expensive, in particular at small rigidities. We have
therefore performed the comparison at a reduced rigidity of
$r_{\text{g}}/L_{\text{c}}=0.089$ much larger than what can be achieved with
the Boris-Push method alone. While we have used the same realisations of the
turbulent magnetic field for both methods, the accumulation of numerical
errors leads to the divergence of individual particle trajectories with time
even if the same initial conditions are chosen. In order to delay this
divergence of trajectories, we have also chosen a rather low turbulence level
of $\eta=0.1$.

In Fig. 15, we show the results for the running parallel and perpendicular
diffusion coefficients, $d_{\parallel}$ and $d_{\perp}$ as a function of time.
The colours blue, orange and green refer to three different random
realisations of the turbulent magnetic field. The results for the Boris-Push
method are shown with solid lines and dark colours, those for the Runge-Kutta
method with dashed lines and light colours. We have chosen a linear ordinate
axis in order to better compare the results of both methods for individual
realisations of the turbulent field.

The running diffusion coefficients in both methods largely show the same
behaviour. While individual trajectories differ, as discussed above, the
running diffusion coefficients agree rather well. Also, the running diffusion
coefficients in _different_ realisations of the turbulent magnetic field agree
at early times (in the ballistic phase) and at late times (modulo the
numerical noise). At intermediate times, they exhibit the largest sensitivity
to the particular realisation of the turbulent magnetic field.

Besides the computational expense, the Runge-Kutta method also does not
conserve particle energy. We know, however, that in the magnetostatic setup,
particle energy is conserved. We have checked the loss of particle energy in
the Runge-Kutta method by tracking the average velocity magnitude in units of
the speed of light, $|\bm{v}|/c$. For low reduced rigidity of
$r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-4}$ and for $\eta=0.5$ we need to
run for several hundred thousands of gyro times in order for particle
transport to become diffusive. For these input values, we have found that
$|\bm{v}|/c$ has on average decreased by more than $15\,\%$ for $\Omega
t=10^{5}$ already. Note that in the Boris-Push method, particle energy is
conserved by construction.

In conclusion, we have confirmed that the Boris-Push and the Runge-Kutta
methods give approximately the same running diffusion coefficients at
intermediate particle rigidities. As the Runge-Kutta method is computationally
more expensive and and does not conserve particle energy, we have adopted the
Boris-Push method for all the results in this paper.

###  IV.3 Derivative vs fractional definition of running diffusion
coefficients

Figure 16: Comparison of the running diffusion coefficients obtained with the
derivative and fractional definitions, see eqs. (33) and (41).

Finally, another systematic error in the running diffusion coefficients is due
to the use of different definitions. We have used the derivative definitions,
see eqs. (32) and (33) throughout. However, in the literature, oftentimes the
fractional definition, see eqs. (40) and (41) is employed.

In Fig. 16, we compare the results for one combination of input parameters,
that is $r_{\text{g}}/L_{\text{c}}=8.9\times 10^{-3}$ and $\eta=0.5$. At early
times, that is in the ballistic phase, the running diffusion coefficients
following from the derivative and fractional definitions differ by about a
factor two as expected for mean-square displacements that are $\propto t^{2}$.
At late times, that is in the diffusive limit where the mean-square
displacements are $\propto t$, the definitions agree, of course. However, note
that the convergence to the diffusive limit takes place later and more
gradually for the fractional definition compared to the derivative definition.
In practice, this means that the code needs to be run for longer times to
achieve convergence of the running diffusion coefficients. While the nuisance
numerical noise at late times is suppressed for the fractional definition,
some of the interesting features at intermediate times are also less
pronounced. Finally, the slope in the running perpendicular diffusion
coefficient at intermediate times also differ.

In summary, the derivative definition of the diffusion coefficient, besides
being the one that allows to relate more closely to the TGK formalism also
offers faster convergence and allows identifying features at intermediate
times.

##  V Discussion

In the following, we will consider some examples where the difference in
scaling of the perpendicular mean free path has important phenomenological
consequences.

In the heliosphere, the large-scale field is often times modelled as a Parker
spiral. If transport would take place only along field lines, spatially
localised events like coronal mass ejections (CMEs) could only be observed in
charged particles if the acceleration site was magnetically connected to the
observer. The observation of solar energetic particles in cases where it is
not connected therefore implies some level of cross-field diffusion. We
refrain from detailed comparison with data, simulations or analytical
theories, since the isotropic turbulence model assumed here is not the one
preferred by the space science community. For those models, that is for a
combination of 2D and slab (composite) turbulence, it is usually claimed that
analytical theories and simulations agree well. Dundovic et al. Dundovic  _et
al._ (2020) suggest however that such an agreement might be limited by the
choice of the 2D to slab power ratio. In fact, analytical theories for
composite turbulence seem to overestimate the ratio
$\lambda_{\perp}/\lambda_{\parallel}$ observed from simulations for the case
where the 2D to slab power ratio is 50%/50% and the predicted rigidity
dependence of $\lambda_{\perp}/\lambda_{\parallel}$ is also different from the
simulated one, especially in the intermediate reduced rigidity range. Since a
more intense slab turbulence would imply stronger fluctuations in the magnetic
field component perpendicular to the background field, one could expect that,
in this case, the effect of the field-line transport might be important in
determining the ratio $\lambda_{\perp}/\lambda_{\parallel}$. In this respect,
subdiffusive episodes in field-line transport could be crucial for a more
comprehensive picture of particle transport in the heliosphere.

Perpendicular transport is also central in shock acceleration at quasi-
perpendicular shocks. For supernova remnants (SNRs), there is consensus that
magnetic fields need to be amplified, in order for shock-accelerated particles
to reach the knee energy. In that situation the turbulence spectrum would
likely differ from the Kolmogorov form assumed here. However, this does not
need to be the case for all SNRs, see e.g. Ferrand  _et al._ (2014).
Interplanetary shocks, e.g. in CMEs, would be another acceleration site where
the acceleration rate and hence the maximum energy is regulated by
perpendicular transport Dosch and Shalchi (2010).

In the remainder of this section, we will confine ourselves to considering the
implications of our findings for the transport of Galactic CRs in more detail,
specifically of the difference in scaling of the perpendicular diffusion
coefficient $\kappa_{\perp}$ with reduced rigidity
$r_{\text{g}}/L_{\text{c}}$. Given the relatively small difference of $\sim
0.2$ in spectral index, one might wonder what such implications could be. For
example, for $\eta=0.5$, if they agree at $r_{\text{g}}/L_{\text{c}}=3\times
10^{-3}$ (which corresponds to a rigidity of $300\,\text{TV}$ for our
benchmark values), they differ by about an order of magnitude at
$r_{\text{g}}/L_{\text{c}}=10^{-8}$ (which corresponds to $1\,\text{GV}$). In
the following, we will consider the competition of perpendicular escape from a
region of size $L$ with three processes: parallel escape, advection and
cooling losses.

##### Diffusive escape

First, we consider competition between parallel and perpendicular escape. To
quantify this, we define the rates $t_{\parallel}\equiv
L_{\parallel}^{2}/(2\kappa_{\parallel})$ and $t_{\perp}\equiv
L_{\perp}^{2}/(2\kappa_{\perp})$. Whether parallel or perpendicular loss
dominates depends on the aspect ratio of the system,
$L_{\perp}/L_{\parallel}$. Perpendicular escape is dominant if
$\kappa_{\perp}/\kappa_{\parallel}>(L_{\perp}/L_{\parallel})^{2}$.

For instance for the Galaxy, we estimate
$L_{\perp}/L_{\parallel}=0.1\mathellipsis 0.3$. Therefore, if the relative
scaling had been
$\kappa_{\perp}/\kappa_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{0.2}$ at
all rigidities, perpendicular transport could have been important only at the
highest rigidities, if at all. However, we have found that
$\kappa_{\perp}/\kappa_{\parallel}=\text{const.}$ Specifically, for
$\eta=0.5$, $\kappa_{\perp}/\kappa_{\parallel}=2\times 10^{-2}$ below
$r_{\text{g}}/L_{\text{c}}=3\times 10^{-3}$. Therefore, perpendicular escape
can be important at all energies below $300\,\text{TV}$.

##### Advective escape

Next, we compare perpendicular diffusive escape with the time scale
$t_{\perp}\equiv L_{\perp}^{2}/(2\kappa_{\perp})$ with escape due to advection
at the speed $U$. We characterise this process with the time scale
$t_{\text{adv}}\equiv L_{\perp}/U$. As $t_{\perp}$ increases with decreasing
rigidity $r_{\text{g}}/L_{\text{c}}$, advection will dominate over diffusion
below some critical rigidity $r_{\text{g,a}}/L_{\text{c}}$ which is defined
through $t_{\text{adv}}=t_{\perp}(r_{\text{g,a}}/L_{\text{c}})$. Note that
$r_{\text{g,a}}/L_{\text{c}}$ will sensitively depend on whether
$\kappa_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ or
$(r_{\text{g}}/L_{\text{c}})^{1/2}$.

For a numerical estimate, we fix the perpendicular diffusion coefficient at
$r_{\text{g}}/L_{\text{c}}=3\times 10^{-3}$ to
$\kappa_{\perp}=(v/3)\lambda_{\perp}=(v/3)10^{-2}\,L_{\text{c}}$. We further
adopt an advection speed $U=2\,\text{km}\,\text{s}^{-1}\simeq(2/3)10^{5}\,v$
with the particle speed $v$ set to $c$ and set $L_{\perp}=100\,L_{\text{c}}$.
Solving for $r_{\text{g,a}}/L_{\text{c}}$, we find
$r_{\text{g,a}}/L_{\text{c}}=3\times 10^{-6}$. Had we assumed
$\kappa_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{1/2}$, this would be
$r_{\text{g,a}}/L_{\text{c}}=3\times 10^{-5}$ instead, that is an order of
magnitude larger.

##### Cooling losses

Finally, we consider the transition between transport dominated by
perpendicular losses with the loss time $t_{\perp}$ and cooling by synchrotron
or Inverse Compton scattering. This is an important phenomenological question
when considering the transport of CR electrons and positrons Mertsch (2018);
Lipari (2019). In the Thomson regime and assuming relativistic particles, the
cooling time scales like
$t_{\text{cool}}\sim(r_{\text{g}}/L_{\text{c}})^{-1}$, hence diffusive losses
will dominate a low rigidities, cooling losses will dominate at high
rigidities. We define the transition to take place at
$r_{\text{g,c}}/L_{\text{c}}$ defined through
$t_{\text{cool}}(r_{\text{g,c}}/L_{\text{c}})=t_{\perp}(r_{\text{g,c}}/L_{\text{c}})$.
We adopt some benchmark numbers for a numerical estimate, that is we demand
$t_{\text{cool}}(\mathcal{R}=1\,\text{GV})=3\times 10^{16}\,\text{s}$ that is
$t_{\text{cool}}(r_{\text{g}}/L_{\text{c}})=10^{7}L_{\text{c}}/c\left[(r_{\text{g}}/L_{\text{c}})/10^{-8}\right]^{-1}$
for our benchmark magnetic field parameters. Thus for the
$\kappa_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{1/2}$ scaling, the
transition takes place at $r_{\text{g,c}}/L_{\text{c}}=1.5\times 10^{-12}$
whereas $r_{\text{g,c}}/L_{\text{c}}=3.1\times 10^{-10}$ for a
$\kappa_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ scaling, a difference
of a factor $\sim 200$.

##  VI Conclusion

Perpendicular transport of high-energy particles is important in a number of
environments and a sound theoretical understanding is important when
interpreting observations, be it of in-situ observations in the heliosphere,
studies of Galactic CRs or non-thermal emission from sources. Our simulations
have shown that the rigidity-dependence of the perpendicular mean free path
$\lambda_{\perp}$ differs from that of the parallel mean free path
$\lambda_{\parallel}$ for reduced rigidities
$r_{\text{g}}/L_{\text{c}}\lesssim 1$; at even lower reduced rigidities
$r_{\text{g}}/L_{\text{c}}\ll 1$, however, the perpendicular diffusion
coefficient returns back to the same scaling. Specifically, for Kolmogorov
turbulence, $\lambda_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{0.5}$ for
$r_{\text{g}}/L_{\text{c}}\lesssim 1$ and
$\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ for $r_{\text{g}}/L_{\text{c}}\ll
1$ whereas $\lambda_{\parallel}\propto(r_{\text{g}}/L_{\text{c}})^{1/3}$ for
all $r_{\text{g}}/L_{\text{c}}\lesssim 1$. Previous analyses had instead
speculated about the $\lambda_{\perp}\propto(r_{\text{g}}/L_{\text{c}})^{0.5}$
behaviour extending to the lowest rigidities. In a companion paper Kuhlen  _et
al._ (2022b), we have provided an analytical model that is able to reproduce
this scaling as long as the subdiffusive phase in the running field line
diffusion coefficient is taken into account.

###### Acknowledgements.

The authors would like to thank Carmelo Evoli for helpful discussions. This
project was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – project number 426614101.

## References

  * Gabici  _et al._ (2019) S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, International Journal of Modern Physics D 28, 1930022-339 (2019), arXiv:1903.11584 [astro-ph.HE] .
  * Jokipii (1966) J. R. Jokipii, ApJ 146, 480 (1966).
  * Kennel and Engelmann (1966) C. F. Kennel and F. Engelmann, Physics of Fluids 9, 2377 (1966).
  * Hall and Sturrock (1967) D. E. Hall and P. A. Sturrock, Physics of Fluids 10, 2620 (1967).
  * Hasselmann and Wibberenz (1970) K. Hasselmann and G. Wibberenz, ApJ 162, 1049 (1970).
  * Note (1) Rigidity $\mathcal{R}$ is defined as the ratio of particle momentum $p$ and charge $q$, that is $\mathcal{R}\equiv pc/q$ with the speed of light $c$.
  * Kolmogorov (1931) A. Kolmogorov, Math. Ann. 104, 415 (1931).
  * Voelk (1975) H. J. Voelk, Reviews of Geophysics and Space Physics 13, 547 (1975).
  * Shalchi (2009) A. Shalchi,  _Nonlinear Cosmic Ray Diffusion Theories_, Vol. 362 (2009).
  * Tautz  _et al._ (2006a) R. C. Tautz, A. Shalchi, and R. Schlickeiser, Journal of Physics G Nuclear Physics 32, 809 (2006a).
  * Shalchi (2020) A. Shalchi, Space Sci. Rev. 216, 23 (2020).
  * Tautz  _et al._ (2006b) R. C. Tautz, A. Shalchi, and R. Schlickeiser, Journal of Physics G Nuclear Physics 32, 1045 (2006b).
  * Bieber and Matthaeus (1997) J. W. Bieber and W. H. Matthaeus, ApJ 485, 655 (1997).
  * Matthaeus  _et al._ (2003) W. H. Matthaeus, G. Qin, J. W. Bieber, and G. P. Zank, ApJ 590, L53 (2003).
  * Shalchi  _et al._ (2004) A. Shalchi, J. W. Bieber, W. H. Matthaeus, and G. Qin, ApJ 616, 617 (2004).
  * Beresnyak  _et al._ (2011) A. Beresnyak, H. Yan, and A. Lazarian, ApJ 728, 60 (2011), arXiv:1002.2646 [astro-ph.GA] .
  * Xu and Yan (2013) S. Xu and H. Yan, ApJ 779, 140 (2013), arXiv:1307.1346 [astro-ph.HE] .
  * Cohet and Marcowith (2016) R. Cohet and A. Marcowith, A&A 588, A73 (2016), arXiv:1601.04971 [astro-ph.HE] .
  * Giacalone and Jokipii (1994) J. Giacalone and J. R. Jokipii, ApJ 430, L137 (1994).
  * Michalek and Ostrowski (1997) G. Michalek and M. Ostrowski, A&A 326, 793 (1997).
  * Michalek and Ostrowski (1998) G. Michalek and M. Ostrowski, A&A 337, 558 (1998).
  * Giacalone and Jokipii (1999) J. Giacalone and J. R. Jokipii, ApJ 520, 204 (1999).
  * Casse  _et al._ (2001) F. Casse, M. Lemoine, and G. Pelletier, Phys. Rev. D 65, 023002 (2001), arXiv:astro-ph/0109223 [astro-ph] .
  * Qin  _et al._ (2002) G. Qin, W. H. Matthaeus, and J. W. Bieber, Geophys. Res. Lett. 29, 1048 (2002).
  * DeMarco  _et al._ (2007) D. DeMarco, P. Blasi, and T. Stanev, J. Cosmology Astropart. Phys 2007, 027 (2007), arXiv:0705.1972 [astro-ph] .
  * Hussein  _et al._ (2015) M. Hussein, R. C. Tautz, and A. Shalchi, Journal of Geophysical Research (Space Physics) 120, 4095 (2015), arXiv:1505.05099 [astro-ph.SR] .
  * Giacinti  _et al._ (2018) G. Giacinti, M. Kachelrieß, and D. V. Semikoz, J. Cosmology Astropart. Phys 2018, 051 (2018), arXiv:1710.08205 [astro-ph.HE] .
  * Reichherzer  _et al._ (2020) P. Reichherzer, J. Becker Tjus, E. G. Zweibel, L. Merten, and M. J. Pueschel, MNRAS 498, 5051 (2020), arXiv:1910.07528 [astro-ph.HE] .
  * Mertsch (2020) P. Mertsch, Ap&SS 365, 135 (2020), arXiv:1910.01172 [astro-ph.HE] .
  * Melrose (1980) D. B. Melrose,  _Plasma astrophysics: Nonthermal processes in diffuse magnetized plasmas. Volume 1 - The emission, absorption and transfer of waves in plasmas_ (1980).
  * Snodin  _et al._ (2016) A. P. Snodin, A. Shukurov, G. R. Sarson, P. J. Bushby, and L. F. S. Rodrigues, MNRAS 457, 3975 (2016), arXiv:1509.03766 [astro-ph.HE] .
  * Dundovic  _et al._ (2020) A. Dundovic, O. Pezzi, P. Blasi, C. Evoli, and W. H. Matthaeus, Phys. Rev. D 102, 103016 (2020), arXiv:2007.09142 [astro-ph.HE] .
  * Evoli  _et al._ (2012) C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Phys. Rev. Lett. 108, 211102 (2012), arXiv:1203.0570 [astro-ph.HE] .
  * Cerri  _et al._ (2017) S. S. Cerri, D. Gaggero, A. Vittino, C. Evoli, and D. Grasso, J. Cosmology Astropart. Phys 2017, 019 (2017), arXiv:1707.07694 [astro-ph.HE] .
  * Tautz (2016) R. C. Tautz, New A 45, 1 (2016).
  * Debye and Hückel (1923) P. Debye and E. Hückel, Physikalische Zeitschrift 24, 185 (1923).
  * Ferrière (2001) K. M. Ferrière, Reviews of Modern Physics 73, 1031 (2001), arXiv:astro-ph/0106359 [astro-ph] .
  * Giacinti  _et al._ (2014) G. Giacinti, M. Kachelrieß, and D. V. Semikoz, Phys. Rev. D 90, 041302 (2014), arXiv:1403.3380 [astro-ph.HE] .
  * Savchenko  _et al._ (2015) V. Savchenko, M. Kachelrieß, and D. V. Semikoz, Astrophys. J. Lett. 809, L23 (2015), arXiv:1505.02720 [astro-ph.HE] .
  * Génolini and Ahlers (2021) Y. Génolini and M. Ahlers, PoS ICRC2021, 160 (2021).
  * Giacinti and Sigl (2012) G. Giacinti and G. Sigl, Phys. Rev. Lett. 109, 071101 (2012), arXiv:1111.2536 [astro-ph.HE] .
  * Ahlers and Mertsch (2015) M. Ahlers and P. Mertsch, ApJ 815, L2 (2015), arXiv:1506.05488 [astro-ph.HE] .
  * López-Barquero  _et al._ (2016) V. López-Barquero, R. Farber, S. Xu, P. Desiati, and A. Lazarian, ApJ 830, 19 (2016), arXiv:1509.00892 [astro-ph.HE] .
  * Ahlers and Mertsch (2017) M. Ahlers and P. Mertsch, Progress in Particle and Nuclear Physics 94, 184 (2017), arXiv:1612.01873 [astro-ph.HE] .
  * Mertsch and Ahlers (2019) P. Mertsch and M. Ahlers, J. Cosmology Astropart. Phys 2019, 048 (2019), arXiv:1909.09052 [astro-ph.HE] .
  * Kuhlen  _et al._ (2022a) M. Kuhlen, V. H. M. Phan, and P. Mertsch, ApJ 927, 110 (2022a), arXiv:2111.06419 [astro-ph.HE] .
  * Schekochihin (2020) A. A. Schekochihin, arXiv e-prints , arXiv:2010.00699 (2020), arXiv:2010.00699 [physics.plasm-ph] .
  * Batchelor and Press (1953) G. Batchelor and C. U. Press,  _The Theory of Homogeneous Turbulence_, Cambridge Science Classics (Cambridge University Press, 1953).
  * Matthaeus and Smith (1981) W. H. Matthaeus and C. Smith, Phys. Rev. A 24, 2135 (1981).
  * Kolmogorov (1941) A. Kolmogorov, Akademiia Nauk SSSR Doklady 30, 301 (1941).
  * Harari  _et al._ (2002) D. Harari, S. Mollerach, E. Roulet, and F. Sánchez, Journal of High Energy Physics 2002, 045 (2002), arXiv:astro-ph/0202362 [astro-ph] .
  * Chepurnov  _et al._ (2010) A. Chepurnov, A. Lazarian, S. Stanimirović, C. Heiles, and J. E. G. Peek, ApJ 714, 1398 (2010), arXiv:astro-ph/0611462 [astro-ph] .
  * Haverkorn  _et al._ (2008) M. Haverkorn, J. C. Brown, B. M. Gaensler, and N. M. McClure-Griffiths, ApJ 680, 362 (2008), arXiv:0802.2740 [astro-ph] .
  * Tautz and Dosch (2013) R. C. Tautz and A. Dosch, Physics of Plasmas 20, 022302 (2013), arXiv:1301.7164 [astro-ph.HE] .
  * Dongarra and Sullivan (2000) J. Dongarra and F. Sullivan, Computing in Science and Engineering 2, 22 (2000).
  * Johnson and Frigo (2008a) S. G. Johnson and M. Frigo, in  _Fast Fourier Transforms_, edited by C. S. Burrus (Connexions, Rice University, Houston TX, 2008) Chap. 11.
  * Giacinti  _et al._ (2012) G. Giacinti, M. Kachelrieß, D. V. Semikoz, and G. Sigl, J. Cosmology Astropart. Phys 2012, 031 (2012), arXiv:1112.5599 [astro-ph.HE] .
  * Cooley and Tukey (1965) J. W. Cooley and J. W. Tukey, Math. Comput. 19, 297 (1965).
  * Johnson and Frigo (2008b) S. G. Johnson and M. Frigo, in  _Fast Fourier Transforms_ , edited by C. S. Burrus (Rice University, Houston TX, 2008) [Online]. Available: http://cnx.org/content/m16336/.
  * Nyquist (1928) H. Nyquist, Transactions of the American Institute of Electrical Engineers 47, 617 (1928).
  * Condon and Ransom (2016) J. J. Condon and S. M. Ransom,  _Essential Radio Astronomy_ (2016).
  * Press  _et al._ (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in C_ , 2nd ed. (Cambridge University Press, Cambridge, USA, 1992).
  * Runge (1895) C. Runge, Math. Ann. 46, 167 (1895).
  * Kutta (1901) M. Kutta, Zeitschrift für Mathematik und Physik 46, 435 (1901).
  * Ahnert and Mulansky (2011) K. Ahnert and M. Mulansky, in  _Numerical Analysis and Applied Mathematics ICNAAM 2011: International Conference on Numerical Analysis and Applied Mathematics_, American Institute of Physics Conference Series, Vol. 1389, edited by T. E. Simos, G. Psihoyios, C. Tsitouras, and Z. Anastassi (2011) pp. 1586–1589, arXiv:1110.3397 [cs.MS] .
  * Griffiths (1995) D. J. Griffiths,  _Introduction to Quantum Mechanics_ (1995).
  * Boris (1970) J. P. Boris, Proceeding of Fourth Conference on Numerical Simulations of Plasmas (1970).
  * Birdsall and Langdon (2005) C. K. Birdsall and A. B. Langdon,  _Plasma physics via computer simulation_ (Taylor and Francis, New York, 2005).
  * Qin  _et al._ (2013) H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M. Tang, Physics of Plasmas 20, 084503 (2013).
  * Higuera and Cary (2018) A. Higuera and J. Cary, in  _APS Division of Plasma Physics Meeting Abstracts_ , APS Meeting Abstracts, Vol. 2018 (2018) p. TO8.003.
  * Ripperda  _et al._ (2018) B. Ripperda, F. Bacchini, J. Teunissen, C. Xia, O. Porth, L. Sironi, G. Lapenta, and R. Keppens, ApJS 235, 21 (2018), arXiv:1710.09164 [astro-ph.IM] .
  * Taylor (1922) G. I. Taylor, Proceedings of the London Mathematical Society s2-20, 196 (1922), https://londmathsoc.onlinelibrary.wiley.com/doi/pdf/10.1112/plms/s2-20.1.196 .
  * Green (1954) M. S. Green, J. Chem. Phys. 22, 398 (1954).
  * Kubo (1957) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Fisher (1953) R. Fisher, Proceedings of the Royal Society of London Series A 217, 295 (1953).
  * Webb  _et al._ (2006) G. M. Webb, G. P. Zank, E. K. Kaghashvili, and J. A. le Roux, ApJ 651, 211 (2006).
  * Giacinti and Kirk (2017) G. Giacinti and J. G. Kirk, ApJ 835, 258 (2017), arXiv:1610.06134 [astro-ph.HE] .
  * López-Coto and Giacinti (2018) R. López-Coto and G. Giacinti, Mon. Not. Roy. Astron. Soc. 479, 4526 (2018), arXiv:1712.04373 [astro-ph.HE] .
  * Gruzinov (2018) A. Gruzinov, (2018), arXiv:1808.00041 [astro-ph.HE] .
  * Sonsrettee  _et al._ (2016) W. Sonsrettee, P. Subedi, D. Ruffolo, W. H. Matthaeus, A. P. Snodin, P. Wongpan, P. Chuychai, G. Rowlands, and S. Vyas, ApJS 225, 20 (2016).
  * Ferrand  _et al._ (2014) G. Ferrand, R. J. Danos, A. Shalchi, S. Safi-Harb, P. Edmon, and P. Mendygral, ApJ 792, 133 (2014), arXiv:1407.6728 [astro-ph.HE] .
  * Dosch and Shalchi (2010) A. Dosch and A. Shalchi, Advances in Space Research 46, 1208 (2010).
  * Mertsch (2018) P. Mertsch, JCAP 11, 045 (2018), arXiv:1809.05104 [astro-ph.HE] .
  * Lipari (2019) P. Lipari, Phys. Rev. D 99, 043005 (2019), arXiv:1810.03195 [astro-ph.HE] .
  * Kuhlen  _et al._ (2022b) M. Kuhlen, V. H. M. Phan, and P. Mertsch, arXiv e-prints (2022b).