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

# Influence of Ion-Neutral Damping on the Cosmic-Ray Streaming Instability: Magnetohydrodynamic Particle-in-cell Simulations

Illya Plotnikov¹ Ⓓ, Eve C. Ostriker² Ⓓ, and Xue-Ning Bai³,⁴ Ⓓ

¹IRAP, Université de Toulouse III-Paul Sabatier, OMP, Toulouse, France; illya plotnikov@irap.omp.eu

²Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton, NJ 08544, USA; eco@astro.princeton.edu

³Institute for Advanced Study, Tsinghua University, Beijing 100084, People's Republic of China

⁴Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China; xbai@mail.tsinghua.edu.cn

Received 2021 February 22; revised 2021 April 12; accepted 2021 April 12; published 2021 June 8

## Abstract

We explore the physics of the gyroresonant cosmic-ray streaming instability (CRSI) including the effects of ion-neutral (IN) damping. This is the main damping mechanism in (partially ionized) atomic and molecular gas, which are the primary components of the interstellar medium (ISM) by mass. Limitation of CRSI by IN damping is important in setting the amplitude of Alfvén waves that scatter cosmic rays (CRs) and control galactic-scale transport. Our study employs the magnetohydrodynamic (MHD)–particle-in-cell hybrid fluid-kinetic numerical technique to follow linear growth as well as post-linear and saturation phases. During the linear phase of the instability—where simulations and analytical theory are in good agreement—IN damping prevents wave growth at small and large wavelengths, with the unstable bandwidth lower for higher IN collision rates $v_{\text{in}}$. Purely MHD effects during the post-linear phase extend the wave spectrum toward larger $k$. In the saturated state, the CR distribution evolves toward greater isotropy (lower streaming velocity) by scattering off of Alfvén waves excited by the instability. In the absence of low-$k$ waves, CRs with sufficiently high momentum are not isotropized. The maximum wave amplitude and rate of isotropization of the distribution function decrease at higher $v_{\text{in}}$. When the IN damping rate approaches the maximum growth rate of CRSI, wave growth and isotropization are suppressed. Implications of our results for CR transport in partially ionized ISM phases are discussed.

*Unified Astronomy Thesaurus concepts: Cosmic rays (329); Magnetic fields (994); Computational methods (1965)*

## 1. Introduction

Cosmic rays (CRs) are a key component of the interstellar medium (ISM). They are potentially important to ISM dynamics, including support against gravity, as their energy density is typically comparable to that of the magnetic field as well as the thermal and turbulent energy of the ISM gas (e.g., Spitzer 1978; Ferrière 2001; Draine 2011; Grenier et al. 2015). Perhaps even more important dynamically, given the large scale height of the CR distribution, is the potential for CRs to contribute in driving galactic winds (e.g., Ipavich 1975; Breitschwerdt et al. 1991; Zirakashvili et al. 1996; Everett et al. 2008; Mao & Ostriker 2018). CRs are also crucial to ISM microphysics, driving ionization and dissociation as well as providing the primary heating of gas in regions shielded against UV (e.g., Draine 2011; Glassgold et al. 2012).

Coupling between CRs and thermal ISM gas occurs through the scattering of CRs off of magnetic waves, preexisting in the turbulent ISM or self-generated by streaming CRs (Kulsrud & Pearce 1969; Wentzel 1969; Skilling 1971; Kulsrud 2005; Amato & Blasi 2018). Higher-energy CRs may interact primarily with externally generated waves (Blasi et al. 2012a). However, particles of energy at or below the gigaelectronvolt level dominate by number and energy the overall CR content. At microparsec scales comparable to the gyroradius of gigaelectronvolt protons, the energy density of ISM turbulence is too low to provide efficient scattering (under the assumption of a cascade from directly observed turbulence at larger scales), and these particles are believed to be mainly scattered by self-excited waves. Alfvén waves generated by streaming therefore bear the primary responsibility for

scattering the dominant portion of the CR distribution in our own and other galaxies (e.g., Zweibel 2013).

There is an increasing necessity to provide satisfactory microphysical understanding of the CR streaming instability (CRSI) and its implications for CR transport in different phases of the ISM. The possible dynamical role of CRs in driving galactic winds has recently led to reconsideration of the traditional assumption that streaming occurs at the Alfvén speed in several analytical studies (Wiener et al. 2013; Recchia et al. 2016; Zweibel 2017). In parallel, a number of groups have implemented CRs as a fluid in numerical magnetohydrodynamic (MHD) codes (e.g., Yang et al. 2012; Dubois & Commerçon 2016; Pakmor et al. 2016; Pfrommer et al. 2017; Jiang & Oh 2018; Thomas & Pfrommer 2019; Hopkins et al. 2020) in order to study the effect of CRs on galactic wind generation (e.g., Hanasz et al. 2013; Girichidis et al. 2016, 2018; Pfrommer et al. 2017; Ruszkowski et al. 2017; Wiener et al. 2017; Butsky & Quinn 2018; Dashyan & Dubois 2020), on the global evolution of supernova remnants (Pais et al. 2018; Dubois et al. 2019), or on the multiphase ISM structure (Bustard & Zweibel 2020). Global simulations testing different models of CR transport coefficients were recently conducted by Hopkins et al. (2021a, 2021b), suggesting that standard models of CR transport are inconsistent with observed constraints on large scales. At the same time, observations that probe ionization and chemistry in ISM clouds at small scales place constraints on CR transport (see review of Padovani et al. 2020, and references therein), e.g., suggesting diffusive rather than free-streaming behavior in the outer layers of molecular clouds (Silsbee & Ivlev 2019).

Fluid implementations of CRs rely on subgrid prescriptions for diffusion and streaming. To date, these subgrid treatments
have relied on empirical estimates of diffusivities, idealized
streaming treatments, or simple models based on analytic
predictions for growth and dissipation of waves. Considering
the important role of microphysics in CR transport, it is
valuable to pursue a deeper investigation that addresses the
processes of wave growth, damping, and particle-wave
interactions directly. The most important parameter to be
obtained from microphysical studies of CR-ISM interactions is
the scattering rate of CRs off of Alfvén waves, since this
influences the determination of diffusivity in traditional
treatments (e.g., Skilling 1971; Hanasz & Lesch 2003) or the
rate of change of the CR flux in two-moment methods (e.g.,
Jiang & Oh 2018; Thomas & Pfrommer 2019).

The energy-dependent scattering frequency of CRs, $\nu_s$, is proportional to the wave energy density in the resonant range. The latter quantity depends on the outcome of the CRSI in the local medium. There are two different classes of the CRSI, resonant (Lerche 1967; Kulsrud & Pearce 1969; Skilling 1971, 1975; Ginzburg et al. 1973; Wentzel 1974; Berezinskii et al. 1990) and nonresonant (Bell 2004; Amato & Blasi 2009; Bykov et al. 2013). The latter is important only in the near environment around accelerators (such as supernova remnant shocks), where the electrical current of CRs is large, i.e., $(4\pi/c)J_{CR}R_{L,0}/B_0 \gg 1$, where $J_{CR} = en_{CR}V_D$ is the electrical current of CRs streaming at speed $V_D$ and with a typical gyroradius $R_{L,0}$ (Amato & Blasi 2009). In the general ISM, the resonant interaction prevails since the density of CRs is very low ($n_{CR} \sim 10^{-10}-10^{-8} \text{cm}^{-3}$). Because we are interested in ambient ISM conditions, here we focus exclusively on the resonant instability.

Wave growth, triggered by the anisotropy of the distribution
function (DF) or CR density gradient, is generally in
competition with different damping mechanisms depending
on the phase of the ISM. In a partially ionized medium the most
important wave damping mechanism is from collisions
between ions and neutrals. The competition of CRSI with
ion-neutral (IN) damping has been studied in the context of
molecular clouds (e.g., Kulsrud & Cesarsky 1971; Zweibel &
Shull 1982; Everett & Zweibel 2011; Morlino & Gabici 2015;
Ivlev et al. 2018) and in the context of particle acceleration at
shocks (e.g., O'C Drury et al. 1996; Bykov & Toptygin 2005;
Reville et al. 2007; Blasi et al. 2012b; Nava et al. 2016;
Brahimi et al. 2020).

Of course, damping mechanisms other than IN collisions can
be important in hotter and more diffuse phases of the ISM.
Other mechanisms that have been discussed include nonlinear
Landau damping (Lee & Völk 1973; Kulsrud 1978) and
damping by interactions with turbulence (Farmer & Gold-
reich 2004; Yan & Lazarian 2011; Lazarian 2016). We also
note the recently studied damping by charged dust grains that
can be important for CRs with energies <300 GeV (Squire
et al. 2021). As discussed in Nava et al. (2016), Xu et al.
(2016), and Brahimi et al. (2020), in neutral-dominated phases
of the ISM, where the ionization fraction is small, IN damping
is dominant for waves that are resonant with CRs of energies
$E < 10$ TeV. Here, we consider only these primarily neutral
phases of the ISM; investigation of effects of alternative
damping mechanisms is deferred to future work.

In the present work, we investigate the interplay between
CRSI and IN damping during the linear phase of the CRSI, and
also assess the behavior of the system in post-linear and late-
time saturated stages of evolution. We are motivated to obtain

better insight into CR transport in ISM conditions that are
largely neutral (neutrals make up most of the ISM mass), where
the IN fraction is controlled by photoionization or ionization by
low-energy CRs.

Our study of CRSI relies on the MHD particle-in-cell (PIC)
numerical method described in Bai et al. (2015) (see also
Reville & Bell 2012; Amano 2018; Lebiga et al. 2018;
Mignone et al. 2018; van Marle et al. 2018). This method is
well adapted to capture the resonant nature of the CR-fluid
interaction while allowing for a relatively large space and
timescale evolution of the system and the fact that $n_{CR}/n_i$ is
extremely low in the parameter range of interest. The
gyroradius of CRs must be resolved, but smaller scales (such
as the ion skin depth and electron scales) do not need to be
resolved. We note that several recent studies have alternatively
investigated certain aspects of the CRSI with fully kinetic
(Holcomb & Spitkovsky 2019; Shalaby et al. 2021) or hybrid
kinetic approaches (Haggerty et al. 2019; Schroer et al. 2020)
that use PIC methods. These approaches have the advantage of
resolving small-scale phenomena on electron (full PIC) or
proton skin depth (hybrid PIC) scales, but are in practice
limited in the range of $n_{CR}/n_i$ and other parameters that can be
studied because the thermal ions are treated via PIC.

In this work, we follow up on Bai et al. (2019, hereafter
Paper I), where the CRSI was studied using the MHD-PIC
approach. There, the main findings include the following: (i)
the linear phase of the instability for both the resonant and
nonresonant branches can be accurately reproduced with our $\delta f$
MHD-PIC method; (ii) the quasi-linear diffusion (QLD)
formalism accurately describes temporal changes in the CR
DF, except near a pitch angle of 90° (the "$\mu=0$ crossing"
problem); (iii) crossing of $\mu=0$ is mainly due to nonlinear
wave-particle interactions; and (iv) the Alfvén wave amplitude
in saturation reflects the expected transfer of net momentum
from the originally anisotropic CRs to forward-propagating
Alfvén waves. Here, we adopt the same numerical methods and
conduct simulations in a similar parameter regime to Paper I,
but now we additionally consider the effects of IN interactions
that damp waves, competing with CR-induced wave genera-
tion. In a separate work, Bambic et al. (2021) have also used
MHD-PIC to investigate CR propagation through the multi-
phase ISM. That study extends Paper I and the present work by
simulating a drifting CR population across a two-component
ISM (with and without damping, modeling neutral and ionized
regions).

This article is organized as follows. In Section 2, we review
the analytical derivation of the CRSI linear growth rate,
including the contribution from IN damping. In Section 3 we
outline our numerical methods, including the implementation
of IN wave damping. The simulation results are presented in
Section 4. We further discuss our findings and their
astrophysical context in Section 5, and summarize our main
conclusions in Section 6.

## 2. CRSI Including IN Interaction

In this section, we discuss the linear growth rate of the CRSI
without and with IN damping. While analytic studies of CRSI
commonly adopt a power-law distribution, here we will instead
consider a $\kappa$-distribution, as adopted in our simulations. This
choice is motivated by the use of the $\delta f$ method, in which the
unperturbed part of the CR distribution must be specified for all
$p$, and should be a smooth function (see Bai et al. 2019).
Because standard expressions for CRSI with IN collisions are based on power-law distributions, here we re-derive growth rates for the case of the $\kappa$-distribution. These analytic results are useful for comparison with the simulations that follow.

Some discussion of the high-frequency regime is given in the last part of this section.

## 2.1. No Damping

### 2.1.1. Resonant Contribution Only

We consider the instability arising from resonant interaction in the low-frequency, $\omega \ll \Omega_0$, and nonrelativistic drift ($V_D \ll c$) limit.<fnref n="5" /> Here, $\omega = kV_A$ is the Alfvén wave frequency and $\Omega_0 = eB_0/(m_p c)$ is the nonrelativistic proton gyrofrequency. We adopt standard notation, where Alfvén velocity $V_A = B_0/(4\pi\rho_i)^{1/2}$ ($B_0$ is the local magnetic field strength and $\rho_i$ is the background ion density), $e$ is the proton charge, $m_p$ is the proton mass, and $c$ is the speed of light. We also assume that the mass and charge of the background fluid ions are identical to those of CR particles<fnref n="6" />: $m_i = m_p$ and $q_i = e$.

We start with the growth rate expression given by Kulsrud (2005; Equation (69) of Ch. 12) for wavenumber $k$:

$$ \Gamma_{\text{CR}}(k) = -\pi^2 e^2 \left(\frac{V_A}{c}\right)^2 \left(\frac{V_D}{V_A} - 1\right) \times \int \frac{\partial F}{\partial p} \frac{p_{\perp}^2}{p} \delta(kp_{\parallel} - m_p \Omega_0) d^3p, \quad (1) $$

where $V_D$ is the initial drift speed between the background gas and the frame in which the CRs have an isotropic distribution (taken along the external magnetic field in the x-direction). The $\delta$ function in the integral accounts for the gyroresonant interaction at the fundamental harmonic, where resonance occurs at $p_{\parallel} = m_p \Omega_0/k \equiv p_{\text{res}}$.

For the (isotropic) distribution $F(p)$ of CRs in the drift frame, we adopt a $\kappa$-distribution:

$$ F(p) = A \left[ 1 + \frac{1}{\kappa} \left( \frac{p}{p_0} \right)^2 \right]^{-(\kappa+1)} . \quad (2) $$

The normalization factor $A$ is given by

$$ A = \frac{n_{\text{CR}} \Gamma(\kappa + 1)}{(\pi \kappa p_0^2)^{3/2} \Gamma(\kappa - 1/2)}, \quad (3) $$

where $n_{\text{CR}}$ is the number density of CRs, and $\Gamma(x)$ is the Euler $\Gamma$-function. We generally adopt $\kappa = 1.25$, corresponding to $F(p) \propto p^{-4.5}$ for $p \gg p_0$. We note that a $\kappa$-distribution reduces to a Maxwellian in the limit $\kappa \to \infty$.

Equation (1) becomes

$$ \Gamma_{\text{CR}} = 2\pi^2 e^2 A \frac{\kappa + 1}{\kappa p_0^2} \left(\frac{V_A}{c}\right)^2 \left(\frac{V_D}{V_A} - 1\right) \times \frac{\int \left[1 + \frac{1}{\kappa}\left(\frac{p}{p_0}\right)^2\right]^{-(\kappa+2)} \delta(kp_{\parallel} - m_p \Omega_0) p_{\perp}^2 d^3p}{I} \quad (4) $$

With $d^3p = d\theta dp_{\parallel} dp_{\perp}$, the integral $I$ in the previous equation becomes

$$ I = \frac{2\pi}{k} \int_{0}^{\infty} \left[ 1 + \frac{1}{\kappa} \frac{p_{\text{res}}^2 + p_{\perp}^2}{p_0^2} \right]^{-(\kappa+2)} p_{\perp}^3 dp_{\perp}. \quad (5) $$

Noting that $p_{\perp}^2 = p^2 - p_{\text{res}}^2$ and that $p_{\perp} dp_{\perp} = pdp$, we change the variable of integration from $p_{\perp}$ to $p$. This leads to

$$ \begin{aligned} I &= \frac{2\pi}{k} \int_{p_{\text{res}}}^{\infty} \left[ 1 + \frac{1}{\kappa} \frac{p^2}{p_0^2} \right]^{-(\kappa+2)} p(p^2 - p_{\text{res}}^2) dp \\ &= \frac{\pi}{k} p_0^4 \frac{\kappa}{\kappa + 1} \left[ 1 + \frac{1}{\kappa} \frac{p_{\text{res}}^2}{p_0^2} \right]^{-\kappa}. \end{aligned} \quad (6) $$

Using $\rho_i = m_p n_i$, $R_{L,0} \equiv p_0/(m_p \Omega_0)$, and $p_{\text{res}}/p_0 = 1/(kR_{L,0})$ and inserting the expression from Equation (3) the growth rate can be written as

$$ \Gamma_{\text{CR}}(k) = \frac{\pi^{1/2} \Gamma(\kappa + 1) n_{\text{CR}} \Omega_0}{2\kappa^{3/2} \Gamma(\kappa - 1/2) n_i k R_{L,0}} \times \left(\frac{V_D}{V_A} - 1\right) \left[1 + \frac{1}{\kappa(kR_{L,0})^2}\right]^{-\kappa}. \quad (7) $$

The growth rate is maximal at $kR_{L,0} = \sqrt{2 - 1/\kappa}$. For the adopted value in this study $\kappa = 1.25$, it is

$$ \Gamma_{\text{max},0} \approx 0.28 \Omega_0 \frac{n_{\text{CR}}}{n_i} \left( \frac{V_D}{V_A} - 1 \right). \quad (8) $$

The subscript “max, 0” stands for the maximum growth rate of the instability without IN damping. The numerical prefactor varies from 0.25 for $\kappa = 1$ to 0.38 for $\kappa \to \infty$. In the long-wavelength limit ($kR_{L,0} \ll 1$) the growth rate varies as $\propto k^{2\kappa-1}$, and in the short-wavelength limit ($kR_{L,0} \gg 1$) it varies as $\propto k^{-1}$.

### 2.1.2. Full Dispersion Relation

In more general cases when CR-induced current can be large, the contribution from nonresonant instability can be significant (see, e.g., Bell 2004; Amato & Blasi 2009; Bykov et al. 2013). The full dispersion relation using the $\kappa$-distribution was derived in Paper I. Here, we reproduce it for completeness:

$$ \omega^2 = k^2 V_A^2 \mp \frac{n_{\text{CR}}}{n_i} \Omega_0 (\omega - k V_D) \left( (1 - Q_1) \pm i Q_2 \right); \quad (9) $$

here the $(1 - Q_1)$ and $Q_2$ terms are due to nonresonant and resonant responses of CRs to waves, respectively (see Equation (9) in Paper I for the definition of these terms). The resonant term given in Equation (40) of Paper I is

$$ Q_2 = \frac{\sqrt{\pi}}{\kappa^{3/2}} \frac{\Gamma(\kappa + 1)}{\Gamma(\kappa - 1/2)} \frac{1}{k R_{L,0}} \left[ 1 + \frac{1}{(\kappa k^2 R_{L,0}^2)} \right]^{-\kappa}. \quad (10) $$

<footnote n="5">⁵ In this study we adopt the notation $\Omega_0$ for the quantity denoted by $\Omega_c$ in Paper I, i.e., $\Omega_0 = \Omega_c$.</footnote>

<footnote n="6">⁶ If this condition is not satisfied, one has to replace $\Omega_0$ by $\Omega_i = q_i B_0/(m_i c)$ in the final expression of the growth rate, i.e., in Equations (7), (8), (9), (11), and (23). We note, however, that if mass densities are used instead of number densities then the characteristic frequency appearing in the growth rate would correspond to the gyrofrequency of CRs, owing to the equality $(n_{\text{CR}}/n_i)\Omega_i = (\rho_{\text{CR}}/\rho_i)\Omega_0$.</footnote>
In the case with negligible contribution from the nonresonant term ((1 − Q$_1$) = 0), the imaginary part of ω obtained by solving Equation (9) is then

$$
\omega_I = \frac{1}{2} \frac{n_{CR}}{n_i} \Omega_0 \left( \frac{V_D}{V_A} - 1 \right) Q_2, \quad (11)
$$

which is identical to the growth rate given in Equation (7).

2.2. With IN Damping

2.2.1. Damping Term from Fluid Equations

In this section, we discuss the Alfvén wave propagation
properties in the presence of neutrals, temporarily ignoring the
destabilizing contribution from CR streaming. The derived
damping rate is then subtracted from the CRSI growth rate
given above to obtain the net growth rate.<fnref n="7" />

Consider coupled two-fluid momentum equations for ions and neutrals, without the contribution of CRs, where ohmic dissipation and the Hall effect are neglected:

$$
\rho_i \frac{Du_i}{Dt} = - \nabla P_i + \frac{J}{c} \times \mathbf{B} - \rho_i \nu_{\text{in}} (\mathbf{u}_i - \mathbf{u}_n) \quad (12)
$$

$$
\rho_n \frac{Du_n}{Dt} = - \nabla P_n - \rho_n \nu_{\text{ni}} (u_n - u_i), \quad (13)
$$

where ρᵢ (ρₙ), Pᵢ (Pₙ), and νᵢₙ (νₙᵢ) are the ion (neutral) mass density, the ion (neutral) pressure, and the frequency of collision of ions with neutrals (neutrals with ions), respectively. The quantities **B** and **J** correspond to the magnetic field and current carried by the plasma. Ions are subject to both Lorentz force and friction with neutrals, while neutrals are only subject to friction with ions.

From momentum conservation, $\rho_i\nu_{in} = \rho_n\nu_{ni} = n_n n_i \mu \langle \sigma v \rangle$, where $\mu = m_i m_n / (m_i + m_n)$ is the reduced mass, and $\langle \sigma v \rangle$ is the momentum transfer rate coefficient. Most relevant to wave damping is the IN collision frequency

$$
\nu_{\text{in}} = \frac{\rho_n}{m_i + m_n} \langle \sigma v \rangle . \qquad (14)
$$

Draine (2011; see Table 2.1) provides values of the coefficient $\langle \sigma v \rangle \approx 3.3 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}$ for neutral atomic gas, where the main collision partner is H with $\mu \approx 0.5 m_p$ for H$^+$ ions, and $\langle \sigma v \rangle \approx 1.9 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}$ in molecular regions, where the collision partner is H$_2$ with $\mu \approx 2 m_p$ for C$^+$ ions in diffuse gas or HCO$^+$ ions in dense gas. The magnetic field in Equations (12)–(13) evolves subject to the induction equation:

$$
\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u}_i \times \mathbf{B}). \quad (15)
$$

For simplicity, we consider a static, homogeneous, and
incompressible background medium. The mean magnetic field
is oriented along the x-direction, **B**₀ = **B**₀**e**ₓ, and the perturba-
tion components are perpendicular to **B**₀: δ**B** = δ**B**ᵧ**e**ᵧ + δ**B**z**e**z.
The velocity and magnetic field perturbations are
∝exp[i(kx − ωt)]. Linearizing Equations (12), (13), and (15)
leads ions to

$$
\omega^2 \delta u_i = k^2 V_{A,i}^2 \delta u_i - i \omega \nu_{\text{in}} (\delta u_i - \delta u_n), \quad (16)
$$

and neutrals to

$$
\omega^2 \delta u_n = -i \omega \nu_{ni} (\delta u_n - \delta u_i). \quad (17)
$$

Here, the Alfvén velocity includes the ion density only, i.e., $V_{A,i} = B_0 / \sqrt{4\pi\rho_i}$. Combining these two equations gives the dispersion relation:

$$
\omega^2 \left( 1 + \frac{i \nu_{\text{in}}}{\omega + i \nu_{\text{in}} \rho_i / \rho_n} \right) = k^2 V_{A,i}^2. \quad (18)
$$

One identifies here the Alfvén wave dispersion relation (in the ion fluid) modified by the presence of neutrals (second term inside brackets on the left-hand side). Equation (18) is a cubic equation for ω. Although available, general solutions are quite cumbersome; hence we do not provide them here. More physical insight can be gained by considering different limits, previously discussed in the literature (e.g., Kulsrud & Pearce 1969; Zweibel & Shull 1982; Tagger et al. 1995; Nava et al. 2016; Xu et al. 2016).

Here we reproduce the limits for completeness. We define
Z = ρᵢ/ρₙ as the ratio of ionized to neutral mass density. It is
related to the ionization fraction as xᵢ ≡ nᵢ/nₙ = Z/q, where
q = mᵢ/mₙ. Three limiting cases are then well defined:

1. Ion-dominated, Z ≫ 1. Here, the Alfvén wave velocity is unmodified by neutrals, ω = kV<sub>A,i</sub>. The damping rate is, for any k and ω,

$$
\Gamma_d = \frac{\nu_{\text{in}}}{2} \frac{k^2 V_{A,i}^2}{k^2 V_{A,i}^2 + \nu_{\text{in}}^2 Z^2}. \quad (19)
$$

For ω<sub>k</sub> = kV<sub>A,i</sub> ≫ ν<sub>ni</sub> = ν<sub>in</sub>Z it is equal to ν<sub>in</sub>/2.

2. Neutral-dominated (Z ≪ 1) and low-frequency (kV<sub>A,i</sub>/ν<sub>in</sub> ≪ √Z). The fluids are strongly coupled; hence the Alfvén wave velocity includes the sum of the ion and neutral densities: ω = kV<sub>A,tot</sub> = kB₀/√(4π(ρᵢ + ρₙ)). The damping rate is

$$
\Gamma_d \approx \frac{k^2 V_{A,i}^2}{2\nu_{\text{in}}}. \qquad (20)
$$

3. Neutral-dominated (Z ≪ 1) and high-frequency (kV<sub>A,i</sub>/ν<sub>in</sub> ≫ 1). Ions and neutrals are decoupled. The Alfvén wave velocity is unmodified by neutrals, ω = kV<sub>A,i</sub>. The damping rate is

$$
\Gamma_d \approx \frac{\nu_{\text{in}}}{2}. \qquad (21)
$$

In neutral-dominated media there is a range of frequencies where there is no propagation (the solution of the dispersion relation is purely imaginary): $2\sqrt{Z} < kV_{A,i}/\nu_{in} < 1/2$.

Figure 1 presents the numerical solutions of Equation (18) and approximate analytical solutions for the ion-dominated case (top panel) and neutral-dominated case (bottom panel).

In the present study, we will adopt a one-fluid numerical approach for the ionized fluid only. In this case, the dispersion relation (ignoring the dynamics of neutrals) is

$$
\omega^2 + i\nu_{\text{in}}\omega = k^2 V_{A,i}^2 \qquad (22)
$$

This implicitly means that we will study the regime where
neutrals are fully decoupled from ions, that is, the case in the
high-frequency limit where the Alfvén wave is not modified by

<footnote n="7">⁷ Such a procedure is typically adopted in the literature (Kulsrud & Cesarsky 1971; Ginzburg et al. 1973; Zweibel & Shull 1982; Nava et al. 2016).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Dispersion relation of Alfvén waves and damping rate in the presence of neutrals. The two panels present the real and imaginary parts of the wave frequency (blue and red solid lines, respectively) in the ion-dominated case (top panel; the analytical and numerical solutions are identical) and in the neutral-dominated case (bottom panel). The black dotted line in the bottom panel corresponds to an approximate function that fits the damping rate reasonably well for any *k* where a solution exists.</figcaption></figure>

neutrals $\omega_R = kV_{A,i}$ and the damping rate is asymptotic at its highest value: $-\omega_I = \Gamma_d = \nu_{in}/2$ (see Figure 1). For realistic parameters in the neutral ISM, this limit generally applies (see Section 2.3).

### 2.2.2. General Dispersion Relation with CR Streaming Contribution

It is straightforward to generalize the dispersion relation of the CRSI with IN interactions (see, e.g., Reille et al. 2007, 2021). Following the procedure of Paper I but with the addition of IN momentum exchange terms in fluid equations, we obtain

$$ \omega^2 \left( 1 + \frac{i\nu_{in}}{\omega + i\nu_{in}\rho_i/\rho_n} \right) = k^2 V_{A,i}^2 \mp \frac{n_{CR}}{n_i} \Omega_0(\omega - kV_D) \times ((1 - Q_1) \pm iQ_2), \quad (23) $$

where $(1 - Q_1)$ and $Q_2$ are from the nonresonant and resonant responses of CRs to waves as discussed above, and $Q_2$ is given for a $\kappa$-DF in Equation (10).

Comparing this to Equation (9), the supplementary term due to IN interaction is the second term inside the parentheses on the left-hand side of Equation (23), while the second term on

<table><caption>**Table 1**

List of Typical Parameter Values for Different Phases of the ISM</caption><thead><tr><th>Phase</th><th>n<br>(cm<sup>-3</sup>)</th><th>m<sub>i</sub><br>(m<sub>p</sub>)</th><th>m<sub>n</sub><br>(m<sub>p</sub>)</th><th>x<sub>i</sub></th><th>&nu;<sub>in</sub><br>(&omega;<sub>max</sub>)</th><th>&zeta;</th></tr></thead><tbody><tr><td>WNM</td><td>0.3</td><td>1</td><td>1.4</td><td>10<sup>-2</sup></td><td>1.6 &times; 10<sup>-5</sup></td><td>15.5</td></tr><tr><td>CNM</td><td>30</td><td>12</td><td>1.4</td><td>10<sup>-4</sup></td><td>1.0 &times; 10<sup>-3</sup></td><td>0.07</td></tr><tr><td>MG</td><td>10<sup>2</sup></td><td>29</td><td>2.3</td><td>10<sup>-6</sup></td><td>3.9 &times; 10<sup>-4</sup></td><td>0.7</td></tr><tr><td>DMG</td><td>10<sup>4</sup></td><td>29</td><td>2.3</td><td>10<sup>-7</sup></td><td>1.2 &times; 10<sup>-1</sup></td><td>6.6 &times; 10<sup>-4</sup></td></tr></tbody></table>

**Notes.** WNM: warm neutral medium; CNM: cold neutral medium; MG: molecular gas; DMG: dense molecular gas. The reported quantities are the total density *n*, the ion mass *m*<sub>i</sub>, the ionization fraction *x*<sub>i</sub>, the corresponding *ω*<sub>HF</sub> above which the high-frequency regime is satisfied, and the drift-normalized instability parameter *ζ*, which is defined in Equation (25). We adopted *B*<sub>0</sub> = 5 μG, *n*<sub>CR</sub> = 10<sup>-9</sup> cm<sup>-3</sup>, and *p*<sub>0</sub> = *m*<sub>p</sub>*c*, corresponding to CR energy *E* ≃ GeV, as the fiducial values.

the right-hand side is due to CR streaming. For sufficiently low $n_{CR}/n_i$ (as is applicable in general), the solution of this dispersion relation is the same as if the damping (IN) and growing (CRSI) imaginary terms were derived separately, with the two terms added to obtain the net growth rate. The growth rate from resonant interaction only is presented in Section 2.1.1 and the damping rate in Section 2.2.1. Taking the high-frequency limit, the net growth rate is

$$ \Gamma_{\text{tot}}(k) = \Gamma_{\text{CR}}(k) - \frac{\nu_{\text{in}}}{2}, \quad (24) $$

where $\Gamma_{CR}(k)$ is given by Equation (7) for a $\kappa$-DF. This equation is valid for any $kV_{A,i} \gg \max[\nu_{in}, \nu_{ni}]$. Note that in the low-frequency (long-wavelength) regime, one needs to distinguish between high and low ionization fraction cases, which we do not consider here.

### 2.3. ISM Parameters and Implications

As discussed in Section 2.2.1, the one-fluid formulation does not capture all the details of wave damping. For this reason we discuss here the conditions for which the one-fluid approach is accurate.

We assume typical ISM values $B = 5$ μG and $n_{CR} = 10^{-9}$ cm<sup>-3</sup> and a typical energy of CRs around $E_0 \sim$ GeV. The frequency of the fastest-growing wave mode driven by gigaelectronvolt protons is defined as $\omega_{max} = 1.1V_{A,i}/R_{L,0}$, since the fastest-growing wavenumber is $k_{max} = 1.1R_{L,0}^{-1}$ (see text below Equation (8)). The typical values for the IN momentum exchange rate $\nu_{in}$ in different phases of the ISM were discussed previously in Section 2.2.1. The high-frequency regime is valid for $\omega_{max}/\nu_{in} \gg 1$. Comparing the two frequencies in different phases of the ISM (as done in Table 1, where the values for the different phases are taken from Draine 2011, Chap. 16) we find that this inequality is largely satisfied in all neutral-dominated phases.

Another important consideration concerns the possibility of wave growth in the presence of IN damping. CRSI-driven waves can only grow if the maximum growth rate of the CRSI is higher than the IN damping rate: $\Gamma_{max,0}/\Gamma_d > 1$. Otherwise, even the fastest-growing mode is damped. For this purpose, we define a $\zeta$ parameter as

$$ \zeta = \frac{\Gamma_{\text{max},0}}{\Gamma_{\text{d}}} \frac{1}{V_{\text{D}}/V_{\text{A,i}} - 1}. \quad (25) $$
Note that by dividing Equation (8) (or more generally Equation (1)) by the factor $V_D/V_{A,i} - 1$, $\zeta$ is independent of $V_D$. CRSI can be triggered at a range of wavelengths provided that $V_D/V_{A,i} - 1 > 1/\zeta$. Thus, if $\zeta > 1$, the CRSI can be triggered for a super-Alfvénic drift velocity $V_D/V_{A,i} - 1$ that is order-unity or lower. From Table 1, we see that $\zeta$ is typically in the range of 0.1–20 (except in the DMG).<fnref n="8" /> Thus in these environments, and considering only IN damping, the CRSI can be triggered for a streaming velocity of $0.1 \lesssim V_D/V_{A,i} - 1 \lesssim 10$. But by the same token, wave damping will significantly limit the range over which the instability can be directly excited if $V_D/V_{A,i} - 1$ is not large compared to $1/\zeta$.

## 2.4. Wave Saturation and Steady-state Streaming Velocity

When CRSI is efficiently triggered, the scattering of CRs off of self-generated waves gradually reduces the anisotropy of the CR DF, lowering the drift velocity and the growth rate. The limitation of CR transport by self-generated Alfvén waves is commonly referred to as self-confinement of CRs.

In the absence of wave damping, the DF anisotropy (in the wave frame) is eventually erased and the wave amplitude saturates. Considering that the momentum lost by CRs is transferred to waves, the saturation level can be estimated approximately as (e.g., Kulsrud 2005; Bai et al. 2019; Holcomb & Spitkovsky 2019)

$$ \left(\frac{\delta B}{B_0}\right)^2 \sim \frac{n_{\text{CR}}}{n_i} \left(\frac{V_D}{V_{A,i}} - 1\right). \qquad (26) $$

The wave intensity at saturation is linearly proportional to the number density of CRs and to the first-order anisotropy of the DF.

In the presence of IN damping (or other damping), a commonly adopted assumption is that the streaming velocity is such that the linear growth rate of CRSI balances the wave damping rate (e.g., Kulsrud & Cesarsky 1971; Krumholz et al. 2020; Hopkins et al. 2021a, 2021b). If we consider the fastest-growing waves and equate to zero the net growth rate with IN damping, $\Gamma_{\text{max},0} - \nu_{\text{in}}/2$, the corresponding streaming velocity (for our fiducial $\kappa = 1.25$) is

$$ \frac{V_{\text{st}}}{V_{A,i}} - 1 \sim 1.8 \frac{n_i}{n_{\text{CR}}} \frac{\nu_{\text{in}}}{\Omega_i}. \qquad (27) $$

Here, the use of $\Omega_i$ instead of $\Omega_0$ in Section 2.1.1 is required because the ionized fluids have $m_i > m_p$ in many of the relevant environments (see Table 1).

The value in Equation (27) sets, in principle, a lower limit on CR streaming velocity when IN damping is present. However, even if waves near the fastest-growing wavenumber can grow, waves at much shorter or longer wavelengths would be damped. As a result, particles with very large (very small) momentum that are resonant with smaller (larger) $k$ would not experience much scattering. Thus, the momentum-weighted value of $V_{\text{st}}$ will always be larger than the one given in Equation (27).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Time evolution of the magnetic wave energy (top panel) and of the streaming velocity of bulk CRs (bottom panel) for different values of $\nu_{\text{in}}$ ranging from 0 (blue line) to $1.9\Gamma_{\text{max},0}$ (black line).</figcaption></figure>

We note that the conventional approach described above considers only linear growth and damping of waves. In reality, both linear growth and pitch angle scattering must be considered, as opposed to wave damping. For a given momentum, resonant waves must be excited and survive IN damping over the timescale required for pitch angle diffusion, in order for isotropization to occur. Quantifying this requires comparison of the isotropization time (after the linear phase), $1/\nu_s$, with the damping time, $1/\nu_{\text{in}}$.

In the absence of damping, the saturation amplitude is given by Equation (26), and since the scattering rate under QLD is $\nu_s \sim (\sim/8)\Omega(\delta B/B_0)^2$ (Kulsrud 2005), comparison to Equation (8) shows that $\nu_s \sim \Gamma_{\text{max},0}$. If we were to assume that full redistribution in $\mu$ requires a time $\sim \nu_s^{-1}$, the total time for isotropization (considering wave growth and scattering) would be at least $2\Gamma_{\text{max},0}^{-1}$; in practice the numerical results described below (see Figure 2) correspond to a prefactor $\sim 10$ for the isotropization time. With IN damping, the growth rate is reduced, and (as we shall show) the maximum level of $(\delta B/B_0)^2$ is also reduced, which slows scattering. One might therefore expect that for isotropization to be successful, the ratio $\Gamma_{\text{max},0}/(\nu_{\text{in}}/2)$ must be above some critical value that is larger (perhaps much larger) than 1, in order to allow for both wave growth and particle diffusion. We return to this issue in Section 4.4.

<footnote n="8">⁸ For the DMG phase, the required initial drift velocity of CRs is extremely high, $V_D > 1.5 \times 10^3 V_{A,i} = 0.3c$, implying that wave excitation by CRSI is extremely challenging in this medium unless a powerful accelerator is present at close distance, increasing locally the CR density well above $n_{\text{CR}} = 10^{-9} \text{ cm}^{-3}$.</footnote>
### 3. Numerical Simulations Using MHD-PIC Approach

We use the MHD-PIC method introduced by Bai et al. (2015) and adopted for the study of CR streaming in Paper I. The focus of the present study is the inclusion of IN wave damping. Let us recall key elements of this method:

1. The background thermal plasma dynamics is governed by the ideal MHD equations without explicit viscosity or resistivity, adopting an adiabatic equation of state with $\gamma = 5/3$ and initial sound speed $c_s = V_{A,i}$. The momentum update includes an explicit damping term, implementation of which is described in Section 3.1. The dynamics of CRs are solved by a PIC method, where CRs are treated as charged macroparticles. The Lorentz equation is solved using a relativistic Boris pusher, the density and electric current are deposited on a grid using second-order shape factors (TSC scheme), and coupling between CRs and the background fluid is achieved by introducing CR source terms into the MHD equations and modifying Ohm's law.

2. The $\delta f$ method is employed to significantly reduce the Poisson noise of macroparticles.

3. Phase scrambling is applied when particles cross the system boundary and reenter on the other side of the periodic box; this effectively mimics a larger numerical box.

4. Seed waves are initialized at $t=0$. This provides better control of the k-by-k growth rate as compared to growth from numerical noise (typically done in full-PIC simulations).

#### 3.1. Numerical Treatment of IN Damping

The presence of neutrals introduces momentum exchange between ions and neutrals through collisions. The corresponding terms in the fluid equations were given previously in Equations (12) and (13). Instead of solving the coupled system of two-fluid equations, we only account for the effect of neutrals on the ions while the dynamical equations for neutrals are not evolved. This procedure is consistent in the limit where the two fluids are decoupled (high-frequency regime). The procedure is as follows. At each time step we reduce transverse momentum fluctuations according to Equation (12). After each numerical time step $\Delta t$ the transverse momentum is updated as $p_{\perp, new} = p_{\perp} \exp(-\nu_{in} \Delta t)$, in order to account for the effect of damping. While the longitudinal velocity should be subject to the same “damping,” it is not incorporated as it is decoupled from the Alfvén waves (which have no longitudinal motion) and we have verified that it does not affect the overall simulation results.

#### 3.2. Numerical Setup

In the present work, there are some differences from the Paper I fiducial setup:

1. We use fewer CR macroparticles per cell, as the $\delta f$ method effectively controls the Poissonian noise level.

2. The physical size of the grid is slightly smaller. Yet, the box size is equal to $\simeq 47$ times the most unstable wavelengths.

3. We include the wave damping term in the MHD momentum equation.

4. The fiducial drift velocity is $10V_{A,i}$ (instead of $2V_{A,i}$ used in Paper I). This choice is motivated by the need to efficiently reach a fully saturated state when there is no IN damping. This case will be used as reference when measuring the effect of different levels of damping.

Table 2 summarizes the numerical parameters and model values adopted in our simulations. Our fiducial choice is $n_{CR}/n_i \sim 10^{-4}$. This is somewhat higher than realistic values in the neutral ISM ($n_{CR}/n_i \sim 3 \times 10^{-7}-10^{-5}$; see Table 1), for numerical expediency; lower values would have very low CRSI growth rates (requiring very long simulation duration) and low saturation amplitudes (exacerbating the numerical issue of crossing $\mu=0$ for a practical resolution). This choice does not affect any of our theoretical conclusions. The box length is chosen to be much larger than the fastest-growing wavelength of the instability. Lengths are given in units of ion skin depth $d_i \equiv C/\omega_{p,i} = V_{A,i}/\Omega_0$, which is 1 in code units. We adopt a reduced speed of light $C = 300V_{A,i}$. Seed waves are initialized with equal amplitudes of right-handed and left-handed polarization and the same amplitude at all $k$. The initial total power in the waves is $\delta B/B_0 = 10^{-4}$. We work in the frame where CRs are initially at rest, and gas moves in the $-e_x$ direction.

It is convenient to parameterize damping rates $\nu_{in}$ relative to the peak growth rate for CRSI; this is what is reported in Table 2. Using $V_D = 10V_{A,i}$ and $n_{CR}/n_i = 10^{-4}$ in Equation (8), this initial peak growth rate (for $\nu_{in}=0$) is $\Gamma_{max,0}/\Omega_0 = 2.5 \times 10^{-4}$.

### 4. Simulation Results

In this section we describe the results of the MHD-PIC simulations. We start with a presentation of the overall time evolution, then show the linear growth rates of the instability. We then discuss the transition into the saturated phase and some aspects of the saturated state of the streaming CR–background fluid system.

#### 4.1. Overall Evolution

In this section we compare the evolution between no-damping and damping cases (models Fid and Fid-Damp reported in Table 2). For the cases with damping we explore values of $\nu_{in}$ between $0.03\Gamma_{max,0}$ and $1.9\Gamma_{max,0}$; the theoretical critical value for no growth is $\nu_{in} = 2\Gamma_{max,0}$. In the following, we use the term “moderate damping” for cases where $\nu_{in}/\Gamma_{max,0}$ is non-negligible but still lower than unity. The representative case is $\nu_{in}/\Gamma_{max,0} = 0.5$.

Figure 2 presents the time evolution of the magnetic wave energy (upper panel) and of the relative velocity between the bulk CRs and the background fluid (bottom panel) for different values of $\nu_{in}$. The undamped case (blue curves) reaches saturation at $t\Omega_0 \simeq 10^5$ with $\delta B^2/B_0^2 \simeq 10^{-3}$. The CR streaming velocity, initially $V_D = 10V_{A,i}$, begins to decrease when $\delta B^2/B_0^2 \simeq 10^{-4}$ and reaches $\nu_{rel} \simeq V_{A,i}$ by the end of the simulation. Here, the CR bulk velocity is defined (in the simulation frame) by

$$ \nu_{rel} \equiv \frac{\int f(p, \mu) v(p) p^2 \mu d\mu dp}{\int f(p, \mu) p^2 d\mu dp} - V_{bg}, \quad (28) $$
<table><caption>**Table 2**
List of Main Simulation Runs</caption><thead><tr><th>Run</th><th>V<sub>D</sub>/V<sub>A,i</sub></th><th>n<sub>CR</sub>/n<sub>i</sub></th><th>ν<sub>in</sub>/Γ<sub>max,0</sub></th><th>Domain Size<br>L<sub>x</sub> (d<sub>i</sub>)</th><th>Domain Size<br>L<sub>x</sub>/λ<sub>m</sub></th><th>Resolution<br>Δx (d<sub>i</sub>)</th><th>N<sub>p</sub><br>(per cell)</th><th>Runtime<br>(Ω<sub>0</sub><sup>-1</sup>)</th></tr></thead><tbody><tr><td>Fid</td><td>10.0</td><td>1.0 × 10<sup>-4</sup></td><td>0</td><td>8 × 10<sup>4</sup></td><td>42.4</td><td>10</td><td>128</td><td>10<sup>6</sup></td></tr><tr><td>Fid-Damp</td><td>10.0</td><td>1.0 × 10<sup>-4</sup></td><td>∈ [0.03, 1.9]</td><td>8 × 10<sup>4</sup></td><td>42.4</td><td>10</td><td>128</td><td>10<sup>6</sup></td></tr><tr><td>HiRes</td><td>10.0</td><td>1.0 × 10<sup>-4</sup></td><td>{0, 0.1, 0.5}</td><td>8 × 10<sup>4</sup></td><td>42.4</td><td>2.5</td><td>128</td><td>&lt;10<sup>6</sup></td></tr></tbody></table>

**Notes.** Fixed parameters: $C/V_{A,i} = 300$, $p_0/(mV_{A,i}) = 300$, $κ = 1.25$, and initial wave amplitude $A = 10^{-4}$. In all models the most unstable wavelength is 0.91$λ_m$ for $λ_m = 2πp_0/(mΩ_0) ≈ 1885d_i$. The ion skin depth $d_i ≡ C/ω_i = V_{A,i}/Ω_0$. The quantity $Γ_{max,0}$ is defined in Equation (8) as the maximum growth rate of the CRSI without IN damping (i.e., $ν_{in} = 0$).

where $V_{bg}$ is the background gas speed. Unlike the initial conditions, the CR distribution is not perfectly isotropic in the frame moving at $ν_{rel}$.

Cases with nonzero damping reach lower saturation levels, with the maximum in $δB^2/B_0^2$ decreasing as $ν_{in}$ increases. After saturation (at late simulation times, $tΩ_0 > 10^5$) the wave intensity stabilizes at an intensity that decreases with increasing $ν_{in}$.<fnref n="9" /> Streaming velocities (bottom panel) decline at lower rates for the stronger-damping models, because the lower wave amplitudes scatter CRs less effectively. The late-time (after $tΩ_0 \sim 10^6$) evolution is asymptotically slow.

In general, we define three representative phases of the time evolution of the instability: linear, post-linear, and saturated (late-time). There is an adjustment phase between the time when the fastest modes stop growing (end of the linear phase) and the saturation of the instability when the anisotropy of the DF is (eventually) erased (start of the saturated phase). We refer to this transitory phase as the post-linear phase.

In Figure 3 we show the magnetic field profile $B_z$ at four different simulation times for the cases with no damping (left column) and moderate damping (right column). For the case with no damping the wave amplitude grows to $δB_z/B_0 \approx 0.02$ at the end of the linear phase ($tΩ_0 = 2 \times 10^4$) and further increases to $0.05$ during the post-linear evolution (between $tΩ_0 \approx 3 \times 10^4$ and $tΩ_0 = 10^5$), during which amplitudes increase at $k$ below the peak (see Section 4.3). For the case with moderate damping, the wave amplitude is reduced, as expected. During the post-linear phase we observe a decrease of the wave amplitude from $δB_z/B_0 \approx 0.02$ to $4 \times 10^{-3}$, while there is less contribution from $k$ below the peak.

The effect of CRSI and subsequent QLD on the CR DF is presented in Figure 4. We show the DF in the wave frame (i.e., moving to the right at $V_{A,i}$ with respect to the gas; see Paper I for its transformation from the simulation frame), $δf_w(p_w, \cos θ_w)/f_0$ for the cases with no damping (left column) and with moderate damping (right column), at the same simulation times as those in Figure 3. By inspecting the left column, we observe the gradual suppression of the initial anisotropy with time (evolution from top to bottom). The anisotropy is globally maintained for the case with damping, although at moderate $p_w$ the distribution becomes relatively flat on each side of $μ_w = 0$. The effect of particle accumulation at $μ_w = \cos θ_w = 0$ is also evident in the two bottom panels on the right.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Profile of the wave component $B_z$ at different simulation times for the cases with no damping (left column) and with moderate damping $ν_{in} = (1/2)Γ_{max,0}$ (right column). Four different times are presented. The top set corresponds to the early linear phase. The two middle panel sets show the evolution during the post-linear phase, and the bottom set shows the final simulation time $tΩ_0 = 10^6$.</figcaption></figure>

## 4.2. Linear Phase

### 4.2.1. No Damping versus Moderate Damping

Paper I demonstrated that our MHD-PIC numerical approach accurately reproduces the linear growth rate of the CRSI in the absence of explicit wave damping. Here, we extend the previous investigation by including IN damping. In Figure 5, we show the linear growth rate $Γ_{tot}$ as a function of the wavenumber $k$ for two cases: no damping (top panel) and moderate damping (bottom panel). There is very good agreement with theory, despite some stochastic noise. In particular, the most unstable mode is very well captured at $k_{max} = \sqrt{2 - 1/κ} m_i Ω_0 / p_0 \approx 1.1 R_{L,0}^{-1}$.

Comparing the top and bottom panels of Figure 5, two immediate effects due to IN damping can be identified. The first is the decrease of the maximum growth rate near $kR_{L,0} = 1$. It is equal to $2.5 \times 10^{-4}Ω_0$ without damping, and equal to

<footnote n="9">⁹ We noticed that the level of the wave intensity “plateau” at the late-time state can evolve somewhat differently depending on the numerical resolution and the number of particles per cell, but does not affect qualitatively the system evolution.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** 2D DF in the wave frame, $\delta f_w(p_w, \cos \theta_w)/f_0$, at different simulation snapshots for the cases with no damping (left column) and with moderate damping $\nu_{in} = (1/2)\Gamma_{max,0}$ (right column). The same simulation times in Figure 3 are presented from top to bottom.</figcaption></figure>

$\simeq 1.8 \times 10^{-4}\Omega_0$ when $\nu_{in} = 1.25 \times 10^{-4}\Omega_0 = \Gamma_{max,0}/2$. The second effect is the suppression of wave growth in the low-end and high-end parts of the spectrum, where $\Gamma_0(k) - \nu_{in}/2 < 0$. In other words, the bandwidth of linearly growing waves is reduced with increasing $\nu_{in}$.

### 4.2.2. Dependence on $\nu_{in}$

In Figure 6 we report the dependence of the maximum growth rate of the instability on $\nu_{in}$. Simulation results are generally in excellent agreement with analytical predictions, namely, the linear growth rate is only weakly reduced for $\nu_{in}/\Gamma_{max,0} < 0.2$, and rapidly decreases for $\nu_{in}/\Gamma_{max,0} > 0.5$. As expected from Equation (24), there can be no linear growth for $\nu_{in}/\Gamma_{max,0} > 2$, as even the fastest-growing mode is damped.

## 4.3. From Linear to Saturated Phase

The instability does not transition directly from linear growth to a fully saturated state. At some point in time the fastest-growing modes at $kR_{L,0} \simeq 1$ cease exponential growth. However, at that time the other modes with lower growth rates can continue to grow because the anisotropy in the DF—which drives all modes—is not erased. At the same time, wave-wave interaction can redistribute some wave energy from $kR_{L,0} \simeq 1$ to other wave modes.

The post-linear phase can be easily identified in Figure 2 for $\nu_{in} = 0$ and $\nu_{in} \lesssim 0.1$. The initial exponential growth slows down roughly at $t\Omega_0 \simeq 3 \times 10^4$ but additional wave growth continues until $t\Omega_0 \sim 10^5$. This phase also corresponds to the fastest rate of decrease in the streaming velocity of CRs, as seen in the lower panel of Figure 2.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Linear growth rate as a function of the wavenumber $k$ for the cases without damping (top panel) and with moderate damping (bottom panel). Blue and red lines correspond to the measured growth rate of the right- and left-handed modes, respectively. Dashed yellow and magenta lines correspond to the analytical expectation from Equation (24). In the bottom panel the analytical expectations for zero damping are also plotted (dotted–dashed) for comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Dependence of the maximum linear growth rate of the instability as a function of the IN momentum exchange rate $\nu_{in}$ for $n_{CR}/n_i = 10^{-4}$ and $V_D/V_{A,i} = 10$.</figcaption></figure>

More detailed insight into the post-linear phase can be gained by considering the wave spectrum. In Figure 7 we present the evolution of the wave spectrum for the cases with no damping (top panel) and with $\nu_{in}/\Gamma_{max,0} = 0.5$ (bottom panel). We note that the spectrum is narrow at the end of the linear phase, as shown by the blue line. It is peaked at $kR_{L,0} \simeq 1$, where the growth rate is fastest. During the post-
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Spectra of forward right-handed modes at different times for the cases with no damping (top panel) and moderate damping (bottom panel). The evolution for forward left-handed modes is identical. In the lower panel, the gray region marks where the initial linear growth is suppressed by IN damping.</figcaption></figure>

linear phase (blue to red line), while the growth at $kR_{L,0} \simeq 1$ becomes slow, the modes at $kR_{L,0} > 1$ start growing at rates comparable to $\Gamma_{\text{max}}$, while the modes with $kR_{L,0} \le 1$ continue to grow approximately at the rate expected from the linear calculation. If there is no damping, during the late evolution the peak in the wave spectrum gradually shifts to larger wavelengths (smaller $k$) while the global level of wave intensity is roughly unchanged. This effect can be clearly seen by comparing the red, yellow, and purple lines in the top panel: the low-$k$ cutoff shifts from $kR_{L,0} \simeq 0.3$ at the end of the linear phase to $kR_{L,0} \simeq 0.02$ at the end of the simulation. The growth of high-$k$ modes (i.e., $kR_{L,0} > 5$) is observed only during the post-linear phase, while the increase of spectral energy in $kR_{L,0} \ll 1$ modes continues well beyond, into the saturated phase.

For the case with moderate damping ($\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.5$), the early evolution of the spectrum follows the same trend as in the case with no damping (blue and red lines in the bottom panel). Interestingly, there is a similar fast rise of modes at $kR_{L,0} \ge 1$ in the post-linear phase, despite the fact that some of these modes are not supposed to grow linearly (there would be a cutoff for $kR_{L,0} > 8$, marked gray in the figure). The late-time evolution of the damped case is significantly different from that of the undamped case. There is a noticeable overall decrease in the wave intensity at all $k$ (difference between red, yellow, and purple lines). The peak of the spectrum shifts slightly toward

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Wave spectrum (forward right-handed mode) for the HiRes simulation with no damping at different simulation times. Blue and yellow lines show the standard CR-driven system at $t\Omega_0 = 2 \times 10^4$ and $t\Omega_0 = 8.2 \times 10^4$, respectively. The red line shows the case of a model where CRs are turned off at the end of the linear phase ($t\Omega_0 = 2 \times 10^4$) and which then freely evolves until $t\Omega_0 = 8.2 \times 10^4$. The difference between the blue and red lines illustrates the growth of high-$k$ modes during the post-linear phase without driving by CRs. The inset in the bottom shows the time evolution of the magnetic wave energy: the blue line follows the standard evolution until $t\Omega_0 = 2 \times 10^4$, the yellow line continues the evolution with CRs after that time, and the red line follows the case where CRs are switched off at $t\Omega_0 = 2 \times 10^4$. The green dot marks the time of the CR switch-off.</figcaption></figure>

smaller $k$ at the end of the post-linear phase. Also, a spectral bump appears at $t\Omega_0 = 10^5$ around $kR_{L,0} = 0.2$, which is due to driving from the part of the CR DF that remains strongly anisotropic: $f(p > 5p_0)$. By the end of the simulation the spectrum stabilizes. It is flat and narrow: $0.1 < kR_{L,0} < 10$ and $kl(k) \sim 2 \times 10^{-7}$. Simulations with different values of $\nu_{\text{in}}$ follow a similar time evolution.

### 4.3.1. Growth of High-k (Small-wavelength) Modes: HiRes Simulations

An important aspect of the post-linear phase is the rapid growth of modes with $kR_{L,0} \ge 5$, even when these modes are not unstable to CRSI because of IN damping. To better study these short-wavelength modes we have conducted additional simulations with higher numerical resolution: $\Delta x = 2.5d_i$ in the HiRes runs instead of $\Delta x = 10d_i$ in the Fid simulations. To understand the mechanism driving this growth, we also performed a numerical experiment where the CRs were removed from the system after a given simulation time, at the end of the linear phase.

In Figure 8 we show the wave spectrum of forward right-handed modes in the HiRes simulation with no damping. Three lines are plotted: (i) the spectrum at the end of the linear phase (blue line), (ii) the spectrum at the end of the post-linear phase if CRs are switched off at the end of the linear phase (red line), and (iii) the spectrum at the end of the linear phase in the standard case with no CRs switched off (yellow line). In the inset, we see that the wave energy stays constant when the CRs are switched off, whereas it continues to grow otherwise. The difference between the blue and red lines illustrates the growth
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Same as Figure 8 but with $\nu_{in}/\Gamma_{max,0} = 0.5$. The additional dotted-dashed black lines compensate analytically for the effect of IN damping on the wave spectrum. The gray-shaded areas delimit the regions where wave growth is not allowed in linear theory (Equation (24)).</figcaption></figure>

of high-k modes during the post-linear phase without driving by CRs. This difference is only seen in the region $kR_{L,0} \ge 5$. For comparison, the difference between the blue and yellow lines shows the post-linear evolution when the driving by CRs is maintained. Here, an overall increase of the wave energy is observed, together with a shift of the spectrum to lower k.

The same effect is present in simulations with IN damping. Figure 9 is the same as Figure 8 but with $\nu_{in}/\Gamma_{max,0} = 0.5$. The additional dotted-dashed black lines analytically compensate for the effect of IN damping on the wave spectrum. During the post-linear phase, a decrease of wave energy at any k is imposed by the IN damping (difference between blue and red lines). If one compensates for the effect of damping, we observe the same excess of high-k modes at the end of the post-linear phase (dotted-dashed black line) as for the undamped case.

Because these high-k modes can grow even in the absence of CRSI, we conclude that the driving mechanism is purely an MHD effect (see further discussion in Section 5).

## 4.4. Saturated Phase

At the end of the post-linear phase, the exponential growth of the instability at all wavelengths is completed. By this time the CRs have also fully experienced the back-reaction from interacting with the waves they generated. In this section, we characterize the particle and wave properties based on analysis of our MHD-PIC simulations at late stages.

### 4.4.1. Late-stage Wave Amplitudes and Particle Distributions

In Figure 10 we present the maximum wave intensity (normalized to $B_0^2$) reached in the simulations as a function of $\nu_{in}$ (plotted using blue circles and red squares for the Fid and HiRes simulations, respectively). For undamped or weakly damped cases we expect the momentum flux associated with the original anisotropy in the CRs to be transferred to forward-propagating Alfvén waves (see Kulsrud 2005). As found in

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Dependence of the maximum value of the magnetic wave energy, $\delta B_{\text{max}}^2 / B_0^2$, as a function of $\nu_{in}$. Blue circles are from Fid simulations and red squares from HiRes simulations. The green solid curve corresponds to the model solving two coupled ordinary differential equations (ODEs) described in the text. The horizontal black line delimits the maximal allowed value at $\nu_{in} = 0$, based on the measured saturation amplitude. For $\nu_{in}/\Gamma_{max,0} > 2$ there is no linear instability possible; this limit is shown by the gray-shaded region on the right side of the figure.</figcaption></figure>

previous simulations (e.g., Bai et al. 2019; Holcomb & Spitkovsky 2019), this leads to magnetic wave energy at saturation given by $\delta B_{\text{sat}}^2 / B_0^2 \approx 1.5(n_{\text{CR}}/n_i)(V_D/V_{A,i} - 1)$. This value is plotted as a horizontal black line in Figure 10. The gray-shaded region at $\nu_{in} > 2\Gamma_{\text{max},0}$ delimits the region where the linear instability becomes impossible, according to Equation (24).<fnref n="10" /> The simulations show a gradual decrease in $(\delta B_{\text{max}} / B_0)^2$ with increasing $\nu_{in}$, as might be expected. The decrease becomes abrupt when approaching $\nu_{in} = 2\Gamma_{\text{max},0}$, resembling an exponential cutoff.

In order to obtain some insight into the dependence of the saturated wave intensity on $\nu_{in}$, we cast our knowledge of dominant dynamical processes into a simple model. This model consists of a system of two coupled ODEs. Let $A = (\delta B / B_0)^2$ be the wave amplitude squared, and $\Gamma$ be the wave growth rate (without damping). In our simple model, the evolution of A in time is determined by

$$ \frac{d(\ln A)}{dt} = 2\left(\Gamma - \frac{\nu_{in}}{2}\right), \qquad (29) $$

$$ \frac{d(\ln \Gamma)}{dt} = -\nu_s \approx -A \frac{\pi}{8} g, \qquad (30) $$

where g is some factor of order unity. The first equation describes linear wave growth partially limited by IN damping. The second equation serves as a proxy for the dynamical adjustment of the growth rate to the changing CR distribution, which is becoming more isotropic under the effect of QLD in the bath of (growing) waves. For initial conditions, we set $A(0) = 10^{-10}$ and $\Gamma(0) = \Gamma_{\text{max},0}$. We further parameterize the IN damping rate as $\nu_{in} = \alpha \Gamma_{\text{max},0}$. In our simulations, we have $\Gamma_{\text{max},0} = 2.5 \times 10^{-4} \Omega_0$, and $\alpha$ ranging from 0 to 2. We

<footnote n="10">¹⁰ We verified that there is no wave growth for $\nu_{in} > 2\Gamma_{\text{max},0}$ with a dedicated simulation, not presented here.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Final spectrum (top panels), DF in the simulation frame $\delta f(p, \cos\theta)/f_0$ (middle panels), and DF in the wave frame $\delta f_w(p_w, \cos\theta_w)/f_0$ (bottom panels) in three simulations with different $\nu_{\text{in}}$. The left, middle, and right columns have $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.1, 0.5, \text{and } 1.5$, respectively. The gray-shaded regions in the upper panels delimit the regions where waves cannot grow during the linear phase, according to Equation (24). Thick dashed lines in the middle and bottom panels are the contours in momentum space resonant with the limiting $k$ (i.e., the inner border of the gray-shaded areas in the upper panels). The gray dotted-dashed lines illustrate resonant contours for a few different wavenumbers.</figcaption></figure>

integrate these equations until $\Gamma/\Gamma_{\text{max},0} < \alpha/2$ so that $\mathcal{A}$ reaches a maximum. In the case without damping, the wave amplitude should reach the value measured in the simulations; from this constraint we find $g \approx 1.2$.

The solutions for $\mathcal{A}$ as a function of $\nu_{\text{in}}$ are plotted with a green solid line in Figure 10. The model qualitatively agrees with the simulation results and shows a good match with our high-resolution simulation results when damping is weak ($\nu_{\text{in}} \lesssim 0.1\Gamma_{\text{max},0}$), but overpredicts $\delta B_{\text{max}}$ in the case of moderate-to-strong damping ($0.2 < \nu_{\text{in}}/\Gamma_{\text{max},0} < 1$).

While a better match can be achieved with additional parameters and fine-tuning, this does not necessarily add further insight into the physical processes involved. In particular, we note that an important effect not captured in our toy model is the $\mu=0$ barrier (which is partly physical and partly numerical). When the isotropization process becomes stuck as particles accumulate at the barrier, further wave growth is suppressed, even though the free energy from the CR anisotropy has not been fully utilized.

Indeed, in Figure 14, it is evident that the HiRes model with $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.1$ (which shows good agreement with the toy model in Figure 10) does not have appreciable particle buildup at $\mu=0$, but the $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.5$ model (which falls below the model prediction) does have particles built up at $\mu=0$.

Figure 11 presents the final state (at $t\Omega_0 = 10^6$) of three representative simulations with different values of $\nu_{\text{in}}$: the columns from left to right show models using $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.1$,

0.5, and 1.5, respectively. The gray-shaded regions in the top panels mark the regions where waves cannot grow during the linear phase, according to Equation (24).

Some interesting features of the saturated phase are apparent in Figure 11. First, there is a reduction in the global wave intensity with an increasing IN damping rate (compare from left to right the top panels). Second, the saturated wave spectrum width exceeds the range imposed by the linear growth, with significant wave amplitude in the gray-shaded regions (best evidenced in the top right panel). Also, the forward-propagating waves (blue and red lines in the top panels) are at a similar level to the backward-propagating waves (yellow and purple lines in the top panels), and a marked reduction from the peak values reached during the post-linear phase brings all modes to comparable intensity levels. High-momentum particles (regions with $p > 10p_0$ in the middle and bottom panels) clearly isotropize less efficiently with increasing $\nu_{\text{in}}$: from left to right in the bottom panels, there is an increasing area that is unaffected by the instability. This corresponds to the lack of waves in the $kR_{L,0} \ll 1$ region that becomes more and more prominent when $\nu_{\text{in}}$ approaches the critical value $2\Gamma_{\text{max},0}$. Finally, we observe an accumulation of particles near $\mu=0$ in all cases with IN damping ("hot spots" in the regions close to $\cos\theta=0$ in the middle panels and close to $\cos\theta_w=0$ in the bottom panels). The issue of crossing the $\mu = \cos\theta = 0$ barrier becomes crucial when non-negligible damping is present. Not only is the global level of waves reduced, but
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Effect of numerical grid resolution on the time evolution of the magnetic wave intensity (top panel) and of the CR streaming velocity (bottom panel). Blue lines show the case with no damping, and red lines show the case with $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.5$. Solid lines correspond to the fiducial resolution and dashed lines correspond to the high-resolution simulations (see Table 2).</figcaption></figure>

also the whole spectrum becomes narrower when $\nu_{\text{in}}$ approaches $2\Gamma_{\text{max},0}$. Both effects contribute in making it difficult for particles to scatter across $\mu = 0$.

4.5. Dependence on Spatial Resolution

By testing different grid resolutions we confirm that there is no noticeable effect on the linear phase of the instability: no difference in the growth rate, wave energy, streaming velocity, or spectrum during this phase. However, the resolution can have an important impact on the post-linear phase. Specifically, resolving a larger dynamic range of wave modes with $kR_{L,0} \gg 1$ allows us to better capture the fast post-linear growth of these modes. This can be seen by comparing the high-$k$ part of the power spectra at $t\Omega_0 = 10^5$ in Figure 7 for the fiducial resolution to those at similar times ($t\Omega_0 = 8.2 \times 10^4$, with CRs kept on) in Figures 8 and 9 at high resolution.

In Figure 12 we compare the Fid and HiRes simulations for the time evolution of the magnetic wave energy (top panel) and of the CR streaming velocity (bottom panel). Two representative cases are shown: no damping (blue lines) and moderate damping $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.5$ (red lines). There is no difference in early evolution, until the post-linear phase, starting at $t\Omega_0 \sim (3-4) \times 10^4$. From the solid (Fid) and dashed (HiRes) blue lines, increased resolution has only a minor effect on the

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Dependence on spatial resolution of the saturated wave spectrum (forward right-handed modes) for models with $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.1$ (top panel) and $\nu_{\text{in}}/\Gamma_{\text{max},0} = 0.5$ (bottom panel). Blue lines correspond to Fid simulations and red lines correspond to HiRes simulations.</figcaption></figure>

case with no damping. More noticeable differences appear for the case with moderate damping. During the post-linear phase, the magnetic energy decreases more slowly in the HiRes simulation and the streaming velocity decreases at a faster rate. This effect is due to the higher level of wave intensity stored in high-$k$ modes during the post-linear phase.

As found for the post-linear stage, the numerical resolution also has an effect on the saturated state of the simulations. We expect, and indeed find, that higher-resolution simulations achieve a higher level of small-scale waves. In Figure 13 we present a comparison of the wave spectrum at $t\Omega_0 = 3 \times 10^5$ between the Fid and HiRes simulations for two values of $\nu_{\text{in}} = 0.1$ (top) and $0.5$ (bottom). For both values of $\nu_{\text{in}}$ there are some differences around $kR_{L,0} = 1$. But more importantly, the HiRes simulations exhibit higher levels of wave intensity at $kR_{L,0} \gg 1$, without regard to the limit of linear growth (which is equal to 0 in gray-shaded areas).

The differences in the wave spectrum between the Fid and HiRes simulations also lead to differences in the final state of the CR DF. In particular, higher amplitudes at large $k$ increase the scattering rate at resonances close to $\mu = 0$, considering the resonance condition $k = (1/\mu)m_p\Omega/p$. Figure 14 presents the DF at $t\Omega_0 = 3 \times 10^5$ of the simulations for the Fid runs (upper panels) and HiRes runs (lower panels). We show cases with both weak (left) and moderate (right) damping, as in Figure 13. For both, the HiRes simulations show a better level of isotropization of $\delta f_w$ at $p/p_0 \in [0.5, 10]$, corresponding to the higher level of high-$k$ modes in the HiRes simulations. No
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Dependence on spatial resolution of the DF $\delta f_w (p_w, \cos \theta_w)/f_0$, at $t\Omega_0 = 3 \times 10^5$. The top panel corresponds to Fid simulations and the bottom panel corresponds to HiRes simulations. Left and right show weak and moderate damping cases.</figcaption></figure>

noticeable difference is observed at $p > 10p_0$, corresponding to quite similar spectra at $kR_{L,0} < 0.1$ (see Figure 13).

## 4.6. Streaming Velocity

We now consider the evolution of the CR bulk or streaming velocity, a key parameter since this characterizes the net CR flux. As mentioned in Section 2.4, in the astrophysical literature the most commonly adopted assumption is that the wave amplitude and streaming velocity are consistent with a state in which the linear growth rate balances the damping rate (Kulsrud & Cesarsky 1971). From Equations (8) and (24), this would imply that if the CRs remain isotropic in a frame moving at $V_{st}$ relative to the gas, this velocity will decline until Equation (27) is satisfied.

Conceptually, the asymptotic streaming velocity obtained through balancing damping and growth relates to the situation where the instability is driven, i.e., through a CR pressure gradient. In contrast, our numerical setup with periodic boundary conditions corresponds more to a transient situation in which instability dies out. This affects the astrophysical relevance of the streaming speed derived in our simulations. Nevertheless, it is interesting to follow the measured time evolution of CR streaming from the simulations.

In the bottom panel of Figure 15 we show the evolution of the measured CR streaming velocity $v_{rel}$ (defined in Equation (28)), including its dependence on the numerical resolution, for the weak damping case $v_{in}/\Gamma_{max,0} = 0.1$. In the Fid simulation (blue line) the final streaming velocity of CRs is $3.3V_{A,i}$. The HiRes simulation isotropizes more efficiently, showing a decrease of the streaming velocity to $v_{rel} = 1.18V_{A,i}$ (but still above full isotropization $v_{rel} = V_{A,i}$). We attribute this difference to better capturing of high-k modes in the HiRes simulations.

If we were to apply Equation (27) to the case with $v_{in} = 0.1\Gamma_{max,0}$, the predicted asymptotic streaming velocity would be $V_{st}/V_{A,i} = 1.45$. This level is marked by the horizontal dashed line in Figure 15. Interestingly, for the HiRes simulation the measured final time $v_{rel}$ is below this. This can be understood from Figure 13: as long as a sufficient

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Dependence on spatial resolution of the time-evolving wave energy (panel a) and streaming velocity (panel b), for $v_{in}/\Gamma_{max,0} = 1/10$. The blue and red lines correspond to the Fid and HiRes simulations, respectively. The horizontal dashed line delimits the final streaming velocity as expected by balancing the fastest growth rate with the damping rate: $\Gamma_{max,0} = v_{in}/2$.</figcaption></figure>

level of waves is present the distribution will continue to isotropize and $v_{rel}$ will continue to decline. In contrast, for the HiRes model with $v_{in} = 0.5\Gamma_{max,0}$ (see Figure 12, red lines), the late-time $v_{rel}$ remains well above the value $3.25V_{A,i}$ that would be predicted by Equation (27), presumably due to the lower wave amplitudes that reduce the scattering rate (compare upper and lower panels of Figure 13, and note also the buildup of particles to the right of $\mu = 0$ in the lower-right corner of Figure 14). Taken together, these results make clear that the simple approach of setting $\Gamma_{max,0} = v_{in}/2$ and solving for streaming velocity is not necessarily applicable. The evolution and final value of $v_{rel}$ instead depends on the evolution of $(\delta B/B_0)^2$. This is illustrated for different values of $v_{in}$ in Figure 12, with more complete isotropization in cases of lower $v_{in}$. In the real astrophysical case, an additional factor affecting evolution would be the history of (anisotropic) energy sources.

Finally, we recall (as pointed out in Section 2.4) that balancing the linear growth time and damping time of waves does not quantitatively take into account the additional time required for QLD to isotropize the distribution, suggesting that the criterion $\Gamma_{max,0} > v_{in}/2$ is necessary but not sufficient. Because the pitch angle diffusion timescale near saturation is approximately proportional to the linear growth timescale, there may still be some critical ratio of growth to damping that allows the distribution to isotropize. This question could be addressed in future work that allows for more realistic driving.
<table><caption>**Table 3**
Maximum Momentum of CRs Able to Self-confine for Different Phases of the ISM</caption><thead><tr><th>Phase</th><th>(<i>p</i>μ)<sub>max</sub>/p<sub>0</sub>, V<sub>D</sub> = 2V<sub>A,i</sub></th><th>(<i>p</i>μ)<sub>max</sub>/p<sub>0</sub>, V<sub>D</sub> = 0.1c</th></tr></thead><tbody><tr><td>WNM</td><td>12.2</td><td>345</td></tr><tr><td>CNM</td><td>0.34</td><td>22.1</td></tr><tr><td>MG</td><td>1.5</td><td>41.7</td></tr><tr><td>DMG</td><td>≈1</td><td>&lt;1</td></tr></tbody></table>

Note. Other parameter values are as in Table 1.

## 5. Discussion

Perhaps one of the most intriguing results of the present work is the evidence of high-k growth during the post-linear phase, presented in Section 4.3. In Figures 8 and 9 we have shown that these modes do not require driving by CRs to be amplified. At present, we do not have a satisfactory explanation for this effect, but hypothesize that it is due to some form of mode coupling or wave steepening into rotational discontinuities (e.g., Cohen & Kulsrud 1974). One could argue that this effect is due to our choice of parameters, which leads to relatively high wave amplitude at saturation, up to $δB/B_0 \simeq 0.05$. This is likely considerably larger than what is present in the general ISM. However, we also performed simulations with $n_{CR}/n_i$ as low as $5 \times 10^{-5}$ and $V_D/V_{A,i} = 2$, observing the same effect of fast high-k mode growth during the post-linear phase. Potentially, the nonlinearity or three-wave interactions could be diagnosed by bispectrum and bicoherence analysis, but we defer this exploration to future work.

We also found that large-wavelength modes with $kR_{L,0} \ll 1$ are not amplified if CR driving is switched off at the end of the linear phase of the instability. This shows that these modes are mainly amplified by resonant interaction with CRs, with growth rate $\Gamma(k) \propto k^{2\kappa-1}$, and that IN damping removes the low-k end of the wave spectrum with $k < k_{min}$, where $\Gamma(k < k_{min}) < \nu_{in}/2$. This can be seen in Figures 7, 9, and 11. This implies, assuming only resonant wave-particle interaction, that high-momentum CRs with $p > p_{max}$ are not isotropized and continue to freely stream. Here, $p_{max}$ is deduced by using Equation (24) with Equation (7) using the resonance condition $kR_{L,0} = p_0/(p\mu)$:

$$ \frac{(p\mu)_{\max}}{p_0} \approx \left[ 1.55 \frac{n_{CR}}{n_i} \frac{\Omega_i}{\nu_{in}} \left( \frac{V_D}{V_{A,i}} - 1 \right) \right]^{\frac{2}{3}}, \quad (31) $$

where we fix $\kappa = 1.25$.

Like Kulsrud & Cesarsky (1971) we estimate the maximum momentum of CRs that could potentially be isotropized (i.e., self-confined) in different phases of the ISM. Using the typical values given in Table 1 and Equation (31), the derived values of $(p\mu)_{max}$ are given in Table 3, adopting $V_D/V_{A,i} = 2$ and $V_D = 0.1c$ as two extreme cases. Assuming $p_0 = m_p c$ for CRs with energy $E \simeq \text{GeV}$, we draw some general conclusions:

1. **DMG:** Strongly suppressed CRSI; no CR isotropization at energies higher than the gigaelectronvolt level. Streaming could become extremely large.

2. **MG:** Marginal instability and isotropy at the gigaelectronvolt level but not beyond 100 GeV. The surface layers of

molecular clouds may therefore be subject to CRSI if the initial anisotropy of CRs is not too small.

3. **CNM:** The conditions for triggering the CRSI are only satisfied for highly super-Alfvénic streaming. The CR distribution could become quite anisotropic in this phase.

4. **WNM:** Most favorable environment for CRSI. Even with tiny anisotropy ($V_D/V_{A,i}$ of order unity) the instability is possible for gigaelectronvolt CRs. CRs with $E \gg \text{GeV}$ could also be isotropized if the initial drift velocity is much higher than $V_{A,i}$.

In the Milky Way and similar galaxies, the WNM is the dominant component of the ISM by mass; the above estimate affirms that CRSI is astrophysically quite important. We note, however, that in the WNM other wave damping mechanisms can compete with IN damping that could prevent isotropization of CRs with energies $\gg \text{TeV}$ (see, e.g., Figure 1 and corresponding text in Brahimi et al. 2020), but our values of $E_{max}$ are in good agreement with similar calculations by Xu et al. (2016, their Section 6.4 and Figure 15).

The estimates above are in broad agreement with those of Kulsrud & Cesarsky (1971)—being based on the same argument—but are updated here with representative parameters for different neutral-dominated media.

The high-frequency limit is adopted in the present study, i.e., $\omega_A = kV_{A,i} \gg \max[\nu_{in}, \nu_{ni}]$. In this regime, the neutral fluid and ionized fluid are decoupled, and the IN wave damping rate does not depend on the wavelength, $\Gamma_d = \nu_{in}/2$. We expect this approximation to hold for any $kR_{L,0} \ge 0.1-1$, in general, and to be satisfied for CRs at gigaelectronvolt energies and below, which make up most of the CR energy density and are responsible for most of the ionization. This justifies the one-fluid approach adopted in the present study. However, in reality very high energy CRs with low resonant frequencies are present in the ISM as well. If the whole CR energy spectrum were fully represented, the low-k part of the wave spectrum (modes with $kR_{L,0} \ll \max[\nu_{in}, \nu_{ni}$]) would be damped at a slower rate than $\nu_{in}/2$. This could have an interesting effect on the amplification and survival of long-wavelength modes. However, we defer investigation of k-dependent damping to future studies.

Similar to Paper I, many of our simulations show particle accumulation at $\mu = 0$. Naively, this would be expected if waves were only present where CRSI is undamped (at intermediate $k$ near $k_{max} = \sqrt{2 - 1/\kappa R_{L,0}^{-1}}$, since there would be no large-k waves that are able to scatter small-$\mu$ particles subject to the resonance condition $k\mu = \Omega m/p$. However, in practice nonlinear MHD effects populate the large-k region regardless of whether CRSI is damped or active, and therefore $\mu = 0$ crossing is possible. Nevertheless, in simulations the high-k regime is subject to numerical dissipation, and we find that higher than standard resolution is required to limit particle buildup. At higher values of $\nu_{in}$, particles still build up near $\mu = 0$ even with higher resolution (see Figure 14).

As we have previously emphasized, the problem of relating the streaming rate and wave amplitude to macroscopic ambient properties clearly merits further study. The traditional “detailed balance” approach of equating the CRSI growth rate to the damping rate has recently been adopted in galaxy formation simulations and other studies as a procedure for setting the scattering rate coefficient, leading to a diffusion coefficient that varies proportionally to $\nu_{in}$ (see, e.g., Hopkins et al. 2021b). To assess and quantitatively improve this kind of prescription,
however, studies similar to the present one but allowing for macroscopic driving (via an imposed CR flux or energy gradient) will be needed.

## 6. Conclusion

Motivated by the recent progress in coupled fluid-kinetic (MHD-PIC) numerical techniques in Paper I, in this work we studied the influence of ambipolar diffusion (IN damping of Alfvén waves) on the CRSI. We adopted parameters for a reference model that in the absence of IN damping leads to full isotropization of the CR DF in the wave frame after the instability saturates. This state corresponds to an asymptotic state in which the CRs stream at Alfvén speed relative to the background fluid. We then explored different values of IN damping rate to study the influence on the outcome of the CRSI.

Our main conclusions are summarized below:

1. The predicted exponential growth rate of the instability including IN damping, $\Gamma(k)$, is well reproduced by the MHD-PIC technique. Thus the linear theory is in good agreement with analytical expectations.

2. Evolution in a post-linear phase is crucial for isotropization of low- and moderate-energy CRs. During this phase, growth of high-$k$ (wavelength $\lambda \ll R_{L,0}$) modes occurs. This growth is not driven by CR anisotropy (the waves are outside of the unstable range) but rather appears to be a result of an MHD wave cascade or wave steepening into rotational discontinuities. The high-$k$ waves are important for scattering when $\mu$ is close to 0. The development of the high-$k$ spectrum is best captured in high-resolution simulations.

3. Systematic comparison between the reference case (no damping) and simulations with IN damping reveals that the width of the wave spectrum decreases with increasing $v_{in}$. The absence of low-$k$ (compared to $R_{L,0}^{-1}$) waves prevents isotropization of high-energy CRs. There is no wave growth at all if $v_{in} > 2\Gamma_{max,0}$, consistent with analytic theory.

4. With stronger IN damping, the maximum amplitude of waves systematically decreases. When $v_{in} = 0.5\Gamma_{max,0}$, the peak wave amplitude is an order of magnitude below the no-damping case. Lower wave amplitudes reduce the rate of QLD and hence slow isotropization.

This work is the second study in a series exploring the physics of CRSI by means of MHD-PIC simulations. This numerical technique has proved to be valuable because it captures the relevant microphysical CR gyroscales and allows long-term simulations at low CR density more affordably than full kinetic simulations. We envision that further improvements of the method and inclusion of additional physics will lead to greater understanding of the microphysics of CR-ISM interactions in different phases, also enabling numerical calibrations of transport coefficients required for large-scale ISM/galactic studies.

We are grateful to the anonymous referee for a careful reading of our manuscript and an insightful report. I.P. was supported by National Science Foundation grants PHY-1804048 and PHY-1523261 and by the Max-Planck/Princeton Center for Plasma Physics. The work of E.C.O. was supported by grant 510940 from the Simons Foundation. X.-N.B.

acknowledges support by NSFC grant 11873033. Computations were conducted on resources at the PICSciE-OIT High Performance Computing Center and Visualization Laboratory at Princeton University, on CALMIP supercomputing resources at Université de Toulouse III (France) under allocations 2016-p1504 and P20028, and on the Orion supercomputer through the Department of Astronomy at Tsinghua University.

## ORCID iDs

Ilya Plotnikov https://orcid.org/0000-0002-0074-4048

Eve C. Ostriker https://orcid.org/0000-0002-0509-9113

Xue-Ning Bai https://orcid.org/0000-0001-6906-9549

## References

Amano, T. 2018, JCoPh, 366, 366

Amato, E., & Blasi, P. 2009, MNRAS, 392, 1591

Amato, E., & Blasi, P. 2018, AdSpR, 62, 2731

Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55

Bai, X.-N., Ostriker, E. C., Plotnikov, I., & Stone, J. M. 2019, ApJ, 876, 60, [Paper I]

Bambic, C. J., Bai, X.-N., & Ostriker, E. C. 2021, arXiv:2102.11877

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

Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of Cosmic Rays (Amsterdam: North-Holland)

Blasi, P., Amato, E., & Serpico, P. D. 2012a, PhRvL, 109, 061101

Blasi, P., Morlino, G., Bandiera, R., Amato, E., & Caprioli, D. 2012b, ApJ, 755, 121

Brahimi, L., Marcowith, A., & Ptuskin, V. S. 2020, A&A, 633, A72

Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79

Bustard, C., & Zweibel, E. G. 2020, arXiv:2012.06585

Butsky, I. S., & Quinn, T. R. 2018, ApJ, 868, 108

Bykov, A. M., Brandenburg, A., Malkov, M. A., & Osipov, S. M. 2013, SSRv, 178, 201

Bykov, A. M., & Toptygin, I. N. 2005, AstL, 31, 748

Cohen, R. H., & Kulsrud, R. M. 1974, PhFl, 17, 2215

Dashyan, G., & Dubois, Y. 2020, A&A, 638, A123

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Dubois, Y., & Commerçon, B. 2016, A&A, 585, A138

Dubois, Y., Commerçon, B., Marcowith, A. r., & Brahimi, L. 2019, A&A, 631, A121

Everett, J. E., & Zweibel, E. G. 2011, ApJ, 739, 60

Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, ApJ, 674, 258

Farmer, A. J., & Goldreich, P. 2004, ApJ, 604, 671

Ferrière, K. M. 2001, RvMP, 73, 1031

Ginzburg, V. L., Ptuskin, V. S., & Tsytovich, V. N. 1973, Ap&SS, 21, 13

Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042

Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJL, 816, L19

Glassgold, A. E., Galli, D., & Padovani, M. 2012, ApJ, 756, 157

Grenier, I. A., Black, J. H., & Strong, A. W. 2015, ARA&A, 53, 199

Haggerty, C., Caprioli, D., & Zweibel, E. 2019, in ICRC (Madison, WI), 36, 279

Hanasz, M., & Lesch, H. 2003, A&A, 412, 331

Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJL, 777, L38

Holcomb, C., & Spitkovsky, A. 2019, ApJ, 882, 3

Hopkins, P. F., Chan, T. K., Garrison-Kimmel, S., et al. 2020, MNRAS, 492, 3465

Hopkins, P. F., Chan, T. K., Squire, J., et al. 2021a, MNRAS, 501, 3663

Hopkins, P. F., Squire, J., Chan, T. K., et al. 2021b, MNRAS, 501, 4184

Ipavich, F. M. 1975, ApJ, 196, 107

Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O., et al. 2018, ApJ, 855, 23

Jiang, Y.-F., & Oh, S. P. 2018, ApJ, 854, 5

Krumholz, M. R., Crocker, R. M., Xu, S., et al. 2020, MNRAS, 493, 2817

Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445

Kulsrud, R. M. 1978, in Astronomical Papers Dedicated to Bengt Stromgren, ed. A. Reiz & T. Andersen (Copenhagen: Copenhagen Univ. Observatory), 317

Kulsrud, R. M. 2005, Plasma Physics for Astrophysics

Kulsrud, R. M., & Cesarsky, C. J. 1971, ApL, 8, 189

Lazarian, A. 2016, ApJ, 833, 131
Lebiga, O., Santos-Lima, R., & Yan, H. 2018, MNRAS, 476, 2779

Lee, M. A., & Völk, H. J. 1973, Ap&SS, 24, 31

Lerche, I. 1967, ApJ, 147, 689

Mao, S. A., & Ostriker, E. C. 2018, ApJ, 854, 89

Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13

Morlino, G., & Gabici, S. 2015, MNRAS, 451, L100

Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552

O'C Drury, L., Duffy, P., & Kirk, J. G. 1996, A&A, 309, 1002

Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, SSRv, 216, 29

Pais, M., Pfrommer, C., Ehlert, K., & Pakmor, R. 2018, MNRAS, 478, 5278

Pakmor, R., Pfrommer, C., Simpson, C. M., Kannan, R., & Springel, V. 2016, MNRAS, 462, 2603

Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017, MNRAS, 465, 4500

Recchia, S., Blasi, P., & Morlino, G. 2016, MNRAS, 462, 4227

Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433

Reville, B., Giacinti, G., & Scott, R. 2021, MNRAS, 502, 4137

Reville, B., Kirk, J. G., Duffy, P., & O'Sullivan, S. 2007, A&A, 475, 435

Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, ApJ, 834, 208

Schroer, B., Pezzi, O., Caprioli, D., Haggerty, C., & Blasi, P. 2020, arXiv:2011.02238

Shalaby, M., Thomas, T., & Pfrommer, C. 2021, ApJ, 908, 206

Silsbee, K., & Ivlev, A. V. 2019, ApJ, 879, 14

Skilling, J. 1971, ApJ, 170, 265

Skilling, J. 1975, MNRAS, 172, 557

Spitzer, L. 1978, Physical Processes in the Interstellar Medium (New York: Wiley) doi:10.1002/9783527617722

Squire, J., Hopkins, P. F., Quataert, E., & Kempski, P. 2021, MNRAS, 502, 2630

Tagger, M., Falgarone, E., & Shukurov, A. 1995, A&A, 299, 940

Thomas, T., & Pfrommer, C. 2019, MNRAS, 485, 2977

van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394

Wentzel, D. G. 1969, ApJ, 156, 303

Wentzel, D. G. 1974, ARA&A, 12, 71

Wiener, J., Oh, S. P., & Guo, F. 2013, MNRAS, 434, 2209

Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906

Xu, S., Yan, H., & Lazarian, A. 2016, ApJ, 826, 166

Yan, H., & Lazarian, A. 2011, ApJ, 731, 35

Yang, H. Y. K., Ruszkowski, M., Ricker, P. M., Zweibel, E., & Lee, D. 2012, ApJ, 761, 185

Zirakashvili, V. N., Breitschwerdt, D., Ptuskin, V. S., & Voelk, H. J. 1996, A&A, 311, 113

Zweibel, E. G. 2013, PhPl, 20, 055501

Zweibel, E. G. 2017, PhPl, 24, 055402

Zweibel, E. G., & Shull, J. M. 1982, ApJ, 259, 859