Magnetohydrodynamic-Particle-in-Cell Simulations of the Cosmic-Ray Streaming Instability: Linear Growth and Quasi-linear Evolution
Xue-Ning Bai Eve C. Ostriker Illya Plotnikov James M. Stone
















^1Institute for Advanced Study, Tsinghua University, Beijing 100084, China
^2Tsinghua Center for Astrophysics, Tsinghua University, Beijing 100084, China
^3Department of Astrophysical Sciences, Peyton Hall, Princeton
University, Princeton, NJ 08544
^4 Université de Toulouse, UPS-OMP, IRAP, 9 av. Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France 
xbai@tsinghua.edu.cn















The gyro-resonant cosmic-ray (CR) streaming instability is believed to play a crucial role
in CR transport, leading to growth of Alfvén waves at small scales that scatter CRs,
and impacts the interaction of CRs with the ISM on large scales.


However, extreme scale separation  (λ≪ pc), low cosmic ray number density
(n_ CR/n_ ISM∼ 10^-9), and weak CR anisotropy (∼ v_A/c) pose
strong challenges for proper numerical studies of this instability  on a microphysical level.
Employing the recently developed magnetohydrodynamic-particle-in-cell (MHD-PIC) method,
which has unique advantages to alleviate these issues, we conduct one-dimensional
simulations that quantitatively demonstrate the growth and saturation of the instability in the
parameter regime consistent with realistic CR streaming in the large-scale ISM.   
Our implementation of the δ f method dramatically reduces Poisson noise and enables
us to accurately capture wave growth over a broad spectrum, equally shared between left
and right handed Alfvén modes. We are also able to 
accurately follow the quasi-linear diffusion of  CRs subsequent to wave growth, which

is achieved by employing phase randomization across periodic boundaries. 
Full isotropization of the CRs in the wave frame requires pitch angles of most CRs to
efficiently cross 90^∘, and can be captured in simulations with
relatively high wave amplitude and/or high spatial resolution.
We attribute this crossing to non-linear wave-particle interaction (rather than mirror reflection) by
investigating individual CR trajectories. 
We anticipate our methodology will open up opportunities for future
investigations that incorporate additional physics.












§ INTRODUCTION


Cosmic-rays (CRs) are an essential constituent of the Galaxy. With energy density
comparable to or exceeding other components of the interstellar medium (ISM), CRs likely
play an important role ISM heating, chemistry, and  dynamics 
over timescales of ∼ 10^3 -10^8  yr, as well as galaxy formation and
evolution over Gyr timescales
<cit.>.  


The dominant CR population is protons, with ∼GeV particles representing
the peak of the CR energy distribution; 
the typical 
cosmic ray energy density of ∼ 1  eV cm^-3 corresponds to a number density 
n_ CR∼ 10^-9 cm^-3
<cit.>.

CRs are expected to be primarily produced in collisionless shocks, most likely from
young supernova remnants (e.g., ), via diffusive shock acceleration
(e.g., ).
They are then scattered across the Galaxy,
primarily following magnetic field lines in the turbulent ISM,
and eventually escaping after a few million years <cit.>.

Being a non-thermal high-energy particle population, the CRs 
only infrequently collide
with the ISM particles (this is greatest at low energies, causing ionization and
heating), but the CRs and bulk ISM gas indirectly interact as they both couple to
the interstellar magnetic field.
Through this indirect coupling, CRs can exert pressure forces on the ISM gas which
are dynamically important. Moreover, being a (trans-)relativistic plasma component,
CRs in the ISM are highly buoyant, with a natural tendency to escape from the Galaxy.
As a result, the CR component can possibly 
contribute to large-scale instabilities and dynamo activity <cit.>,
and to driving galactic winds <cit.>.

At microscopic level, the physics of the 
interaction between CRs and background thermal ISM gas
is extremely rich.

Considered as passive test particles, the CRs can
be scattered diffusively by background ISM turbulence at the gyro-radius scale (e.g., 
), which is important for understanding
CR transport. For typical ∼ GeV CR particles, the gyro-radius for an ISM magnetic field of a few
μ G is on the order of 10^12cm – i.e. a micro-pc.
Considering the huge dynamic range from the energy-containing 
scales of ISM turbulence to the relevant wavelength 
(a factor of 10^8 in length), as well as the anisotropy of the MHD cascade at sub-pc scales <cit.>,
it is considered unlikely for background turbulence to be the
dominant source of CR scattering for GeV particles (e.g., ). 
There are even some suggestions from spectral signatures that only higher-energy
(≳ 200  GeV)
particles are affected by background ISM turbulence <cit.>.

Anisotropy in the distribution of CRs can destabilize Alfvén waves.
Of particular interest
and relevance is the cosmic-ray streaming instability (CRSI,
).  Instability occurs when the bulk drift
speed of the CRs exceeds the Alfvén speed in the background plasma.

The instability drives the growth of Alfvén waves at the cost of free energy from CR anisotropy,

and in the absence of wave damping  
the resulting scattering of particles off the waves
will effectively reduce the CR drift speed to Alfvén speed. This process limits the bulk flow 
speed of the CRs, and simultaneously 
leads to energy and momentum exchange between the CRs and background gas.
As a consequence, CRSI plays a fundamental role in the interaction between the CRs and ISM gas,
and is the key element behind the picture of CR self-confinement as well as CR-driven winds 
(see, e.g.,  for reviews).  On much larger scales in galaxy clusters, the CRSI may also be important in regulating heating <cit.>.  Here, however, we shall focus on the parameter regime appropriate for cosmic ray interactions with the ISM.  

Studies of CR interactions with the ISM on large scales, whether analytic or numerical, treat the CRs as fluid-like
(e.g., ), generally adopting the assumption of streaming, diffusion, or some combination.
This is necessary given the huge scale separation between the macroscopic scale of
astrophysical systems and the tiny gyro-scale of the CRs, but this approach is subject to the
prescriptions adopted to treat the microphysical effects of CR diffusion and streaming.

In recent years, there has been a resurgence of interest in  CR-driven galactic winds.  While most numerical studies to date have adopted diffusion prescriptions 
<cit.>, there have also been studies that compare wind outcomes when streaming vs. diffusive prescriptions are adopted for the CRs
<cit.>.  
The dynamics of winds in which CRs stream  
at the Alfv'en speed have distinctive features, since the effective sound speed of gas increases as its density decreases, enabling steady-state CR-driven acceleration of cool gas to high velocities in typical galactic potentials <cit.>. 
Most recently, several studies of wind driving have considered even more sophisticated treatment of the CR fluid, implementing prescriptions for wave damping and for decoupling of CRs from cold ISM components <cit.>.   
While all of the above studies highlight the important role played by the CRs, it is also
evident that different microphysical prescriptions can yield dramatically different
outcomes for wind properties. Currently, however, microphysical coefficients that are adopted in simulations are motivated on phenomenological grounds or based on classical analytic theory of the CRSI, diffusion, and wave damping.   
Development of a modern theoretical foundation for CR fluid treatments, with coefficients calibrated from numerical studies that directly follow the microphysics, will clearly aid progress in this field.  


At microscopic level,
our current understanding of the CRSI is limited to quasi-linear theory.
We note that on the one hand, the instability can become non-linear when
CR streaming is sufficiently strong, and eventually transitions to the Bell
instability <cit.> in more extreme environments such as super-Alfvénic
shocks in supernova remnants. On the other hand, there are additional
uncertainties related to mechanisms that lead to wave damping (see discussions
in Section <ref>), many of which again require numerical studies. 


Numerical studies of the CRSI at microphysical level 
per force must involve
particle-in-cell (PIC) methods, and are inherently challenging 

because of the extreme ratios between the density of thermal and nonthermal particles and the scales relevant for the background thermal plasma compared to CR gyroradii.
Another difficulty is that the weak level of CR anisotropy (of order a few times v_A/c)
must be accurately represented, demanding a huge number of computational
particles (see Section <ref>). Both of these issues make it very
challenging for conventional particle-in-cell (PIC) methods to study the CRSI, 
especially in the regime relevant for the bulk of the ISM away from CR injection sites
(but see  and further discussion 
in Section <ref>).

Recently, we developed 
a new hybrid MHD-PIC method <cit.>, 
implemented within the Athena code package <cit.>.
By treating the background thermal plasma as a fluid described by magnetohydrodynamics
(MHD), 
we greatly alleviate the issue of scale separation between the thermal gas and CRs. 
Using this method, we 
are able to properly capture the basic properties
of the CRSI, particularly in astrophysically relevant/realistic parameter regimes.
For this purpose, 
in our first study we simplify by ignoring wave damping and concentrating on
one-dimensional simulations that provide the highest possible resolution. 
By employing several novel techniques, we
demonstrate that our method can accurately reproduce the linear growth of CSRI, and
subsequently  follow the quasi-linear diffusion of particles by scattering off of
waves generated by the CRSI. We anticipate this effort will pave the  way for future
studies that incorporate additional physics under a variety of environments.



This paper is organized as follows. We describe the basic theory of the CRSI
in Section <ref>. In Section <ref>, we present the techniques
and methodology employed in our CRSI simulations using the MHD-PIC approach.
Our simulation setup is described in Section <ref>. In Section <ref>,
we present simulation results from two fiducial runs. This is followed by a parameter
study in Section <ref>. Our results are discussed in a  broader context in
Section <ref>, before we conclude in Section <ref>.
Additional numerical aspects of our simulations are described in the Appendices.



§ BASIC THEORY


In this section, we briefly review the theory of the CRSI, which largely follows the
description in <cit.>. We focus on the growth and saturation of the
instability. While wave damping is another important ingredient for the CRSI, it will
be addressed in our follow-up work and is not considered here.



 §.§ Physical Mechanism


Without rigorous mathematical derivation, the basic physics of the CRSI can be
understood schematically as illustrated in Figure <ref>.
Consider a static background thermal plasma with constant density ρ_0 and
constant magnetic field B=B_0x̂ with B_0>0. The Alfvén speed is thus
v_A=B_0/√(4πρ_0). Coexisting with the background plasma is a population
of CRs. Let the initial CR velocity distribution be isotropic in a frame that drifts relative
to the background gas at velocity v_D, with the drift velocity along the background field
x̂. Without loss of generality, let v_D>0 (along the direction of magnetic field).
Let f_0( v) or f_0( p) denote the initial particle velocity or momentum
distribution function, monotonically decreasing with increasing p.


Note that individual CR particle speeds greatly exceed v_D, hence approximately half
of the particles travel forward/backward with respect to the gas.
For particles with mass m, charge q and Lorentz factor γ, the cyclotron
frequency Ω_c and gyro frequency Ω read

    Ω_c=qB_0/mc ,   Ω=Ω_c/γ .


The CR particles can experience gyro-resonance with circularly polarized Alfvén waves
when the pattern of particle gyration matches that of the waves, and the frequency of the
waves experienced by the particles matches their gyro-frequency. To satisfy the first
condition, forward traveling CRs can resonate with left-handed, forward-propagating
Alfvén waves, whereas backward traveling CRs can resonate with
right-handed, forward-propagating Alfvén waves (see Appendix <ref>
for our adopted conventions for left/right handedness). To satisfy the second condition,
the parallel velocity of a particle v_∥ must satisfy

    ω-kv_∥=±Ω≈ -kv_∥ ,

where ω and k are the frequency and wave vector of the Alfvén wave
(we only consider wave vector k along the magnetic field), and the +/- sign
corresponds to right/left polarized modes, respectively. Note that ω=kv_A,
while typically the CR velocity (and its projection) v_∥≫ v_A. We
can thus omit ω in the above equation. In other words, the CR
particles travel so rapidly that Alfvén waves simply appear as a static pattern.

Now consider CR particles in velocity space decomposed into components parallel
and perpendicular to the background magnetic field (v_∥, v_⊥), and
pick a representative particle located at point “O" shown in Figure <ref>.
To facilitate the illustration, we now assume particles are non-relativistic (but it can be
easily generalized to relativistic CRs), and hence f_0( v) is constant along a
circle centered on (v_D, 0) (the blue solid circle in Figure <ref>). Clearly,
f_0 is larger (smaller) inside (outside) the circle.
For linear Alfvén waves, their associated electric field vanishes in the wave frame.
Because magnetic fields do no work on particles, each particle trajectory in
velocity space must follow a circle centered in (v_A, 0) (black dashed
circle in Figure <ref>). Wave-particle interaction generally leads to a
quasi-linear diffusion (QLD) of particle pitch angle, towards removing the gradient of
f_0( v) across the dashed circle near point O. We see that when
v_D>v_A, for forward traveling CR particles near O (left panel), the net outcome
should be that more particles with higher v_∥ on the right hand
side of point O diffuse along the dashed circle towards smaller v_∥. 
In net, when v_D>v_A the diffusion (via resonant interactions) of particles irreversibly
gives momentum to the left-handed Alfvén waves, feeding their growth.
Similarly, for backward traveling CR particles (right panel), more particles on the right
side of O diffuse along the dashed circle to the left side, with the change in
v_∥ feeding the growth of right-handed Alfvén waves.
We can also easily see that if v_D<v_A, instead of wave excitation, the Alfvén
waves will be damped.

This discussion illustrates how the CRSI is driven by the anisotropy of the CRs.
At the same time, it suggests that CRSI growth is sustained by 
QLD of the CRs 
which produces secular 
and irreversible changes in the background plasma.
This means that for a numerical method to capture linear
growth of CRSI, it must properly capture QLD as well, which is an important fact to bear in
mind when we discuss numerical methods in Section <ref>.



 §.§ Linear Instability of Waves


We now proceed to a more quantitative description of the CRSI.
Its derivation follows from standard procedures, as detailed in <cit.>.
Here, we only summarize the necessary results needed for setting up and analyzing
the simulations.

For convenience, in the discussions below and also throughout this paper, we
assume that the mass of CR ions and background ions are identical (being m).
As will be seen in Section <ref>, the value of m has no
physical significance and can always be absorbed into a factor for the
charge-to-mass ratio. We thus drop it and express the momentum of and energy of individual CR
particles as

    p=γ v ,   E_p=γ c^2 ,

where v is a particle's
velocity and γ is its Lorentz factor. In other words, we effectively take
m=1. However, we do retain the symbol m in certain places for better clarity.

As defined earlier, we use f_0( p)=f_0(p) to denote the initial CR distribution
function in the frame moving at v_D,

where the distribution is isotropic, with total CR number density given by

n_ CR=∫ 4π p^2f_0(p)dp.


For Alfvén waves with perturbations 

proportional to exp[ i(kx-ω t)], the
CR-modified dispersion relation reads

    ω^2=k^2v_A^2∓n_ CR/n_iΩ_c(ω-kv_D)
    [(1-Q_1)± iQ_2] ,

where n_i is the density of background ions. The role of the CRs is reflected in
the two dimensionless factors

    Q_1(k)   ≡∫ dp4π p^2f_0(p)/n_ CR(p_ res/2plog|1+p/p_ res/1-p/p_ res|) ,
    
    Q_2(k)   ≡π/2∫_p_ res^∞dp4π p^2f_0(p)/n_ CR(p_ res/p) ,

where the dependence of Q_1 and Q_2 on k is reflected in the resonant
momentum 

    p_ res(k)=Ω_c/k ,

which is the minimum CR momentum to resonate with an Alfvén wave
with wave number k (at zero pitch angle).[The dispersion relation
shown here is equivalent to Equations (11)-(13) in <cit.> with slight
differences in sign conventions. The terms Q_1 and Q_2 arise from integrals
containing df/dp (see also Equation (69) in Chapter 12 of <cit.>); this
is the mathematical embodiment of the situation depicted in Figure <ref>,
in which a gradient in f leads to wave growth. We have performed integration
by parts to simplify the results, and df/dp no longer appears.]

Physically, the term associated with 1 in the (1-Q_1) factor is due to the
background CR current along B_0. The Q_1 term is from the
non-resonant response of CR particles to the waves. It is straightforward to
show that Q_1 approaches 1 in long-wavelength limit, and monotonically
decreases to 0 in short wavelength limit. Thus, the real part of the CR
contribution is bounded. The imaginary part, namely, the Q_2 term, results
from the resonant response of CR particles to the waves. Note that at a given
k, only particles with p>p_ res contribute to Q_2, as expected.

In the limit n_ CR/n_i≪1, the real part of the dispersion relation
(<ref>) is clearly dominated by the k^2v_A^2 term, and hence the
term proportional to (1-Q_1) can be neglected. In this case, the dispersion
relation largely gives normal Alfvén waves, which can slowly grow or damp
owing to the imaginary part of the dispersion relation (i.e., the Q_2 term
from resonant response). This is the general CRSI originally derived
<cit.>, where the forward propagating modes, both left and
right handed, grow at the same rates, given by

    Γ(k)=1/2n_ CR/n_iΩ_c(v_D/v_A-1)Q_2(k) ,

and hence growing modes are essentially linearly polarized.
Clearly, instability occurs only when v_D>v_A. Moreover, the maximum growth
rate is of the order (n_ CR/n_i)Ω_c at wavelengths that maximizes
Q_2 (to reach order unity), which typically correspond to the resonant scale with
lowest energy CRs. Backward propagating modes are damped.

More specifically, the initial distribution function is commonly considered to be
a power law with

    f_0(p)∝ p^-4+α ,   p_0≤ p≤ p_ max

and f_0(p)=0 otherwise. Fermi acceleration gives α=0 (e.g.,
). In typical ISM conditions, the Galactic CR
population corresponds to α≈0.7 for CR particle energy beyond
∼GeV, and n_ CR/n_i∼10^-9. The Q_2 factor is given by

    Q_2(k)=π/2(1+α/2+α)p_ resp_ min^-(2+α)-p_ max^-(2+α)/p_0^-(1+α)-p_ max^-(1+α) ,

where p_ min= Min[ Max(p_ res, p_0), p_ max].
The maximum growth rate is achieved at p_ res=p_0, or

    k=k_0≡Ω_c/p_0 .

This corresponds to the resonant wavelength for particles with p=p_0 (and zero pitch
angle). In the limit p_ max→∞, we have

    Q_2(k)=π2(1+α2+α)
    (kp_0Ω_c)^1+α ,    
    k<k_0 , 
    π2(1+α2+α)
    (Ω_ckp_0) ,    
    k>k_0 .




Taking into account the term from (1-Q_1) breaks the degeneracy between
left and right handed modes, especially when n_ CR/n_i increases.
In particular, it makes the right-handed modes grow faster and more dominant
towards short wavelength (see Figure <ref>), and eventually, it smoothly
transitions to the non-resonant Bell mode at k≳ k_0 <cit.>. 
This transition occurs when <cit.>

    n_ CR/n_i≳v_A^2/p_0v_D .


For simulations in this work, we fiducially adopt v_D/v_A=2, p_0/v_A=300,
α=0.5, and n_ CR/n_i=10^-4 to 10^-3. The resulting
dispersion relation is shown in Figure <ref>, which displays all features
discussed above. Note that  
for application to realistic large-scale ISM conditions, a very low value of low
n_ CR/n_i∼ 10^-9 is appropriate, such that the left/right handed modes are
degenerate.  However, in practice, the values of n_ CR/n_i we use in the
simulations have to be much higher 
in order to follow wave growth and particle
quasi-linear diffusion within a reasonable computational time. We see that
for our fiducial parameters, the transition to non-resonant modes
occurs when n_ CR/n_i≳10^-3.



 §.§ Quasi-Linear Diffusion of Particles


The growth of the CRSI creates a spectrum of Alfvén waves. The intensity of the
waves at wave number k is described by I(k), so that

    ∫ I(k)dk=⟨δ B^2/B_0^2⟩ ,

where angle bracket represents spatial averages.
In Appendix <ref>, we describe the procedures to decompose parallel
propagating Alfvén waves into different modes (left/right handed, forward/backward
propagation), from which we obtain I(k) for each mode.

In a spatially uniform media (as in our simulations), quasilinear diffusion (QLD) leads to
pitch angle scattering in the frame of the Alfvén waves, with pitch angle θ
more commonly replaced by μ≡cosθ. We use subscript `_w' to
denote quantities measured in this frame. The evolution of f_w, the CR distribution
function in the wave frame in a homogeneous system, is described by
(e.g., )

    f_w/ t=
    /μ_w[1-μ_w^2/2ν(μ_w) f_w/μ_w] ,

where




    ν(μ_w)=πΩ k_ resI(k_ res) ,

and

    k_ res≡Ω_c/(p_wμ_w)

is the wave number of the resonant waves.[In calculating ν from  Equation 
(<ref>), we use I(|k_ res|) and choose the intensity for left/right handed
waves for positive/negative μ_w. Also note our normalization of I(k) only accounts
for the magnetic energy spectrum (will be doubled for total energy spectrum). Because
of these normalization choices, our expression for ν does not have a factor 1/4.]

Our simulations are conducted in the frame where the initial CR distribution is isotropic,
whose relative velocity to the wave frame is Δ v≡ v_D-v_A. Using subscripts
`_∥' and `_⊥' to denote components parallel and perpendicular to
background field B_0, and to first order of Δ v≪ c, we have for particle
momenta

    p_w,⊥=p_⊥ ,   p_w,∥=p_∥+γΔ v .

With conservation of phase-space volume, it follows that the distribution functions
f=f_w.



The conversion relations between the two frames read



    p_w=p(1+μΔ v/v) ,  μ_w=μ+(1-μ^2)Δ v/v ,

where v∼ c is the velocity corresponding to momentum p.
In the quasi-linear regime, the distribution function only deviates slightly from f_0, and
hence we can write f(p,μ)=f_0(p)+δ f(p,μ). Transforming to the wave frame,
and to the leading order of δ f/f_0 and Δ v/c, we obtain

    f_w(p_w,μ_w)=f(p,μ)
    ≈   f_0(p_w)+δ f(p_w,μ_w)-μ_w f_0/ln pΔ v/v_w .


In the wave frame, we have

    f_w/μ_w≈δ f/μ_w- f_0/ln pΔ v/v_w .

In the fully saturated state, we expect f_w/μ_w=0, and hence the pitch angle
distribution should satisfy

    δ f_ sat/μ_w≈ f_0/ln pΔ v/v 
    
    
     ⇒ δ f_ sat≈μ_w f_0/ln pΔ v/v .




The amplitude of the Alfvén waves in the saturated state can also be estimated.
The 
net momentum density acquired by the CRs after saturation is reached is

    Δ P_ CR=∫μ p(δ f_ sat) d^3p
    =4π/3∫ f_0/ln pΔ v/vp^3dp .

Note that with ∂ f_0/∂ p <0, this  implies that if initially Δ v=v_D -v_A >0,  in saturation the CR distribution must have acquired a momentum deficit, Δ P_ CR <0.
In the limit that all CRs are relativistic, this yields

    Δ P_ CR
    ≈ - 4/3n_ CRΔ v/c⟨ p⟩ ,

where ⟨ p⟩≡∫ f_0pd^3p/∫ f_0d^3p. For the non-relativistic CR limit,
the momentum reduction is
Δ P_ CR≈ - mn_ CRΔ v.


The deficit in CR momentum must have been transferred to 
forward-propagating Alfvén waves. Note that there is an equipartition
in kinetic and magnetic energies in Alfvén waves, with effective momentum density in waves
(the Poynting flux divided by v_A^2; cf. )
given by

    P_ wave=ρ v_A⟨δ B^2/B_0^2⟩ .

Following the discussion in Kulsrud <cit.>, this  momentum density represents the growth of Maxwell and Reynolds stresses in waves needed to transfer the CR momentum to the bulk gas,
and is exhibited as a small fractional change in the momentum of the bulk gas.

Equating |Δ P_ CR| with P_ wave, we can obtain the expected saturation level for magnetic fluctuations.
For ultra-relativistic and non-relativistic cases, we simply have

    δ B^2/B_0^2≈4/3⟨ p⟩/mcn_ CR/n_iΔ v/v_A ,   or  ≈n_ CR/n_iΔ v/v_A .


The real situation (as we consider) is of course somewhere in between 
nonrelativistic and ultrarelativistic, but the dependence of wave amplitude on the basic problem parameters is clear.
These relations show that the saturation level of the instability 

is limited by the initial free energy in CR streaming. The above 

upper limit for magnetic field amplitudes would be reduced in the presence of damping.

The timescale for the CRs to relax 
to a state of isotropy in the wave frame 

is characterized by the scattering frequency
ν, which depends on both the pitch angle and the wave amplitudes. This means that
scattering across different pitch angle ranges can take different amounts of time. Still, one may
define a characteristic scattering frequency for the bulk CR population as

    ν_ QLD∼ν(p,μ)∼Ωδ B^2/B_0^2 .

Note that wave growth is eventually quenched as the gradient in distribution function is smeared out through QLD. This suggests that another way to estimate the wave saturation amplitude is by equating the above ν_ QLD with Γ from Equation (<ref>), which is indeed similar (to order-of-magnitude) to the more rigorous estimate discussed above.



 §.§ Reflection Across 90^∘ Pitch Angle


In Equation (<ref>), when μ_m→0 (i.e., pitch angle approaches
90^∘), k_ res→∞, i.e., the resonant wavelength becomes
infinitely small. There is little energy in the waves generated from the CRSI towards
the shortest wavelength, leading to a bottleneck for particles to be scattered across
this 90^∘ pitch angle. This is a well-known problem, and without crossing the
90^∘ pitch angle, the CRSI would saturate prematurely far from fully utilizing the
free-energy in the system.

This problem has been discussed extensively in the literature. Within quasi-linear
theory, relaxing the magnetostatic approximation [Equation (<ref>)] allows for
additional resonant interactions when both forward and backward propagating waves
are present <cit.>, covering a pitch angle of Δθ∼ 2v_A/c
around 90^∘.

Magnetic mirroring is more commonly invoked in overcoming the 90^∘ problem
(e.g., ). With magnetic moment M≡ p_⊥^2/2B being
an adiabatic invariant, a particle will experience a backward mirror force when its guiding
center travels along a positive gradient in total field strength. For particle pitch angle
sufficiently close to 90^∘, this mirror force can lead to reflection, thus directly jumping
over the 90^∘ barrier. With total field strength
B=√(B_0^2+δ B^2)≈ B_0(1+δ B^2/2B_0^2), the critical pitch angle
cosine below which mirror reflection can occur is thus given by

    μ_ mir≈1/√(2)δ B/B_0 .


One thus expects mirror reflection to take over when QLD manages to scatter particles
to pitch angle μ∼μ_ mir.

Note that mirror reflection requires slow changes in field strength over many particle
gyro orbits. 

Given that a broad range of wave spectrum is excited by the CRSI, including those
at small scales (≲ k_ res^-1) for particles of a given momentum, we can also imagine the opposite situation,
where changes in field strength occurs more abruptly. For instance, an abrupt change of
perpendicular field by δ B_dis would change the direction of particle motion by
δ B_dis/B_0, and in the mean time, the change the particle pitch angle is also on
the order of ∼δ B_dis/B_0. Depending on the gyro-phase, particles have a good
chance to be reflected when μ≲δ B_dis/B_0.

The (extreme) scenario outlined above falls into the category of non-linear wave-particle
interaction: considering particle orbits in perturbed fields (rather than unperturbed fields as
in quasi-linear theory) gives resonance broadening 
(e.g., ), and alleviates the 90^∘ problem.
It is more effective towards larger wave amplitudes, and eventually enters the regime of
strong MHD turbulence <cit.>.
As we will show, our analysis suggests that
this effect is primarily responsible for overcoming the 90^∘ barrier; 
for more discussion see Section <ref>.




§ NUMERICAL METHOD



There are several challenges that confront  numerical study of the CRSI.
First, there is substantial physical scale separation, which is exhibited as follows.


  * The separation between microscopic scales of the background plasma, including
the electron and ion skin depths (c/ω_pe or c/ω_pi that full-PIC or
hybrid-PIC simulations must resolve, where ω_pe and ω_pi are the
electron and ion plasma frequencies), and the resonant scale of the CRSI (∼
gyro-radii of lowest-energy, trans-relativistic CRs). 
The ratio of the two scales is at least of the order c/v_A, which is prohibitively large for
conventional PIC simulations if one were to approach realistic conditions.


  * For quasi-linear diffusion of particles, individual CR particles must encounter a sufficient
number of independent wave packets, so as to experience random-walk behavior in pitch angle evolution.


  * The issue with particles scattering across the 90^∘ pitch angle. Particles with pitch
angle near 90^∘ have resonant wavelength much shorter than their gyro-radii,
implying that the most unstable wavelength of the CRSI must be very finely resolved.



The MHD-PIC approach uniquely alleviates the first scale separation issue by bypassing the microscopic plasma
scales (ion skin depth) of the background thermal plasma. However, the remaining two issues
remain challenging. Clearly, very high resolution is needed to capture scattering across
90^∘ pitch angle. In the mean time, when employing periodic boundary conditions, a
sufficiently long simulation domain that cover many resonant scales is essential so that particles
do not experience the same wave packets upon re-entering the domain. This issue will be
discussed further in Section <ref> and Appendix <ref>.

In addition to the above challenges of physical scale separation, another challenge in modeling the regime of interest for the large-scale ISM is the huge ratio between the number density of CR particles and number density of background ISM particles, n_ CR/n_i ∼ 10^-9. Cosmic ray and ISM pressures, however, are generally within an order of magnitude of each other.
For modeling the large-scale ISM with its extremely small n_ CR/n_i, the MHD-PIC method is advantageous in comparison to full- or hybrid-PIC in that there there is no need to represent the background plasma with individual particles.

The third major challenge is that capturing the resonant condition still requires very large number of CR particles. For a
given resonant mode, particles of a given momentum can (if capable) contribute only at
a specific pitch angle. Correspondingly, only a tiny fraction of particles can potentially
be in resonance with a given MHD wave mode. This issue may be seen from another viewpoint.
The level of anisotropy in the CR distribution function is of the order v_D/c≪1. Such weak
level of anisotropy must be accurately represented by the angular distribution of particles,
again demanding for a huge number of particles in each cell.

Even though the MHD-PIC method is able to concentrate all particles in the CRs, 
properly capturing CRSI is still challenging when v_D/c≪1.
To overcome this issue, we employ the δ f method, which dramatically improves the
signal-to-noise and allows us to successfully simulate the CRSI with substantially reduced
number of particles (see Section <ref> and Appendix <ref>).
Below, we describe our methodology in detail.



 §.§ Formulation


We solve CR-modified MHD equations in conservative form[These equations are
simplified from those derived in <cit.> in that terms related to the CR-modified
Hall effect are dropped. As discussed there, these terms are important only when
(n_ CRv_D)/(nv_A)≳ 1,  corresponding to very strong CR streaming.
In our simulations, this ratio is ≲10^-3, and in reality the ratio is several orders
of magnitude lower.]

    ρ/ t+∇·(ρv_g)=0 ,


    (ρv_g)/ t+∇·(   ρ v_g^T v_g
    
    - B^T B+P^* I)
    
    =   -(q_ CRn_ CRE+
     J_ CR/c× B) ,


    E/ t+∇·[(E+P^*) v_g
    
    -( B· v_g) B]
    =- J_ CR·E ,

where P^*≡ P_g+B^2/2, E≡- v_g× B/c is the
electric field, I is the identity tensor, and the total energy density of the gas is defined as

    E=P_g/Γ-1+1/2ρ v_g^2+1/2B^2 .

In the above, ρ, v_g, P_g are gas density, velocity, and pressure,
and Γ is the adiabatic index. 
Note we have adopted the units where magnetic permeability is unity so that factors
of (4π)^-1/2 that would otherwise appear with the magnetic field are eliminated.


The CR number density n_ CR and current
density J_ CR are defined as

    n_ CR(t,  x)   =∫ f(t,  x,  p)d^3 p ,
     J_ CR(t,  x)   =q_ CR∫ vf(t,  x,  p)d^3 p ,

where q_ CR is individual CR particle charge, and p, v are
the momentum and velocity of individual CR particles, with f(t,  x,  p) is the local
CR momentum distribution function.

While the MHD formulation is non-relativistic, the CR particles can be relativistic.
We define an artificial speed of light ℂ for the CR particles, and the overall
formulation is consistent as long as ℂ≫ v_A (typical MHD velocities).
For an individual particle j, we

have p_j=γ_j v_j, with
Lorentz factor given by





    γ_j=√(ℂ^2+p_j^2)/ℂ
    =ℂ/√(ℂ^2-v_j^2) .


The particle equation of motion reads

    dp_j/d t
    =(q/mc)_j(cE
    + v_j× B) .

Note that as mentioned earlier, we have dropped the individual particle mass in the
definition of particle momentum. This mass is absorbed to the factor q/(mc), representing
particle charge-to-mass ratio. Also note that the physical speed of
light c has no significance 
in all equations where it appears: inspection of Equations 
(<ref>), (<ref>), (<ref>), and (<ref>) shows that c
appears only in the combinations J_ CR/c,  c  E, and
q/mc, and for numerical purposes can be absorbed into the parameter (q/mc).


In our simulations, ideal MHD equations for the background thermal plasma are solved
using the Athena MHD code <cit.>, which is a higher-order Godunov code
with constrained transport to enforce the divergence-free condition of the magnetic
field. The corner transport upwind (CTU, )
method is adopted for time integration. We use the Roe Riemann solver <cit.>
and third-order reconstruction in characteristic variables.
CR particles are implemented as Lagrangian particles, and the coupling between CRs
and the background plasma is handled by adding source terms 
in the gas momentum and energy updates, as indicated on the right hand sides of
Equations (<ref>), (<ref>). The CR particle equation of motion
(<ref>) is solved by the standard Boris integrator <cit.>.
A standard triangular-shaped cloud (TSC) scheme <cit.> is
used for interpolating grid quantities to particle locations, and for depositing particle
quantities back to the grid (e.g., to evaluate n_ CR and J_ CR
in Equation (<ref>)).
Details of the CR implementation are described in <cit.>.



 §.§ The δf Method


The original CR implementation in <cit.> interprets individual particles 
as representing the full distribution function f( x,  p), and evaluates
physical quantities such as n_ CR, J_ CR directly according
to Equation (<ref>). We refer to this approach as the “full-f method”, which can be
subject to large Poisson noise.

In the case where the distribution function is close to some equilibrium distribution
f_0( x,  p), one can use the fact 
that f_0 is already known
analytically, and employ individual particles as Lagrangian markers taken to represent the
difference, δ f, between f_0 and the full distribution function f. This is known
as the δ f method (e.g.,
). 
The basis of the δ f method is the Liouville theorem, which requires that the full distribution function f be constant along
particle trajectories in phase space (i.e., characteristics). 

To implement the δ f method, we first record the initial value of f at t=0 for
all particles (which is essentially f_0). Then at every time t, we assign a weight w_j to
each particle j given by

    w_j≡δ f(t,  x_j(t),  p_j(t))/f(t,  x_j(t),  p_j(t))
    =1-f_0( x_j(t),  p_j(t))/f(0,  x_j(0),  p_j(0)) .

Physical quantities such as the CR number density and current density are obtained by

    n_ CR(t,  x)   =n_ CR,0+∫δ f(t,  x,  p)d^3 p
       ≃ n_ CR,0+∑_j=1^N_pw_jS( x- x_j) ,


    J_ CR(t,  x)
    
       = J_ CR,0+q_ CR∫ vδ f(t,  x,  p)d^3 p
       ≃ J_ CR,0+q_ CR∑_j=1^N_pw_j v_jS( x- x_j) ,

where S( r) is the shaping function used in particle interpolation (i.e., the TSC
scheme), and n_ CR, 0 and J_ CR, 0 are obtained analytically
from f_0. By contrast, full-f method corresponds to setting w_j=1, n_ CR, 0=0
and J_ CR, 0=0.
In this way, the δ f method dramatically reduces the Poisson noise of
the background particle distribution, allowing the signal from the δ f part to be
substantially boosted.


In practice, we have found that the δ f method is essential. Without
employing it (i.e., using the full-f method), and using our fiducial simulation parameters,
we barely observe the development of the CRSI even using ∼10^4 particles per cell,
and the system is almost entirely dominated by Poisson noise (see  Appendix
<ref> for more information).


We note that with the δ f method, exact conservation of the total
(gas and CR) momentum and energy, achieved in the full-f method, is
lost. This is inevitable, since the formulation of the δ f method is intrinsically
non-conservative. However, the benefit from low noise is of overwhelming importance
for the CRSI problem, and in practice we find that the error in total momentum and energy
is negligible throughout all our simulations.










 §.§ The κ distribution


While the CRSI is commonly analyzed with f_0 being a simple truncated power law
(<ref>), this is incompatible with the δ f method. The reason is that
the δ f method
requires f_0 to be finite at all p since it serves as the normalization factor in the
weighting function, Equation (<ref>).
To avoid divergence, we modify f_0 to be a κ distribution (e.g., see
 and references therein for motivation)

    f_0(p)=n_ CR/(πκ p_0^2)^3/2Γ(κ+1)/Γ(κ-1/2)[1+1/κ(p/p_0)^2]^-(κ+1) ,



which has the property that f(p) is a constant when p≪ p_0,
and f(p)∝ p^-2(κ+1) as p≫ p_0.

The dispersion relation remains the same except with factors of Q_1(k) and Q_2(k) in
Equation (<ref>) replaced by

    Q_1
    =   Γ(κ+1)/s_0^3(πκ)^3/2Γ(κ-1/2)∫_0^∞2π s· ds/[1+(s/s_0)^2/κ]^κ+1log|1+s/1-s| ,


    Q_2
    
    =√(π)/κ^3/2Γ(κ+1)/Γ(κ-1/2)1/s_0[1+1/(κ s_0^2)]^κ ,

where the dependence on k is expressed in 
s_0≡ p_0/p_ res=kp_0/mΩ_c.
The Q_1 integral can be evaluated numerically, and in the limit n_ CR/n_i≪1,
the growth rate can be obtained analytically by substituting (<ref>) into
(<ref>).

In Figure <ref> (top panel), we show the linear growth rate in the same way as in Figure
<ref>, but for f_0 being a κ distribution. 
We see that in the long wavelength limit (p_ res≫ p_0 or s_0≪1), we have
Q_2∝ s_0^2κ-1. In the short wavelength limit, we have
Q_2∝ s_0^-1. The results are consistent with a power law distribution
function in the same limits, with α=2(κ+1). The growth rate varies smoothly
in between the two limits, with peak growth rate only slightly reduced compared to the
truncated power-law case. It is this analytical growth rate that we aim to test with further study.

We further show in the bottom panel of Figure <ref> the analytical growth rate
for different drift velocities at fixed n_ CR/n_i=10^-4. Increasing the
drift velocity has a similar effect as increasing n_ CR/n_i, which increases the overall
growth rate (roughly linearly). It also makes the right-handed mode stand out, which becomes
more notable for v_D≳10v_A. In this paper, we restrict our explorations to v_D≲8v_A.



 §.§ Simulation Box Size and Phase Randomization



Even equipped with the δ f methods and a κ distribution function, 
it is difficult to accurately recover the linear growth of the CRSI with a numerical simulation.
Over a large suite of tests with wide range of numerical setups, we have found that sustained growth of the
CRSI at the rate predicted by theory 

is only achieved when the simulation box size is extremely long.
Otherwise, the desired growth rates are found only at very early stages, and they
become erroneous shortly afterwards (we demonstrate this in Appendix <ref>).

Further adjusting the numerical speed of light ℂ reveals that 
in order to properly recover growth rates,
the required simulation box length L should be such that the time for the
CR particles to traverse the simulation domain, t_ cross=L/(μℂ), exceeds the growth timescale for CRSI, 1/Γ_ max.  
 Given that
Γ_ max∼ (n_ CR/n_i)Ω_c, and using Equation <ref>, 
this may be further translated to

    L>L_ min≡ℂμ/Γ_ max∼n_i/n_ CRk_ res^-1/γ.



This phenomenon is related to the validity of the random phase approximation for QLD
discussed at the beginning of Section <ref>, which feeds instability growth.
For CRSI, CR particles
should constantly experience different wave packets throughout the growth process.
With wave growth, we may regard the wave packets to be altered over the timescale of
∼Γ_ max^-1. With the fastest particles traveling at ℂ along
the background field, a simulation box longer than L_ min guarantees all particles
experience “new" wave packets after traversing the entire box and reenter from the
other side.

However, for reasonable simulation values of n_i/n_CR, the required L_ min could easily reach a few
hundreds or thousands times the resonant/most-unstable wavelengths 
(using n_i/n_CR∼ 10^9 as in the bulk ISM, the number could be billions).
The condition of L≳ L_ min thus would impose a prohibitive requirement for
simulating the CRSI, even in 1D.

We resolve this issue by realizing that the effect of long simulation box can be
equivalently achieved by directly randomizing particle phases when a particle crosses
the periodic boundary: both approaches allow particles to effectively see different wave
packets at all times to conform to the random phase approximation.
In our simulations, upon boundary crossing, we fix a particle's total velocity, as well as
its parallel velocity along background magnetic field, whereas we randomize its gyro-phase.
By this simple approach, we completely eliminate the requirement of L>L_ min.
Also, under the δ f method, the phase randomization procedure does not
introduce more errors to momentum/energy conservation.



§ SIMULATION SETUP AND DIAGNOSTICS




We set up the simulations in the frame of the “unperturbed” CRs, where their background
distribution function is isotropic f_0( p)=f_0(p). This is also the f_0 for the
δ f method, with background CR current J_ CR, 0=0.
We set background density to ρ_0=1 and pressure P_0=0.6, with an ideal gas
equation of state with adiabatic index Γ=5/3 (so that sound speed
c_s=√(Γ P_0/ρ_0)=1). While thermodynamics has no effect to the
system, we adopt adiabatic rather than isothermal equation of state to better characterize
energy conservation of the system (see Section <ref>).
In 1D (along the x̂ axis), the background magnetic field is set to
B_0=1 along x̂, so that the Alfvén speed v_A=B_0/√(ρ_0)=1.

For all simulations, we initialize the background gas to a velocity v_0=-v_Dx̂.  
Since we are working in the frame where the CRs are initially isotropic, this means that the
CR distribution streams to the right (along the direction of the mean magnetic field) with
speed v_D relative to the gas. 

On top of the background state, we initialize the system with a spectrum of Alfvén
waves propagating along background B_0. For each wavenumber |k|,
there are four modes corresponding to forward and backward propagation, and left
and right polarization. The wavenumber coverage ranges from k=2π/L_x to
k=2π/(2Δ x) (except that the longest wave is initiated with 0 amplitude),
where L_x=N_xΔ x is the simulation domain size, and
N_x and Δ x are the number of grid cells and cell size in x̂.
Correspondingly, a total of 4(N_x/2-1) modes are initialized. We choose the wave
amplitudes to be such that

    I(k)=A^2/k ,

so that wave energy is equally distributed in logarithmic k-space, and the total
wave energy integrated over one dex in k is a fraction 2A^2ln10 of the
background field energy 
(cf. Equation <ref>;  the factor of two accounts for velocity perturbations).
Each mode is initialized with a random phase.


We choose code units such that individual CR particles  have cyclotron frequency
Ω_c=1 in the background field B_0. This corresponds to the factor q/mc=1 in
Equation (<ref>). Thus, our code unit of length is d_i = v_A/Ω_c=1.
We note that  while irrelevant to MHD, d_i is equivalent to the ion skin depth
c/ω_pi for the background plasma. Our MHD-PIC framework typically applies to
scales greater than d_i. 


The CR particles are injected with f_0(p) being a κ distribution.

We choose κ=1.25 throughout this work, corresponding to f_0(p)∝ p^-4.5,
at large p or α=0.5 in Equation (<ref>).
In practice, we divide the momentum space into 8 bins ranging from 0.01p_0 to
100p_0, with half a dex per bin. We inject equal number (N_p) of CR particles
per cell per bin so as to guarantee sufficient number of particles in each bin
(otherwise there would be too few particles towards lower and higher energies).
Within each bin in each cell, particles are sampled according to f_0(p) in that
momentum range, and are symmetrically distributed in pitch angle and phase.
More specifically, once a set of random momentum p, pitch
angle θ and gyro-phase ϕ is generated, four particles are injected at
cell center whose momenta are
(pcosθ, psinθcosϕ, psinθsinϕ) and its permutations with
θ→π-θ and ϕ→π+ϕ.

We note that particles with p<0.01p_0 and p>100p_0 are not included in the
simulation whereas we still use a background distribution function f_0 that covers
the entire momentum space. This will introduce some inconsistencies. However,
since these momenta are sufficiently far from the peak of f_0, waves that
they drive are typically beyond what can be accommodated in our simulation box
(and moreover the simulation time is typically not long enough for these modes to grow).



 §.§ Choice of Parameters


The interstellar medium in our Galaxy spans a wide range of conditions in various
phases. In a averaged sense, n_ CR/n_i is on the order of ∼10^-9,
and v_A/c∼10^-5-10^-4, and the bulk of the energy  is in GeV CR particles that are mildly
relativistic. While the parameters that we adopt are by no means fully realistic 
(this would be unachievable), we aim to
achieve as much scale separation as possible to mimic conditions in the large-scale ISM.

We choose the artificial speed of light ℂ/v_A=300, and
choose the characteristic momentum 
in the CR distribution to be 
p_0/v_A=300=ℂ/v_A.  
With this choice, Equation (<ref>) yields γ_0 = √(2) so as to mimic the fact that the dominant CR population
in the ISM is trans-relativistic ∼GeV particles.

The fact that ℂ=300 is larger than v_A=1 by more than two orders of
magnitude allows for sufficient separation, making the transformation between
different reference frames well in the non-relativistic limit (e.g., a drift velocity
as large as ∼10v_A only amounts to about 3% of ℂ).
With our choices, the wavelength for most unstable mode is given by
λ_m≈ 2π/k_0=2π p_0/Ω_c≈1885 in code units.

Based on Equation <ref> for the growth rate of CRSI, the two main physical parameters are the CR number density ratio n_ CR/n_i,
and the CR drift velocity v_D/v_A. As already seen in Figure <ref>, given
our choices of ℂ and p_0, the density ratio n_ CR/n_i has to
be well below 10^-3 to make left and right handed modes grow at about equal
rates (as in standard CRSI). We thus choose n_ CR/n_i=10^-4 as fiducial,
and vary it within [10^-5, 10^-3]. Note that further smaller values would make the
waves grow too little within a reasonable computational time.
We take v_D=2v_A for our fiducial model, and also consider higher drift velocities up to
8v_A.[In terms of energy density, we have E_ CR/E_B∼45 under the fiducial
parameters,
which is much larger than the order-unity value that applies through most of the ISM.
We note, however, that while the total energy density of CRs matters for large-scale ISM
dynamics, 

it is not directly relevant for the CRSI because the free energy of the CRSI comes from
CR anisotropy (i.e., streaming) rather than total CR kinetic energy. 
Nonetheless, with  E_ CR/E_B ∼ (n_ CR/n_i)(ℂ/v_A)^2, our run
with n_ CR/n_i=10^-5 better approaches the realistic limit. One can also
reduce ℂ/v_A to 100 or smaller to reduce this ratio, although
this would limit the range of v_D that can be more reliably covered by simulations.]

The remainder are numerical parameters. Fiducially, we choose simulation box size
L/d_i=9.6×10^4, resolved by 9600 cells (Δ x/d_i=10).

The fiducial box 
is about ∼50 times the most unstable wavelength, which is more than

sufficient for particles with p∼ p_0, but it better accommodates more energetic
particles with p≳10p_0. Note that without implementing phase randomization,
from Equation (<ref>) the box size L_ min would have to be 
∼ 300 times longer than our fiducial choice.
We also conduct a simulation with shorter box size for comparison.

The level of noise in the simulation depends on the number of particles. Fiducially, we
choose N_p=256 particles per momentum bin per cell (in total 2048 per cell),
but we also consider a case with fewer particles.  The initial wave amplitude is fiducially
chosen to be A=10^-4.
Note that A should not be too small, otherwise the initial evolution would be noise
dominated. It should neither be too large, which would trigger artificial QLD of particles
from the beginning. A more detailed discussion about noise in our simulations is
given in Appendix <ref>.

In Table <ref>, we list all simulation runs presented in this paper. We will
primarily focus on fiducial run Fid and M3 (higher n_ CR/n_i). For other simulations,
we generally vary just one parameter from run Fid at a time, and examine the role of
both physical and numerical parameters. Note that in the simulations, almost all
computational cost is spent on the CR parts given the large number of particles per
cell. For our fiducial run, a typical timestep is 6×10^-2Ω_c^-1, and
the entire run to t=10^6Ω_c^-1 costs about 4×10^4 CPU hours on
intel Xeon E5-2690 v4 CPUs.







 §.§ Diagnostics


We first measure the linear growth rate of the CRSI. Over constant time intervals,
we decompose the perpendicular magnetic field and velocity profiles into four
Alfvén modes and obtain I(k, t) for each mode (see Appendix <ref>).
A linear fit is then performed to ln[I(k,t)], with the slope being 2Γ(k).
This is done within early phases of evolution before QLD starts to modify the CR distribution
function. 

Another diagnostic is to directly measure the CR distribution function.
Thanks to the δ f method, this again can be done with great precision.
In the simulation frame, we simply compute f(p,μ)=f_0(p)+δ f(p,μ),
where f_0 is analytic and δ f can be obtained in the same way binning particles
onto a p-μ grid using the weighting scheme in Equation (<ref>). 

We can further combine the above information aiming to test QLD directly.
This involves transforming f(p,μ) into the wave frame using Equation (<ref>),
and then solving Equation (<ref>), with I(k) feeding into the calculation.
The results can be compared with the measured CR distribution function.

Finally, we have tracked a (small) subset of particles in certain runs over certain
time intervals; we combine this information  with the field data to investigate how particles overcome
the 90^∘ barrier.



§ SIMULATION RESULTS: RUN FID AND M3


In Figure <ref>, we show the evolutionary history of total wave energy
for the first six runs listed in Table <ref>. 
In this section, we first focus on the two main simulation runs: Fid and M3,
and will also briefly discuss run M5.
With wave amplitude A=10^-4, total wave energy integrated over the wave
spectrum in k- space, and summed over all four modes, is about
3.2×10^-7 (compared to 0.5 for background magnetic field energy
density). Wave energy slowly decreases at first, which is due to 
numerical damping of high-k modes at grid scale, as well as damping of
backward propagating modes. Shortly afterwards, exponential growth of
the CRSI 

becomes clear,
leading to rapid increase of wave energy. This
continues to time t∼10^5 for run Fid and t∼10^4 for run M3,
after which wave growth slows down as QLD starts to substantially modify the particle
distribution function. Subsequent growth is slower, and is characterized by a
combination of saturation of the fastest growing modes excited by low-energy
particles, and the linear growth of longer wavelength modes resonant with
(less abundant) higher-energy particles. Eventually, for run M3, the wave energy
plateaus, marking system saturation. The wave energy density in the
saturated state for run M3 is comparable to the estimates in Equation
(<ref>). However, for run Fid, in which the saturated energy
density should simply be one order of magnitude smaller, this is not yet reached.

Clearly, the wave energy evolution is much better discussed on the
k-by-k basis, and should be accompanied by the evolution of particle
distribution function.These will be addressed in the following subsections,
for individual evolutionary phases.









 §.§ Linear Wave Growth


In Figure <ref>, we show the Alfvén wave spectrum (in wave
intensity kI(k) for the four modes (forward/backward propagation, left/right
polarization) at some early times during the linear phases for the two runs.
Initial wave intensity of all modes has kI(k)=10^-8 being constant, with
δ v/v_A∼δ B/B_0 at the level of a few times 10^-4.
Modes with k≳ a few times 10^-2 are numerically damped.
With very low numerical noise, linear growth of all other forward-propagating
modes is observed from the beginning, whereas backward propagating modes
always damp. The fastest growing mode, resonant with particles near p=p_0,
is found to be near k=k_m≈3.3×10^-3=Ω_c/p_0 as expected.

More quantitatively, we measure the linear growth rate of forward-propagating modes
on the k-by-k bases. This is done within t=3×10^4Ω_c^-1 for
run Fid and t=10^4Ω_c^-1 for run M3. We also consider run M5, and
measure its growth rate within t=5×10^5Ω_c^-1. The results are
shown in Figure <ref>, and are compared with analytical growth rates.


Overall, even though the modes are very densely populated, the numerical growth rate
matches analytical results remarkably well. For runs Fid and M5, both left and right
handed modes grow at about the same rate, whereas for run M3, right-handed
mode grows slightly faster, all being consistent with theoretical expectations.
We stress that phase randomization is essential in achieving these results,
as is further discussed in Appendix <ref>.

The growth of the CRSI is better captured towards low k, whereas for high k,
instability growth competes with wave damping by numerical dissipation.
As a result, the measured growth rates cut off at some k∼ k_ cut, which
decreases from run M3 to run M5 as the instability growth rate decreases by about two
orders of magnitude. In fact, because of the small CRSI growth rate, k_ cut
corresponds to about 20, 30, and 45 grid cells for runs M3, Fid, and M5.
This again indicates much higher numerical resolution is needed to study the 
CRSI 

compared with typical pure MHD instabilities.

Note that for run M3, although the difference in growth rate between left and
right handed modes is small, the resulting difference in wave amplitudes can
be significant after a few e-folding times. This is already evident in the right
bottom panel of Figure <ref>: within one e-folding time,
the intensity of the right handed mode is already about twice that of the left
handed mode. This will affect subsequent QLD of particles to be discussed later.

As a side note, we have also conducted simulations with drift speed v_D
between 0 and v_A. We find that all wave modes are damped, and
damping is the fastest for at k∼ k_m. In the case of v_D=0, 
forward/backward propagating modes damp at the same rate, whereas for
0<v_D<v_A, damping of forward propagating modes are slower than that
for backward propagating modes. They all agree well with the dispersion relations.



 §.§ Quasi-Linear Evolution for run Fid


Past the linear stage, QLD starts to modify the particle distribution functions, and
drives the particle distribution towards isotropy in the wave frame. In Figure
<ref>, we show the time evolution of the particle distribution function
(δ f/f_0) via four successive snapshots of the our fiducial run. For each snapshot,
the color scale represents occupation at a given p and μ. We have also
marked with dashed lines the loci in phase space resonating with the
same wavelength

    pμ=Ω_c/k= const .

From Figure <ref>, it is clear that over time a deficit of particles develops near μ∼ 1 while an excess develops near μ∼ -1.  This change in δ f reflects the tendency of QLD to scatter CRs towards a new distribution that is isotropic in a frame that is moving to the left with respect to the initial CR distribution.  Figure <ref> also shows that an excess and deficit of particles develops respectively to the right and left of μ=0; we discuss this further below.

In Figure <ref> we show the wave spectrum for the same
snapshots as Figure <ref>. In these plots, we mark with dashed vertical lines the k values for
the corresponding dashed lines in Figure <ref>.
Over time, the spectrum of forward-propagating waves grows at all but the lowest and highest k, with comparable amplitudes for left- and right- polarizations.



  §.§.§ Overall Evolution


At t∼3 to 10×10^4Ω_c^-1, QLD has started to make appreciable
changes to the particle distribution function. The initial QLD is the fastest around
particles in resonance with the fastest growing mode at
k∼ k_m=3.3×10^-3, as can be clearly identified in the second
panel of Figure <ref>.[The initial wave spectrum also
leads to some QLD that slightly affects the appearance of the distribution function (owing
to contributions from high-k modes) at early times t≲5×10^4Ω_c^-1.]
As discussed in Section <ref>, the way QLD acts to isotropize the CR
distribution is by scattering forward-traveling CRs to from higher to lower (positive) parallel
velocities, and by scattering backward-traveling CRs from lower to higher (negative)
parallel velocities. This is exactly what we see in Figure <ref>, with excess
and deficits on the left and right sides of the thick black dashed line.

With some CRs resonant with these waves giving up their free energy through
QLD, the growth of the most unstable mode slows down. 
More slowly growing modes at longer wavelengths have not yet grown to sufficient
amplitudes to yield appreciable QLD to particles in resonance with them at this time,
thus still grow at the same rate. As a result, the peak of the wave spectrum shifts towards
longer wavelength, as seen from the time sequence in Figure <ref>.

Such wave growth then leads to further QLD for more energetic particles, as can be
seen in the later two panels of Figure <ref>. By the end of our simulation
at t=10^6Ω_c^-1, all particles with p≲10 p_0 have undergone
substantial QLD, whereas for particles with p>10p_0, only the ones with smaller
pitch angle have undergone some QLD. 
While the trend should continue over time, we are also limited by the simulation box
size, with corresponds to 50 times the most unstable wavelength. Modes available for scattering
particles with p>10p_0 are fairly limited and discretized, precluding following QLD
accurately (the discreteness can already be seen in the last panel of Figure <ref>).
We thus do not run this simulation for longer.

Across μ=0, there is substantial excess (deficit) of particles with μ>0 (μ<0)
for p_0≲ p≲ 10p_0, which clearly illustrates the 90^∘ pitch angle
problem. Evidently, redistribution of particle pitch angle occurs only on either side of
μ=0 but particles do not scatter across μ=0. 
In part, the failure to scatter across μ=0 owes to lack of wave power at large
k_x; particles with small |μ| would only be able to resonate and scatter off of
short-wavelength waves.  In our simulation, k_x > 0.1 waves have very low amplitude
due to damping (see Figure <ref>; waves at k_x >  π N_x/L_x are
not resolved at all). In addition, for run Fid the amplitude δ B/B_0 is
insufficient to trigger reflection near m=0 for the bulk CR particles.  We further
discuss the circumstances that lead to “pre-mature” saturation in the next subsection.

Particles with p≲ p_0 have also undergone QLD, which is caused by
shorter-wavelength modes with k>k_m. These modes initially grow slower than the
fastest growing mode at k=k_m, but they catch up towards later time, and
eventually saturate at similar amplitudes, as can be seen from Figure
<ref>. Correspondingly, the evolution in the distribution function
also propagates towards lower p side, as seen in Figures <ref>.
However, this is eventually limited by numerical dissipation. Given our resolution, we
already see from Figure (<ref>) that modes with k≳ 6k_m grow much
slower than theoretical expectation due to numerical dissipation. Considering scattering
across the entire pitch angle range, and given that even particles with p∼ p_0
are subject to pre-mature saturation, we do not discuss particles with p<p_0 in general.

We have also examined run M5. The overall evolution is largely similar to run Fid:
the evolutionary stage at time t in run M5 is similar to that at time t/10 in run Fid,
which is also visible from Figure <ref>. This is reasonable, given the exactly
factor of 10 difference in wave growth rates. QLD starts to modify the CR distribution
function after t∼3×10^5Ω_c^-1, and by the end of the simulation at t=1.5×10^6Ω_c^-1, QLD is already close to being completed for particles
with p∼1-3p_0, except that they are stuck at μ=0. Given even smaller
wave amplitudes in run M5 than in run Fid, we do not expect further evolution would help
with this 90^∘ problem.






  §.§.§ Quantitative Analysis


To see how well our simulation results compare with the quasilinear theory, we
transform the distribution function to the wave frame using Equation (<ref>).
For each particle bin (p_w, μ_w) at a given snapshot, we can find its resonant
wave number k_ res and its corresponding I(k_ res), and compute the
scattering frequency ν from Equation (<ref>). Then, we solve the QLD
Equation (<ref>) over time numerically, with I(k) being updated over a time
interval of Δ t=500Ω_c^-1 (which is the frequency of simulation data output).

In Figure <ref>, we show four snapshots in this process. Note that for better
comparison, we have chosen snapshots 1, 3 and 4 to correspond to snapshots 2, 3, 4 in
Figures <ref> and <ref>. The scattering frequency gradually
increases to the level of ν_ max∼10^-4Ω_c as the waves grow, and
its appearance in p-μ space again shows the characteristic features corresponding
to pμ=const from the resonance condition (modulo that ν decreases towards
larger p as particles become relativistic). 

With short wavelength modes being damped by numerical dissipation, the scattering
frequency for small values of |pμ| has a cutoff, leaving a gap around μ=0
that is broader towards smaller p. This causes the 90^∘ barrier for QLD, where
no particles can cross μ=0, as is the case when solving Equation (<ref>).

In the bottom panel, we show the particle distribution function in the wave frame.
Initially, f_w-f_0 is simply given by Equation (<ref>) taking δ f=0,
which is linear in μ_w.
Note that for particles with different momenta, both
ln f_0/ln p and v_w can be different, yielding different slopes. With the
κ distribution, the slope is smaller for smaller p, and approach a constant
2(κ+1)Δ v/c for relativistic particles (thus the red and green dotted lines
overlap). Subsequent evolution has already been discussed in Section <ref>
in the simulation frame, which is straightforward to be translated to the wave frame.
Several important results and features are worth noticing.

First, at early time (t≲2×10^5Ω_c^-1), evolution of the numerical distributions follow the QLD prediction 
remarkably well. This gives us confidence that our simulations properly capture QLD,
thanks to the phase randomization technique that we employed.

Second, QLD drives the local slope of f_w-f_0 towards zero 
on either side of μ=0 [see Equations (<ref>), (<ref>)];
piecewise flattening is
achieved at later stages of evolution, reached by lower energy particles first. 
This marks the saturation of the CRSI in the local momentum space.

Third, particles undergoing QLD tend to get stuck at μ=0. 
Eventually, 
the CRSI saturates because f_w=const on both
sides of μ=0 for each p, although there is a discontinuous jump at μ=0.
This is the hallmark of pre-mature saturation of QLD.

Fourth, 
solving the QLD equation allows us to well reproduce the saturated state
for lower-energy particles
(p≲3p_0), but this is not true for higher energy particles (p≳10p_0). 
Given that the only missing ingredient in the QLD equation (<ref>) is the process that governs the crossing of 90^∘ pitch angle, this indicates
that the μ=0 crossing does happen (slowly) for higher-energy particles in our simulations.

Given the above results, we now check the requirements for QLD to 
scatter particles to the
critical pitch angle such that the condition for mirror reflection
(<ref>) can be satisfied. Seen from Figure <ref>, we
have the rms δ B/B_0∼9×10^-3, giving μ_ mir∼6×10^-3.
Within our simulation time of 10^6Ω_c^-1, from Equation (<ref>) we need kI(k)≳10^-6
for effective QLD. Reading from Figure <ref>, this gives a maximum
wave number k_ cut∼8k_m. Based on resonance conditions, the
minimum momentum for which QLD can scatter the particle to reach
μ=μ_ mir is then 

    p_ min∼k_m/μ_ mirk_ cutp_0∼20p_0 .

In other words, waves in our simulation only have sufficient power to scatter
particles with p≳ p_ min to reach pitch angle μ_ mir, potentially
allowing the mirror effect to yield reflection (and deviate from pure QLD result).
However,
from the bottom right panel of Figure <ref> representing the
last snapshot of the simulation, we see that even for p as small as ∼3p_0, deviation from QLD
can be identified, though just by a very small amount. Reading off the figure, 
deviation from QLD appear to begin for p_ min between 3p_0 and 10p_0, which is in tension with the expectation  from
mirror reflection (Equation <ref>). Further discussion of this issue is deferred to Section <ref>,
where a more detailed study based on run M3 is presented.  





  §.§.§ Evolution of CR Drift Velocity


While it is generally convenient to use the reduction of CR bulk drift velocity to describe how much
CR particles have undergone scattering through the CRSI, this drift velocity is clearly energy-dependent.
We can evaluate the energy-dependent drift velocity v_d(p) (to distinguish from the initial CR bulk
drift velocity v_D), in the simulation frame, as 

    v_d(p)=∫ f(p,μ)v(p)μ dμ/∫ f(p,μ)dμ ,

where v(p) is the velocity that corresponds to momentum p, and f=f_0+δ f. Values of v_d(p) can be simply
obtained integrating over a horizontal line in Figure <ref>. 
We can further integrate the distribution function over p, to obtain the full drift velocity v_d, full.

In the left panel of Figure <ref>,
we show the results for four different representative particle momenta from p=p_0 to p=33p_0.
Based on the previous discussion, the results are straightforward to interpret. The initial drift velocity in the simulation
frame is by definition 0. Low-energy CRs, which are responsible for driving faster wave growth, undergo more rapid initial
QLD and a corresponding reduction their drift speed; higher-energy CRs take much longer time to catch up. However, the CRSI and drift reduction of low-energy CRs saturates pre-maturely, leading to only modest reduction of v_d, whereas
higher-energy CRs (p≳10p_0) can overcome the μ=0 barrier and are on the way to achieve complete isotropy. The overall average drift v_d, full decreases over time. Contribution to the overall reduction is mostly from lower-energy particles with momenta between ∼ p_0 and ∼3p_0 at early
times, and from higher-energy particles with p∼10p_0 at late times.

The drift velocity is frame-dependent. Similar to the above procedure, we can measure the drift velocity in the wave
frame v_d,w(p), with the results are shown in the right panel of Figure <ref>. While the velocity
difference between the two frames is just v_A, initial particle drift speeds within individual momentum bins do
not equal v_A. This is because after a frame transformation, particles with a given momentum p
come from a range of momenta near p from the original frame, making the results depend on dln f_0/dln p
(see Equation (<ref>)). Only the full drift speed v_d,w, full precisely equals v_A, as
expected. On the other hand, 
in a fully saturated state, we would expect  v_d,w(p)=0 for all p. 
However, we see clearly from the Figure that low-energy particles are stuck
at nonzero v_d,w(p) whereas higher-energy particles require a longer time to
achieve saturation.







 §.§ Quasi-Linear Evolution for run M3


With n_ CR/n_i a factor of 10 larger, run M3 evolves much more rapidly than run Fid. We
choose four representative snapshots at Ω_ct=10^4, 3×10^4, 10^5
and 9×10^5 and analyze the results.



  §.§.§ Overall Evolution


In Figure <ref>, we show the particle distribution function at these snapshots,
while Figure <ref> shows the corresponding wave intensity spectra.
As discussed earlier, right-handed modes grow faster near k∼ k_m. This continues
to t≳10^4Ω_c^-1. Since these modes are resonant with backward-traveling
CRs, we see that at early time, QLD proceeds much more rapidly for μ<0 particles
resonant with these waves, leading to an initial asymmetry in the particle distribution function.

With asymmetric wave growth and QLD, backward traveling CRs have used up more
free-energy, thus slowing down the growth of right-handed waves. The left-handed modes,
on the other hand, continue to grow normally. This allows left-handed modes to eventually
catch up with the right-handed modes, and we see that by the time
t=3×10^4Ω_c^-1, left and right handed modes at all wavelengths 
have similar amplitudes.  
At the same time, the distribution functions for forward and backward traveling CRs also
tends to become more symmetric about μ=0. However, the μ=0 barrier remains
at this stage for low-energy particles.

In between Ω_ct=3×10^4 and 10^5, we see that crossing of the μ=0
barrier occurs, paving the way for full saturation. The release of more free energy has led to
some further wave growth. This process continues until the end of our simulation, with
QLD continuing to extend to higher and higher energies.









  §.§.§ More Detailed Analysis Towards Saturation


As before, for run M3 we compute the CR drift velocity in four different momentum bins for both
the simulation frame and the wave frame, and show the results in Figure <ref>.
We see that in the simulation frame, particles of all four momentum bins eventually flatten
at different levels of v_d(p).  In the wave frame, they all reach v_d(p) equal to zero[The reason that the
saturated drift velocity in the wave frame is not exactly zero in Figure <ref> is purely
numerical, depending on the resolution employed in binning particle momentum during
simulation output (due to the sensitive dependence of f_0 on p). Finer binning makes the
results closer to zero but with more noise, and vice versa.], indicating that they are
fully isotropized and have  achieved complete saturation.

This saturation is more evident from Figure <ref>, showing a time sequence of
the full distribution function in the wave frame, with the same snapshots as in previous Figures
<ref> and <ref>. We can clearly identify the isotropization of
particles in the later two snapshots, extending from the range of p_0≲ p≲10p_0
at t=10^5Ω_c^-1 to covering most of the particle population at the end of the simulation. Because of this, and with most of the particle population having momentum close to p_0, we see that the full drift speed v_d, full
averaged over all particle population shown is reduced almost exactly
by v_A (Figure <ref>).

For run M3, we have also computed the scattering frequency according to Equation
(<ref>) and solve the QLD equation (<ref>) based on output data, and
show the results at four different snapshots in Figure <ref>. The asymmetry
in scattering of forward/backward traveling CRs at early times (t∼10^4Ω_c^-1)
can be clearly identified in the leftmost upper panel. At later times, we can see by comparing with Figure <ref>
that the extent to which QLD has proceeded in momentum space is strongly related to the region
where scattering frequency is high. 

The bottom panels of Figure <ref> give more
quantitative information about how quasi-linear evolution has progressed over time.
For lower-energy particles (with p∼ p_0-3p_0), evolution proceeds rapidly and can initially
be well reproduced by solving the QLD equation. During subsequent evolution, particle
distribution functions in fixed momentum bins are largely flat (due to rapid QLD) on both sides with a jump at μ=0.
The jump gradually decreases, as particles get scattered across μ=0, which is the
rate-limiting process for particle isotropization.
For higher-energy particles (p≳10p_0), on the other hand, crossing of μ=0 is
more rapid (leading to substantial deviation from directly solving the QLD equation). The
rate-limiting process for isotropization then becomes QLD itself, especially at large μ.

A key result from our run M3 is efficient particle crossing of the μ=0 barrier.

We again first examine whether this is consistent with mirror reflection.
From Figure <ref>, we see that by the time of t=6×10^4Ω_c^-1, the
magnetic energy density approaches saturated value (δ B^2/B_0^2∼2×10^-3),
corresponding to a δ B/B_0∼0.045. This gives the mirror reflection threshold of
μ_ mir∼0.032. Within a simulation time of ∼10^5Ω_c^-1, we need
kI(k)∼10^-5 for effective QLD. This is achieved for k up to k_ cut∼10k_m,
by looking at Figure <ref> at t=10^5Ω_c^-1.
According to the same calculation as in Equation (<ref>), we find mirror reflection can
be achieved for particles with p≳ p_ min=(k_m/k_ cut)μ_ mir^-1∼3p_0.
However, we see that by this time, particles with p≳ p_0 are already fully isotropized,
and lower-energy particles have also undergone partial isotropization. Again, we conclude

mirror reflection appears insufficient to explain the detailed evolution of the  momentum distributions.





 §.§ Overcoming the 90^∘ Barrier


To reveal the mechanism behind particle crossing the 90^∘ pitch angle in run M3, we have
randomly selected a subsample of particles (1280 per energy bin), and closely followed their
trajectories for a time interval of 200Ω_c^-1 starting from t_0=6×10^4Ω_c^-1.


We count the number of particles that have undergone reflection across 90^∘ over this
interval. Here, we consider particle pitch angle in the wave frame that is relative to the full
magnetic field vector (rather than background field), which is marked by μ'_w.
To be counted, we require that μ'_w  before and after this time interval to be beyond
± v_A/ℂ∼0.0033 (chosen somewhat arbitrarily, but this does not qualitatively affect the
counting result). Among the 4 central momentum bins from p=0.1p_0 to 10p_0, we find 9,
13, 20, 23 particles are reflected. We looked at the trajectories of all these particles, and choose to
show four particles with characteristic behaviors from each momentum bin in Figure <ref>.



  §.§.§ Representative Cases


For the first particle with p≈0.2p_0 that undergoes a reflection, we see that surprisingly,
the field strength it experiences is almost constant. Quantitatively, the perpendicular
field strength over this time interval is δ B∈(0.0398, 0.0417). The range it spans translates to
μ_ mir=1.3×10^-3. While the particle does evolve quasi-adiabatically with many
gyrations (which is typical for lower energy particles given their small gyro-radii and lack of power at
small scales), its μ'_w changes by a factor of several larger than the expected μ_ mir.
Checking the magnetic field profile in the bottom panel, we see that the particle is traveling against
a shallow gradient in B_z. The range of B_z experienced by the particle ∈(-0.081, -0.041),
which changes the magnetic field orientation by ∼0.04. While this is no proof that the particle
must be reflected, it can well accommodate the change in particle μ'_w, which is less than
5×10^-3.

For the second particle with p≈0.9p_0, its higher velocity allows it to traverse larger distances
over the course of reflection, and experiences more rapid magnetic fluctuations. We have again
checked that variations in total field strength are insufficient for mirror reflection. Reflection occurs
again when the particle travels against the gradient of individual field components, and in this time,
we have B_y increasing rapidly and B_z decreasing more gradually. These variations change the 
magnetic field orientation by ∼0.05 within just a few particle gyrations, which is again more than
sufficient to account for the change in particle pitch angle by ∼0.02.

In the remaining two examples, reflection occurs within two gyro-orbits. These cases are more common
among relatively energetic particles: because of their larger gyro-radii, they more often experience
abrupt changes in magnetic field.

These changes average out to the leading order in quasi-linear theory, except when they
lead to reflection. We see that the threshold pitch angle for reflection can be higher for these particles,
and can reach as large as μ∼0.1. This is accompanied by abrupt changes in B_y and
B_z at comparable levels within the distance that the particle travels in one gyration.




  §.§.§ Non-linear Wave-Particle Interaction


Overall, the four examples roughly cover the range of behaviors of particles during reflection. They
can be interpreted under the framework of non-linear wave-particle interaction <cit.>.

Reflection is most effective (almost instant) for the last two cases in Figure <ref>,
for which abrupt changes in fields occur on a scale shorter than the distance that particle travels in one
gyro time. In this regard, we can discuss the critical pitch angle below which non-linear wave-particle
interaction takes over to reflect particles as follows. For a particle with momentum p (and velocity
v), let μ_ crit be this critical pitch angle. It travels by
L_ crit∼ 2π vμ_ crit/Ω in one gyro time, corresponding to a wave number
k_ crit=Ω/(vμ_ crit)=Ω_c/(pμ_ crit). Reflection occurs when an
abrupt change in magnetic field orientation is comparable to μ_ crit. Note
this depends on the spectrum of waves at scales below some ζ(k_ crit)^-1,
we arrive at a relation

    μ_ crit^2∼ C∫_k_ crit/ζ^∞ I(k)dk ,

where there is a square in the left hand side because μ_ crit∼δ B/B_0, while
the right hand side corresponds to its square. We have also include dimensionless factors ζ
and C that encapsulate the roughness of the model.



The solution to this equation gives μ_ crit, which exhibits the following features


  * Higher wave intensity leads to smaller μ_ crit, because more wave
power at smaller scales makes reflection easier.

  * Higher particle momentum leads to larger μ_ crit, because larger gyro scale
accommodates more wave power to build up abrupt field profiles.

  * If there is a lack of power at small scales, there can be no solution, meaning that direct
reflection in a couple gyro orbits is unlikely.

Clearly, all these features qualitatively agree with our findings. They are also justified
through our discussions in Section <ref>. 

To be more quantitative, we note that while the integral of I(k) over k gives the rms
fluctuations, actual fluctuations can easily exceed the rms value, and deviations by up to
3σ are not rare. We thus take C∼3^2 to account for the upper end of
fluctuations. We may further take ζ∼4 considering that waves with wavelengths
up to 4L_ crit can contribute significantly to the local slope within this length scale of
L_ crit (i.e., sinusoidal variations have roughly constant maximum slope in just 1/4 of
its periodic domain consecutively). With these numbers, and taking p=4.5p_0 for run M3,
we can derive μ_ crit∼0.1 (which gives
k_ crit/ζ∼1.8×10^-3, and the integral gives 0.0011).

Equation (<ref>) closely resembles Equation (54) of <cit.>, with their
v_ osc∼ 0.25-0.3(vμ_ crit), interpreted as the rms velocity oscillation
amplitude of a particle in “the fluctuation field that contributes to its diffusion". The
effective diffusion coefficient at μ=0 is then given by
∼3Ω v_ osc^2.

The discussions above apply to instant reflection upon encountering abrupt field gradients.
This corresponds to the n=1 term in Equation (48) of <cit.>.
As one lowers particle momentum, μ_ crit decreases, and eventually no
solution can be found. This situation applies to lower-energy particles in Figure
<ref>. More generally, higher order effects enter (corresponding to the
n>1 terms), allowing particles to be reflected over multiple gyrations in a fluctuating
field. Being less effective, it is consistent with fewer crossing particles found in our
particle counting in lower-energy bins.



§ SIMULATION RESULTS: PARAMETER STUDY


Our fiducial run selects a parameter range that is physically more
realistic but numerically challenging, with small n_ CR/n_i and low CR
drift velocity v_D, leading to relatively small growth rate and relatively low
wave saturation level. In this Section, we vary parameters around our fiducial run,
and conduct a two-part parameter study.
In Section <ref>, we examine
the saturation of the CRSI when v_D is increased to higher values.
In Section <ref>, we discuss results of 
varying solely numerical parameters for our run Fid. 



 §.§ Dependence on Drift Speeds


With v_D increased to 4v_A and 8v_A, we find that we recover the
linear growth rates very well (not shown), with left and right handed modes growing at
comparable rates (as expected).
Figure <ref> includes the time evolution of wave energy for runs vD4 and
vD8 in dotted lines. They initially follow very similar evolutionary paths as our run Fid,
except that everything (linear growth and QLD) proceeds faster.
More notably, later evolution of these two runs are more similar to run M3, where
wave energy eventually plateaus, indicating saturation is close to completion.

In Figure <ref>, we show results related to quasi-linear evolution of runs vD4 and vD8. The overall outcomes are in between those of run Fid and run M3. Crossing
across 90^∘ becomes easier as v_D increases. By the end of the simulations,
particles with 4p_0≲ p≲20p_0 are roughly fully isotropized in run vD4.
For run vD8, the range enlarges to p_0≲ p≲40p_0.
[Two numerical artifacts can be identified in Figure <ref>.
The banded horizontal feature in the top panels is due to the division into 8 particle
momentum bins. Within each bin, there is more noise towards higher energies because
there are fewer particles. One may also find that particle mean drift velocities at various
energies do not reach zero even they are theoretically expected to be fully isotropized
(bottom panel). This is related to coarse particle binning during outputs (mentioned
earlier in Section <ref>).]
Lower-energy particles are in the process of (but have not yet completed)
isotropization. Isotropization is limited by the rate of 90^∘ crossing, which is much slower for these lower-energy particles (see discussion in Section <ref>).





  §.§.§ The Saturation Level


At this point, we can gather all simulation runs discussed so far to address the saturation
level of the CRSI as a function of v_D. The saturation level is measured in terms of
wave energy density. With equipartition of magnetic and kinetic energies, it can simply
be represented by ⟨δ B^2/B_0^2⟩ in dimensionless form.
Theoretically, the predicted wave energy saturation level is obtained by equating Equations (<ref>) and (<ref>),
and is proportional to (n_ CR/n_i)Δ v. Integrating the derivative of the
distribution function numerically, the result is shown as the dashed line in Figure
<ref>.

The final wave energy density obtained from our runs Fid, vD4, vD8, M3, and
Fid-hires (discussed in the next subsection) are also shown in Figure <ref>.
Even though runs M3 and vD8 are close to complete saturation, their final
wave energy densities appear to be far from the expected level. This is because 
numerical dissipation, while very low in the Athena MHD code, is still non-negligible
for the extremely long simulation time necessary to reach saturation.
[ Dissipation of low-amplitude waves imposed in the initial conditions is also
evident at early times in Figure <ref>.]
With an adiabatic equation of state, we can capture the change in gas internal energy
over time (thanks to the excellent energy conservation properties in the gas component),
and find that it increases steadily as result of numerical dissipation, and reaches a level
that is comparable to the overall wave energy density by the end of the runs. When we add
this energy lost to dissipation to the direct wave energy, we find much more reasonable
saturation levels, as shown in thick/larger symbols in Figure <ref>. In particular, 
run Fid remains at a level below the full saturation prediction because particles become
“stuck” at 90^∘, whereas runs vD4 and vD8 are more and more close to the expected
saturation level, as more particles are fully isotropized. Run M3, which based on its distribution
function is clearly isotropized, is excellent agreement with theoretical expectations for the
wave amplitude.  












 §.§ Numerical Parameters


Here, we examine our runs Fid-Hires, Fid-short, Fid-Np64 and Fid-Np1024, in which we keep physical parameters the same and vary only  
numerical parameters (note that for run Fid-Hires, we have slightly reduced
the box size to save computational cost). In Figure <ref>, we show the time
evolution of the wave energy density of these runs. Overall, these three runs proceed very
similarly to run Fid.[All these runs have the same initial wave amplitude A=10^-4 (see Equation <ref>). Run Fid-Hires thus has more initial wave energy because of more
spectral range it covers.] QLD proceeds slightly slower in run Fid-Short, but it catches up to the
other runs shortly after. Run Fid-Hires eventually grows to reach larger wave amplitudes, which we
will show later is due to more efficient crossing of 90^∘ pitch angles.

In Figure <ref>, we show the wave spectrum from these runs at time
t=2×10^5Ω_c^-1. The overall shape of the spectra are comparable to each other.
The high-resolution run has less numerical dissipation, allowing the wave spectra to extend to
larger k. Run Fid-short has a box size that is only 1/3 of the fiducial box (but still long
enough to fit ∼17 most unstable modes), and hence modes are a lot more sparsely sampled,
especially at low-k. This is likely related to its initially slower quasi-linear evolution. We have
also tried simulations with 1/8 the fiducial box size (not shown); while evolution is smooth, the
lack of power at larger scales start to adversely affect the evolution of higher-energy particles.

Figure <ref> shows that there are backward-propagating modes excited in all these runs. The amplitude
of these modes reflects the level of noise. Clearly, inserting more particles helps reduce the noise.
The bottom line for properly choosing particle number is that the amplitude of such
backward-propagating modes should always be much smaller than the forward-propagating
modes driven by the CRSI. In this regard, our run Fid-Np64 is acceptable, but more caution
should be exercised if one were to use even fewer particles.

In Figure <ref>, we show the particle distribution functions in the wave frame at
time t=10^5. These may be compared with Figure <ref>. Overall, runs
Fid-short and Fid-Np64 are more noisy, whereas run Fid-Np1024 appears exceptionally smooth.
Our run Fid lies in between. Despite these  differences in noise levels, 
 subsequent QLD does not appear to be affected.   

Of particular interest is the high-resolution run Fid-Hires, where waves grow to substantially
higher amplitudes. In Figure <ref>, we show three snapshots of the particle
distribution function in the wave frame for four representative momentum bins. These frames are to
be contrasted with Figure <ref> for the same snapshots. Clearly, in this run, a
substantial fraction of particles with p≳3p_0 have managed to cross the 90^∘
barrier. This fraction steadily increases with time and, were we to run the simulation longer
(e.g., to t=2×10^6), particles in these momentum bins would essentially be fully
isotropized. This trend can also be seen from the last panel, which shows the time evolution
of the mean CR drift velocity: v_d for particles with p≳3p_0 are steadily decreasing
towards zero. Clearly, crossing the 90^∘ barrier is easier compared to run Fid. This is
natural because higher resolution allows more small-scale modes to be better preserved,
thus enhancing the chance of non-linear wave-particle interaction to reflect the particles
(see Section <ref>). In Figure <ref>, we see that the wave energy
density in our run Fid-Hires is also a factor of several closer to the expected saturation
level.  We expect that increasing the resolution further (which would be much more computationally
costly) would eventually lead to full isotropization of the bulk CR populations.

In sum, our exploration of numerical parameters suggests that full saturation of the CRSI
(ideal case, no damping) can be achieved as long as there is sufficient numerical resolution.
The required resolution is mainly limited by the demands imposed by crossing the 90^∘ pitch angle,
which is sensitive to wave amplitudes at small scales. Similar results to our fiducial model can be achieved
when one reduces the number of particles and/or the box size at some modest level (∼ a
factor of ∼3). 






§ DISCUSSION




 §.§ CR Feedback to Background Gas


Through the CRSI and QLD, CRs transfer momentum to background gas, which serves as
an important dynamical feedback mechanism in broad range of astrophysical 
environments.
In the following, we take run vD8 as an example and discuss CR feedback in terms
of momentum and energy.

Momentum feedback is manifested in the top panel of Figure <ref>, which
shows changes in the mean gas and CR momentum density 
(Δ P_g and Δ P_ CR).
In our simulation frame, the initial CR momentum is zero, and initial
gas momentum is negative. Over time, the CR momentum density decreases, accompanied by an
increase in gas momentum density of exactly the same amount.
Physically, the backreaction from the CR isotropization process effectively gives rise to a “parallel force"
F_∥. 
This parallel force is mediated by forward-propagating Alfvén waves,
and transfers CR momentum to the gas (e.g., ).
The amount of momentum transferred to the gas is in fact the same as the “effective” wave
momentum defined in Equation (<ref>), which is the momentum needed to excite such
forward-propagating waves <cit.>.
By the end of our run vD8, although the evolution of gas and CR momenta are not yet
fully saturated, both are approaching the expected saturation level indicated by the
dashed lines, obtained from integrating Equation (<ref>). In the case of
run M3 (not shown here), we do observe that saturation of gas and CR momenta is
achieved.  
That is, the mean CR streaming velocity relative to the gas is reduced from 2 v_A to v_A (Figure <ref> shows that the final mean CR velocity is zero in the wave frame), and the gas momentum density is increased to balance the loss from CRs.

Energy feedback can be discussed by decomposing the energy budget into various
components. These include gas bulk kinetic energy
E_k, bulk=⟨ρ v_x^2/2⟩, representing gas motion parallel to
background field, wave energy E_ wave, CR kinetic energy E_ CR,
and internal energy of the gas E_ internal=P_g/(Γ-1). 
The time evolution of these quantities is shown in the bottom panel of Figure <ref>, where
we only plot their changes relative to the initial condition.
The overall energy conservation is expressed as

    E_k, bulk+E_ CR+E_ wave+E_ internal= const .

Note that E_k, bulk and E_ CR are frame-dependent.
Here it suffices to report results in our simulation frame.

In the simulation frame, the gas bulk energy density decreases over time (as the bulk gas velocity becomes less negative). This is associated by the work done by
the same parallel force F_∥ that changes the gas momentum. In fact, we can
exactly reproduce the black solid line in the bottom panel of Figure <ref>
using F_∥ inferred from acceleration in the gas (from the top panel), with
Δ E_k, bulk≈ v_gΔ P_g.
(Since v_g only changes by a tiny fraction over the process, it is approximately taken to be constant here.)
This same force also does opposite work on the CRs, increasing their total kinetic energy.
Some of this energy is used to excite and amplify the waves, leading to their continued
growth. In the mean time, the waves are slowly dissipated numerically, increasing the
internal energy of the gas.[While damping in our simulations is purely numerical,
we can see that it is negligible in early phases of linear and quasi-linear evolution. Only
towards later stages where saturation is largely achieved (and hence there is little wave growth)
does numerical dissipation become more noticible. Thus, our results are largely unaffected
by numerical dissipation.]

We emphasize that the CRSI itself does not lead to heating of background gas. Heating
only occurs through internal dissipation in the gas, which is generally accompanied by
wave damping. This fact is reflected in the more recent formulation of
CR-(magneto)hydrodynamics of <cit.>, though not necessarily
correctly modeled in some earlier studies. 
If wave growth is balanced by wave damping in a quasi-steady state (not achieved in
our present simulations), then the heating rate is simply determined by the rate of wave
growth/damping. 



 §.§ Towards More Realistic Simulations


There are several directions to further improve the realism of our simulations.
In particular, there are several wave damping mechanisms to be incorporated,
including ion-neutral damping, non-linear Landau damping, and turbulent damping of waves. 

Ion-neutral damping, or ampibolar damping, is well known (e.g.,
). It is present when background gas is partially
ionized. The ions are directly coupled and respond to magnetic fields, but the
neutrals hardly feel the magnetic fields and can possess different bulk velocity from the
ions. Ion-neutral collisions give friction, and lead to damping of Alfvén waves
in the ion fluid, especially at small scales. Given that crossing the μ=0 barrier
is sensitive to the spectrum of waves at small scales, ambipolar damping may cause
severe challenge to the isotropization of the CRs.
This effect will be incorporated and studied in greater detail in a follow-up paper.

Non-linear Landau damping is a result of wave-particle interaction at Landau resonance
between background thermal particles and a beat wave formed by two propagating
waves <cit.>, and is more effective in high-β plasmas. Our MHD-PIC
formulation do not capture the kinetic physics of the background plasma, and hence can
not capture non-linear Landau damping. However, this is not necessarily a limitation: we
would lose the advantage of MHD-PIC by resolving the kinetic physics, and in principle,
this damping mechanism can be incorporated in the form of sub-grid physics.

Turbulent damping is due to the distortion of CR-generated Alfvén waves when they
collide with oppositely directed wave packets originating from background turbulence.
As a result, the wave energy cascades into smaller scales and is ultimately dissipated
<cit.>. This damping rate has only been calculated approximately
(e.g., see more recent work of ). 
Incorporating this effect is fully compatible with our MHD-PIC
framework, and can be studied in future multi-dimensional simulations.

For this first study, we focused on the classical theory with 1D simulations with uniform
background density and magnetic field, which enable us to confirm resolution and box
size requirements. Multi-dimensional extensions are straightforward and are underway,
although are much more computationally expensive. Such extensions would allow modes
propagating obliquely with respect to background field to be investigated, and open up
more parameter space, as the relative importance of thermal, magnetic and CR pressures,
and background turbulence can play certain roles, 

Closely related to the CRSI is another type of gyroresonance instability, which is
driven by CR pressure anisotropy (e.g., ),
and can be easily triggered in compressible MHD flows (when CR pressure anisotropy 
exceeds ∼ v_A/c). It excites circularly polarized Alfvén waves propagating to
opposite directions, which then isotropize the CRs via QLD. Recently, <cit.>
studied this instability numerically based on the MHD-PIC framework in full f
approach (with relatively strong anisotropy), and found good agreement with quasi-linear
theory. In the absence of CR streaming, the CR distribution function is always symmetric
about μ=0, and hence crossing 90^∘ pitch angle is not required for CR
isotropization. In reality, CR streaming and pressure anisotropy likely co-exist. This
would excite a wide ranges of waves propagating in both directions, calling for further
investigations.

Finally, in any simulation with periodic boundary conditions (even with arbitrarily high resolution such that crossing 90^∘ pitch angle is not a limitation), there is a finite free energy of the system based on the initial energy and momentum densities of gas and CRs.  In the absence of damping (numerical or physical), growth of Alfvén waves would halt when QLD has reduced the mean velocity difference between the gas and the CRs to v_A.  

 In reality, at a given location the CRSI is typically driven by a large-scale CR pressure gradient. For instance, CRs are continuously generated from the galactic midplane region, and diffusively propagate downhill through the gradient to escape from the galaxy.  Energy and momentum are continually added to the system.  
This configuration likely yields continued wave growth, balanced by damping to
reach a steady state, a situation that has yet to be addressed with more realistic
boundary conditions.



 §.§ Comparison with conventional PIC simulations


Our MHD-PIC approach to study the CRSI strongly contrasts with conventional full-PIC or hybrid-PIC approaches that are widely used in modeling highly non-linear plasma problems.
The full-PIC approach treats charged particles of all kinds as particles. It resolves all
scales down to the electron skin depth (eventually the Debye length), and hence all plasma
wavemodes are self-consistently captured along with relevant wave-particle interactions.
The hybrid-PIC approach

treats electrons as a massless conducting fluid. By compromising the electron-scale physics, it alleviates the stringent resolution requirement of the full-PIC method, while
it can still capture the ion-scale physics well.
For the CRSI problem, employing full-PIC or hybrid-PIC methods means that background ISM
fluid must be represented by particles, hence resolving the corresponding microscopic
scales. This causes two major drawbacks: substantially enhanced numerical cost, and PIC
noise due to the finite number of particles representing the background plasma. 

As already discussed in Section <ref>, there is substantial scale separation in
the CRSI problem between the CR gyro-radii and the ion skin depth (and below). Conducting
a hybrid-PIC simulation using the same parameter settings as our fiducial simulations
would require more than an order of magnitude more cells, and another order of magnitude
more cells would be needed in the full-PIC case.
The numerical cost in the full-PIC case is further exaggerated by the ion-to-electron mass
ratio, requiring a factor of m_i/m_e more timesteps per CR gyro period.

Moreover, the number of background particles per cell needs to be sufficiently high in
order to reduce the noise floor as compared to the signal level of growing waves and
avoid artificial QLD of the CRs by such noise.


Despite these difficulties, <cit.> have recently presented a 1D study
of the CRSI using full-PIC simulations. Typical parameters of their simulations are: m_i/m_e=100, υ_A/c=0.1 (note in full-PIC, c=ℂ in our notation),
10/100 cells per electron/ion skin depth, N_ ppc=50-1000 (all plasma particles
equally divided among background and CR species),
n_ CR/n_i ∈ [2× 10^-4, 2 × 10^-2] and  υ_D/υ_A from
1.4 (low-anisotropy) to 7.9 (high-anisotropy) [Let us recall that in present study
(Fid run) the parameters are:
m_i/m_e →∞, υ_A/ℂ=3.3 × 10^-3,

10 ion skin depths per cell, N_ ppc=2048 (CRs only),
n_ CR/n_i=10^-4 and  υ_D/υ_A=2.].
They considered two different CR distribution functions: ring and power-law, with the former
representing a mono-energetic CR distribution in the drift frame.
They have successfully reproduced the peak growth rate in runs with the ring distribution,
and similar to our results, achieved full isotropization in runs with relatively large
n_ CR/n_i and/or CR drift velocity.

We note that the large ratio of v_A/c adopted in <cit.> makes the
resonant condition (<ref>) asymmetric with respect to forward and backward traveling
particles, thus modifying the standard CRSI dispersion relation for left/right handed modes
in different ways. Together with the choices of large n_ CR/n_i ratio, drift
velocity being mildly relativistic, and the use of a ring distribution, these parameters
represent substantial compromises to alleviate the issue of scale separation and particle
noise. As a result, reasonable agreement with linear theory was found mainly in the special
case of the ring distribution, and power-law distribution with exaggerated high anisotropy
and high n_ CR/n_i values. These more extreme run parameters (approaching the Bell
regime), together with the asymmetry in left/right handed modes caused by large v_A/c, make
the left handed modes in these simulations grow at much reduced rates compared with right
handed modes, leading to asymmetric scattering for forward/backward traveling CRs. 


While this situation might apply near strong shocks where the Bell instability operates, the
parameter regime is not representative of typical conditions in the ISM.

Additionally, they have suggested mirror reflection as being responsible for crossing
90^∘ pitch angle. However, the reflection event shown in their Figure 10 in
fact appears more consistent with non-linear wave-particle interaction, as we have
demonstrated in this work.
 
Perhaps the chief advantage of full-PIC and hybrid-PIC methods over the MHD-PIC approach
is the ability to capture non-linear Landau damping, as the background plasma is treated
kinetically. However, the physics of non-linear Landau damping is highly subtle, and
damping relies on second-order longitudinal electric fields resulting from a beat wave
formed from two circularly polarized Alfvén waves <cit.>.
It is conceivable that capturing such effect must require extremely low level of
background noise, and hence exceedingly large number of background plasma particles:
another major challenge yet to be resolved.

Finally, we comment that several of the difficulties for simulating the CRSI outlined in
Section <ref> are common for all variations of PIC methods. The δ f
approach coupled with a κ distribution that we have adopted could also be
incorporated in full-PIC and hybrid-PIC simulations, likely leading to substantial
improvements in future studies of the CRSI.




§ SUMMARY AND CONCLUSIONS


In this paper, we have presented the first numerical studies of the CRSI using the MHD-PIC
method. We use 1D simulations in a periodic box and an idealized setup (no damping, 
a reduced speed of light ℂ), for a range of model and numerical parameters.  We have
been able to accurately reproduce the theoretically predicted linear growth rate, follow quasi-linear evolution as CR particles interact with CSRI-driven Alfvén waves, and
identify the mechanism responsible for CR particles to cross the 90^∘ pitch angle. In the course of this study, we have developed techniques to investigate CRSI and QLD with substantial precision,
paving way for future studies with more realistic setups and incorporating more physics.

Unlike conventional PIC methods, the MHD-PIC method has unique advantages in
bypassing the microscopic scales of the background thermal plasma, allowing one to focus on
the interaction between the CRs and background gas around the  gyro-resonant scale. However,
the CRSI still involves substantial scale separation, ranging from scales much smaller than,
to scales much larger than, the CR gyro-resonant scale. Covering this large range of scales is necessary
to properly capture the crossing of the 90^∘ pitch angle and QLD.
Moreover, the extremely weak level of anisotropy (∼ v_A/c) in the CR distribution function
over a wide energy range makes it challenging to properly represent with a sufficiently large number of particles. 

To address these challenges, we have extended the MHD-PIC formulation to use a δ f
(instead of full f) representation of the CR particles on top of a known initial equilibrium
distribution function f_0, which dramatically reduces the Poisson noise. While our fiducial
simulations employ as many as 2048 particles per cell to achieve
better precision, similar results can be achieved with fewer particles, and we recommend
testing the level of Poisson noise in a non-CR-streaming simulation to determine the optimal
number of particles. 
Using the conventional full f representation, while
the signature of the CRSI can be identified, the high noise level precludes any quantitative
measurement of its overall properties.

The δ f method requires the background distribution function f_0 to be continuous,
which is incompatible with f_0 being a truncated power-law. We have therefore chosen f_0 to be a
κ distribution, and the corresponding dispersion relation of the CRSI is analytically
derived. We note that under realistic large-scale ISM conditions (n_ CR/n_i extremely small and moderate v_D/v_A), the left and right handed Alfvén waves
grow at equal rates. The Bell-type modes, in which the right-handed modes grow
faster, emerge when n_ CR/n_i and/or v_D/v_A become large. We have carefully chosen
parameters to ensure sufficient scale separation v_A/ℂ=1/300, and that the left
and right handed modes growth at about the same rates (except for run M3 where there is
a small difference), so that the conditions are as close to ISM as possible. Our particle
population also covers a wide range of CR momenta from 0.01p_0 to 100p_0, with
p_0=ℂ being the characteristic momentum of the κ distribution.

We have found that when employing a standard periodic boundary condition,
properly capturing the linear growth rate requires an excessively long simulation box, corresponding to the
box-crossing time of a typical particle approaching the wave growth time. 
We attribute this
to the fact that the growth of the CRSI is driven by QLD, which requires particles to
encounter sufficient number of independent wave packets, and is related to the
random-phase approximation commonly employed in describing QLD. 
To address this issue, we implemented a phase-randomization technique for particles that cross the periodic
boundaries, which allows us to achieve even better results using much smaller
simulation boxes (size ≲50 instead of ≳1000 most unstable wavelength).

Equipped with these techniques, we are able  not only to accurately reproduce predicted CRSI linear
growth rates over a broad range of wave spectrum, but also to accurately follow the
quasi-linear evolution of the system. We observe that QLD proceeds to isotropize the
CR distribution function f_w in the frame of forward-propagating Alfvén waves, so
that f_w/μ_w approaches 0 on either side of the pitch angle cosine
μ_w≈μ=0. With an increasingly flat distribution of particle pitch angles, wave growth slows down and eventually stops.
As longer-wavelength modes that resonate with higher-energy CRs grow more slowly,
wave growth proceeds from the most unstable to longer wavelengths (and in principle also
to shorter wavelengths, though we lack resolution for these).  In parallel,  quasi-linear evolution
proceeds from low-energy to higher-energy CRs as waves become available.

Whether full isotropization of CRs in the wave frame is achieved depends on whether
particles manage to overcome the μ=0 barrier (crossing the 90^∘ pitch angle).  The mechanisms involved go beyond quasi-linear theory. While mirror reflection is commonly invoked
as the main mechanism enabling crossing μ=0, we have found that reflection generally occurs
via non-linear wave-particle interaction, with particles encountering abrupt changes in
field orientation. Reflection is more easily achieved for higher-energy particles, and is
facilitated by the presence of waves at wavelength much smaller than the
gyro-resonance scale (see Section <ref> for details). 

In our study, efficient crossing of the μ=0 barrier is achieved in simulation runs
with strong wave growth (runs M3 and vD8). For these runs, full saturation is reached, with the mean CR drift velocity reduced to Alfvén speed.  The final particle distribution function,
as well as saturation amplitudes of the waves, are in excellent agreement with
theoretical expectations.
Our fiducial run, with a relatively low n_ CR/n_i=10^-4 and small initial drift
speed v_D=2v_A, becomes  stuck in a state of pre-mature saturation: while the distribution is flat on each side of μ=0 within most momentum bins, the bulk of the CR
population fails to overcome the μ=0 barrier. The inability of particles to breach the μ=0 barrier is primarily due 
to numerical dissipation of low-amplitude waves at small scales, together with absence of even shorter-wavelength waves due to limited resolution.   We have further found that enhancing
numerical resolution can substantially alleviate the situation. Thus, proper choice of
numerical resolution for given physical parameters is also crucial for saturation.

The CRs provide direct momentum feedback to background gas, which accelerates
in accordance with momentum conservation of the gas-CR system. Wave power is directly 

proportional to the momentum transferred from the CRs to the gas.
In terms of energy feedback, the work
done by the CRs leads to gas acceleration and wave growth, but no direct
heating. Gas heating is only achieved through wave damping, which is
purely numerical in our simulations.

As a first study, we have restricted ourselves to 1D simulations in a periodic box in
the ideal MHD limit. Our explorations 
have identified and clarified the

numerical methods necessary to properly capture the essential 
mechanisms of the CRSI, and have
paved the way for future investigations to incorporate more realistic physics. These
may include implementing various damping mechanisms, extensions to
multi-dimensions, and employing more realistic boundary conditions.



We thank Cole Holcomb, Matt Kunz, Peng Oh, and Anatoly Spitkovsky for useful discussions.
XNB acknowledges support by the Youth Thousand Talent program.  The work of ECO  was supported by grant 510940 from the Simons Foundation. 
The work of IP was supported by NSF grants PHY-1804048 and PHY-1523261 and facilitated by the Max-Planck/Princeton Center for Plasma Physics.
Computation of this work has been conducted at computational facilities at Princeton
Institute for Computational Science and Engineering (PICSciE), and at the YZ cluster
from GyroTech Co. Ltd.





§ A. DECOMPOSITION OF ALFVÉN WAVES


Given that modifications from the CRs to the Alfvén waves are tiny in this context, it
suffices to consider pure linear Alfvén waves. At fixed wavelength, there are four independent
Alfvén modes, depending on left/right-handed polarization, and direction of propagation relative
to fluid rest frame.
We here take the perspective that the left/right handedness is defined from the receiver's
point of view. That is, for a wave propagating toward an observer (relative to the fluid
rest frame), the perturbed electric (and magnetic) field vector rotates in a anti-clockwise (clockwise)
sense for a left (right) handed wave.

Let δ v and δ B be perturbed transverse velocities and magnetic
fields, and take the perturbations in the form of exp[ i(kx-ω t)]. Note that k can
have both signs indicating direction of propagation, whereas ω always has positive real part.
The eigenvectors for Alfvén waves satisfy

    δ v/v_A=- sgn(k)δ B/B_0 .

The sign of B_0 is implicitly contained in the above, whereas v_A is considered to be
positive. In the following, we always take B_0>0 without loss of generality (for B_0<0, sign chages
must be employed in many of the subsequent descriptions), and consider forward propagation to
be along the direction of B_0 (positive x̂). Under this convention, 
the +/- sign of k corresponds to  forward/backward propagation. The left/right-handed modes have
δ B_y=∓ iδ B_z and δ v_y=∓ iδ v_z for forward-propagating waves.
The opposite sign should be adopted in the above for backward-propagating waves.

Based on the above description, the natural way to decompose circularly polarized
Alfvén modes is to consider the following combination:

    C^±(x)=1/2[δ v_y± iδ v_z/v_A
    - sgn(k)δ B_y± iδ B_z/B_0] .

For instance, for forward propagating (k>0), right-handed (δ B_y= iδ B_z)
waves, we choose the plus sign. A discrete Fourier transform of C^±(x),

    W^±(k_i)=1/N∑_n=0^N-1C^±(x_n)e^- ik_ix_n ,
         where k_i=2π/L_xi ,

gives the wave amplitude and phase (which is what is being plotted in the Figures).
The wave intensity is then given by

    I^±(k)=|W^±(k)|^2L_x/2π .

To initialize a wave spectrum with I^±(k)=A^2/k, we can simply choose
|W^±(k)|=A√(Δ k/|k|), where Δ k≡ 2π/L_x.






§ B. LEVEL OF NOISE IN THE Δ F AND FULL F METHODS


In this Appendix, we study the noise properties in our simulations. We do so by conducting
simulations with almost identical setup to our fiducial runs, except that we change the CR
drift velocity to v_D=0, and that we start from an initial amplitude that is much smaller
with A=10^-12. We also use a shorter simulation box L=24000 (which we verify has no
influence on the result).
In this case, we only expect wave damping to occur, but Poisson noise from particle sampling
will stand out and causes fluctuations that overwhelm the initially low-amplitude waves.

In Figure <ref>, we show the time evolution of mean fluctuation energy density in the
system, which is given by (ρ v_⊥^2+B_⊥^2)/2. This is done for both the δ f
method (left), and the full f method (right). Note that we have chosen n_ CR/n_i=10^-3
in tests using the δ f method (a smaller value would take longer time to reach steady state),
whereas we use n_ CR/n_i=10^-4 in the full f case (larger value would make the noise
somewhat too large). Clearly, we see that even with a higher n_ CR/n_i, the noise with
the δ f method is dramatically smaller than that in the full f method by more than five
orders of magnitude (this would be more than six orders of magnitude if they had the same
n_ CR/n_i).

The energy density of the waves in this test results from the Poisson noise in the particles. We
can see that in both cases, reducing the number of particles by a factor of 4 leads to a factor
of ∼4 increase in wave energy density (or a factor √(4) in wave amplitude). The wave
spectrum (not shown) is consistent with white noise (with I(k) being flat) for all modes with 
cutoff at high k due to numerical damping. The spectrum at intermediate k grows and
saturates faster, whereas modes at the low-k end grow slower, explaining the secular trend
at late times.

Interestingly, the noise level also depends on resolution, but the dependence is different for the
full f and δ f methods. For the δ f method, we find a sensitive dependence, which
scales roughly as Δ x^4, and higher spatial resolution can greatly reduce Poisson noise.
On the other hand, for the full f method, the dependence is much weaker and is sub-linear.
In addition, we find that using a larger p_0 while keeping other parameters fixed also gives
very different noise levels. Enlarging p_0 by a factor of 10, noise is dramatically reduced by
a factor of ∼300 when using the δ f method. 

Noise is enhanced
by a factor of ∼3-4 when using the full-f method. We do not aim at further studying and
explaining these trends, but simply report them here for reference.







§ C. RANDOM PHASE APPROXIMATION AND SIMULATION BOX SIZE


In this appendix, we show that simulations without implementing phase randomization
across periodic boundaries fail to properly capture the linear grow rate of the CRSI unless
the simulation box size is extremely long. 
This failure reflects the fact that the random phase
approximation must be satisfied for particles to properly follow quasi-linear diffusion,
which drives the CRSI.

We first show in Figure <ref> a comparison between our fiducial run with
and without phase randomization. At some early time (less than one e-folding time), we
plot the growth factors across the wave spectrum. Despite some slight deviation, our fiducial run
with phase randomization nicely matches theoretical expectation (except for numerical damping
at large k). Without phase randomization, we see that while the fastest growth rate is not
different from theoretical expectation by much, the wave spectrum has significant offset, with
most unstable wavelength being about a factor of three shorter.

Further looking at the time history of the wave intensity of three representative modes
(with the fastest growing mode k=(1/3)×10^-2 shown in black) in the bottom two
panels, we see that with phase randomization, wave growth is clearly exponential. Without
randomizing the phases, the exponential growth is only maintained for up to
t∼5000Ω_c^-1, after which the wave amplitudes wildly vary.



We then conducted a serious of simulations using fiducial parameters but in different box
sizes and without phase randomization. The box sizes are L_x=2.88×10^5,
7.68×10^5 and 1.92×10^6, respectively. They are 3, 8 and 20 times
the length of the fiducial box that we use (0.96×10^5), while they are run for shorter
amount of time (about 2.5-4×10^4Ω_c^-1) just within the expected linear
growth phase. Note that even for the longest box that we have attempted, the box crossing
time for a relativistic particle is 1.92×10^6/300=6.4×10^3Ω_c^-1, which
is about a factor ∼4-5 smaller than the growth time of the CRSI.

In Figure <ref>, we show again the spectrum of growth factors for these runs.
Clearly, the offset in the fastest growing wavelength from theoretical expectation gets smaller
with increasing box size, and eventually vanishes in the run with longest box size. There is
still some deviation in the absolute growth rate, though only by a small factor (the numerical
growth rate is accurate within about ∼10%). We thus consider this longest box size run
to be marginally appropriate to study the CRSI.





§ D. COMPARISON WITH THE FULL-F METHOD


To demonstrate the advantage of the δ f method, we show in this Appendix
results obtained from the full f approach. The simulation setup is identical to our
run Fid except that the δ f weighting scheme is turned off and that the box size is halved. Phase randomization
is also included. To further reduce
Poisson noise, we use N_p=2048 particles per bin per cell (in total 16384 per cell).
According to the study in Appendix <ref>, the level of noise is reduced by
a factor of 8 in wave intensity, with wave energy density ∼10^-4, or
δ B^2/B_0^2≲5×10^-5, in the form of Alfvén waves propagating
in both directions. This level of noise would yield a rate of QLD, according to
(<ref>), that is comparable to the maximum wave growth rate of
∼3×10^-5Ω_c. The simulation is run to time t=2.25×10^5Ω_c^-1.

In Figure <ref>, we show the two snapshots the wave pattern and spectrum.
We see that noise level indeed quickly reaches the level of δ B/B≲10^-2,
causing appreciable QLD of particles. 
Note that since Alfvén waves in both directions
are present due to noise, CR energy would no longer be conserved in the (forward-propagating)
wave frame, compromising the physics of the CRSI. As a result, we can only barely identify
the signature of the CRSI in the wave spectrum at early time of t=2×10^4. At
later time, forward-propagating waves are more easily seen to be amplified. However, these
waves are at longer wavelengths (with k∼10^-3 instead of
Ω_c/p_0=3.3×10^-3) where intrinsic numerical noise is smaller, allowing
the CRSI to better stand out. The growth rates at these longer wavelength are somewhat
too noisy to measure with sufficient precision to compare with theoretical values, and
through our analysis, particles responsible for wave growth must have already undergone 
some QLD from pure noise.



We have further inspected the evolution of the distribution function, and show in Figure
<ref> this distribution function for the last snapshot of the run. Note that because
the binning is based on the full-f prescription (i.e., direct particle counting in each bin),
it is much noisier than we had before under the δ f prescription. On the other hand,
we can find signatures of isotropization. In fact, particles at p∼3p_0 are already nearly
isotropized, and particles with p≳10p_0 are still in early stages of isotropization.
For particles with p≲ p_0, accumulation near μ∼0 is also clearly
visible, as they are stuck and can not cross the 90^∘ pitch angle. While these features
are qualitatively consistent with expectations, QLD proceeds much faster than that obtained
using the δ f method owing to substantial enhancement of scattering through noise.
In particular, recall from discussions in Section <ref>, isotropization is never achieved
when using the δ f method. 

Overall, we conclude that while the full-f method can capture some basic properties of
the CRSI, numerical noise adversely affects its evolution at all stages, leading to severe
discrepancies at a quantitative (and sometimes even qualitative) level, which precludes making any
accurate/reliable measurement.
The situation can only get worse under more realistic conditions (low n_ CR/n_i
and small v_D).


apj
59
natexlab#1#1

[Achterberg(1981)]Achterberg81
Achterberg, A. 1981, , 98, 161

[Amato & Blasi(2009)]Amato09 Amato, E., & Blasi, P. 2009, , 392, 1591

[Amato & Blasi(2018)]Amato18 Amato, E., & Blasi, P. 2018, Advances in Space Research, 62, 2731 

[Baade & Zwicky(1934)]BaadeZwicky34
Baade, W. & Zwicky, F. 1934, Proceedings of the National Academy of
  Science, 20, 259

[Bai et al.(2015)Bai, Caprioli, Sironi, &
  Spitkovsky]Bai_etal15
Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ,
  809, 55

[Bell(1978)]Bell78
Bell, A. R. 1978, , 182, 147

[Bell(2004)]Bell04
—. 2004, , 353, 550

[Birdsall & Langdon(2005)]BirdsallLangdon05
Birdsall, C. K. & Langdon, A. B. 2005, Plasma Physics Via Computer
  Simulation (Taylor & Francis Group, 2005)

[Blandford & Ostriker(1978)]BlandfordOstriker78
Blandford, R. D. & Ostriker, J. P. 1978, , 221, L29

[Blasi et al.(2012)]Blasi12 Blasi, P., Amato, E., & Serpico, P. D. 2012, Physical Review Letters, 109, 061101 


[Booth et al.(2013)Booth, Agertz, Kravtsov, &
  Gnedin]Booth_etal13
Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013,
  , 777, L16

[Boris(1970)]Boris70
Boris, J. P. 1970, in Proceedings of the Fourth Conference on Numerical
  Simulation Plasmas (Navel Research Laboratory, Washington, D.C.), pp. 3–67

[Breitschwerdt et al.(1991)Breitschwerdt, McKenzie, &
  Voelk]Breitschwerdt_etal91
Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, , 245, 79

[Breitschwerdt et al.(1993)]Breitschwerdt93 Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1993, , 269, 54 



[Denton & Kotschenreuther(1995)]DentonKotschenreuther95
Denton, R. E. & Kotschenreuther, M. 1995, Journal of Computational
  Physics, 119, 283

[Dimits & Lee(1993)]DimitsLee93
Dimits, A. M. & Lee, W. W. 1993, Journal of Computational Physics, 107,
  309
  
[Draine(2011)]Draine11 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium by Bruce T. Draine. Princeton University Press, 2011. ISBN: 978-0-691-12214-4,  



[Drury(1983)]Drury83
Drury, L. O. 1983, Reports on Progress in Physics, 46, 973

[Dupree(1966)]Dupree66
Dupree, T. H. 1966, Physics of Fluids, 9, 1773

[Enßlin et al.(2007)Enßlin, Pfrommer, Springel, &
  Jubelgas]Enblin_etal07
Enßlin, T. A., Pfrommer, C., Springel, V., & Jubelgas, M. 2007,
  , 473, 41
  
  
[Enßlin et al.(2011)]Ensslin11 Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, , 527, A99 
  

[Everett et al.(2008)Everett, Zweibel, Benjamin, McCammon,
  Rocks, & Gallagher]Everett_etal08
Everett, J. E., Zweibel, E. G., Benjamin, R. A., McCammon, D., Rocks,
  L., & Gallagher, III, J. S. 2008, , 674, 258

[Farber et al.(2018)]Farber18 Farber, R., Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. G. 2018, , 856, 112 

[Farmer & Goldreich(2004)]FarmerGoldreich04
Farmer, A. J. & Goldreich, P. 2004, , 604, 671

[Felice & Kulsrud(2001)]FeliceKulsrud01
Felice, G. M. & Kulsrud, R. M. 2001, , 553, 198

[Ferrière(2001)]Ferriere01 Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031 

[Fujita & Ohira(2011)]Fujita11 Fujita, Y., & Ohira, Y. 2011, , 738, 182 

[Gardiner & Stone(2005)]GardinerStone05
Gardiner, T. A. & Stone, J. M. 2005, Journal of Computational Physics,
  205, 509

[Gardiner & Stone(2008)]GardinerStone08
—. 2008, Journal of Computational Physics, 227, 4123

[Ginzburg & Syrovatskii(1964)]GinzburgSyrovatskii64
Ginzburg, V. L. & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays

[Girichidis et al.(2016)Girichidis, Naab, Walch, Hanasz,
  Mac Low, Ostriker, Gatto, Peters, Wünsch, Glover, Klessen,
  Clark, & Baczynski]Girichidis_etal16
Girichidis, P., Naab, T., Walch, S., Hanasz, M., Mac Low, M.-M.,
  Ostriker, J. P., Gatto, A., Peters, T., Wünsch, R., Glover,
  S. C. O., Klessen, R. S., Clark, P. C., & Baczynski, C. 2016, ,
  816, L19

[Girichidis et al.(2018)]Girichidis18 Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, , 479, 3042 


[Goldreich & Sridhar(1995)]GoldreichSridhar95
Goldreich, P. & Sridhar, S. 1995, , 438, 763

[Guo & Oh(2008)]GuoOh08
Guo, F. & Oh, S. P. 2008, , 384, 251

[Grenier et al.(2015)]Grenier15 Grenier, I. A., Black, J. H., & Strong, A. W. 2015, , 53, 199 


[Hanasz et al.(2013)Hanasz, Lesch, Naab, Gawryszczak,
  Kowalik, & Wóltański]Hanasz_etal13
Hanasz, M., Lesch, H., Naab, T., Gawryszczak, A., Kowalik, K., &
  Wóltański, D. 2013, , 777, L38

[Holcomb & Spitkovsky(2018)]HolcombSpitkovsky19
Holcomb, C. & Spitkovsky, A. 2018, arXiv e-prints

[Holguin et al.(2018)]Holguin18 Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H.-Y. K. 2018, arXiv:1807.05494 


[Hu & Krommes(1994)]HuKrommes94
Hu, G. & Krommes, J. A. 1994, Physics of Plasmas, 1, 863
[Ipavich(1975)]Ipavich75 Ipavich, F. M. 1975, , 196, 107 

[Jacob & Pfrommer(2017)]Jacob17 Jacob, S., & Pfrommer, C. 2017, , 467, 1449 


[Jacob et al.(2018)]Jacob18 Jacob, S., Pakmor, R., Simpson, C. M., Springel, V., & Pfrommer, C. 2018, , 475, 570 


[Jiang & Oh(2018)]JiangOh18
Jiang, Y.-F. & Oh, S. P. 2018, , 854, 5

[Jokipii(1966)]Jokipii66
Jokipii, J. R. 1966, , 146, 480

[Jubelgas et al.(2008)Jubelgas, Springel, Enßlin, &
  Pfrommer]Jubelgas_etal08
Jubelgas, M., Springel, V., Enßlin, T., & Pfrommer, C. 2008, ,
  481, 33

[Krymskii(1977)]Krymsky77 Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, 234, 1306

[Kulsrud & Pearce(1969)]KulsrudPearce69
Kulsrud, R. & Pearce, W. P. 1969, , 156, 445

[Kulsrud(2005)]Kulsrud05
Kulsrud, R. M. 2005, Plasma physics for astrophysics (Princeton University
  Press)

[Kunz et al.(2014)Kunz, Stone, & Bai]Kunz_etal14
Kunz, M. W., Stone, J. M., & Bai, X.-N. 2014, Journal of Computational
  Physics, 259, 154
  
[Lazarian(2016)]Lazarian16
Lazarian, A. 2016, , 833, 131

[Lazarian & Beresnyak(2006)]LazarianBeresnyak06 Lazarian, A., & Beresnyak, A. 2006, , 373, 1195

[Lebiga et al.(2018)]Lebiga18 Lebiga, O., Santos-Lima, R., & Yan, H. 2018, , 476, 2779

[Lee & Völk(1973)]LeeVolk73
Lee, M. A. & Völk, H. J. 1973, , 24, 31

[Mao & Ostriker(2018)]Mao18 Mao, S. A., & Ostriker, E. C. 2018, , 854, 89 

[Naab & Ostriker(2017)]Naab17 Naab, T., & Ostriker, J. P. 2017, , 55, 59 

[Pakmor et al.(2016)Pakmor, Pfrommer, Simpson, &
  Springel]Pakmor_etal16
Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ,
  824, L30

[Parker(1966)]Parker66
Parker, E. N. 1966, , 145, 811

[Parker(1992)]Parker92
—. 1992, , 401, 137

[Parker & Lee(1993)]ParkerLee93
Parker, S. E. & Lee, W. W. 1993, Physics of Fluids B, 5, 77





[Pfrommer et al.(2017)Pfrommer, Pakmor, Schaal, Simpson,
  & Springel]Pfrommer_etal17
Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V.
  2017, , 465, 4500

[Recchia et al.(2016)Recchia, Blasi, &
  Morlino]Recchia_etal16
Recchia, S., Blasi, P., & Morlino, G. 2016, , 462, 4227

[Roe(1981)]Roe81
Roe, P. L. 1981, Journal of Computational Physics, 43, 357

[Ruszkowski et al.(2017)Ruszkowski, Yang, &
  Zweibel]Ruszkowski_etal17
Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, , 834, 208

[Salem & Bryan(2014)]SalemBryan14
Salem, M. & Bryan, G. L. 2014, , 437, 3312

[Salem et al.(2014)]Salem_etal14 Salem, M., Bryan, G. L., & Hummels, C. 2014, , 797, L18 


[Schlickeiser(1989)]Schlickeiser89
Schlickeiser, R. 1989, , 336, 264

[Schlickeiser & Miller(1998)]SchlickeiserMiller98
Schlickeiser, R. & Miller, J. A. 1998, , 492, 352

[Sharma et al.(2010)]Sharma10 Sharma, P., Colella, P., & Martin, D. F. 2010, SIAM J Sci. Comp. 32, 3564 

[Simpson et al.(2016)]Simpson16 Simpson, C. M., Pakmor, R., Marinacci, F., et al. 2016, , 827, L29 



[Skilling(1975a)]Skilling75a
Skilling, J. 1975a, , 172, 557

[Skilling(1975b)]Skilling75b
—. 1975b, , 173, 245

[Skilling(1975c)]Skilling75c
—. 1975c, , 173, 255



[Soler et al.(2016)Soler, Terradas, Oliver, &
  Ballester]Soler_etal16
Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2016, ,
  592, A28
  
[Stone et al.(2008)Stone, Gardiner, Teuben, Hawley, &
  Simon]Stone_etal08
Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon,
  J. B. 2008, , 178, 137

[Summers & Thorne(1991)]SummersThorne91
Summers, D. & Thorne, R. M. 1991, Physics of Fluids B, 3, 1835

[Thomas & Pfrommer(2018)]ThomasPfrommer19
Thomas, T. & Pfrommer, C. 2018, arXiv e-prints

[Uhlig et al.(2012)Uhlig, Pfrommer, Sharma, Nath,
  Enßlin, & Springel]Uhlig_etal12
Uhlig, M., Pfrommer, C., Sharma, M., Nath, B. B., Enßlin, T. A.,
  & Springel, V. 2012, , 423, 2374

[Völk(1973)]Volk73
Völk, H. J. 1973, , 25, 471

[Wentzel(1974)]Wentzel74 Wentzel, D. G. 1974, , 12, 71 

[Wiener et al.(2013)Wiener, Oh, & Guo]Wiener_etal13
Wiener, J., Oh, S. P., & Guo, F. 2013, , 434, 2209

[Wiener et al.(2017)Wiener, Pfrommer, & Oh]Wiener_etal17
Wiener, J., Pfrommer, C., & Oh, S. P. 2017, , 467, 906

[Wiener et al.(2018)]Wiener18a Wiener, J., Zweibel, E. G., & Oh, S. P. 2018, , 473, 3095 

[Wiener & Zweibel(2018)]Wiener18b Wiener, J., & Zweibel, E. G. 2018, arXiv:1812.02179 


[Yan & Lazarian(2002)]YanLazarian02
Yan, H. & Lazarian, A. 2002, Physical Review Letters, 89, 281102

[Yan & Lazarian(2008)]YanLazarian08
—. 2008, , 673, 942

[Zirakashvili et al.(1996)Zirakashvili, Breitschwerdt,
  Ptuskin, & Voelk]Zirakashvili_etal96
Zirakashvili, V. N., Breitschwerdt, D., Ptuskin, V. S., & Voelk, H. J.
  1996, , 311, 113

[Zweibel(2013)]Zweibel13
Zweibel, E. G. 2013, Physics of Plasmas, 20, 055501

[Zweibel(2017)]Zweibel17 Zweibel, E. G. 2017, Physics of Plasmas, 24, 055402 







