# Diffusion of cosmic rays in MHD turbulence with magnetic mirrors

A. Lazarian  Department of Astronomy, University of Wisconsin, 475 North
Charter Street, Madison, WI 53706, USA  Center for Computation Astrophysics,
Flatiron Institute, 162 5th Ave, New York, NY 10010  Siyao Xu  Institute for
Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA111Hubble Fellow

###### Abstract

As the fundamental physical process with many astrophysical implications, the
diffusion of cosmic rays (CRs) is determined by their interaction with
magnetohydrodynamic (MHD) turbulence. We consider the magnetic mirroring
effect arising from MHD turbulence on the diffusion of CRs. Due to the
intrinsic superdiffusion of turbulent magnetic fields, CRs with large pitch
angles that undergo mirror reflection, i.e., bouncing CRs, are not trapped
between magnetic mirrors, but move diffusively along the turbulent magnetic
field, leading to a new type of parallel diffusion, i.e., mirror diffusion.
This mirror diffusion is in general slower than the diffusion of non-bouncing
CRs with small pitch angles that undergo gyroresonant scattering. The critical
pitch angle at the balance between magnetic mirroring and pitch-angle
scattering is important for determining the diffusion coefficients of both
bouncing and non-bouncing CRs and their scalings with the CR energy. We find
non-universal energy scalings of diffusion coefficients, depending on the
properties of MHD turbulence.

††software: MATLAB (MATLAB, 2018)

##  1 Introduction

Charged energetic particles or cosmic rays (CRs) are an important ingredient
in the physical processes in space and astrophysical environments. It is
customary to use the term “energetic particles” in space physics. The
theoretical understanding on their acceleration and diffusion in the Solar
atmosphere, solar wind, Earth magnetosphere, and heliosphere is important for
studying the properties of the interplanetary magnetic field, solar modulation
of Galactic CRs, and space weather forecasting (Parker, 1965; Jokipii, 1971;
Singer et al., 2001).

The energetic particles with higher energies outside our direct neighborhood,
i.e., of Galactic and extragalactic origin, are usually referred to as CRs.
The knowledge on the acceleration and diffusion of CRs is essential for
probing their sources, explaining their chemical composition, studying their
roles in ionizing molecular gas and circumstellar discs (e.g., Schlickeiser et
al. 2016; Padovani et al. 2018), driving galactic winds (e.g., Ipavich 1975;
Holguin et al. 2019), and feedback heating in clusters of galaxies (e.g., Guo
& Oh 2008; Brunetti & Jones 2014), as well as modeling the synchrotron
foreground emission for cosmic microwave background (CMB) radiation and
redshifted $21$ cm radiation (e.g., Cho & Lazarian 2002a; Cho et al. 2012). In
this work, we focus on the diffusion physics that is generally applicable to
energetic particles of Solar origin and CRs. Thus we do not distinguish
between them and will only use the term “CRs”.

The diffusion of CRs is governed by their interaction with turbulent magnetic
fields, as well as the magnetic fluctuations arising from the instabilities
induced by CRs, i.e., the streaming instability (Kulsrud & Pearce, 1969). In
this work we focus on the former type of interaction. The CRs interaction with
magnetic turbulence has been a subject of intensive research for decades
(Jokipii, 1966; Kulsrud & Pearce, 1969; Schlickeiser & Miller, 1998; Giacalone
& Jokipii, 1999). However, ad hoc models for magnetohydrodynamic (MHD)
turbulence were adopted in early studies. 222The model of isotropic MHD
turbulence (see Schlickeiser 2002) and 2D + slab superposition model of solar
wind turbulence (Matthaeus et al., 1990) adopted in early studies are in
contradiction with MHD turbulence simulations and are also challenged by solar
wind observations (e.g., Horbury et al. 2008).

The recently developed theories of MHD turbulence (Goldreich & Sridhar, 1995;
Lazarian & Vishniac, 1999; Lithwick & Goldreich, 2001) have been numerically
tested (Cho & Vishniac 2000; Maron & Goldreich 2001; Cho et al. 2002; Cho &
Lazarian 2003; Beresnyak 2014; Kowal & Lazarian 2010; Kowal et al. 2017, see
also the book by Beresnyak & Lazarian 2019), which are also supported by solar
wind observations (e.g., Horbury et al. 2008; Luo & Wu 2010; Forman et al.
2011). The development of MHD turbulence theories allows important advances in
studying the pitch-angle scattering, stochastic acceleration, spatial
diffusion of CRs along the magnetic field (Chandran, 2000a; Yan & Lazarian,
2002, 2004; Brunetti & Lazarian, 2007; Yan & Lazarian, 2008; Lynn et al.,
2012; Xu & Yan, 2013; Xu & Lazarian, 2018; Sioulas et al., 2020; Lemoine &
Malkov, 2020), superdiffusion and diffusion perpendicular to the mean magnetic
field (Yan & Lazarian, 2008; Lazarian & Yan, 2014), propagation of CRs in
weakly ionized astrophysical media (Xu et al., 2016), and interactions of CRs
with relativistic MHD turbulence (Demidem et al., 2020). These studies bring
significant changes to the standard diffusion models of CRs based on the ad
hoc models of MHD turbulence (e.g., Matthaeus et al. 1990; Kóta & Jokipii
2000; Qin et al. 2002) and shed light on some long-standing problems and
observational puzzles (e.g., Palmer 1982; Evoli & Yan 2014; López-Barquero et
al. 2016; Krumholz et al. 2020). Naturally, the propagation of CRs should be
modeled using the tested MHD turbulence theories in order to explain
multifrequency observations and direct CR measurements (e.g., Nava & Gabici
2013; Cohet & Marcowith 2016; Orlando 2018; Gabici et al. 2019; Amato &
Casanova 2021; Fornieri et al. 2021). We are still far from fully
understanding the diffusion of CRs to interpret the observations near the
Earth and in the vicinity of CR sources and their differences. This motivates
us to further study the fundamental physics of CR diffusion in this work.

For CRs interacting with magnetic irregularities, in addition to the pitch-
angle scattering, it has long been known that the CR propagation can also be
affected by magnetic mirror reflection (Fermi, 1949; Noerdlinger, 1968;
Cesarsky & Kulsrud, 1973; Klepach & Ptuskin, 1995; Chandran, 2000b). The
magnetic mirroring effect was explored, for instance, for solving the
$90^{\circ}$ problem of the quasilinear theory (QLT, Jokipii 1966) for pitch-
angle scattering. In these early studies, the mirroring effect was invoked for
trapping CRs that bounce back and forth between two mirror points, but it has
not been considered in the context of the numerically tested modern theories
of MHD turbulence.

In MHD turbulence, compressions of magnetic fields, which are generated by
pseudo-Alfvén modes in an incompressible medium and slow and fast modes in a
compressible medium, naturally give rise to the mirroring effect over a range
of length scales following the energy cascade of turbulence. The interaction
of CRs with magnetic mirrors is regulated by the dynamics of turbulent
magnetic fields. In particular, this work will demonstrate that the intrinsic
perpendicular superdiffusion of turbulent magnetic fields (Lazarian &
Vishniac, 1999; Eyink et al., 2013; Lazarian & Yan, 2014) and its interaction
with the parallel diffusion should be taken into account when studying the CR
diffusion subject to the mirroring effect.

By using the numerically tested scalings of MHD turbulence (Cho & Lazarian,
2003), Xu & Lazarian (2020) (hereafter XL20) investigated the scattering of
CRs with the mirroring effect included. There we confirmed the dominant role
of fast modes in gyroresonant scattering, and we found that the resulting
diffusion coefficient can be significantly smaller than that in the absence of
magnetic mirroring.

In this work, by taking into account the intrinsic dynamics of MHD turbulence,
i.e., magnetic field perpendicular superdiffusion, we investigate the effect
of magnetic mirroring on the parallel diffusion of CRs. This new diffusion
mechanism arising from the mirroring effect in MHD turbulence should be
considered together with other diffusion processes related to scattering and
streaming of CRs for a more comprehensive description of CR propagation. Here
we focus on the fundamental physics of the mirror diffusion mechanism and the
formulation of its diffusion coefficient. Its applicability to various
astrophysical media and its confrontation with observations will be studied in
our future work (see e.g., Xu 2021).

In what follows, in Section 2, we introduce the perpendicular superdiffusion
of both turbulent magnetic fields and CRs and its effect on the parallel
diffusion of bouncing CRs. In Sections 3 and 4, we formulate the diffusion
coefficients of bouncing CRs in compressible and incompressible MHD
turbulence, respectively. In Section 5, we consider the exchange between
bouncing and non-bouncing CRs and discuss the averaged diffusion coefficient
on scales much larger than their mean free paths. Finally, the discussion and
the summary of our main results are given in Sections 6 and 7.

##  2 Spatial diffusion of bouncing particles

The mirroring effect in MHD turbulence causes diffusion of CRs along turbulent
magnetic fields. This diffusion parallel to the local magnetic field, i.e.,
”parallel diffusion”, is affected by the superdiffusion of turbulent magnetic
fields in the direction perpendicular to the mean magnetic field. Below we
will discuss in detail these two types of diffusion and their relation.

###  2.1 Magnetic mirrors in MHD turbulence

MHD turbulence can be decomposed into Alfvén, slow, and fast modes (Goldreich
& Sridhar, 1995; Lithwick & Goldreich, 2001; Cho & Lazarian, 2002b, 2003).
Alfvén modes induce magnetic field wandering, which accounts for the
superdiffusion of CRs perpendicular to the mean magnetic field (Lazarian &
Vishniac, 1999; Eyink et al., 2013; Lazarian & Yan, 2014). Slow modes are
passively mixed by Alfvén modes and have the same anisotropic scaling as
Alfvén modes (Lithwick & Goldreich, 2001; Cho & Lazarian, 2003). Fast modes
have an independent energy cascade and isotropic energy distribution (Cho &
Lazarian, 2002b, 2003; Kowal et al., 2009). Magnetic compressions arising from
slow and fast modes can act as magnetic mirrors that result in bouncing of
particle among the mirror points. Their detailed statistical properties are
presented in Appendix.

The magnetic mirroring effect caused by static magnetic bottles is well known
in plasma physics (e.g., Post 1958; Budker 1959; Noerdlinger 1968; Kulsrud &
Pearce 1969). A particle with the Larmor radius smaller than the variation
scale of the magnetic field preserves its adiabatic invariant, i.e.
$p_{\bot}^{2}/B=$const. It implies

$\frac{p_{\bot}^{2}}{B_{0}}=\frac{p^{2}}{B_{0}+\delta b}$ (1)

for the condition of magnetic mirroring, where $B_{0}$ and $B_{0}+\delta b$
are the magnetic field strengths in the weak and strong magnetic field
regions, and $p$ is the total momentum of the particle. As the total momentum
is preserved, it is easy to see that the particles with $\mu<\mu_{lc}$, where
$\mu$ is the pitch-angle cosine and $\mu_{lc}$ satisfies

$\mu_{lc}^{2}=\cos^{2}\theta_{lc}=\frac{\delta b}{B_{0}+\delta b},$ (2)

are subject to magnetic mirroring. Here $\theta_{lc}$ is angular size of the
loss cone for escaping particles with smaller pitch angles. For slow and fast
modes in MHD turbulence with a spectrum of magnetic fluctuations, there is
$\delta b=b_{k}$. When $b_{k}$ at wavenumber $k$ is significantly smaller than
the mean magnetic field strength $B_{0}$, the above expression can be
approximated by

$\mu_{lc}^{2}\approx\frac{b_{k}}{B_{0}}.$ (3)

The amplitude of $b_{k}$ depends on the driven magnetic perturbations $\delta
B_{s}$ and $\delta B_{f}$ of slow and fast modes. Their relation depends on
the compressibility of the medium and plasma $\beta$, where
$\beta=P_{\text{gas}}/P_{\text{mag}}>1$, , and $P_{\text{gas}}$ and
$P_{\text{mag}}$ are gas/plasma and magnetic pressures. In a in high-$\beta$
medium, with the sonic Mach number $M_{s}=V_{L}/c_{s}<1$, where $V_{L}$ is the
driven turbulent velocity, $c_{s}$ is the sound speed, $\delta B_{s}$ is
expected to be larger than $\delta B_{f}$. In the opposite limit of
magnetically dominated medium with $\beta<1$, there is $\delta B_{s}<\delta
B_{f}$.

In the literature, particles with $\mu<\mu_{lc}$ are considered “trapped” in
magnetic bottles and thus they cannot diffuse. However, as we will show below
(Section 2.3), this is an erroneous notion. In fact, for 3D propagation in MHD
turbulence, magnetic mirroring results in a new type of parallel diffusion of
bouncing particles.

###  2.2 Perpendicular superdiffusion

CRs following magnetic field lines 333In the presence of turbulence, CRs
follow the magnetic flux tubes averaged within their gyro-orbit. We
nevertheless adopt the generally accepted way to describe magnetic fields as
“magnetic field lines”. undergo pitch-angle scattering via the interaction
with small-scale magnetic fluctuations and bouncing among magnetic mirrors
(see above). In the meantime, the dispersion of their trajectories in the
direction perpendicular to the mean magnetic field is determined by the
dispersion of magnetic field lines. This dispersion increases with the
propagation distance of CRs along the magnetic field and accounts for the
perpendicular diffusion of CRs, as pointed out by Jokipii (1966). The
superdiffusive nature of this dispersion within the range of length scales of
strong MHD turbulence was later found by Lazarian & Vishniac (1999) when they
analytically quantified the magnetic field wandering induced by the Alfvénic
component444The separation of MHD turbulence into Alfvenic, slow and fast
modes was demonstrated numerically in Cho & Lazarian (2002, 2003). The
dominance of Alfvenic modes in terms of inducing magnetic field wandering was
shown analytically in Lazarian & Vishniac (1999). of MHD turbulence for both
super-Alfvénic ($M_{A}=V_{L}/V_{A}>1$) and sub-Alfvénic ($M_{A}<1$)
turbulence, where $V_{A}$ is the Alfvén speed. We note that the perpendicular
superdiffusion is with respect to the mean magnetic field in sub- and trans-
Alfvénic ($M_{A}=1$) turbulence, and with respect to the local mean field,
i.e., the magnetic field averaged over scales less than $l_{A}=LM_{A}^{-3}$ in
super-Alfvénic turbulence (Lazarian, 2006), where $L$ is the turbulence
driving scale, and $l_{A}$ is the scale where the local turbulent velocity
becomes equal to $V_{A}$. For the turbulent motions on scales larger than
$l_{A}$, the effect of magnetic field is subdominant, and they gradually get
isotropic, similar to the hydrodynamic Kolmogorov turbulence.

For sub-Alfvénic turbulence, there exists a scale $l_{\text{tran}}=LM_{A}^{2}$
for the transition from the weak turbulence over the scales
$[l_{\text{tran}},L]$ to the strong turbulence on smaller length scales. In
the weak turbulence regime, the magnetic fields exhibit normal diffusion, i.e.
the mean squared separation between magnetic field lines increases with $s$,
where $s$ is the distance measured along the magnetic field lines (Lazarian &
Vishniac, 1999).

We next consider the divergence of magnetic field lines over scales less than
$l_{\text{tran}}$. Naturally, when one follows the magnetic field line over
the parallel scale $l_{\|}$ of a turbulent eddy, the magnetic field line has
its perpendicular displacement equal to the transverse size $l_{\bot}$ of the
eddy. As $l_{\bot}$ can be both positive and negative, the dispersion $\langle
y^{2}\rangle$ of magnetic field lines increases as

$d\langle y^{2}\rangle\approx l_{\bot}^{2}\frac{ds}{l_{\|}},$ (4)

where $ds$ is the distance measured along the magnetic field line, and the
bracket denotes an ensemble average. Using the above relation and the scaling
relation between $l_{\|}$ and $l_{\bot}$ for sub-Alfvénic turbulence (Lazarian
& Vishniac (1999), see also Appendix A), one can find

$d\langle y^{2}\rangle\approx\langle y^{2}\rangle^{2/3}M_{A}^{4/3}L^{-1/3}ds.$
(5)

If the expression on the right-hand side were “constant$\times ds$”, Eq. (5)
would describe the random walk of field lines. However, it has a power-law
dependence on the field line separation with a positive power index, i.e.
$\propto\langle y^{2}\rangle^{2/3}$, leading to an accelerating superdiffusion
as magnetic field lines spread out. The numerical demonstration of the
superdiffusion of turbulent magnetic fields can be found in, e.g., Lazarian et
al. (2004); Beresnyak (2013).

The physical explanation of the effect is the following. As we follow the
magnetic field line the distance $\sim s$, the divergence rate of field lines
increases with larger and larger turbulent eddies contributing to the
dispersion of field line separations. A similar effect in hydrodynamic
turbulence is related to the accelerated separation of a pair of test
particles with time, which is known as Richardson dispersion (Richardson,
1926).555In the presence of magnetic field, this results in the Richardson
dispersion of magnetic field lines with time (Eyink et al., 2011, 2013), as a
time-dependent analog of the superdiffusion of magnetic field lines in space.
The accelerated divergence of magnetic field lines results in the
superdiffusive field line separations. The analogy between the Richardson
dispersion and the superdiffusion of magnetic fields introduced in Lazarian &
Vishniac (1999) was discussed in detail in Eyink et al. (2011).

Eq. (5) also applies to the dispersion of separations of CRs that follow
magnetic field lines. The resulting superdiffusive perpendicular divergence of
CR trajectories was studied in Yan & Lazarian (2008); Lazarian & Yan (2014) in
two regimes, one describing the ballistic motion of CRs along the magnetic
field, the other describing the diffusion of CRs along the magnetic field.
Both regimes were numerically tested in Xu & Yan (2013). For the ballistic
propagation, the corresponding dispersion of separations of particles in the
perpendicular direction is

$\langle y^{2}\rangle\sim
l_{\bot}^{2}\approx\frac{s^{3}}{L}M_{A}^{4},~{}~{}~{}M_{A}<1,~{}s<\lambda_{\|},~{}s<L.$
(6)

Here $s$ is the distance of particles measured along the magnetic field line,
and $\lambda_{\|}$ is the CR parallel mean free path. $y$ can be identified
with the displacement of particles perpendicular to the mean magnetic field.
Note that the dependence $\langle y^{2}\rangle\propto s^{3}$ in the above
expression is also applicable to super-Alfvénic turbulence at scales less than
$l_{A}$, i.e., $s<l_{A}$.

In the presence of the parallel diffusion, i.e. over the distance
$s>\lambda_{\|}$, the perpendicular dispersion of CRs is modified to scale as
$\langle y^{2}\rangle\propto s^{3/2}$. Naturally, this is still in a
superdiffusion regime that arises from the superdiffusive magnetic field line
dispersion. Therefore, in strong MHD turbulence, i.e., on scales smaller than
$l_{A}$ in super-Alfvénic turbulence and $l_{\text{tran}}$ in sub-Alfvénic
turbulence, for both ballistic and diffusive CR propagation along the
turbulent magnetic field, CR trajectories get separated superdiffusively in
the direction perpendicular to the mean magnetic field, although the law of
superdiffusion differs in the two cases. A detailed discussion on the
perpendicular (super)diffusion of CRs in different turbulence regimes is
provided in Lazarian & Yan (2014) (see Table 1 in Lazarian & Yan 2019).

We note that the concept of perpendicular superdiffusion of CRs contradicts to
some existing views of CR transport, e.g. the Non-Linear Guiding Center Theory
(NLGCT) (Matthaeus et al., 2003) that is formulated using the 2D/slab model of
MHD turbulence. As the (super)diffusion behavior of CRs strongly depends on
the properties of MHD turbulence, CR (super)diffusion should be studied using
the tested model of MHD turbulence, instead of synthetic models.

###  2.3 Parallel mirror diffusion

The effect of magnetic mirrors on trapping particles was extensively studied
in the literature (e.g., Noerdlinger 1968; Kulsrud & Pearce 1969; Cesarsky &
Kulsrud 1973). There it was assumed that the trapped particles undergo
oscillatory motions between two mirror points of a magnetic bottle, without
cumulative diffusion along the magnetic field. As the main difference between
this work and earlier studies, here we consider the existence of Alfvénic
component of turbulence and the resulting perpendicular superdiffusion of
magnetic fields and CRs. As we explain below, in the presence of Alfvénic
turbulence, the superdiffusion of magnetic fields does not allow the particles
to be trapped in the same magnetic bottle.

For every crossing of a magnetic bottle induced by the compressive component
of MHD turbulence with a size $s$, the particle experiences perpendicular
superdiffusion with the perpendicular displacement $y$ given in Section 2.2.
As a result, the particle escapes the magnetic bottle within one crossing time
and then encounters another magnetic bottle within a perpendicular distance
$y$ from the previous magnetic bottle. As a result, the bouncing with the
mirror points of different magnetic bottles leads to the diffusive motion of
particles along turbulent magnetic field lines. We term this diffusion
parallel to local magnetic field arising from the mirroring effect as
“parallel mirror diffusion” to distinguish it from the traditional CR parallel
diffusion associated with the resonance scattering (see Schlickeiser 2002). We
illustrate the parallel mirror diffusion in Fig. 1. CRs that follow diffusing
magnetic field lines are not trapped, but move diffusively parallel to the
turbulent magnetic field when bouncing with different magnetic mirrors.
Through this paper, “bouncing” and “mirroring” are equivalent.

Due to the 3D character of CR motions and the complexity of magnetic field
structure at different $M_{A}$, it is important to clarify that the parallel
diffusion of CRs in MHD turbulence is the diffusion with respect to the local
magnetic field sampled by CRs. At $M_{A}<1$, the local magnetic field has the
direction close to that of the global mean magnetic field, while at $M_{A}>1$,
the magnetic field direction changes significantly on scales larger than
$l_{A}$. In the latter case, we consider the parallel mirror diffusion on
scales smaller than $l_{A}$. The parallel mirror diffusion of CRs is
accompanied by the dispersion of their trajectories in the direction
perpendicular to the mean magnetic field (see Fig. 1). The perpendicular
superdiffusion of CRs can be observed when their initial separations in space
are small, as shown in simulations (Xu & Yan, 2013). The exact scaling
relation between $\langle y^{2}\rangle$ and $s$ depends on $\lambda_{\|}$ of
bouncing CRs (Section 2.2). $\lambda_{\|}$ is related to the parallel
diffusion coefficient $D_{\|}$ by $D_{\|}=1/3v\lambda_{\|}$, where $v$ is
particle velocity. $D_{\|}$ of bouncing CRs will be derived in the next
sections.

Figure 1:  Because of the perpendicular superdiffusion of particles following
turbulent magnetic field lines, particles that undergo bouncing move
diffusively along the local magnetic field. Thin lines represent magnetic
field lines. Thick lines represent the trajectories of two particles whose
initial separation is small.

##  3 Parallel diffusion of bouncing particles induced by fast modes in
compressible MHD turbulence

###  3.1 Bouncing and non-bouncing particles

In compressible MHD turbulence, if fast modes carry a significant fraction of
the injected turbulent energy, they act as the main agent for scattering
particles (Yan & Lazarian, 2002, 2004). More recently, XL20 identified the
important role of fast modes in both bouncing and scattering CRs. XL20 studied
the parallel diffusion associated with the gyroresonant scattering of non-
bouncing particles. Here we will focus on the parallel mirror diffusion of
bouncing particles.

The scattering causes diffusion in pitch angle. In the quasilinear
approximation (Jokipii, 1966), for the gyroresonant scattering by fast modes,
the pitch-angle diffusion coefficient is (Voelk, 1975)

$D_{\mu\mu,\text{QLT},f}=C_{\mu}\int
d^{3}k\frac{k_{\|}^{2}}{k^{2}}[J_{1}^{\prime}(x)]^{2}I_{f}(k)R(k).$ (7)

Here

$C_{\mu}=(1-\mu^{2})\frac{\Omega^{2}}{B_{0}^{2}},~{}~{}x=\frac{k_{\perp}v_{\perp}}{\Omega}=\frac{k_{\perp}}{r_{g}^{-1}},$
(8)

$B_{0}$ is the mean magnetic field strength, $v_{\perp}$ is the perpendicular
component of $v$, $\Omega$ is the gyrofrequency, $r_{g}$ is the particle
gyroradius, and $k_{\|}$ and $k_{\perp}$ are parallel and perpendicular
components of wavenumber $k$. Besides,

$I_{f}(k)=C_{f}k^{-\frac{7}{2}}$ (9)

is the energy spectrum of fast modes (Cho & Lazarian, 2002b), with

$C_{f}=\frac{1}{16\pi}\delta B_{f}^{2}L^{-\frac{1}{2}},$ (10)

where $\delta B_{f}$ is the magnetic perturbation of fast modes at $L$, and

$R=\pi\delta(\omega_{k}-v_{\|}k_{\|}+\Omega)$ (11)

is the resonance function for gyroresonance in the quasilinear approximation,
where $\omega_{k}$ is the wave frequency. We see that when the pitch angle
approaches $90^{\circ}$, the above resonance condition cannot be satisfied
with a limited range of $k$ and decreasing magnetic fluctuation amplitude with
increasing $k$, leading to the $90^{\circ}$ problem666The vanishing scattering
close to $90^{\circ}$ and infinite mean free path in the quasi-linear theory
is known as the $90^{\circ}$ problem (Fisk et al., 1974).. For our analytical
estimate, we can use the approximate expression of $D_{\mu\mu,\text{QLT},f}$
(Xu et al., 2016; Xu & Lazarian, 2018)

$\displaystyle D_{\mu\mu,\text{QLT},f}$
$\displaystyle\approx\frac{\pi}{56}\frac{\delta
B_{f}^{2}}{B_{0}^{2}}\Big{(}\frac{v}{L\Omega}\Big{)}^{\frac{1}{2}}\Omega(1-\mu^{2})\mu^{\frac{1}{2}}.$
(12)

We define the rate of change in $\mu$ due to scattering as the scattering rate
(XL20),

$\Gamma_{s,f}=\frac{2D_{\mu\mu,\text{QLT},f}}{\mu^{2}}\approx\frac{\pi}{28}\frac{\delta
B_{f}^{2}}{B_{0}^{2}}\Big{(}\frac{v}{L\Omega}\Big{)}^{\frac{1}{2}}\Omega(1-\mu^{2})\mu^{-\frac{3}{2}}.$
(13)

In the presence of the stochastic magnetic mirrors induced by fast modes, the
particles that satisfy the bouncing condition (see below) undergo the bouncing
motions among different magnetic mirrors as discussed in Section 2.3. For
particles with (see Eq. (3))

$\mu\approx\sqrt{\frac{b_{fk}}{B_{0}}},$ (14)

the bouncing is dominated by the magnetic mirrors at $k$ with the magnetic
perturbation $b_{fk}$ (Cesarsky & Kulsrud, 1973), where

$b_{fk}=\delta B_{f}(kL)^{-\frac{1}{4}}$ (15)

according to the scaling of isotropic fast modes (Cho & Lazarian, 2002b). The
rate of the adiabatic change in $\mu$ due to the spatial variation of magnetic
field is (Cesarsky & Kulsrud 1973,XL20)

$\displaystyle\Gamma_{b,f}=\Big{|}\frac{1}{\mu}\frac{d\mu}{dt}\Big{|}$
$\displaystyle=\frac{v}{2B_{0}}\frac{1-\mu^{2}}{\mu}b_{fk}k$ (16)
$\displaystyle=\frac{v}{2L}\frac{\delta
B_{f}^{4}}{B_{0}^{4}}\frac{1-\mu^{2}}{\mu^{7}}$

at $\mu>\mu_{\text{min,f}}$, where

$\mu_{\text{min,f}}=\sqrt{\frac{b_{fk}(r_{g})}{B_{0}}}=\sqrt{\frac{\delta
B_{f}}{B_{0}}}\Big{(}\frac{r_{g}}{L}\Big{)}^{\frac{1}{8}},$ (17)

and $b_{fk}(r_{g})$ is the magnetic perturbation at $r_{g}$. At
$\mu<\mu_{\text{min,f}}$, the bouncing rate is

$\displaystyle\Gamma_{b,f}$
$\displaystyle=\frac{v}{2B_{0}}\frac{1-\mu^{2}}{\mu}b_{fk}(r_{g})r_{g}^{-1}$
(18) $\displaystyle=\frac{v}{2r_{g}}\frac{\delta
B_{f}}{B_{0}}\Big{(}\frac{r_{g}}{L}\Big{)}^{\frac{1}{4}}\frac{1-\mu^{2}}{\mu}.$

At the balance between scattering and bouncing, i.e.,
$\Gamma_{s,f}=\Gamma_{b,f}$, one can find the cutoff $\mu$ (Eqs. (13) and
(16),

$\mu_{c}\approx\bigg{[}\frac{14}{\pi}\frac{\delta
B_{f}^{2}}{B_{0}^{2}}\Big{(}\frac{v}{L\Omega}\Big{)}^{\frac{1}{2}}\bigg{]}^{\frac{2}{11}},$
(19)

in agreement with the result in XL20. We define the particles with
$\mu<\mu_{c}$ as “bouncing particles” and those with $\mu>\mu_{c}$ as “non-
bouncing” particles. Their motions are dominated by bouncing and scattering,
respectively. As shown in Fig. 2, $\mu_{c}$ increases with the CR energy
$E_{\text{CR}}$ until reaching its maximum value

$\mu_{c,\text{max}}=\sqrt{\frac{\delta B_{f}}{B_{0}+\delta B_{f}}}.$ (20)

Here as an illustration, we consider CR protons, the magnetic field strength
$\delta B_{f}=B_{0}=3\mu$G, and $L=30$ pc as the driving scale of interstellar
turbulence.

Figure 2: $\mu_{c}$ as a function of $E_{\text{CR}}$ for fast modes. The
analytical approximation is given by Eq. (19). The numerical result is
obtained based on the numerical evaluation of Eq. (7).

###  3.2 Parallel diffusion of bouncing particles

When the mirroring condition $\mu<\mu_{c}$ is satisfied, the parallel
diffusion of particles with $\mu$ is regulated by the mirroring effect. Based
on the relations in Eqs (14) and (15), bouncing particles diffuse along the
magnetic field with a step size

$k^{-1}=L\Big{(}\frac{\delta B_{f}}{B_{0}}\Big{)}^{-4}\mu^{8}$ (21)

for $k^{-1}>r_{g}$, or equivalently $\mu>\mu_{\text{min,f}}$. At
$\mu<\mu_{\text{min,f}}$, the step size of parallel diffusion is $r_{g}$,
corresponding to the Bohm diffusion. The maximal step size of bouncing
particles corresponds to $\mu_{c}$,

$k_{c}^{-1}=L\Big{(}\frac{\delta B_{f}}{B_{0}}\Big{)}^{-4}\mu_{c}^{8}.$ (22)

It decreases with the increase of $\delta B_{f}/B_{0}$ (see Eq. (19)) and
should not exceed $L$. The bouncing of the particles with $\mu>\mu_{c}$ by the
magnetic fluctuations on scales larger than $k_{c}^{-1}$ is less efficient
than their gyroresonant scattering by small-scale magnetic fluctuations.

The corresponding parallel diffusion coefficient arising from bouncing is

$\displaystyle D_{\|,f,b}(\mu)=$ $\displaystyle v\mu
k^{-1}=vL\Big{(}\frac{\delta B_{f}}{B_{0}}\Big{)}^{-4}\mu^{9},$ $\displaystyle
D_{\|,f,b}(\mu)=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mu_{\text{min,f}}<\mu<\mu_{c},$
(23a) $\displaystyle D_{\|,f,b}(\mu)=$ $\displaystyle v\mu
r_{g},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mu<\mu_{\text{min,f}}.$
(23b)

The strong dependence on $\mu$ comes from the $\mu$-dependence of the size of
the dominant magnetic mirror field (Eq. (21)). The bouncing makes the parallel
diffusion of CRs with a smaller $\mu$ more inefficient. This is opposite to
the expectation from the quasi-linear scattering theory.

We first assume an isotropic pitch angle distribution for a simple analytical
estimate of the parallel diffusion coefficient. By integrating
$D_{\|,f,b}(\mu)$ over $\mu$, we find the parallel diffusion coefficient as a
sum of the contributions from the multi-scale magnetic mirror field for
bouncing particles with different $\mu$,

$\displaystyle D_{\|,f,b}$
$\displaystyle=\int_{0}^{\mu_{c}}D_{\|,f,b}(\mu)d\mu$ (24)
$\displaystyle\approx\frac{1}{10}vL\Big{(}\frac{\delta
B_{f}}{B_{0}}\Big{)}^{-4}\mu_{c}^{10},$

where we consider that the contribution from $D_{\|,f,b}(\mu)$ at
$\mu<\mu_{\text{min,f}}$ is small.

In the situation with an anisotropic pitch angle distribution resulting from
the anisotropic scattering (see Eq. (12)), we focus on the pitch-angle
diffusion and use the spatially averaged Fokker-Planck equation

$\frac{\partial f}{\partial
t}=\frac{\partial}{\partial\mu}\bigg{[}D_{\mu\mu}\frac{\partial
f}{\partial\mu}\bigg{]}$ (25)

to derive the particle distribution function $f$ under the effect of
scattering. If the steady state, i.e., $\partial f/\partial t=0$, can be
reached, we have

$D_{\mu\mu}\frac{df}{d\mu}=\text{C},$ (26)

where $-$C is the constant flux of particles in $\mu$ space. The steady-state
solution is

$f(\mu)=-\text{C}\int_{\mu}^{\mu_{c}}\frac{d\mu^{\prime}}{D_{\mu\mu}(\mu^{\prime})}+f(\mu_{c}).$
(27)

To simplify the evaluation of $f(\mu)$, we further adopt the boundary
condition

$f(\mu_{\text{c}})=0$ (28)

by assuming that the diffusion of non-bouncing particles is relatively fast
and the non-bouncing particles at $\mu>\mu_{c}$ instantly escape from the
system. This assumption implies that the diffusion coefficient of the non-
bouncing particles is much larger than that of the bouncing particles. Its
validity will be examined later. We then normalize the distribution function
of the bouncing particles

$f^{\prime}=\frac{f}{\int_{0}^{\mu_{c}}fd\mu}$ (29)

to have unit integral,

$\int_{0}^{\mu_{c}}f^{\prime}d\mu=1.$ (30)

The parallel diffusion coefficient of bouncing particles under the
consideration of anisotropic pitch angle distribution is then

$D_{\|,f,b}=\int_{0}^{\mu_{c}}D_{\|,f,b}(\mu)f^{\prime}d\mu,$ (31)

where $D_{\|,f,b}(\mu)$ is given by Eq. (23b). It is evident that the above
calculation of $D_{\|,f,b}$ for bouncing particles depends on $D_{\mu\mu}$
that originates from scattering.

To provide an example for illustrating the parallel diffusion of bouncing
particles, in Fig. 3, we present the numerically calculated $D_{\|,f,b}$ as a
function of $E_{\text{CR}}$ and its analytical approximation (Eq. (24)). The
same parameters as used in Fig. 2 are adopted. We see that the simplification
by using an isotropic pitch angle distribution in Eq. (24) does not
significantly affect the result. This is because for scattering by fast modes,
$D_{\mu\mu,\text{QLT},f}$ and the resulting $f(\mu)$ are not significantly
anisotropic. 777The result derived from isotropic pitch angle distribution is
supposed to be larger than the one derived from anisotropic distribution. Our
analytical approximation underestimates the former due to the underestimation
of $\mu_{c}$ (see Fig. 2). The analytical estimate has the energy dependence
as $D_{\|,f,b}\propto E_{\text{CR}}^{10/11}$ for relativistic particles (Eqs.
(19) and (24)), which is close to $D_{\|,f,b}\propto E_{\text{CR}}$. The
numerical result is a bit shallower and can be fitted by $D_{\|,f,b}\propto
E_{\text{CR}}^{0.7}$. The discrepancy between the analytical and numerical
results at the high-energy end comes from the energy dependence of
$D_{\|,f,b}(\mu)$ at $\mu<\mu_{\text{min,f}}$ (Eq. (23b)), which is not taken
into account in the analytical approximation.

As a comparison, in Fig. 3 we also add the parallel diffusion coefficient of
non-bouncing particles determined by the gyroresonant scattering by fast modes
derived in XL20. The analytical approximation is given by

$\displaystyle D_{\|,f,nb}$
$\displaystyle=\frac{v^{2}}{4}\int_{\mu_{c}}^{1}d\mu\frac{(1-\mu^{2})^{2}}{D_{\mu\mu,\text{QLT},f}(\mu)}$
(32) $\displaystyle\approx\frac{28}{5\pi}\frac{B_{0}^{2}}{\delta
B_{f}^{2}}\Big{(}\frac{v}{L\Omega}\Big{)}^{-\frac{1}{2}}\frac{v^{2}}{\Omega}\big{[}4-\sqrt{\mu_{c}}(5-\mu_{c}^{2})\big{]},$

where the lower bound in integration is given by $\mu_{c}$ (see Cesarsky &
Kulsrud 1973). The resulting $D_{\|,f,nb}$ for non-bouncing particles is
shallower than $\propto E_{\text{CR}}^{0.5}$ except for the high-energy end
with a constant $\mu_{c}$. Its scaling with $E_{\text{CR}}$ is not an exact
power law, but can be approximated by $\propto E_{\text{CR}}^{1/3}$ (see Fig.
3). We see that there is $D_{\|,f,b}<D_{\|,f,nb}$ for all CR energies
considered here. This justifies the assumption on the boundary condition in
Eq. (28). The bouncing particles have a more inefficient diffusion compared
with the non-bouncing particles.

We note that in the above illustration, we assume $\delta B_{f}=B_{0}$. In
this case with $\delta B_{f}<B_{0}$, both bouncing and non-bouncing particles
would have larger diffusion coefficients. Nevertheless, the diffusion of non-
bouncing particles is expected to be faster than that of bouncing particles.

Figure 3: Parallel diffusion coefficients $D_{\|,f,b}$ of bouncing CRs and
$D_{\|,f,nb}$ of non-bouncing CRs as a function of CR energy in compressible
MHD turbulence with fast modes dominating both bouncing and scattering. The
analytical approximations are given by Eq. (24) and Eq. (32), respectively.

##  4 Parallel diffusion of bouncing particles in incompressible MHD
turbulence

###  4.1 Trans-Alfvénic turbulence

CRs propagation in incompressible MHD turbulence was considered in many
earlier studies using isotropic MHD turbulence model or 2D/slab superposition
model (Matthaeus et al., 1990). Later it was found that if turbulent energy is
injected at large scales, as this is the case of interstellar turbulence
(e.g., Chepurnov et al. 2010), the gyroresonant scattering by both Alfvén
modes and pseudo-Alfvén modes is inefficient especially for low-energy CRs due
to the scale-dependent anisotropy of MHD turbulence (Chandran, 2000a; Yan &
Lazarian, 2002). Yan & Lazarian (2002) further argued that fast modes in
compressible MHD turbulence, rather than Alfvén modes, determine the
scattering of GeV CRs in our galaxy (see Xu & Lazarian 2018 for discussion on
importance of slow modes in TTD interaction with CRs). However, considering
the variety of astrophysical conditions, there are situations where fast modes
are subject to severe damping (Yan & Lazarian, 2004; Brunetti & Lazarian,
2007; Xu et al., 2016; Xu & Lazarian, 2018). When fast modes are severely
damped, or plasma $\beta$ is large (see Section 2.1), or the injected energy
of fast modes is small, additional sources for confining CRs are needed. This
motives us to study CR diffusion in idealized incompressible MHD turbulence or
compressible MHD turbulence with a very small energy fraction in fast modes.
In this situation, we focus on the mirroring effect induced by pseudo-Alfvén
(slow) modes. The magnetic mirroring effect has been found in XL20 to be
important for confining CRs in incompressible MHD turbulence. Here we will
formulate the corresponding parallel diffusion coefficient.

The pitch-angle diffusion coefficients for gyroresonant interactions with
Alfvén and pseudo-Alfvén modes in the quasi-linear approximation are (Voelk,
1975),

$D_{\mu\mu,\text{QLT},A}=C_{\mu}\int d^{3}kx^{-2}[J_{1}(x)]^{2}I_{A}(k)R(k),$
(33)

and

$D_{\mu\mu,\text{QLT},s}=C_{\mu}\int
d^{3}k\frac{k_{\|}^{2}}{k^{2}}[J_{1}^{\prime}(x)]^{2}I_{s}(k)R(k).$ (34)

The magnetic energy spectra are (Cho et al., 2002)

$I_{A}(k)=C_{A}k_{\perp}^{-\frac{10}{3}}\exp{\Bigg{(}-L^{\frac{1}{3}}\frac{k_{\|}}{k_{\perp}^{\frac{2}{3}}}\Bigg{)}},~{}C_{A}=\frac{1}{6\pi}\delta
B_{A}^{2}L^{-\frac{1}{3}}$ (35)

for Alfvén modes, and

$I_{s}(k)=C_{s}k_{\perp}^{-\frac{10}{3}}\exp{\Bigg{(}-L^{\frac{1}{3}}\frac{k_{\|}}{k_{\perp}^{\frac{2}{3}}}\Bigg{)}},~{}C_{s}=\frac{1}{6\pi}\delta
B_{s}^{2}L^{-\frac{1}{3}}$ (36)

for pseudo-Alfvén modes, where $\delta B_{A}$ and $\delta B_{s}$ are their
magnetic perturbations at $L$. The corresponding scattering rate with
contributions from both Alfvén and pseudo-Alfvén modes is

$\Gamma_{s,\text{inc}}=\frac{2(D_{\mu\mu,\text{QLT},A}+D_{\mu\mu,\text{QLT},s})}{\mu^{2}}.$
(37)

Besides scattering, pseudo-Alfvén modes also generate magnetic mirrors
accounting for the bouncing of CRs. The particles with

$\mu\approx\sqrt{\frac{b_{sk}}{B_{0}}}$ (38)

mainly undergo the bouncing motion at $k_{\|}$, where the magnetic
perturbation of pseudo-Alfvén modes is

$b_{sk}=\delta B_{s}(k_{\perp}L)^{-\frac{1}{3}}=\delta
B_{s}(k_{\|}L)^{-\frac{1}{2}}$ (39)

based on the anisotropic scaling of MHD turbulence (Cho et al., 2002). The
bouncing rate is (Cesarsky & Kulsrud 1973, XL20)

$\displaystyle\Gamma_{b,s}=\Big{|}\frac{1}{\mu}\frac{d\mu}{dt}\Big{|}$
$\displaystyle=\frac{v}{2B_{0}}\frac{1-\mu^{2}}{\mu}b_{sk}k_{\|}$ (40)
$\displaystyle=\frac{v}{2L}\aleph^{2}_{s}\frac{1-\mu^{2}}{\mu^{3}}$

at $\mu>\mu_{\text{min,s}}$, where

$\aleph_{s}=\frac{\delta B_{s}}{B_{0}},$ (41)
$\mu_{\text{min,s}}=\sqrt{\frac{b_{sk}(r_{g})}{B_{0}}}=\sqrt{\aleph_{s}}\Big{(}\frac{r_{g}}{L}\Big{)}^{\frac{1}{4}},$
(42)

and $b_{sk}(r_{g})$ is the magnetic perturbation at $k_{\|}=1/r_{g}$. At
$\mu<\mu_{\text{min,s}}$, the bouncing rate becomes

$\displaystyle\Gamma_{b,s}$
$\displaystyle=\frac{v}{2B_{0}}\frac{1-\mu^{2}}{\mu}b_{sk}(r_{g})r_{g}^{-1}$
(43)
$\displaystyle=\frac{v}{2r_{g}}\aleph_{s}\Big{(}\frac{r_{g}}{L}\Big{)}^{\frac{1}{2}}\frac{1-\mu^{2}}{\mu}.$

Here we note that as $\Gamma_{b,s}$ and $\Gamma_{b,f}$ (Eqs. (16) and (18))
have different $\mu$ dependence, in the presence of both slow and fast modes
in compressible MHD turbulence, the relative importance between the mirroring
effects induced by slow and fast modes can change with $\mu$ (see Section
5.2).

As found in XL20, because of the inefficient scattering by Alfvén and pseudo-
Alfvén modes, the effect of magnetic mirrors dominates over the gyroresonant
scattering, i.e. $\Gamma_{b,s}>\Gamma_{s,\text{inc}}$, for $\mu<\mu_{c}$, and
$\mu_{c}$ is given by its maximum value 888$\mu_{c}$ in incompressible MHD
turbulence was approximated to be unity in XL20 for simplicity.

$\mu_{c}=\sqrt{\frac{\delta B_{s}}{B_{0}+\delta B_{s}}},$ (44)

which is independent of $E_{\text{CR}}$.

The non-bouncing particles with $\mu>\mu_{c}$ are only poorly constrained by
the inefficient scattering, and their diffusion coefficients are expected to
be large (see Section 5). For the parallel diffusion of bouncing particles,
the step size is

$k_{\|}^{-1}=L\aleph^{-2}_{s}\mu^{4}$ (45)

at $\mu>\mu_{\text{min,s}}$, or equivalently, $k_{\|}^{-1}>r_{g}$, and the
step size is $r_{g}$ at $\mu<\mu_{\text{min,s}}$. Accordingly, the
$\mu$-dependent parallel diffusion coefficient is

$\displaystyle D_{\|,\text{inc},b}(\mu)=$ $\displaystyle v\mu
k_{\|}^{-1}=vL\aleph^{-2}_{s}\mu^{5},$ $\displaystyle
D_{\|,\text{inc},b}(\mu)=$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mu_{\text{min,s}}<\mu<\mu_{c},$
(46a) $\displaystyle D_{\|,\text{inc},b}(\mu)=$ $\displaystyle v\mu
r_{g},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mu<\mu_{\text{min,s}}.$
(46b)

If the pitch angle distribution of bouncing particles is isotropic, then the
parallel diffusion coefficient can be approximated by

$D_{\|,\text{inc},b}\approx\int_{0}^{\mu_{c}}D_{\|,\text{inc},b}(\mu)d\mu=\frac{1}{6}vL\aleph_{s}^{-2}\mu_{c}^{6},$
(47)

which is energy independent with a constant $\mu_{c}$. Here we neglect the
contribution from $D_{\|,\text{inc},b}(\mu)$ at $\mu<\mu_{\text{min,s}}$.

We next consider the situation with an anisotropic pitch angle distribution.
We follow the similar analysis in Section 3.2 to derive $D_{\|,\text{inc},b}$
numerically by using Eqs. (27)-(31). In Eq. (27), we adopt (Eqs. (33) and
(34))

$D_{\mu\mu}(\mu)=D_{\mu\mu,\text{QLT},A}+D_{\mu\mu,\text{QLT},s}$ (48)

to include the gyroresonant scattering by both Alfvén and pseudo-Alfvén modes.
Eq. (31) becomes

$D_{\|,\text{inc},b}=\int_{0}^{\mu_{c}}D_{\|,\text{inc},b}(\mu)f^{\prime}d\mu,$
(49)

where $D_{\|,\text{inc},b}(\mu)$ is given by Eq. (46b).

In Fig. 4, we illustrate the numerically calculated $D_{\|,\text{inc},b}$ with
an anisotropic pitch angle distribution in comparison with the analytical
approximation with an isotropic pitch angle distribution (Eq. (47)). The same
parameters as in Fig. 2 are used, and the turbulence is considered to be
trans-Alfvénic with $\delta B_{A}=B_{0}$. Here we also use $\delta
B_{s}=B_{0}$, i.e., $\aleph_{s}=1$. We see that the anisotropic distribution
leads to a significantly smaller $D_{\|,\text{inc},b}$ than that derived from
an isotropic distribution. Fig. 5 displays the anisotropic $f^{\prime}$ caused
by the anisotropic pitch-angle scattering for different CR energies. The
highly anisotropic scattering originates from the scale-dependent anisotropy
of Alfvén and pseudo-Alfvén modes. At a small $\mu$, the gyroresonance with
many uncorrelated eddies with the perpendicular eddy size much smaller than
$r_{g}$ makes the scattering inefficient (Chandran, 2000a; Yan & Lazarian,
2002). Consequently, most particles are concentrated at a small $\mu$.

The more pronounced anisotropy of $f^{\prime}$ for low-energy CRs is caused by
the inner cutoff at $k_{\text{max}}$ of turbulent energy spectrum. XL20 found
that for the gyroresonance with Alfvén and pseudo-Alfvén modes, the magnetic
fluctuations at

$k_{\perp
p}=\Big{(}\frac{L^{\frac{1}{3}}k_{\|,\text{res}}}{8}\Big{)}^{\frac{3}{2}}=8^{-\frac{3}{2}}k_{\perp,\text{res}}$
(50)

play a dominant role in determining the scattering efficiency. Here the
parallel and perpendicular resonant wavenumbers are
$k_{\|,\text{res}}\approx\Omega/v_{\|}$ and
$k_{\perp,\text{res}}=L^{\frac{1}{2}}k_{\|,\text{res}}^{\frac{3}{2}}$ for
trans-Alfvénic turbulence. The former is derived from the resonance condition
in Eq. (11). When $k_{\perp p}>k_{\text{max}}$ at a small $\mu$ and a low
$E_{\text{CR}}$, the scattering becomes ineffective, leading to a more
anisotropic $f^{\prime}$. As a result, the integral in Eq. (49) is dominated
by $D_{\|,\text{inc},b}(\mu)$ at a small $\mu$ (Eq. (46b)) (see Fig. 5). For
the parameters adopted here, the CR energy corresponding to $k_{\perp
p}=k_{\text{max}}=10^{-8}$ cm-1 and the minimum $\mu=0.01$ is indicated by the
dashed line in Fig. 4, below which the damping of turbulent spectrum has a
significant effect on the diffusion of bouncing particles. The enhanced
anisotropy due to damping gives rise to more suppressed diffusion.

At higher energies, Eq. (49) is dominated by $D_{\|,\text{inc},b}(\mu)$ at a
larger $\mu$, where both $f^{\prime}$ and $D_{\|,\text{inc},b}(\mu)$ remain
unchanged (Figs. 5 and 5), and the resulting $D_{\|,\text{inc},b}$ becomes
energy independent (Fig. 4). For high-energy CRs, Eq. (49) is again dominated
by $D_{\|,\text{inc},b}(\mu)$ at a small $\mu$ similar to low-energy CRs, but
due to the increase of $r_{g}$ with $E_{\text{CR}}$ (Eq. (46b)). At both low-
and high-energy ends, $D_{\|,\text{inc},b}$ has the energy dependence close to
$D_{\|,\text{inc},b}\propto E_{\text{CR}}$ as dictated by Eq. (46b).

In the vicinity of a CR source, the initial pitch angle distribution of the
injected CRs is important for determining the parallel diffusion of bouncing
CRs until they lose the memory about the initial distribution via scattering.
Therefore, the CRs closer to the source can have slower diffusion. We note
that the pitch angle distribution of low-energy CRs can also be affected by
CR-driven instabilities especially near CR sources. Intensive scattering by
Alfvén waves generated by streaming CRs can reduce the anisotropy in pitch
angle distribution. This effect should be taken into account for more
realistic modeling of the paralell diffusion coefficient of low-energy CRs
(see, e.g., Blasi et al. 2012).

We see that even though the pitch-angle scattering is inefficient in
incompressible trans-Alfvén MHD turbulence, the mirroring effect due to the
presence of pseudo-Alfvén modes significantly suppresses the parallel
diffusion of CRs.

Figure 4: Parallel diffusion coefficient $D_{\|,\text{inc},b}$ of bouncing CRs
as a function of CR energy in incompressible MHD turbulence. The numerical and
analytical (Eq. (47)) results are obtained using anisotropic and isotropic
pitch angle distribution, respectively. The dashed line indicates the energy,
below which the diffusion is affected by damping of turbulence.

Figure 5:  (a) Normalized particle distribution function $f^{\prime}$, (b)
$\mu$-dependent parallel diffusion coefficient $D_{\|,\text{inc},b}(\mu)$, and
(c) $D_{\|,\text{inc},b}(\mu)f^{\prime}$ vs. $\mu$ for $4$ GeV (solid line),
$10^{2}$ GeV (dotted line), $10^{3}$ GeV (dashed line), and $10^{5}$ GeV
(dash-dotted line) bouncing CRs in incompressible MHD turbulence.

###  4.2 Super- and sub-Alfvénic turbulence

In the above analysis, we discussed the case of trans-Alfvénic turbulence with
$M_{A}=1$, or equivalently, $\delta B_{A}=B_{0}$. Our analysis can be
generalized for both super- and sub-Alfvénic turbulence.

For super-Alfvénic turbulence ($M_{A}>1$), the perpendicular superdiffusion
(Lazarian & Vishniac, 1999; Yan & Lazarian, 2008; Lazarian, 2006; Lazarian &
Yan, 2014) takes place on scales less than $l_{A}$. 999On scales larger than
$l_{A}$, CRs following magnetic field lines undergo isotropic diffusion with
the step size equal to $l_{A}$ (see Brunetti & Lazarian 2007). For the
analysis in Section 4.1 for anisotropic MHD turbulence to be applicable, we
use $l_{A}$ (Eq. (A4)) as the effective injection scale of turbulence with the
injected turbulent velocity equal to $V_{A}$. Therefore, after replacing $L$
by $l_{A}$, we can still use the expressions in Section 4.1 for describing the
parallel diffusion of bouncing particles on scales smaller than $l_{A}$ in
super-Alfvénic turbulence. We note that $\delta B_{s}$ in this case should be
measured at $l_{A}$ instead of $L$.

In sub-Alfvénic turbulence with $M_{A}<1$, as we discuss in Appendix A, the
so-called weak turbulence (Lazarian & Vishniac, 1999; Galtier et al., 2000)
exists on scales from $L$ down to $l_{\text{tran}}$. Strong MHD turbulence is
only developed on scales below $l_{\text{tran}}$. In the strong turbulence
regime in sub-Alfvénic turbulence, the turbulence scaling is different from
that in trans-Alfvénic turbulence, and the turbulent eddies are more
elongated, with the elongation depending on $M_{A}$ (Lazarian & Vishniac,
1999). We notice in Appendix A, if we introduce an effective injection scale
$L_{\text{eff}}$ for sub-Alfvénic turbulence that is given by Eq. (A9), then
we can describe the sub-Alfvénic turbulence at scales less than
$l_{\text{tran}}$ as arising from the fictitious driving with the injection
scale $L_{\text{eff}}$ and injection velocity $V_{A}$. $L_{\text{eff}}$ is
larger than the actual injection scale $L$ by a factor of $M_{A}^{-4}$. Unlike
$l_{A}$ in super-Alfvénic turbulence, where the transition from isotropic to
anisotropic turbulence occurs, $L_{\text{eff}}$ does not have a well-defined
physical meaning. Its introduction, nevertheless, allows us to generalize the
results for trans-Alfvénic turbulence to the strong turbulence regime of sub-
Alfvénic turbulence by replacing $L$ used in Section 4.1 by $L_{\text{eff}}$.
$\delta B_{s}$ in this case should be measured at $L_{\text{eff}}$, and it is
related to the magnetic perturbation $\delta B_{s,L}$ of pseudo-Alfvén modes
at $L$ by

$\delta B_{s}=\delta B_{s,L}M_{A}^{-2}.$ (51)

In Table 1, we summarize the parameters introduced for generalizing the
results on diffusion of bouncing particles in trans-Alfvénic turbulence to
super- and sub-Alfvénic turbulence.

Table 1:  | $L_{\text{eff}}$ |  $\delta B_{s}$ at $L_{\text{eff}}$ | Range of length scales
---|---|---|---
Super-Alfvénic turbulence ($M_{A}>1$) | $l_{A}=LM_{A}^{-3}$ |  $\delta B_{s}$ at $l_{A}$ | $<l_{A}$
Sub-Alfvénic turbulence ($M_{A}<1$) | $LM_{A}^{-4}$ | $\delta B_{s,L}M_{A}^{-2}$ | $<LM_{A}^{2}$

##  5 Averaged parallel diffusion coefficients on scales $\gg\lambda_{\|}$

The analysis on the parallel diffusion of bouncing particles in Sections 3 and
4 can be applied to the bouncing CRs near a CR source. On scales much larger
than $\lambda_{\|}$ of both bouncing and non-bouncing particles, we need to
consider their exchange and the isotropization of pitch angle distribution due
to scattering. In this case we define an averaged parallel diffusion
coefficient for all particles as

$D_{\|,\text{tot}}\approx\alpha D_{\|,b}+(1-\alpha)D_{\|,nb},$ (52)

where

$\alpha=\frac{\tau_{b}}{\tau_{b}+\tau_{nb}},~{}~{}1-\alpha=\frac{\tau_{nb}}{\tau_{b}+\tau_{nb}},$
(53)

$\tau_{b}$ and $\tau_{nb}$ represent the times for particles to stay in the
bouncing state and non-bouncing state, respectively. As a simple estimate, we
have

$\tau_{b}\approx\frac{\mu_{c}^{2}}{D_{\mu\mu}},~{}~{}\tau_{nb}\approx\frac{1-\mu_{c}^{2}}{D_{\mu\mu}}$
(54)

to account for the diffusion in pitch angle by scattering and the resulting
transition from bouncing (non-bouncing) to non-bouncing (bouncing) particles.
For the weakly anisotropic pitch angle distribution under consideration, we
can use $D_{\mu\mu}$ at $\mu_{c}$ in the above expression as an approximation.
Therefore, Eq. (52) can be rewritten as

$D_{\|,\text{tot}}\approx\mu_{c}^{2}D_{\|,b}+(1-\mu_{c}^{2})D_{\|,nb}.$ (55)

Over a timescale $T$ much longer than $\tau_{b}$ and $\tau_{nb}$, the mean
squared displacement of CRs in the bouncing state is

$\Delta_{b}^{2}=D_{\|,b}\alpha T,$ (56)

and the mean squared displacement of CRs in the non-bouncing state is

$\Delta_{nb}^{2}=D_{\|,nb}(1-\alpha)T.$ (57)

The total mean squared displacement
$\Delta_{\text{tot}}^{2}=\Delta_{b}^{2}+\Delta_{nb}^{2}$ during $T$
corresponds to a total diffusion coefficient

$D_{\|,\text{tot}}=\frac{\Delta_{\text{tot}}^{2}}{T}\approx\alpha
D_{\|,b}+(1-\alpha)D_{\|,nb},$ (58)

which recovers Eq. (52). In the case with $D_{\|,b}\ll D_{\|,nb}$,
$D_{\|,\text{tot}}$ can be dominated by the diffusion of non-bouncing
particles, that is,

$D_{\|,\text{tot}}\approx(1-\alpha)D_{\|,nb}.$ (59)

It is different from the parallel diffusion coefficient in XL20 by a factor of
$1-\alpha$.

###  5.1 Compressible MHD turbulence with fast modes dominating magnetic
fluctuations

In Section 3, we derived the parallel diffusion coefficients of bouncing and
non-bouncing particles separately (see Fig. 3). Based on these calculations
and by using Eq. (55), here we present $D_{\|,f,\text{tot}}$ for all particles
interacting with fast modes in compressible MHD turbulence, as shown in Fig.
6. By comparing Fig. 6 with Fig. 3, we see that as $D_{\|,f,nb}$ is
considerably larger than $D_{\|,f,b}$ for the range of $E_{\text{CR}}$ under
consideration, the resulting $D_{\|,f,\text{tot}}$ is determined by
$D_{\|,f,nb}$ alone and slightly smaller than $D_{\|,f,nb}$ with the increase
of $\mu_{c}$ toward high energies. We note that unlike $D_{\|,f}$ usually
defined for pitch-angle scattering that corresponds to the change of pitch
angle by $90^{\circ}$, $D_{\|,f,nb}$ corresponds to a smaller change of pitch
angle and is smaller than $D_{\|,f}$ with $\mu_{c}=0$ (see Eq. (32)).

Figure 6:  Total parallel diffusion coefficient $D_{\|,f,\text{tot}}$ as a
function of $E_{\text{CR}}$ in compressible MHD turbulence with magnetic
fluctuations dominated by fast modes.

###  5.2 Compressible MHD turbulence with a varying fraction of fast modes

In realistic astrophysical media, the energy fraction of fast modes depends on
the driving condition of turbulence and plasma $\beta$ (see Section 2.1).

In the presence of both fast and slow modes, the ratio between their bouncing
rates (Eqs (16) and (40)) is

$\frac{\Gamma_{b,f}}{\Gamma_{b,s}}=\frac{\delta B_{f}^{4}}{B_{0}^{2}\delta
B_{s}^{2}}\mu^{-4},$ (60)

when $\mu>\mu_{\text{min,f}}$ and $\mu>\mu_{\text{min,s}}$. We use the larger
one, i.e., max$[\Gamma_{b,f},\Gamma_{b,s}]$, as the bouncing rate. In the case
of $\delta B_{f}\sim\delta B_{s}\sim B_{0}$, there is always
$\Gamma_{b,f}>\Gamma_{b,s}$ and fast modes dominate bouncing.

As illustrative examples, in Fig. 7, we consider the parallel diffusion
coefficients of bouncing and non-bouncing particles in trans-Alfvénic
turbulence with a varying fraction of fast modes, and we adopt an isotropic
pitch angle distribution. $\aleph_{s}$ is fixed at $0.5$ to address the role
of fast modes in affecting the diffusion coefficients.

(a) $\delta B_{f}=0$. In the limit case of incompressible MHD turbulence with
$\delta B_{f}=0$ (Fig. 7(a)), scattering is inefficient and $\mu_{c}$ is given
by Eq. (44) (see Fig. 8). $D_{\|,b}$ is consistent with Eq. (47). $D_{\|,nb}$
is

$D_{\|,nb}=\frac{v^{2}}{4}\int_{\mu_{c}}^{1}d\mu\frac{(1-\mu^{2})^{2}}{D_{\mu\mu}(\mu)},$
(61)

where $D_{\mu\mu}(\mu)$ is given by Eq. (48). As both
$D_{\mu\mu,\text{QLT},A}$ and $D_{\mu\mu,\text{QLT},s}$ have energy dependance
as $\propto E_{\text{CR}}^{3/2}$ (XL20), the resulting $D_{\|,nb}$ decreases
with increasing energy following $\propto E_{\text{CR}}^{-3/2}$. Because the
turbulence anisotropy is weaker toward larger length scales, higher-energy CRs
are more efficiently scattered and have smaller $D_{\|,nb}$.

(b) $\delta B_{f}=0.01B_{0}$. In compressible MHD turbulence, we consider the
scattering rate with contributions from all three modes,

$\Gamma_{s,\text{tot}}=\frac{2(D_{\mu\mu,\text{QLT},A}+D_{\mu\mu,\text{QLT},s}+D_{\mu\mu,\text{QLT},f})}{\mu^{2}}.$
(62)

When there is a tiny fraction of fast modes with $\delta B_{f}\ll B_{0}$ (see
Fig. 7(b)), the scattering of low-energy CRs is dominated by fast modes
despite their small energy fraction. With a constant $\mu_{c}$ (Fig. 8),
$D_{\|,nb}$ determined by $D_{\|,f,nb}$ increases with energy as $\propto
E_{\text{CR}}^{0.5}$ (see Eq. (32)). The scattering of higher-energy CRs is
taken over by Alfvén and slow modes. For bouncing CRs, $D_{\|,b}$ is
determined by slow modes and is the same as that in the incompressible
turbulence.

(c) $\delta B_{f}=0.1B_{0}$. When we further increase the fraction of fast
modes (Fig. 7(c)), scattering is dominated by fast modes for most energies,
resulting in the energy dependence of $\mu_{c}$ until it reaches its maximum
value (Fig. 8). With increasing $\mu_{c}$, $D_{\|,nb}$ has a weaker dependence
on $E_{\text{CR}}$ than $\propto E_{\text{CR}}^{0.5}$ (Eq. (32)) for low-
energy CRs. For bouncing of CRs, although $\Gamma_{b,f}$ is larger than
$\Gamma_{b,s}$ at a small $\mu$ for low-energy CRs (Eq. (60)), $D_{\|,b}$ is
still determined by slow modes for most energies and Eq. (47) applies.

(d) $\delta B_{f}=0.5B_{0}$. With a large fraction of fast modes, fast modes
dominate both scattering and bouncing. $D_{\|,nb}$ (see Eq. (32))
approximately follows $\propto E_{\text{CR}}^{1/3}$ (see Fig. 7(d)).
$D_{\|,b}$ in this case is given by Eq. (24) and is close to $\propto
E_{\text{CR}}$.

The thick solid lines in all cases represent the averaged $D_{\|,\text{tot}}$
(Eq. (55)). It is basically determined by $D_{\|,nb}$ as $D_{\|,nb}$ is
considerably larger than $D_{\|,b}$.

Based on the above results, our main findings are:

(i) The difference between $D_{\|,nb}$ dominated by incompressible modes and
fast modes comes from not only the different slopes of turbulent energy
spectra, but also the anisotropy of turbulence. The scale-dependent turbulence
anisotropy in the former case results in a decreasing $D_{\|,nb}$ with
increasing $E_{\text{CR}}$.

(ii) Even with a small fraction of fast modes, scattering becomes much more
efficient than that in incompressible turbulence.

(iii) With an energy-dependent $\mu_{c}$, $D_{\|,nb}$ dominated by fast modes
has a weaker dependence on $E_{\text{CR}}$ than the case with a constant
$\mu_{c}$.

(iv) With a non-zero $\mu_{c}$ and a factor $1-\alpha$, $D_{\|,nb}$ is in
general smaller than the parallel diffusion coefficient for gyroresonant
scattering with $\mu_{c}=0$ and $\alpha=0$.

(v) The maximum $D_{\|,b}$ depends on the maximum $\mu_{c}$, which is
determined by the amplitude of magnetic fluctuations ($\delta B_{s}$ or
$\delta B_{f}$) of the modes that dominate bouncing.

(vi) On scales larger than $\lambda_{\|}$ of both bouncing and non-bouncing
particles, the averaged parallel diffusion coefficient is determined by that
of non-bouncing particles.

(vii) The energy scaling of $D_{\|,\text{tot}}$ is non-universal, depending on
the energy fractions of different turbulence modes. In diverse astrophysical
environments, the energy fractions of turbulence modes vary with the
turbulence parameters, e.g., $M_{s}$, $M_{A}$ (Cho & Lazarian, 2003). For
realistic modeling of CR diffusion, prior knowledge of turbulence parameters
is necessary, which can be measured with the Gradient Techniques (e.g.,
Lazarian et al. 2018). This falls beyond the scope of this work and is left to
future studies.

(a) $\delta B_{f}/B_{0}=0$, $\aleph_{s}=0.5$

(b) $\delta B_{f}/B_{0}=0.01$, $\aleph_{s}=0.5$

(c) $\delta B_{f}/B_{0}=0.1$, $\aleph_{s}=0.5$

(d) $\delta B_{f}/B_{0}=0.5$, $\aleph_{s}=0.5$

Figure 7: Parallel diffusion coefficients of bouncing particles $D_{\|,b}$,
non-bouncing particles $D_{\|,nb}$, and the averaged parallel diffusion
coefficient $D_{\|,\text{tot}}$ in MHD turbulence with different fractions of
fast modes. The pitch angle distribution is assumed to be isotropic.  Figure
8:  $\mu_{c}$ vs. $E_{\text{CR}}$ for different fractions of fast modes.

##  6 Discussion

###  6.1 Toward a more comprehensive description of CR propagation

Important new findings on CR propagation have been reported over the past two
decades. First of all, the fast modes of MHD turbulence were identified as the
major scattering agent for diffusion of CRs with $r_{g}\ll L$ in the
interstellar medium (Yan & Lazarian, 2002, 2003). On the contrary, the Alfvén
and slow modes, which were traditionally considered for CR diffusion, were
found to be very inefficient in scattering (Chandran, 2000a; Yan & Lazarian,
2002). Therefore, it was believed that, when the CR streaming instability is
inefficient, e.g. due to Alfvenic turbulence suppression (see Yan & Lazarian
2002; Farmer & Goldreich 2004; Lazarian 2016), and fast modes are absent, CRs
propagate ballistically along magnetic field lines without scattering.
101010We note that in super-Alfvénic turbulence, CRs that follow the tangled
magnetic fields without scattering have an effective mean free path given by
$l_{A}$ (Brunetti & Lazarian, 2007).

In this work, we introduce a new diffusion mechanism for CRs that encounter
magnetic compressions induced by slow (pseudo-Alfvén) and fast modes. This
parallel mirror diffusion of CRs can take place even in incompressible MHD
turbulence and cause slower diffusion than that induced by scattering in both
incompressible and compressible MHD turbulence. In the idealized situation
with incompressible turbulence and the streaming instability suppressed by
Alfvénic turbulence, bouncing CRs can stay in the system for an extended
period of time until they become non-bouncing CRs due to scattering or
acceleration processes. This finding has important implications for CR
diffusion in high-$\beta$ intracluster media. In addition, the slow diffusion
of bouncing CRs near CR sources may provide a possible explanation for the
small diffusion coefficient of high-energy electrons and positrons around
pulsar wind nebulae suggested by recent observations (e.g., Abeysekara et al.
2017; Huang et al. 2018), as well as the steep high-energy CR spectra
indicated by gamma-ray observations of middle-aged supernova remnants (Xu
2021; see also, e.g., Nava et al. 2016; D’Angelo et al. 2016, 2018; Nava et
al. 2019, for suppressed diffusion of self-confined CRs around Galactic
sources).

In addition, the mirroring of CRs also affects the parallel diffusion
associated with scattering. In the presence of mirroring, the scattering even
in the quasilinear approximation does not face the $90^{\circ}$ problem
(Noerdlinger, 1968; Kulsrud & Pearce, 1969; Felice & Kulsrud, 2001), and also
the resulting mean free path corresponds to the pitch angle change over a
smaller range $[\mu_{c},1]$, as pointed out by Cesarsky & Kulsrud (1973). In
XL20, we derived the parallel diffusion coefficient for scattering by fast
modes in the presence of bouncing. With an energy-dependent $\mu_{c}$, the
diffusion coefficient has a energy dependence close to $\propto E^{1/3}$
instead of $\propto E^{0.5}$. In earlier studies (e.g., Ptuskin et al. 2006),
the energy dependence $\propto E^{1/3}$ of diffusion coefficient is only
expected for the Kolmogorov spectrum of turbulence, but the Alfvén and slow
modes with the Kolmogorov spectrum are incapable of scattering. Our finding
provides a possible solution to this theoretical problem and a possible
explanation for the energy scaling of diffusion coefficient of high-energy CRs
indicated by observations (e.g., Blasi et al. 2012; Fornieri et al. 2021).

###  6.2 Mirroring and acceleration

The mirroring of CRs also induces a new acceleration mechanism. The mirroring
does not change a particle’s pitch angle if the magnetic mirrors are static.
When encountering moving magnetic mirrors in MHD turbulence, CRs undergo the
second order Fermi acceleration with the stochastic increase of parallel
momentum $p_{\|}$ and $\mu$. Together with other acceleration mechanisms,
e.g., TTD that also causes the stochastic increase of $\mu$ (Section 6.4), the
mirroring can enhance the diffusion in $\mu$ and the transition of particles
from bouncing to non-bouncing state. This new acceleration mechanism and its
implications will be studied in our future work.

The diffusive propagation of CRs plays an important role in the first order
Fermi acceleration processes at shocks (e.g.,Schlickeiser & Oppotsch 2017) and
in regions of turbulent magnetic reconnection (de Gouveia dal Pino & Lazarian,
2005; Lazarian, 2005; Kowal et al., 2012). Both the perpendicular
superdiffusion and the parallel diffusion resulting from mirroring and
scattering should be taken into account when studying the acceleration of CRs.

###  6.3 Mirroring and instabilities

Mirroring and resonant scattering by MHD turbulence are not the only processes
that affect the diffusion of CRs. The fluid of low-energy CRs is subject to
instabilities. The streaming instability is the most studied one (Kulsrud,
2005). This instability arises when CRs move preferentially in the same
direction. The Alfvénic fluctuations that are induced are different from MHD
turbulence. These are resonant waves that effectively scatter CRs and prevent
free escape of particles. Streaming instability can be important for
explaining the energy scaling of diffusion coefficient of low-energy CRs
indicated by AMS-02 measurements (e.g., Fornieri et al. 2021).

Apart from the streaming instability, CRs can face other types of
instabilities, the importance of which is less explored. For instance, the
compression of CRs by turbulence induces anisotropy in the distribution of CRs
in momentum space. This can result in the gyroresonance, mirror, and firehose
collisionless instabilities of CRs (Lazarian & Beresnyak, 2006). In addition,
the faster escape of non-bouncing particles from the vicinity of CRs sounce
changes the momentum distribution toward having more CRs with $\mu<\mu_{c}$.
This difference in the diffusion coefficients of bouncing and non-bouncing CRs
is another source of momentum space anisotropy, which can result in an
instability and induce resonant Alfvénic waves. The list above does not
exhaust the list of the CR instabilities that can be present in the vicinity
of CR sources. For instance the current instability (Bell, 2004) was
introduced in the context of shock acceleration. It is also expected to be
present beyond the shock context.

The instabilities, along with the mirroring, can significantly suppress the
diffusion of CRs in the vicinity of their sources. The environmental effects,
e.g. ambient turbulence, gas ionization, have different influence on mirroring
and the scattering caused by instabilities. We note that the CR streaming
instability can be suppressed in galactic disks through both ion-neutral
collisional damping (Kulsrud & Pearce, 1969; Xu et al., 2016; Krumholz et al.,
2020) and the interaction with Alfvénic turbulence (Yan & Lazarian, 2002;
Farmer & Goldreich, 2004; Lazarian, 2016). More severe damping of streaming
instability is expected at a larger $M_{A}$ (Lazarian, 2016), while the
magnetic compressions that cause mirroring increase with $M_{A}$ (Hu,
Lazarian, & Xu, in prep). Therefore, the relative importance between streaming
instability and mirroring effect on diffusion of low-energy CRs depends on
$M_{A}$, which varies in the multi-phase interstellar medium (Lazarian et al.,
2018). The comparison between scattering by instabilities and mirroring, as
well as their interplay, requires further studies.

###  6.4 Mirroring and TTD

The TTD effect is an important component of CR dynamics (Schlickeiser, 2002;
Yan & Lazarian, 2002; Xu & Lazarian, 2018). The TTD interaction arises from
magnetic compressions in turbulent media. This is the common feature of TTD
and mirroring. Nevertheless, they are two different processes and should be
distinguished.

The resonant TTD interaction happens when the parallel phase speed of
compressible waves matches the parallel speed of the particle. Such
interactions can induce the second order Fermi acceleration of $p_{\|}$. The
resulting stochastic increase of $\mu$ associated with TTD acceleration should
be distinguished from the pitch-angle diffusion due to resonant scattering. In
contrast, the mirror diffusion does not require any resonance condition. It
can take place even in the limit case of stationary magnetic bottles, i.e.
$V_{A}\rightarrow 0$, where stochastic acceleration is not involved. As
another difference, mirroring can not happen at small pitch angles in the loss
cone, but there is no such constraint for TTD.

In brief, TTD is a stochastic acceleration process with acceleration-induced
diffusion in $\mu$. Mirror diffusion is a diffusion process in space, which
can be accompanied by stochastic acceleration.

###  6.5 Damping effect

Fast modes can be subject to ion-neutral collisional damping in partially
ionized interstellar phases (Xu et al., 2015, 2016; Xu & Lazarian, 2018) and
collisionless damping in galactic halo. The effect of damping of fast modes on
scattering of CRs has been addressed in earlier studies, e.g., Yan & Lazarian
(2002, 2004); Xu et al. (2016); Xu & Lazarian (2018). The effect of damping on
mirroring by fast modes will be addressed in our future work. In the case when
fast modes are severely damped, slow modes can dominate the bouncing of CRs
(see Sections 4 and 5.2).

##  7 Summary

To achieve a more comprehensive description of CR propagation, different
diffusion mechanisms depending on the properties of MHD turbulence should be
taken into account. In this work we identify the mirror diffusion as an
essential process of CR propagation.

The perpendicular superdiffusion of CRs originates from the superdiffusion of
magnetic fields induced by Alfvénic turbulence. This effect makes the trapping
of CR within magnetic bottles improbable. Instead, CRs bounce with different
magnetic mirrors and move diffusively parallel to the turbulent magnetic
field. As a result, both magnetic mirroring and pitch-angle scattering
contribute to the parallel diffusion. The former governs the parallel
diffusion of CRs at large pitch angles with $\mu<\mu_{c}$, and the latter is
dominant at $\mu>\mu_{c}$, where $\mu_{c}$ corresponds to the balance between
mirroring and scattering. The two diffusion processes interact with each other
through exchanging CRs via the pitch angle diffusion.

In the case when fast modes dominate both scattering and mirroring, $\mu_{c}$
increases with $E_{\text{CR}}$. As a result, the diffusion coefficient for
scattering by fast modes has a weaker dependence on $E_{\text{CR}}$ compared
to the case without mirroring, i.e., $\mu_{c}=0$. It is not an exact power law
function of $E_{\text{CR}}$, but can be approximated by $\propto
E_{\text{CR}}^{1/3}$. The scaling $\propto E_{\text{CR}}^{1/3}$ is expected
for isotropic turbulence with the Kolmogorov spectrum. However, the scattering
by Alfvén and slow modes with the Kolmogorov spectrum is inefficient and the
resulting diffusion coefficient depends on $E_{\text{CR}}$ as $\propto
E_{\text{CR}}^{-3/2}$ due to the scale-dependent turbulence anisotropy. Our
finding provides the physical justification for the commonly used energy
scaling of diffusion coefficient.

With the mirroring of CRs taken into account, the $90^{\circ}$ problem of
quasilinear gyroresonant scattering can be solved. Moreover, as the mean free
path of bouncing CRs is determined by the size of compressive magnetic
fluctuations, which cannot exceed the driving scale of turbulence, the
corresponding diffusion is slow. The energy scaling of the diffusion
coefficient of bouncing CRs depends on the anisotropy of pitch angle
distribution and the energy dependence of $\mu_{c}$.

In the vicinity of a CR source, the injected CRs can have an anisotropic pitch
angle distribution. Note that the damping of turbulence can also affect the
anisotropic distribution. If there are more particles at larger pitch angles,
the diffusion of bouncing CRs can be further suppressed.

For the galactic environment with a nearly uniform and isotropic distribution
of CRs, the average diffusion coefficient of bouncing and non-bouncing CRs on
scales much larger than their mean free paths is usually determined by the
diffusion coefficient of the latter, as it is larger than that of bouncing
CRs. Its energy scaling can be non-universal, depending on the properties of
interstellar turbulence.

S.X. acknowledges the support for this work provided by NASA through the NASA
Hubble Fellowship grant # HST-HF2-51473.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

##  Appendix A Basic scalings of MHD turbulence

###  A.1 Scalings of incompressible MHD turbulence

MHD turbulence is a major agent determining the dynamics of CRs. The
turbulence in magnetized media in typical astrophysical settings is injected
at the scale $L$ with the injection velocity $V_{L}$, and the turbulent energy
then cascades down to smaller scales. If the injection happens with
$V_{L}>V_{A}$, where $V_{A}$ is the Alfvén speed, this is the case of super-
Alfvénic turbulence. If $V_{L}<V_{A}$, the turbulence is sub-Alfvénic. The
ratio $M_{A}=V_{L}/V_{A}$ is the Alfvén Mach number. $M_{A}=1$ corresponds to
the trans-Alfvénic turbulence.

Note that we consider turbulence as a result of energy cascade. Therefore the
Alfvén waves excited by CR instabilities, e.g. streaming instabilities (see
Farmer & Goldreich 2004), gyro-resonance instabilities (see Lazarian &
Beresnyak 2006), are not classified as turbulence.

The theory of trans-Alfvénic turbulence was developed by Goldreich & Sridhar
(1995) (henceforth GS95) in the global system of reference with respect to the
mean magnetic field. Lazarian & Vishniac (1999) (henceforth LV99) presented an
alternative way to derive the GS95 theory based on the theory of turbulent
reconnection of magnetic fields. More importantly, they found that the GS95
theory is only valid in the “local system of reference” with respect to the
local mean magnetic field averaged over the length scale of interest, as
confirmed by numerical simulations (Cho & Vishniac, 2000; Maron & Goldreich,
2001; Cho et al., 2002). According to the LV99 theory, magnetic reconnection
takes place within one eddy turnover time, which allows the turbulent motions
to mix up magnetic fields without bending them. The unconstrained eddy
rotation occurs in the system of reference of the local magnetic field
averaged over the eddy scale, and the eddy is aligned with the local magnetic
field. The direction of the local magnetic field averaged over a small scale
can differ significantly from the global mean magnetic field direction
averaged over a large scale.

The turbulent cascade in the direction perpendicular to the local magnetic
field remains Kolmogorov-like with

$v_{l}\approx V_{A}\left(\frac{l_{\bot}}{L}\right)^{1/3},$ (A1)

where it is taken into account that $V_{L}=V_{A}$ for trans-Alfvénic
turbulence, $v_{l}$ is the turbulent velocity at $l_{\perp}$, and $l_{\perp}$
is the perpendicular size of a turbulent eddy. The eddy rotation perpendicular
to the local magnetic field induces Alfvénic perturbation that propagates
along the magnetic field with the speed $V_{A}$. The timescale of this
perturbation $l_{\|}/V_{A}$ should be equal to the eddy turnover time
$l_{\perp}/v_{l}$. The corresponding relation between the parallel and
perpendicular sizes of the eddy

$\frac{l_{\|}}{V_{A}}\approx\frac{l_{\perp}}{v_{l}},$ (A2)

was termed as the “critical balance” in GS95. By combining Eq. (A1) and Eq.
(A2), one can obtain the scale-dependent anisotropy of trans-Alfvénic
turbulence

$l_{\|}\approx L\left(\frac{l_{\perp}}{L}\right)^{2/3}.$ (A3)

Smaller eddies are more elongated along the local magnetic field. Note that
Eqs. (A1) and (A3) should be understood in the statistical sense. They
represent the scaling relations between the most probable values of the
quantities involved. For instance, on the basis of numerical simulations, Cho
et al. (2002) provided an analytical fit for the detailed distribution
function describing the probability of finding a $l_{\perp}$ at a given
$l_{\|}$. This was used later in Yan & Lazarian (2002, 2004) and subsequent
studies on CR scattering.

For super-Alfvénic turbulence injected with $M_{A}>1$, the magnetic field is
of marginal importance at $L$. The super-Alfvénic turbulence is initially
hydrodynamic-like with the isotropic Kolmogorov energy spectrum. With the
decrease of turbulent velocity along the energy cascade, $v_{l}\sim
V_{L}(l/L)^{1/3}$, where $v_{l}$ is the turbulent velocity at the length scale
$l$, the effect of magnetic field becomes more and more manifested.
Eventually, at the scale (Lazarian, 2006)

$l_{A}=LM_{A}^{-3},$ (A4)

$v_{l}$ becomes equal to $V_{A}$, and the turbulence becomes fully
magnetohydrodynamic. To describe the MHD cascade on scales less than $l_{A}$
in super-Alfvénic turbulence, $L$ in Eqs. (A1) and (A3) should be replaced by
$l_{A}$.

For sub-Alfvénic turbulence with $M_{A}<1$, it was shown in LV99 that at $L$
the turbulence is weak, and the parallel scale of wave packets remains
unchanged, i.e. $l_{\|}=L$. The scaling obtained in LV99 for the weak
turbulence under the assumption of the isotropic turbulence driven at $L$ is

$v_{l}\approx V_{L}\left(\frac{l_{\perp}}{L}\right)^{1/2},$ (A5)

and this result was supported by the subsequent study by Galtier et al.
(2000). With the decrease of $l_{\perp}$, the intensity of interactions of
Alfvénic perturbations increases. At a scale (LV99)

$l_{\text{tran}}\approx LM_{A}^{2},$ (A6)

where $M_{A}<1$, the turbulence gets strong. For the sub-Alfvénic MHD
turbulence at $l<l_{\text{tran}}$, LV99 derived the relations

$v_{l}\approx V_{L}\left(\frac{l_{\bot}}{L}\right)^{1/3}M_{A}^{1/3},$ (A7)

and

$l_{\|}\approx L\left(\frac{l_{\bot}}{L}\right)^{2/3}M_{A}^{-4/3},$ (A8)

which at $M_{A}=1$ can recover the relations in Eqs. (A1) and (A3) for trans-
Alfvénic turbulence.

To generalize our analysis for diffusion of bouncing particles in trans-
Alfvénic turbulence to super- and sub-Alfvénic turbulence, we introduce the
following effective injections scales.

(1) Super-Alfvénic turbulence. As discussed above, the super-Alfvénic
turbulence at $l<l_{A}$ is equivalent to the trans-Alfvénic turbulence
injected at $l_{A}$. So by using $l_{A}$ as the effective injection scale to
replace $L$ in trans-Alfvénic turbulence, the results derived for trans-
Alfvénic are still applicable to super-Alfvénic on scales less than $l_{A}$.

(2) Sub-Alfvénic turbulence. By introducing an effective injection scale

$L_{\text{eff}}=LM_{A}^{-4},$ (A9)

Eqs. (A7) and (A8) can be rewritten as

$v_{l}\approx V_{A}\left(\frac{l_{\bot}}{L_{\text{eff}}}\right)^{1/3},$ (A10)

and

$l_{\|}\approx
L_{\text{eff}}\left(\frac{l_{\bot}}{L_{\text{eff}}}\right)^{2/3},$ (A11)

which take the same forms as Eqs. (A1) and (A3) with $L$ replaced by
$L_{\text{eff}}$. Therefore, by using $L_{\text{eff}}$ instead of $L$, the
expressions for trans-Alfvénic turbulence can be applied to sub-Alfvénic
turbulence on scales less than $l_{\text{tran}}$.

###  A.2 Magnetic compressions in MHD turbulence

The magnetic mirroring effect arise from the magnetic compressions in MHD
turbulence, which are associated with slow and fast modes. We use the term
“modes” rather than “waves”, as the nonlinear interactions intrinsic to the
turbulent cascade make the properties of magnetic fluctuations different from
the properties of linear waves. We stress that the compression of magnetic
field can also occur in incompressible MHD turbulence. The incompressible
limit of slow modes, i.e., pseudo-Alfvén modes (see GS95), is a typical
example for turbulent compression of magnetic field in an incompressible
medium.

It was pointed out in GS95 that Alfvén modes impose their scaling on the
pseudo-Alfvén modes. With a finite compressibility, slow modes also follow the
same scaling as Alfvén modes. This property was discussed in Lithwick &
Goldreich (2001) and confirmed numerically by Cho & Lazarian (2002b, 2003). As
discussed above, the anisotropic scaling of slow modes is defined and should
be measured in the local system of reference. For bouncing CRs that interact
with the magnetic compression generated by slow modes, they only feel the
local magnetic field averaged over the scale of the magnetic compression. When
describing the gyroresonant scattering of non-bouncing CRs, they also only
feel the local magnetic field averaged over the scale comparable to $r_{L}$.
The QLT for scattering is formally only applicable to infinitesimally small
magnetic perturbations. In the local system of reference, CRs only interact
with the small magnetic fluctuations on small scales, while the large
fluctuations on large scales are not involved. This extends the applicability
of the QLT.

Fast modes are different from Alfvén and slow modes. Their evolution is not
related to the local system of reference. Therefore, the traditional wave
representation is applicable to fast modes. The scaling of fast modes is
somewhat less certain. The theoretical considerations in LG01 for high-$\beta$
plasma and in Cho & Lazarian (2002b, 2003) for low-$\beta$ plasma suggested
the “isotropic” energy spectrum of fast modes with $k^{-3/2}$ similar to the
acoustic-type turbulence. However, the simulations in Kowal & Lazarian (2010)
indicated a steeper spectrum close to $k^{-2}$ for fast modes. This
discrepancy may be attributed to the effect of shocks in the latter study.

## References

  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911, doi: 10.1126/science.aan4880
  * Amato & Casanova (2021) Amato, E., & Casanova, S. 2021, Journal of Plasma Physics, 87, 845870101, doi: 10.1017/S0022377821000064
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550, doi: 10.1111/j.1365-2966.2004.08097.x
  * Beresnyak (2013) Beresnyak, A. 2013, ApJ, 767, L39, doi: 10.1088/2041-8205/767/2/L39
  * Beresnyak (2014) —. 2014, ApJ, 784, L20, doi: 10.1088/2041-8205/784/2/L20
  * Beresnyak & Lazarian (2019) Beresnyak, A., & Lazarian, A. 2019, Turbulence in Magnetohydrodynamics
  * Blasi et al. (2012) Blasi, P., Amato, E., & Serpico, P. D. 2012, Phys. Rev. Lett., 109, 061101, doi: 10.1103/PhysRevLett.109.061101
  * Brunetti & Jones (2014) Brunetti, G., & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007, doi: 10.1142/S0218271814300079
  * Brunetti & Lazarian (2007) Brunetti, G., & Lazarian, A. 2007, MNRAS, 378, 245, doi: 10.1111/j.1365-2966.2007.11771.x
  * Budker (1959) Budker, G. I. 1959, in Plasma Physics and the Problem of Controlled Thermonuclear Reactions, Ed. by M. A. Leontovich (Pergamon, New York), Vol. 1
  * Cesarsky & Kulsrud (1973) Cesarsky, C. J., & Kulsrud, R. M. 1973, ApJ, 185, 153
  * Chandran (2000a) Chandran, B. D. G. 2000a, Physical Review Letters, 85, 4656, doi: 10.1103/PhysRevLett.85.4656
  * Chandran (2000b) Chandran, B. D. G. 2000b, ApJ, 529, 513, doi: 10.1086/308232
  * Chepurnov et al. (2010) Chepurnov, A., Lazarian, A., Stanimirović, S., Heiles, C., & Peek, J. E. G. 2010, ApJ, 714, 1398, doi: 10.1088/0004-637X/714/2/1398
  * Cho & Lazarian (2002a) Cho, J., & Lazarian, A. 2002a, ApJ, 575, L63, doi: 10.1086/342722
  * Cho & Lazarian (2002b) —. 2002b, Physical Review Letters, 88, 245001, doi: 10.1103/PhysRevLett.88.245001
  * Cho & Lazarian (2003) —. 2003, MNRAS, 345, 325, doi: 10.1046/j.1365-8711.2003.06941.x
  * Cho et al. (2012) Cho, J., Lazarian, A., & Timbie, P. T. 2012, ApJ, 749, 164, doi: 10.1088/0004-637X/749/2/164
  * Cho et al. (2002) Cho, J., Lazarian, A., & Vishniac, E. T. 2002, ApJ, 564, 291, doi: 10.1086/324186
  * Cho & Vishniac (2000) Cho, J., & Vishniac, E. T. 2000, ApJ, 539, 273, doi: 10.1086/309213
  * Cohet & Marcowith (2016) Cohet, R., & Marcowith, A. 2016, A&A, 588, A73, doi: 10.1051/0004-6361/201527376
  * D’Angelo et al. (2016) D’Angelo, M., Blasi, P., & Amato, E. 2016, Phys. Rev. D, 94, 083003, doi: 10.1103/PhysRevD.94.083003
  * D’Angelo et al. (2018) D’Angelo, M., Morlino, G., Amato, E., & Blasi, P. 2018, MNRAS, 474, 1944, doi: 10.1093/mnras/stx2828
  * de Gouveia dal Pino & Lazarian (2005) de Gouveia dal Pino, E. M., & Lazarian, A. 2005, A&A, 441, 845, doi: 10.1051/0004-6361:20042590
  * Demidem et al. (2020) Demidem, C., Lemoine, M., & Casse, F. 2020, Phys. Rev. D, 102, 023003, doi: 10.1103/PhysRevD.102.023003
  * Evoli & Yan (2014) Evoli, C., & Yan, H. 2014, ApJ, 782, 36, doi: 10.1088/0004-637X/782/1/36
  * Eyink et al. (2013) Eyink, G., Vishniac, E., Lalescu, C., et al. 2013, Nature, 497, 466, doi: 10.1038/nature12128
  * Eyink et al. (2011) Eyink, G. L., Lazarian, A., & Vishniac, E. T. 2011, ApJ, 743, 51, doi: 10.1088/0004-637X/743/1/51
  * Farmer & Goldreich (2004) Farmer, A. J., & Goldreich, P. 2004, ApJ, 604, 671, doi: 10.1086/382040
  * Felice & Kulsrud (2001) Felice, G. M., & Kulsrud, R. M. 2001, ApJ, 553, 198, doi: 10.1086/320651
  * Fermi (1949) Fermi, E. 1949, Physical Review, 75, 1169, doi: 10.1103/PhysRev.75.1169
  * Fisk et al. (1974) Fisk, L. A., Goldstein, M. L., Klimas, A. J., & Sandri, G. 1974, ApJ, 190, 417, doi: 10.1086/152893
  * Forman et al. (2011) Forman, M. A., Wicks, R. T., & Horbury, T. S. 2011, ApJ, 733, 76, doi: 10.1088/0004-637X/733/2/76
  * Fornieri et al. (2021) Fornieri, O., Gaggero, D., Cerri, S. S., De La Torre Luque, P., & Gabici, S. 2021, MNRAS, 502, 5821, doi: 10.1093/mnras/stab355
  * Gabici et al. (2019) Gabici, S., Evoli, C., Gaggero, D., et al. 2019, arXiv e-prints, arXiv:1903.11584.  https://arxiv.org/abs/1903.11584
  * Galtier et al. (2000) Galtier, S., Nazarenko, S. V., Newell, A. C., & Pouquet, A. 2000, Journal of Plasma Physics, 63, 447, doi: 10.1017/S0022377899008284
  * Giacalone & Jokipii (1999) Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204, doi: 10.1086/307452
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763, doi: 10.1086/175121
  * Guo & Oh (2008) Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251, doi: 10.1111/j.1365-2966.2007.12692.x
  * Holguin et al. (2019) Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H. Y. K. 2019, MNRAS, 490, 1271, doi: 10.1093/mnras/stz2568
  * Horbury et al. (2008) Horbury, T. S., Forman, M., & Oughton, S. 2008, Phys. Rev. Lett., 101, 175005, doi: 10.1103/PhysRevLett.101.175005
  * Huang et al. (2018) Huang, Z.-Q., Fang, K., Liu, R.-Y., & Wang, X.-Y. 2018, ApJ, 866, 143, doi: 10.3847/1538-4357/aadfed
  * Ipavich (1975) Ipavich, F. M. 1975, ApJ, 196, 107, doi: 10.1086/153397
  * Jokipii (1966) Jokipii, J. R. 1966, ApJ, 146, 480, doi: 10.1086/148912
  * Jokipii (1971) —. 1971, Reviews of Geophysics and Space Physics, 9, 27, doi: 10.1029/RG009i001p00027
  * Klepach & Ptuskin (1995) Klepach, E. G., & Ptuskin, V. S. 1995, Astronomy Letters, 21, 411
  * Kóta & Jokipii (2000) Kóta, J., & Jokipii, J. R. 2000, ApJ, 531, 1067, doi: 10.1086/308492
  * Kowal et al. (2012) Kowal, G., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, Physical Review Letters, 108, 241102, doi: 10.1103/PhysRevLett.108.241102
  * Kowal et al. (2017) Kowal, G., Falceta-Gonçalves, D. A., Lazarian, A., & Vishniac, E. T. 2017, ApJ, 838, 91, doi: 10.3847/1538-4357/aa6001
  * Kowal & Lazarian (2010) Kowal, G., & Lazarian, A. 2010, ApJ, 720, 742, doi: 10.1088/0004-637X/720/1/742
  * Kowal et al. (2009) Kowal, G., Lazarian, A., Vishniac, E. T., & Otmianowska-Mazur, K. 2009, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 36, Revista Mexicana de Astronomia y Astrofisica Conference Series, 89–96.  https://arxiv.org/abs/0812.2024
  * Krumholz et al. (2020) Krumholz, M. R., Crocker, R. M., Xu, S., et al. 2020, MNRAS, 493, 2817, doi: 10.1093/mnras/staa493
  * Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445, doi: 10.1086/149981
  * Kulsrud (2005) Kulsrud, R. M. 2005, Plasma physics for astrophysics, ed. R. M. Kulsrud
  * Lazarian (2005) Lazarian, A. 2005, in American Institute of Physics Conference Series, Vol. 784, Magnetic Fields in the Universe: From Laboratory and Stars to Primordial Structures., ed. E. M. de Gouveia dal Pino, G. Lugones, & A. Lazarian, 42–53, doi: 10.1063/1.2077170
  * Lazarian (2006) Lazarian, A. 2006, ApJ, 645, L25, doi: 10.1086/505796
  * Lazarian (2016) —. 2016, ApJ, 833, 131, doi: 10.3847/1538-4357/833/2/131
  * Lazarian & Beresnyak (2006) Lazarian, A., & Beresnyak, A. 2006, MNRAS, 373, 1195, doi: 10.1111/j.1365-2966.2006.11093.x
  * Lazarian & Vishniac (1999) Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700, doi: 10.1086/307233
  * Lazarian et al. (2004) Lazarian, A., Vishniac, E. T., & Cho, J. 2004, ApJ, 603, 180, doi: 10.1086/381383
  * Lazarian & Yan (2014) Lazarian, A., & Yan, H. 2014, ApJ, 784, 38, doi: 10.1088/0004-637X/784/1/38
  * Lazarian & Yan (2019) —. 2019, ApJ, 885, 170, doi: 10.3847/1538-4357/ab50ba
  * Lazarian et al. (2018) Lazarian, A., Yuen, K. H., Ho, K. W., et al. 2018, ApJ, 865, 46, doi: 10.3847/1538-4357/aad7ff
  * Lemoine & Malkov (2020) Lemoine, M., & Malkov, M. A. 2020, MNRAS, 499, 4972, doi: 10.1093/mnras/staa3131
  * Lithwick & Goldreich (2001) Lithwick, Y., & Goldreich, P. 2001, ApJ, 562, 279, doi: 10.1086/323470
  * López-Barquero et al. (2016) López-Barquero, V., Farber, R., Xu, S., Desiati, P., & Lazarian, A. 2016, ApJ, 830, 19, doi: 10.3847/0004-637X/830/1/19
  * Luo & Wu (2010) Luo, Q. Y., & Wu, D. J. 2010, ApJ, 714, L138, doi: 10.1088/2041-8205/714/1/L138
  * Lynn et al. (2012) Lynn, J. W., Parrish, I. J., Quataert, E., & Chandran, B. D. G. 2012, ApJ, 758, 78, doi: 10.1088/0004-637X/758/2/78
  * Maron & Goldreich (2001) Maron, J., & Goldreich, P. 2001, ApJ, 554, 1175, doi: 10.1086/321413
  * MATLAB (2018) MATLAB. 2018, 9.7.0.1190202 (R2019b) (Natick, Massachusetts: The MathWorks Inc.)
  * Matthaeus et al. (1990) Matthaeus, W. H., Goldstein, M. L., & Roberts, D. A. 1990, J. Geophys. Res., 95, 20673, doi: 10.1029/JA095iA12p20673
  * Matthaeus et al. (2003) Matthaeus, W. H., Qin, G., Bieber, J. W., & Zank, G. P. 2003, ApJ, 590, L53, doi: 10.1086/376613
  * Nava & Gabici (2013) Nava, L., & Gabici, S. 2013, MNRAS, 429, 1643, doi: 10.1093/mnras/sts450
  * Nava et al. (2016) Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552, doi: 10.1093/mnras/stw1592
  * Nava et al. (2019) Nava, L., Recchia, S., Gabici, S., et al. 2019, MNRAS, 484, 2684, doi: 10.1093/mnras/stz137
  * Noerdlinger (1968) Noerdlinger, P. D. 1968, Phys. Rev. Lett., 20, 1513, doi: 10.1103/PhysRevLett.20.1513
  * Orlando (2018) Orlando, E. 2018, MNRAS, 475, 2724, doi: 10.1093/mnras/stx3280
  * Padovani et al. (2018) Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111, doi: 10.1051/0004-6361/201732202
  * Palmer (1982) Palmer, I. D. 1982, Reviews of Geophysics and Space Physics, 20, 335, doi: 10.1029/RG020i002p00335
  * Parker (1965) Parker, E. N. 1965, Planet. Space Sci., 13, 9, doi: 10.1016/0032-0633(65)90131-5
  * Post (1958) Post, R. F. 1958, Proc. of Second U.N. Int. Conf. on Peaceful Uses of Atomic Energy, Vol. 32, Paper A/Conf. 15/P/377, Geneva, pp. 245-265
  * Ptuskin et al. (2006) Ptuskin, V. S., Moskalenko, I. V., Jones, F. C., Strong, A. W., & Zirakashvili, V. N. 2006, ApJ, 642, 902, doi: 10.1086/501117
  * Qin et al. (2002) Qin, G., Matthaeus, W. H., & Bieber, J. W. 2002, ApJ, 578, L117, doi: 10.1086/344687
  * Richardson (1926) Richardson, L. F. 1926, Proceedings of the Royal Society of London Series A, 110, 709, doi: 10.1098/rspa.1926.0043
  * Schlickeiser (2002) Schlickeiser, R. 2002, Cosmic Ray Astrophysics, ed. R. Schlickeiser
  * Schlickeiser et al. (2016) Schlickeiser, R., Caglar, M., & Lazarian, A. 2016, ApJ, 824, 89, doi: 10.3847/0004-637X/824/2/89
  * Schlickeiser & Miller (1998) Schlickeiser, R., & Miller, J. A. 1998, ApJ, 492, 352, doi: 10.1086/305023
  * Schlickeiser & Oppotsch (2017) Schlickeiser, R., & Oppotsch, J. 2017, ApJ, 850, 160, doi: 10.3847/1538-4357/aa970e
  * Singer et al. (2001) Singer, H. J., Heckman, G. R., & Hirman, J. W. 2001, Washington DC American Geophysical Union Geophysical Monograph Series, 125, 23, doi: 10.1029/GM125p0023
  * Sioulas et al. (2020) Sioulas, N., Isliker, H., Vlahos, L., Koumtzis, A., & Pisokas, T. 2020, MNRAS, 491, 3860, doi: 10.1093/mnras/stz3259
  * Voelk (1975) Voelk, H. J. 1975, Reviews of Geophysics and Space Physics, 13, 547, doi: 10.1029/RG013i004p00547
  * Xu (2021) Xu, S. 2021, submitted
  * Xu & Lazarian (2018) Xu, S., & Lazarian, A. 2018, ApJ, 868, 36, doi: 10.3847/1538-4357/aae840
  * Xu & Lazarian (2020) —. 2020, ApJ, 894, 63, doi: 10.3847/1538-4357/ab8465
  * Xu et al. (2015) Xu, S., Lazarian, A., & Yan, H. 2015, ApJ, 810, 44, doi: 10.1088/0004-637X/810/1/44
  * Xu & Yan (2013) Xu, S., & Yan, H. 2013, ApJ, 779, 140, doi: 10.1088/0004-637X/779/2/140
  * Xu et al. (2016) Xu, S., Yan, H., & Lazarian, A. 2016, ApJ, 826, 166, doi: 10.3847/0004-637X/826/2/166
  * Yan & Lazarian (2002) Yan, H., & Lazarian, A. 2002, Physical Review Letters, 89, B1102+, doi: 10.1103/PhysRevLett.89.281102
  * Yan & Lazarian (2003) —. 2003, ApJ, 592, L33, doi: 10.1086/377487
  * Yan & Lazarian (2004) —. 2004, ApJ, 614, 757, doi: 10.1086/423733
  * Yan & Lazarian (2008) —. 2008, ApJ, 673, 942, doi: 10.1086/524771