Influence of Ion-Neutral Damping on the Cosmic-Ray Streaming Instability: Magnetohydrodynamic Particle-in-cell Simulations
DDraft version February 25, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Influence of Ion-Neutral Damping on the Cosmic Ray StreamingInstability: MHD-PIC Simulations
Illya Plotnikov, Eve C. Ostriker, and Xue-Ning Bai
3, 4 IRAP, Universit´e de Toulouse III - Paul Sabatier, OMP, Toulouse, France Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton, NJ 08544, USA Institute for Advanced Study, Tsinghua University, Beijing 100084, China Department of Astronomy, Tsinghua University, Beijing 100084, China
ABSTRACTWe explore the physics of the gyro-resonant cosmic ray streaming instability (CRSI) including theeffects of ion-neutral (IN) damping. This is the main damping mechanism in (partially-ionized) atomicand 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´en waves that scattercosmic rays and control galactic-scale transport. Our study employs the MHD-PIC hybrid fluid-kineticnumerical technique to follow linear growth as well as post-linear and saturation phases. During thelinear phase of the instability – where simulations and analytical theory are in good agreement – INdamping prevents wave growth at small and large wavelengths, with the unstable bandwidth lower forhigher ion-neutral collision rate ν in . Purely MHD effects during the post-linear phase extend the wavespectrum towards larger k . In the saturated state, the cosmic ray distribution evolves toward greaterisotropy (lower streaming velocity) by scattering off of Alv´en waves excited by the instability. In theabsence of low- k waves, CRs with sufficiently high momentum are not isotropized. The maximumwave amplitude and rate of isotropization of the distribution function decreases at higher ν in . Whenthe IN damping rate approaches the maximum growth rate of CSRI, wave growth and isotropization issuppressed. Implications of our results for CR transport in partially ionized ISM phases are discussed. Keywords: instabilities — waves — scattering — ISM: magnetic fields — ISM: cosmic rays — methods:numerical INTRODUCTIONCosmic Rays (CRs) are a key component of the inter-stellar medium (ISM). They are potentially importantto ISM dynamics including support against gravity, astheir energy density is typically comparable to that ofthe magnetic field as well as the thermal and turbulentenergy of the ISM gas (e.g., Spitzer 1978; Ferri`ere 2001;Draine 2011; Grenier et al. 2015). Perhaps even moreimportant dynamically, given the large scale height ofthe CR distribution, is the potential for CRs to con-tribute in driving galactic winds (e.g. Ipavich 1975; Bre-
Corresponding author: Illya Plotnikov, Eve C. Ostriker, [email protected]@[email protected] itschwerdt et al. 1991; Zirakashvili et al. 1996; Everettet al. 2008; Mao & Ostriker 2018). CRs are also crucialto ISM microphysics, driving ionization and dissociationas well as providing the primary heating of gas in regionsshielded to UV (e.g. Draine 2011; Glassgold et al. 2012).The coupling between CRs and the thermal ISMgas occurs through scattering of CRs off of magneticwaves, preexisting in the turbulent ISM or self-generatedby streaming CRs (Kulsrud & Pearce 1969; Wentzel1969; Skilling 1971; Kulsrud 2005; Amato & Blasi2018). Higher-energy CRs may interact primarily withexternally-generated waves (Blasi et al. 2012a). How-ever, particles of energy (cid:46)
GeV dominate by numberand energy the overall CR content. At the micro-pcscales comparable to the gyro-radius of GeV protons,the energy density of the ISM turbulence is too smallto provide efficient scattering (under the assumption ofa cascade from directly-observed turbulence at largerscales), and these particles are believed to be mainly a r X i v : . [ a s t r o - ph . H E ] F e b I. Plotnikov, E. Ostriker & X.-N. Bai scattered by self-excited waves. Alfv´en waves generatedby streaming therefore bear primary responsibility forscattering the dominant portion of the CR distributionin our own and other galaxies (e.g. Zweibel 2013).There is an increasing necessity to provide satisfac-tory microphysical understanding of the CR streaminginstability and its implications for CR transport in dif-ferent phases of the ISM. The possible dynamical role ofCRs in driving galactic winds has recently led to recon-sideration of the traditional assumption that streamingoccurs at the Alfv´en speed in several analytical studies(Wiener et al. 2013; Recchia et al. 2016; Zweibel 2017).In parallel, a number of groups have implemented CRsas a fluid in numerical MHD codes (e.g., Yang et al.2012; Dubois & Commer¸con 2016; Pakmor et al. 2016;Pfrommer et al. 2017; Thomas & Pfrommer 2019; Jiang& Oh 2018; Hopkins et al. 2020) in order to study theeffect of CRs on galactic wind generation (e.g., Hanaszet al. 2013; Girichidis et al. 2016, 2018; Pfrommer et al.2017; Wiener et al. 2017; Ruszkowski et al. 2017; But-sky & Quinn 2018; Dashyan & Dubois 2020), on theglobal evolution of supernova remnants (Pais et al. 2018;Dubois et al. 2019) or on the multiphase ISM structure(Bustard & Zweibel 2020). Global simulations testingdifferent models of CR transport coefficients were re-cently conducted by Hopkins et al. (2021a,b), suggestingthat standard models of CR transport are inconsistentwith observed constraints on large scales. At the sametime, observations that probe ionization and chemistryin ISM clouds at small scales also place constraints onCR transport (see review of Padovani et al. 2020, andreferences therein), e.g. suggesting diffusive rather thanfree-streaming behavior in the outer layers of molecularclouds (Silsbee & Ivlev 2019).Fluid implementations of CRs rely on subgrid pre-scriptions for diffusion and streaming. To date, thesesubgrid treatments have relied on empirical estimatesof diffusivities, idealized streaming treatments, or sim-ple models based on analytic predictions for growth anddissipation of waves. Considering the important role ofmicrophysics to CR transport, it is valuable to pursuea deeper investigation that addresses the processes ofwave growth, damping, and particle-wave interactionsdirectly. The most important parameter to be obtainedfrom microphysical studies of CR-ISM interactions is thescattering rate of CRs off of Alfv´en waves, since this en-ters in determining the diffusivity in traditional treat-ments (e.g. Skilling 1971; Hanasz & Lesch 2003), or therate 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, ν s , is proportional to the wave energy density in the reso-nant range. The latter quantity depends on the outcomeof the cosmic-ray streaming instability (CRSI) in the lo-cal medium. There are different classes of the CRSI,resonant (Lerche 1967; Kulsrud & Pearce 1969; Skilling1971; Ginzburg et al. 1973; Wentzel 1974; Skilling 1975; Berezinskii et al. 1990) and non-resonant (Bell 2004;Amato & Blasi 2009; Bykov et al. 2013). The latter isimportant only in the near environment around accelera-tors (such as supernova remnant shocks) where the elec-trical current of CRs is large, i.e. (4 π/c ) J CR R L, /B (cid:29)
1, where J CR = en CR V D is the electrical current of CRsstreaming at speed V D and with typical gyroradius R L, (Amato & Blasi 2009). In the general ISM, the resonantinteraction prevails since the density of CRs is very small( n CR ∼ − − − cm − ). Because we are interestedin ambient ISM conditions, here we focus exclusively onthe resonant instability.Wave growth, triggered by the anisotropy of the dis-tribution function or CR density gradient, is generallyin competition with different damping mechanisms de-pending on the phase of the ISM. In a partially ionizedmedium the most important wave damping mechanismis from the collisions between ions and neutrals. Thecompetition of CRSI with ion-neutral (IN) damping wasstudied in context of molecular clouds (e.g., Kulsrud &Cesarsky 1971; Zweibel & Shull 1982; Everett & Zweibel2011; Morlino & Gabici 2015; Ivlev et al. 2018) or inthe context of particle acceleration at shocks (e.g., O’CDrury 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 ion-neutral collisions can be important in hotter and morediffuse phases of the ISM. Other mechanisms thathave been discussed include non-linear Landau damping(NLLD Lee & V¨olk 1973; Kulsrud 1978) and damping byinteractions with turbulence (Farmer & Goldreich 2004;Yan & Lazarian 2011; Lazarian 2016). We also notethe recently studied damping by charged dust grainsthat can be important for CRs with energies < E <
10 TeV. Here, weconsider only these primarily-neutral phases of the ISM;investigation of effects of alternative damping mecha-nisms is deferred to future work.In the present work, we investigate the interplay be-tween CRSI and ion-neutral damping during the linearphase of the CRSI, and also assess the behavior of thesystem in post-linear and late-time saturated stages ofevolution. We are motivated to obtain better insightinto CR transport in ISM conditions which are largelyneutral (this makes up most of the ISM mass), where theion/neutral fraction is controlled by photo-ionization orionization by low-energy CRs.Our study of CRSI relies on the MHD-PIC numericalmethod described in Bai et al. (2015) (see also, Re-ville & Bell 2012; van Marle et al. 2018; Amano 2018;Mignone et al. 2018; Lebiga et al. 2018). This methodis well adapted to capture the resonant nature of the osmic ray streaming with ion-neutral damping n CR /n i is extremely small in the parameter rangeof interest. The gyro-radius of CRs must be resolved,but smaller scales (such as ion skin depth and electronscales) not need to be resolved. We note several recentstudies have alternatively investigated certain aspectsof the CRSI with fully-kinetic (Holcomb & Spitkovsky2019; Shalaby et al. 2020) or hybrid-kinetic approaches(Haggerty et al. 2019; Schroer et al. 2020) that useparticle-in-cell (PIC) methods. These approaches havethe advantage of resolving small-scale phenomena onelectron- (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 thethermal ions are treated via PIC.In this work, we follow up on Bai et al. (2019) (here-after Paper I), where the CRSI was studied using theMHD-PIC approach. There, the main findings include:(i) the linear phase of the instability for both reso-nant and non-resonant branches can be accurately repro-duced with our “ δf ” MHD-PIC method; (ii) the quasi-linear diffusion (QLD) formalism accurately describestemporal changes in the CR DF, except near pitch an-gle of 90 degrees (the “ µ = 0 crossing” problem); (iii)crossing of µ = 0 is mainly due to nonlinear wave-particle interactions; (iv) the Alfv´en wave amplitudein saturation reflects the expected transfer of net mo-mentum from the originally-anisotropic CRs to forward-propagating Alfv´en waves. Here, we adopt the same nu-merical methods and conduct simulations in a similarparameter regime to Paper I, but now we additionallyconsider the effects of ion-neutral interactions that dampwaves, competing with CR-induced wave generation.This article is organized as follows. In section 2, we re-view the analytical derivation of the CRSI linear growthrate, including the contribution from ion-neutral damp-ing. In section 3 we outline our numerical methods in-cluding the implementation of ion-neutral wave damp-ing. Simulation results are presented in section 4. Wefurther discuss our findings and their astrophysical con-text in section 5, and summarize our main conclusionsin section 6. COSMIC RAY STREAMING INSTABILITYINCLUDING ION-NEUTRAL INTERACTIONIn this section, we discuss the linear growth rate of theCRSI without and with ion-neutral damping. While an-alytic studies of CRSI commonly adopt a power-law dis-tribution, here we will instead consider a κ -distribution,as adopted in our simulations. This choice is motivatedby the use of δf method, in which the unperturbed partof the CR distribution must be specified for all p , andshould be a smooth function (see Bai et al. 2019). Be-cause standard expressions for CRSI with ion-neutralcollisions are based on power-law distributions, here were-derive growth rates for the case of the κ − distribution. These analytic results are useful for comparison with thesimulations that follow.Some discussion of the high-frequency regime is givenin the last part of this section.2.1. No damping
Resonant contribution only
We consider the instability arising from resonant in-teraction in the low-frequency, ω (cid:28) Ω , and non-relativistic drift ( V D (cid:28) c ) limit . Here, ω = kV A isthe Alfv´en wave frequency and Ω = eB / ( m p c ) is thenon-relativistic proton gyro-frequency. We adopt stan-dard notation, with Alfv´en velocity V A = B / (4 πρ i ) / for B the local magnetic field strength and ρ i the back-ground ion density; e is the proton charge, m p is theproton mass, and c is the speed of light. We also as-sume that the mass and charge of the background fluidions are identical to those of CR particles : m i = m p and q i = e .We start with the growth rate expression as given by(Kulsrud 2005, Eq. 69 of Ch. 12) for wavenumber k :Γ CR ( k ) = − π e (cid:18) V A c (cid:19) (cid:18) V D V A − (cid:19) × (cid:90) ∂F∂p p ⊥ p δ ( kp (cid:107) − m p Ω )d p , (1)where V D is the initial drift speed between the back-ground gas and the frame in which the cosmic rays havean isotropic distribution (taken along the external mag-netic field in the x -direction). The δ function in theintegral accounts for the gyro-resonant interaction atthe fundamental harmonic, where resonance occurs at p (cid:107) = m p Ω /k ≡ p res .For the (isotropic) distribution F ( p ) of CRs in thedrift frame, we adopt a κ -distribution, F ( p ) = A (cid:34) κ (cid:18) pp (cid:19) (cid:35) − ( κ +1) . (2)The normalization factor A is given by A = n CR ( πκp ) / Γ( κ + 1)Γ( κ − / , (3)where n CR is the number density of cosmic rays, andΓ( x ) is the Euler Γ-function. We generally adopt κ = In this study we adopt the notation Ω for the quantity denotedby Ω c in Paper I, i.e. Ω ≡ Ω c . If this condition is not satisfied, one has to replace Ω byΩ i = q i B / ( m i c ) in the final expression of the growth rate,i.e., in Equations 7, 8, 9, 11, 23. We note however, that ifmass densities are used instead of number densities then thecharacteristic frequency appearing in the growth rate would cor-respond to the gyrofrequency of CRs, owing to the equality( n CR /n i )Ω i = ( ρ CR /ρ i )Ω I. Plotnikov, E. Ostriker & X.-N. Bai .
25, corresponding to F ( p ) ∝ p − . for p (cid:29) p . Wenote that a κ − distribution reduces to a Maxwellian inthe limit κ → ∞ .Equation 1 becomesΓ CR = 2 π e A κ + 1 κp (cid:18) V A c (cid:19) (cid:18) V D V A − (cid:19) × (cid:90) (cid:34) κ (cid:18) pp (cid:19) (cid:35) − ( κ +2) δ ( kp (cid:107) − m p Ω ) p ⊥ d p (cid:124) (cid:123)(cid:122) (cid:125) I . (4)With d p = d θ d p (cid:107) p ⊥ d p ⊥ , the integral I in the previousequation becomes I = 2 πk (cid:90) ∞ (cid:20) κ p + p ⊥ p (cid:21) − ( κ +2) p ⊥ d p ⊥ . (5)Noting that p ⊥ = p − p and that p ⊥ d p ⊥ = p d p , wechange the variable of integration from p ⊥ to p . Thisleads to I = 2 πk (cid:90) ∞ p res (cid:20) κ p p (cid:21) − ( κ +2) p ( p − p )d p = πk p κκ + 1 (cid:20) κ p p (cid:21) − κ . (6)Using ρ i = m p n i , R L, ≡ p / ( m p Ω ), p res /p =1 / ( kR L, ) and inserting the expression from Equation 3the growth rate can be written asΓ CR ( k ) = π / κ / Γ( κ + 1)Γ( κ − / n CR n i Ω kR L, × (cid:18) V D V A − (cid:19) (cid:20) κ kR L, ) (cid:21) − κ . (7)The growth rate is maximal at kR L, = (cid:112) − /κ .For the adopted value in this study κ = 1 .
25, it isΓ max , (cid:39) . n CR n i (cid:18) V D V A − (cid:19) . (8)The subscript “max , (cid:48)(cid:48) stands for the maximum growthrate of the instability without IN damping. The nu-merical prefactor varies from 0 .
25 for κ = 1 to 0 .
38 for κ → ∞ . In the long-wavelength limit ( kR L, (cid:28)
1) thegrowth rate varies as ∝ k κ − , and in short-wavelengthlimit ( kR L, (cid:29)
1) as ∝ k − .2.1.2. Full dispersion relation
In more general cases when CR-induced current canbe large, the contribution from non-resonant instabilitycan be significant (see, e.g., Bell 2004; Amato & Blasi2009; Bykov et al. 2013). The full dispersion relation using the κ − distribution was derived in Paper I. Here,we reproduce it for completeness: ω = k V A ∓ n CR n i Ω ( ω − kV D ) ((1 − Q ) ± iQ ) ;(9)here the (1 − Q ) and Q terms are due to non-resonantand resonant responses of CRs to waves, respectively(see Eq. 9 in Paper I for definition of these terms). Theresonant term given in Eq. 40 of Paper I is: Q = √ πκ / Γ( κ + 1)Γ( κ − /
2) 1 kR L, (cid:2) / ( κk R L, ) (cid:3) − κ . (10)In the case with negligible contribution from the non-resonant term ((1 − Q ) = 0), the imaginary part of ω obtained by solving the Equation 9 is then ω I = 12 n CR n i Ω (cid:18) V D V A − (cid:19) Q , (11)which is identical to the growth rate given in Equation 7.2.2. With ion-neutral damping
Damping term from fluid equations
In this section, we discuss the Alfv´en wave propaga-tion properties in presence of neutrals, temporarily ig-noring the destabilizing contribution from CR stream-ing. The derived damping rate is then subtracted fromthe CRSI growth rate given above to obtain the netgrowth rate. Consider coupled two-fluid momentum equations forions and neutrals, without the contribution of CRs, andwhere Ohmic dissipation and the Hall effect are ne-glected: ρ i D u i Dt = −∇ P i + J c × B − ρ i ν in ( u i − u n ) (12) ρ n D u n Dt = −∇ P n − ρ n ν ni ( u n − u i ) , (13)where ρ i ( ρ n ), P i ( P n ), ν in ( ν ni ) are the mass density,pressure, and collision frequencies of ions with neutral(neutrals with ions), respectively. The quantities B and J correspond to the magnetic field and current carriedby the plasma. The ions are subject to both the Lorentzforce and friction with neutrals, while neutrals are onlysubject to friction with ions.From momentum conservation, ρ i ν in = ρ n ν ni = n n n i µ (cid:104) σv (cid:105) , where µ = m i m n / ( m i + m n ) is the reducedmass, and (cid:104) σv (cid:105) is the momentum transfer rate coeffi-cient. Most relevant to wave damping is the ion-neutralcollision frequency ν in = ρ n m i + m n (cid:104) σv (cid:105) . (14) Such a procedure is typically adopted in the literature (Kulsrud& Cesarsky 1971; Ginzburg et al. 1973; Zweibel & Shull 1982;Nava et al. 2016) osmic ray streaming with ion-neutral damping (cid:104) σv (cid:105) ≈ . × − cm s − for neutral atomic gaswhere the main collision partner is H with µ ≈ . m p for H + ions, and (cid:104) συ (cid:105) ≈ . × − cm s − in molecularregions where the collision partner is H with µ ≈ m p for C + ions in diffuse gas or HCO + ions in dense gas.The magnetic field in Equations 12-13 evolves subjectto the induction equation: ∂ B ∂t = ∇ × ( u i × B ) (15)For simplicity, we consider a static, homogeneous andincompressible background medium. The mean mag-netic field is oriented along the x -direction, B = B e x and the perturbation components are perpendicular to B : δ B = δB y e y + δB z e z . The velocity and magneticfield perturbations are ∝ exp[ i ( kx − ωt )]. LinearizingEquations 12, 13 and 15 leads for ions to ω δ u i = k V A,i δ u i − iων in ( δ u i − δ u n ) , (16)and for neutrals to ω δ u n = − iων ni ( δ u n − δ u i ) . (17)Here, the Alfv´en velocity includes the ion density only,i.e., V A,i = B / √ πρ i . Combining these two equationsgives the dispersion relation: ω (cid:18) iν in ω + iν in ρ i /ρ n (cid:19) = k V A,i . (18)One identifies here the Alfv´en wave dispersion rela-tion (in the ion fluid) modified by the presence of neu-trals (second term inside brackets on the left-hand side).Equation 18 is a cubic equation for ω . Despite existing,general solutions are quite cumbersome, hence we do notprovide them here. More physical insight can be gainedby considering different limits, previously discussed inthe 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. Wedefine Z = ρ i /ρ n as the ratio of ionized to neutral massdensity. It is related to the ionization fraction as x i ≡ n i /n n = Z/q , where q = m i /m n . Three limiting casesare then well defined:1. Ion-dominated, Z (cid:29)
1. Here, the Alfv´en wavevelocity is unmodified by neutrals, ω = kV A,i . Thedamping rate is, for any k and ω ,Γ d = ν in k V A,i k V A,i + ν Z (19)For ω k = kV A,i (cid:29) ν ni = ν in Z it is equal to ν in / Z (cid:28)
1) and low-frequency( kV A,i /ν in (cid:28) √ Z ). The fluids are strongly cou-pled, hence the Alfv´en wave velocity includes the sum of ion and neutral densities: ω = kV A,tot = kB / (cid:112) π ( ρ i + ρ n ). The damping rate isΓ d ≈ k V A,i ν in (20)3. Neutral-dominated ( Z (cid:28)
1) and high-frequency( kV A,i /ν in (cid:29) ω = kV A,i . The damping rate isΓ d ≈ ν in √ Z A,i /ν in < / ω + iν in ω = k V A,i . (22)This implicitly means that we will study the regimewhere neutrals are fully decoupled from ions, that is thecase in high-frequency limit where the Alfv´en wave is notmodified by neutrals ω R = kV A,i and the damping rateis asymptotic at its highest value: − ω I = Γ d = ν in / General dispersion relation with CR streamingcontribution It is straightforward to generalize the dispersion rela-tion of the CRSI with ion-neutral interactions (see, e.g.,Reville et al. 2007, 2021). Following the procedure ofPaper I but with the addition of ion-neutral momentumexchange terms in fluid equations, we obtain: ω (cid:18) iν in ω + iν in ρ i /ρ n (cid:19) = k V A,i ∓ n CR n i Ω ( ω − kV D ) × ((1 − Q ) ± iQ ) , (23)where (1 − Q ) and Q from non-resonant and resonantresponses of CRs to waves are discussed above, and Q is given for a κ -DF in Equation 10.Comparing to Equation 9, the supplementary termdue to ion-neutral interaction is the second term in-side the parentheses on the left-hand side of the Equa-tion 23, while the second term on the right-hand side isdue to the CR streaming. For sufficiently small n CR /n i I. Plotnikov, E. Ostriker & X.-N. Bai -6 -3 -6 -3 Figure 1. Dispersion relation of Alfv´en waves and dampingrate in presence of neutrals. The two panels present the realand imaginary parts of the wave frequency (blue and redsolid lines, respectively) in ion-dominated case (top panel,analytical and numerical solution are identical) and neutraldominated case (bottom panel). The black dotted line inthe bottom panel corresponds to an approximate functionthat fits the damping rate reasonably well for any k where asolution exists. (as is applicable in general), the solution of this disper-sion relation is the same as if the damping (ion-neutral)and growing (CRSI) imaginary terms were derived sep-arately, with the two terms added to obtain the netgrowth rate. The growth rate from resonant interactiononly was presented in the Section 2.1.1 and the dampingrate in Section 2.2.1. Taking the high-frequency limit,the net growth rate isΓ tot ( k ) = Γ CR ( k ) − ν in , (24)where Γ CR ( k ) is given by Equation 7 for a κ -DF. Thisequation is valid for any kV A,i (cid:29) max[ ν in , ν ni ]. Notethat in the low-frequency (long-wavelength) regime, oneneeds to distinguish between high- and low-ionizationfraction cases, which we do not consider here. 2.3. ISM parameters and implications As discussed in Section 2.2.1, the one-fluid formulationdoes not capture all the details of wave damping. Forthis reason we discuss here the conditions for which theone-fluid approach is accurate.We assume typical ISM values B = 5 µ G, n CR =10 − cm − , and typical energy of CRs around E ∼ GeV. The frequency of the fastest growing wavemode driven by GeV protons is defined as ω max =1 . V A,i /R L, , since the fastest growing wavenumber is k max = 1 . R − L, (see text below Equation 8). The typi-cal values for the ion-neutral momentum exchange rate ν in in different phases of the ISM were discussed pre-viously in Section 2.2.1. The high-frequency regime isvalid for ω max /ν in (cid:29) 1. Comparing the two frequen-cies in different phases of the ISM (as done in Table 1,where values for different phases are taken from Draine2011, Chap. 16) we obtain that this inequality is largelysatisfied in all neutral-dominated phases.Another important consideration concerns the pos-sibility of wave growth in presence of IN damping.The CRSI-driven waves can only grow if the maximumgrowth rate of the CRSI is larger than the IN dampingrate: Γ max , / Γ d > 1. Otherwise, even the fastest grow-ing mode is damped. For this purpose, we define a ζ parameter as ζ = Γ max , Γ d V D /V A,i − V D /V A,i − ζ is independentof V D . CRSI can be triggered at a range of wavelengthsprovided that V D /V A,i − > /ζ . Thus, if ζ > 1, theCRSI can be triggered for super-Alfv´enic drift velocity V D /V A,i − ζ is typically in the range of 0.1-20 (exceptin the DMG), Thus in these environments, and con-sidering only IN damping, the CRSI can be triggered forstreaming velocity of 0 . (cid:46) v D /V A,i − (cid:46) 10. But bythe same token, wave damping will significantly limit therange over which the instability can be directly excitedif v D /V A,i − /ζ .2.4. Wave saturation and steady-state streamingvelocity When CRSI is efficiently triggered, scattering ofCRs off of self-generated waves gradually reduces theanisotropy of the CR DF, lowering the drift velocityand the growth rate. The limitation of CR transport For the DMG phase, the required initial drift velocity of CRs isextremely high, V D > . × V A,i = 0 . c , implying that waveexcitation by CRSI is extremely challenging in this medium un-less a powerful accelerator is present at close distance, increasinglocally the CR density well above n CR = 10 − cm − . osmic ray streaming with ion-neutral damping Table 1. List of typical parameter values for different phasesof the ISM.Phase n m i m n x i ν in ζ (cm − ) ( m p ) ( m p ) ( ω max )WNM 0.3 1 1.4 10 − . × − − . × − 29 2.3 10 − . × − 29 2.3 10 − . × − . × − WNM: Warm Neutral Medium, CNM: Cold NeutralMedium, MG: Molecular Gas , DMG: Dense Molecular Gas.Reported quantities are the total density n , ion mass m i ,ionization fraction x i , corresponding ω HF above which thehigh-frequency regime is satisfied, and the drift-normalizedinstability parameter ζ , that is defined in Equation 25. Weadopted B = 5 µ G, n CR = 10 − cm − , and p = m p c corre-sponding to CR energy E (cid:39) GeV, as fiducial values. by self-generated Alfv´en waves is commonly referred toas self-confinement of CRs.In the absence of wave damping, the DF anisotropy(in the wave frame) is eventually erased and the waveamplitude saturates. Considering that the momentumlost by CRs is transferred to waves, the saturation levelcan be estimated approximately as (e.g., Kulsrud 2005;Bai et al. 2019; Holcomb & Spitkovsky 2019) (cid:18) δBB (cid:19) ∼ n CR n i (cid:18) V D V A,i − (cid:19) . (26)The wave intensity at saturation is linearly proportionalto the number density of CRs and to the first-orderanisotropy of the distribution function.In the presence of IN damping (or other damping),a commonly adopted assumption is that the streamingvelocity is such that the linear growth rate of CRSI bal-ances the wave damping rate (e.g., Kulsrud & Cesarsky1971; Krumholz et al. 2020; Hopkins et al. 2021a,b). Ifwe consider the fastest-growing waves and equate to zerothe net growth rate with IN damping, Γ max , − ν in / κ = 1 . 25) is V st V A,i − (cid:39) . n i n CR ν in Ω i . (27)Here, the use of Ω i instead of Ω in Section 2.1.1 isrequired because the ionized fluids have m i > m p inmany of the relevant environments (see Table 1).The value in Equation 27 sets, in principle, a lowerlimit on CR streaming velocity when IN damping ispresent. However, even if waves near the fastest-growingwavenumber can grow, waves at much shorter or longerwavelengths would be damped. As a result, particles with very large (or very small) momentum that are res-onant with smaller (or larger) k would not experiencemuch scattering. Thus, the momentum-weighted valueof V st will always be larger than the one given in Equa-tion 27.We note that the conventional approach describedabove considers only linear growth and damping ofwaves. In reality, both linear growth and pitch anglescattering must be considered, in opposition to wavedamping. For a given momentum, resonant waves mustbe excited and survive IN damping over the timescalerequired for pitch angle diffusion, in order for isotropiza-tion to occur. Quantifying this requires comparison ofthe isotropization time (after the linear phase), 1 /ν s ,with the damping time, 1 /ν in .In the absence of damping, the saturation amplitudeis given by Equation 26, and since the scattering rateunder QLD is ν s ∼ ( π/ ( δB/B ) (Kulsrud 2005),comparing to Equation 8 shows that ν s ∼ Γ max , . Ifwe were to assume that full redistribution in µ requiresa time ∼ ν − s , the total time for isotropization (con-sidering wave growth and scattering) would be at least2Γ − , ; in practice the numerical results described be-low (see Figure 2) correspond to a prefactor ∼ 10 forthe isotropization time. With IN damping, the growthrate is reduced, and (as we shall show) the maximumlevel of ( δB/B ) is also reduced, which slows scatter-ing. One might therefore expect that for isotropizationto be successful, the ratio Γ max , / ( ν in / 2) must be abovesome critical value which is larger (perhaps much larger)than 1, in order to allow for both wave growth and par-ticle diffusion. We return to this issue in subsection 4.4. NUMERICAL SIMULATIONS USING MHD-PICAPPROACHWe use the MHD-PIC method introduced by Bai et al.(2015) and adopted for the study of CR-streaming inPaper I. The focus of the present study is the inclusionof ion-neutral wave damping. Let us recall key elementsof this method: • The background thermal plasma dynamics are de-scribed using MHD equations. The dynamics ofCRs are solved by a PIC method, where CRs aretreated as charged macro-particles. The Lorentzequation is solved using a Boris pusher, the den-sity and electric current are deposited on grid us-ing 2nd order shape factors (TSC scheme), and thecoupling between CRs and the background fluid isachieved by introducing CR source terms into theMHD equations and modifying Ohm’s law. • The δf method is employed to significantly reducethe Poisson noise of macro-particles • Phase scrambling is applied when particles crossthe system boundary and re-enter on the other I. Plotnikov, E. Ostriker & X.-N. Bai side of the periodic box; this effectively mimics alarger numerical box. • Seed waves are initialized at t = 0. This providesbetter control of k -by- k growth rate as comparedto growth from numerical noise (typically done infull-PIC simulations).3.1. Numerical treatment of the ion-neutral damping The presence of neutrals introduces momentum ex-change between ions and neutrals through collisions.The corresponding terms in the fluid equations weregiven previously in Equation 12 and Equation 13. In-stead of solving the coupled system of two-fluid equa-tions, we only account for the effect of neutrals on theions while the dynamical equations for neutrals are notevolved. This procedure is consistent in the limit wherethe two fluids are decoupled (high-frequency regime).The procedure is as follows. At each time step wereduce the transverse momentum fluctuations accord-ing to Equation 12. After each numerical time step∆ t the transverse momentum is updated as δu new = δu exp( − ν in ∆ t ), in order to account for the effect ofdamping. While the longitudinal velocity should be sub-ject to the same “damping”, it is not incorporated as itis decoupled from the Alf´ven waves (which have no lon-gitudinal motion) and we have verified that it does notaffect the overall simulation results.3.2. Numerical setup In the present work, there are some differences com-pared to the Paper I fiducial setup. Namely:(i) We use fewer CR macro-particles per cell, as the δf method effectively controls the Poissonian noise level.(ii) The physical size of the grid is slightly smaller. Yet,the box size is equal to (cid:39) 47 times the most unstablewavelengths.(iii) We include the wave damping term in the MHDmomentum equation.(iv) The fiducial drift velocity is 10 V A,i (instead of 2 V A,i used in Paper I). This choice is motivated by the need toefficiently reach a fully saturated state when there is noion-neutral damping. This case will be used as referencewhen measuring the effect of different levels of damping.Table 2 summarizes numerical parameters and modelvalues adopted in our simulations. Our fiducial choice is n CR /n i = 10 − . This is somewhat larger than realisticvalues in the neutral ISM ( n CR /n i ∼ × − − − , seeTable 1), for numerical expediency; lower values wouldhave very low CRSI growth rate (requiring very longsimulation duration) and low saturation amplitudes (ex-acerbating the numerical issue of crossing µ = 0 for prac-tical resolution). This choice does not affect any of ourtheoretical conclusions. The box length is chosen to bemuch larger than the fastest growing wavelength of theinstability. Lengths are given in units of ion skin depth d i ≡ C /ω p,i = V A,i / Ω , which is 1 in code units. Weadopt a reduced speed of light C = 300 V A,i . Seed wavesare initialized with equal amplitudes of right-handed andleft-handed polarization and the same amplitude at all k . The initial total power in the waves is δB/B = 10 − .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 ν in rel-ative to the peak growth rate for CRSI; this is what isreported in Table 2. Using V D = 10 V A,i and n CR /n i =10 − in Equation 8, this initial peak growth rate (for ν in = 0) is Γ max , / Ω = 2 . × − . SIMULATION RESULTSIn this section we describe the results of MHD-PICsimulations. We start with a presentation of the overalltime-evolution, then show the linear growth rates of theinstability. We then discuss the transition into saturatedphase and some aspects of the saturated state of thestreaming CR-background fluid system.4.1. Overall evolution In this section we compare the evolution betweenno-damping and damping cases (models Fid and Fid-damp reported in Table 2). For the cases with damp-ing we explored values of ν in between 0 . max , and1 . max , ; the theoretical critical value for no growthis ν in = 2Γ max , . In the following, we use the term“moderate damping” for cases where ν in / Γ max , is non-negligible but still smaller than unity. The representa-tive case is ν in / Γ max , = 0 . ν in . The undampedcase (blue curves) reaches saturation at t Ω (cid:39) with δB /B (cid:39) − . The CR streaming velocity, initially V D = 10 V A,i , begins to decrease when δB /B (cid:39) − and reaches v rel (cid:39) V A,i by the end of the simulation.Here, the CR bulk velocity is defined (in the simulationframe) by v rel ≡ (cid:82) f ( p, µ ) v ( p ) p µ d µ dp (cid:82) f ( p, µ ) p d µ dp − V bg , (28)where V bg the background gas speed. Unlike the initialconditions, the CR distribution is not perfectly isotropicin the frame moving at v rel .Cases with non-zero damping reach lower saturationlevels, with the maximum in δB /B decreasing as ν in increases. After saturation (at late simulation times, t Ω > ) the wave intensity stabilizes at intensitythat decreases with increasing ν in 5 . Streaming ve- We noticed that the level of the wave intensity ‘plateau’ at thelate-time state can evolve somewhat differently depending on nu-merical resolution and number of particles per cell, but does notaffect qualitatively the system evolution. osmic ray streaming with ion-neutral damping Table 2. List of main simulation runsRun V D /V A,i n CR /n i ν in / Γ max , Domain size Domain size resolution N p Runtime L x ( d i ) L x /λ m ∆ x ( d i ) (per cell) (Ω − )Fid 10.0 1 . × − × Fid-Damp 10.0 1 . × − ∈ [0 . , . 9] 8 × HiRes 10.0 1 . × − { , . , . } × . < Fixed parameters: C /V A,i = 300, p / ( mV A,i ) = 300, κ = 1 . 25, and initial wave amplitude A = 10 − . In all models the mostunstable wavelength is 0 . λ m for λ m = 2 πp / ( m Ω ) ≈ d i . The ion skin depth d i ≡ C /ω i = V A,i / Ω . The quantity Γ max , is defined in Equation 8 as the maximum growth rate of the CRSI without IN damping (i.e. ν in = 0). -8 -7 -6 -5 -4 -3 -2 Figure 2. Time evolution of the magnetic wave energy (toppanel) and of the streaming velocity of bulk CRs (bottompanel) for different values of ν in ranging from 0 (blue line)to 1 . max , (black line). locities (bottom panel) decline at lower rates for thestronger-damping models, because the lower wave am-plitudes scatter CRs less effectively. The late-time (after t Ω ∼ ) evolution is asymptotically slow.In general, we define three representative phases oftime-evolution of the instability: linear, post-linear, andsaturated (late-time). There is an adjustment phase be-tween the time when the fastest modes stop growing -0.0500.05-0.0500.05-0.0500.05 25 30 35 40-0.0500.05 t = 2.0e+04t = 3.6e+04t = 1.0e+05 25 30 35 40 t = 1.0e+06 Figure 3. Profile of the wave component B z at differentsimulation times for the case with no damping (left column)and with moderate damping ν in = (1 / max , (right col-umn). Four different times are presented. The top set corre-sponds to the early linear phase. The two middle panel setsshow the evolution during post-linear phase, and the bottomset shows the final simulation time t Ω = 10 . (end of linear phase) and the saturation of the instabil-ity when the anisotropy of the distribution function is(eventually) erased (start of saturated phase). We referto this transitory phase as post-linear .In Figure 3 we show the magnetic field profile B z at four different simulation times for the case with nodamping (left column) and moderate damping (right col-umn). For the case with no damping the wave ampli-0 I. Plotnikov, E. Ostriker & X.-N. Bai Figure 4. 2D distribution function in the wave frame, δf w ( p w , cos θ w ) /f , at different simulation snapshots for thecase with no damping (left column) and with moderatedamping ν in = (1 / max , (right column). The same simu-lation times are presented as in Figure 3 from top to bottom. tude grows to δB z /B (cid:39) . 02 at the end of the linearphase ( t Ω = 2 × ) and further increases to 0 . 05 dur-ing the post-linear evolution (between t Ω (cid:39) × and t Ω = 10 ), during which amplitudes increase at k belowthe peak (see Section 4.3). For the case with moderatedamping, the wave amplitude is reduced, as expected.During the post-linear phase we observe a decrease ofthe wave amplitude from δB z /B (cid:39) . 02 to 4 × − ,while there is less contribution from k below the peak.The effect of CRSI and subsequent QLD on the CRdistribution function is presented in Figure 4. We showthe distribution function in the wave frame (i.e. mov-ing to the right at V A , i with respect to the gas, seePaper I for its transformation from simulation frame), δf w ( p w , cos θ w ) /f for the case with no damping (leftcolumn) and with moderate damping (right column), atthe same simulation times as in Figure 3. By inspectingthe left column, we observe the gradual suppression ofthe initial anisotropy with time (evolution from top tobottom). The anisotropy is globally maintained for thecase with damping, although at moderate p w the distri-bution becomes relatively flat on each side of µ w = 0. -1 -5 -4 -5 -4 Simulation: forw rightSimulation: forw leftAnalytic: forw rightAnalytic: forw left Figure 5. Linear growth rate as function of the wavenumber k for the case without damping (top panel) and with moder-ate damping (bottom panel). Blue and red lines correspondto the measured growth rate of right and left handed modes,respectively. Dashed yellow and magenta lines correspond tothe analytical expectation from Equation 24. In the bottompanel the analytical expectations for zero damping are alsoplotted (dot-dashed) for comparison. The effect of particle accumulation at µ w = cos θ w = 0is also evident in the two bottom panels on the right.4.2. Linear phase No damping vs moderate damping Paper I demonstrated that our MHD-PIC numericalapproach accurately reproduces the linear growth rateof the CRSI in the absence of explicit wave damping.Here, we extend the previous investigation by includ-ing ion-neutral damping. In Figure 5, we show the lin-ear growth rate Γ tot as function of the wavenumber k for two cases: no damping (top panel) and moderatedamping (bottom panel). There is very good agreementwith theory, despite some stochastic noise. In partic-ular, the most unstable mode is very well captured at k max = (cid:112) − /κm i Ω /p (cid:39) . R − L, .Comparing the top and bottom panels of Figure 5,two immediate effects due to the IN damping can beidentified. The first is the decrease of the maximumgrowth rate near kR L, = 1. It is equal to 2 . × − Ω without damping, and equal to (cid:39) . × − Ω when osmic ray streaming with ion-neutral damping -2 -1 -5 -4 -3 Figure 6. Dependence of the maximum linear growth rateof the instability as function of the ion-neutral momentumexchange rate ν in for n CR /n i = 10 − and V D /V A,i = 10. ν in = 1 . × − Ω = Γ max , / 2. The second effect isthe suppression of wave growth in the low-end and high-end parts of the spectrum, where Γ ( k ) − ν in / < 0. Inother words, the bandwidth of linearly growing waves isreduced with increasing ν in .4.2.2. Dependence on ν in In Figure 6 we report the dependence of the maximumgrowth rate of the instability on ν in . Simulation resultsare generally in excellent agreement with analytical pre-dictions, namely, the linear growth rate is only weaklyreduced for ν in / Γ max , < . 2, and rapidly decreases for ν in / Γ max , > . 5. As expected from Equation 24, therecan be no linear growth for ν in / Γ max , > 2, as even thefastest growing mode is damped.4.3. From linear to saturated phase The instability does not transition directly from lin-ear growth to a fully saturated state. At some pointin time the fastest growing modes at kR L, (cid:39) kR L, (cid:39) ν in = 0 and ν in ≤ . 1. The initial exponentialgrowth slows down roughly at t Ω (cid:39) × but ad-ditional wave growth continues until t Ω ∼ . Thisphase also corresponds to the fastest rate of decrease inthe streaming velocity of CRs, as seen in the lower panelof Figure 2.More detailed insight into the post-linear phase can begained by considering the wave spectrum. In Figure 7 wepresent the evolution of the wave spectrum for the case Figure 7. Spectra of forward right-handed modes at dif-ferent times for the case with no damping (top panel) andmoderate damping (bottom panel). Evolution for forwardleft-handed modes is identical. In the lower panel, the greyregion marks where initial linear growth is suppressed by INdamping. with no damping (top panel) and with ν in / Γ max , . kR L, (cid:39) 1, where the growth rate isthe fastest. During the post-linear phase (blue to redline), while the growth at kR L, (cid:39) kR L, > max , while the modes with kR L, ≤ k ) while the global level ofwave intensity is roughly unchanged. This effect can beclearly seen by comparing red, yellow and purple lines inthe top panel: the low- k cut-off shifts from kR L, (cid:39) . kR L, (cid:39) . 02 at theend of the simulation. The growth of high- k modes (i.e., kR L, > 5) is observed only during the post linear phase,while the increase of spectral energy in kR L, (cid:28) I. Plotnikov, E. Ostriker & X.-N. Bai For the case with moderate damping ( ν in / Γ max , . kR L, ≥ kR L, > 8, marked in grey in the figure). Thelate-time evolution of the damped case is significantlydifferent from the undamped case. There is a noticeableoverall decrease in the wave intensity at all k (differ-ence between red, yellow and purple lines). The peakof the spectrum shifts slightly towards smaller k at theend of the post-linear phase. Also, a spectral bump ap-pears at t Ω = 10 around kR L, = 0 . 2, which is dueto driving from the part of the CR distribution functionthat remains strongly anisotropic: f ( p > p ). By theend of the simulation the spectrum stabilizes. It is flatand narrow: 0 . < kR L, < 10 and kI ( k ) ∼ × − .Simulations with different values of ν in follow similartime-evolution.4.3.1. Growth of high- k (small wavelength) modes: HiRessimulations An important aspect of the post-linear phase is therapid growth of modes with kR L, ≥ even when thesemodes are not unstable to CRSI because of ion-neutraldamping . To better study these short-wavelength modeswe have conducted additional simulations with highernumerical resolution: ∆ x = 2 . d i in the HiRes runsinstead of ∆ x = 10 d i in Fid simulations. To understandthe mechanism driving this growth, we also performeda numerical experiment where the CRs were removedfrom the system after a given simulation time, at theend of the linear phase.In Figure 8 we show the wave spectrum of forwardright-handed modes in the HiRes simulation with nodamping. Three lines are plotted: (i) the spectrum atthe end of the linear phase (blue line), (ii) the spectrumat the end of the post-linear phase if CRs are switchedoff at the end of the linear phase (red line), and (iii) thespectrum at the end of the linear phase in the standardcase with no CRs switch off (yellow line). In the inset,we see that wave energy stays constant after switchingoff the CRs, while continue to grow otherwise. The dif-ference between blue and red lines illustrates the growthof high- k modes during the post-linear phase without driving by CRs. This difference is only seen in the re-gion kR L, ≥ 5. For comparison, the difference betweenthe blue and yellow lines shows the post linear evolu-tion when the driving by CRs is maintained. Here, theoverall increase of the wave energy is observed, togetherwith the shift of the spectrum to lower k .The same effect is present in simulations with ion-neutral damping. Figure 9 is the same as Figure 8 butwith ν in / Γ max , = 0 . 5. The additional dot-dashed blacklines analytically compensate for the effect of ion-neutral -2 -1 -14 -12 -10 -8 -6 -4 -2 t -8 -5 B / B Figure 8. Wave spectrum (forward right-handed mode) forthe HiRes simulation with no damping at different simulationtimes. Blue and yellow lines show the standard CR-drivensystem at t Ω = 2 × and t Ω = 8 . × , respectively.The red line shows the case of a model where CRs are turnedoff at the end of the linear phase ( t Ω = 2 × ) and thenfreely evolves until t Ω = 8 . × . The difference betweenblue and red lines illustrates the growth of high- k modesduring the post-linear phase without driving by CRs. Theinset in the bottom shows the time evolution of the magneticwave energy: the blue line follows the standard evolutionuntil t Ω = 2 × , the yellow line continues the evolutionwith CRs after that time, and the red line follows the casewhere CRs are switched off at t Ω = 2 × . The green dotmarks the time of the CRs switch-off. damping on the wave spectrum. During the post-linearphase, a decrease of wave energy at any k is imposed bythe ion-neutral damping (difference between blue andred lines). If one compensates the effect of damping, weobserve the same excess of high- k modes at the end ofthe post-linear phase (dot-dashed black line) as for theundamped case.Because these high- k modes can grow even in the ab-sence of CRSI, we conclude that the driving mechanismis purely an MHD effect (see further discussion in Sec-tion 5). 4.4. Saturated phase At the end of the post-linear phase exponential growthof the instability at all wavelengths is completed. Bythis time the CRs also fully experience the back-reactionfrom interacting with the waves they generated. In thissection, we characterize the particle and wave propertiesin based on analysis of our MHD-PIC simulations at latestages. osmic ray streaming with ion-neutral damping Figure 9. Same as Figure 8 but with ν in / Γ max , = 0 . 5. Theadditional dot-dashed black lines compensate analytically forthe effect of ion-neutral damping on the wave spectrum. Thegray-shaded areas delimit the regions where wave growth isnot allowed in linear theory, according to Equation 24. Late-stage wave amplitudes and particle distributions In Figure 10 we present the maximum wave intensity(normalized to B ) reached in the simulations as func-tion of ν in (plotted using blue circles and red squaresfor Fid and HiRes simulations, respectively). For theundamped or weakly damped cases we expect the mo-mentum flux associated with the original anisotropy inthe CRs to be transferred to forward-propagating Alfv´enwaves (see Kulsrud 2005). As found in previous simula-tions (e.g., Bai et al. 2019; Holcomb & Spitkovsky 2019),this leads to magnetic wave energy at saturation givenby δB /B ≈ . n CR /n i ) ( V D /V A,i − ν in > max , delimits the regionwhere the linear instability becomes impossible, accord-ing to Equation 24. The simulations show a gradualdecrease in ( δB max /B ) with increasing ν in , as mightbe expected. The decrease becomes abrupt when ap-proaching ν in = 2Γ max , , resembling an exponential cut-off.In order to get some insight into the dependence ofsaturated wave intensity on ν in , we cast our knowledgeof dominant dynamical processes into a simple model.This model consists of a system of two coupled ODEs.Let A = ( δB/B ) be the wave amplitude squared, andΓ be the wave growth rate (without damping). In our We verified that there is no wave growth for ν in > max , witha dedicated simulation, not presented here. Figure 10. Dependence of the maximum value of the mag-netic wave energy, δB /B , as function of ν in . Blue circlesare from Fid simulations and red squares from HiRes sim-ulations. The solid green curve corresponds to the modelsolving two coupled ODEs, described in the text. The hor-izontal black line delimits the maximal allowed value at ν in = 0, based on the measured saturation amplitude. For ν in / Γ max , > simple model, evolution of A in time is determined byd(ln A )d t = 2 ∗ (Γ − ν in , (29)d(ln Γ)d t = − ν s (cid:39) −A Ω π g , (30)where g is some factor of order unity. The first equa-tion describes linear wave growth partially limited byIN damping. The second equation serves as a proxyfor the dynamical adjustment of the growth rate tothe changing CR distribution, which is becoming moreisotropic under the effect of QLD in the bath of (grow-ing) waves. For initial condition, we set A (0) = 10 − ,Γ(0) = Γ max , . We further parameterize IN dampingrate as ν in = α Γ max , . In our simulations, we haveΓ max , = 2 . × − Ω , and α ranging from 0 to 2. Weintegrate these equations until Γ / Γ max , < α/ A reaches a maximum. In the case without damping,the wave amplitude should reach the value measured insimulations; from this constraint we find g ≈ . A as function of ν in are plottedwith a green solid line in Figure 10. The model quali-tatively agrees with the simulation results, but in theregion 0 . < ν in / Γ max , < I. Plotnikov, E. Ostriker & X.-N. Bai points remain well below the model prediction. Whilethis discrepancy reflects the simplifications behind ourtoy model, another important contributing factor in thesimulations is the µ = 0 barrier. When the isotropiza-tion process becomes stuck at the barrier, further wavegrowth is suppressed, even though the free energy fromthe CR anisotropy has not been fully utilized. This ef-fect is not easily accounted for in our toy model.Figure 11 presents the final state (at t Ω = 10 ) ofthree representative simulations with different valuesof ν in : columns from left to right show models using ν in / Γ max , = 0 . 1, 0 . . 5, respectively. The gray-shaded regions in the top panels mark the regions wherewaves cannot grow during the linear phase, according toEquation 24.Some interesting features of the saturated phase areapparent in Figure 11. Firstly, there is a reduction inthe global wave intensity with increasing IN dampingrate (compare from left to right top panels). Secondly,the saturated wave spectrum width exceeds the rangeimposed by the linear growth, with significant wave am-plitude in gray-shaded regions (best evidenced in topright panel). Also, the forward-propagating waves (blueand red lines in top panels) are at a similar level to thebackward-propagating waves (yellow and purple linesin top panels), a marked reduction from peak valuesreached during the post-linear phase brings all modesto comparable intensity level. High-momentum parti-cles (regions with p > p in middle and bottom pan-els) clearly isotropize less efficiently with increasing ν in :from left to right in bottom panels, there is an increas-ing area that is unaffected by the instability. This cor-responds to the lack of waves in kR L, (cid:28) ν in approachesthe critical value 2Γ max , . Finally, we observe an accu-mulation of particles near µ = 0 in all cases with INdamping ( “hot spots” in the region close to cos θ = 0in the middle panels and close to cos θ w = 0 in the bot-tom panels). The issue of crossing µ = cos θ = 0 barrierbecomes crucial when non-negligible damping is present.Not only is the global level of waves reduced, but also thewhole spectrum becomes narrower when ν in approaches2Γ max , . Both effects contribute in making it difficult toscatter across µ = 0.4.5. Dependence on spatial resolution By testing different grid resolutions we confirm thatthere is no noticeable effect on the linear phase of theinstability: no difference in the growth rate, wave en-ergy, streaming velocity, nor spectrum during this phase.However, the resolution can have an important impacton the post-linear phase. Specifically, resolving a largerdynamic range of wave modes with kR L, (cid:29) k part ofthe power spectra at t Ω = 10 in Figure 7 for the fidu-cial resolution to those at similar times ( t Ω = 8 . × , with CRs kept on) in Figure 8 and Figure 9 at high res-olution.In Figure 12 we compare the Fid and HiRes simula-tions for the time-evolution of the magnetic wave energy(top panel) and of the CR streaming velocity (bottompanel). Two representative cases are shown: no damp-ing (blue lines) and moderate damping ν in / Γ max , . t Ω ∼ − × .From the solid (Fid) and dashed (HiRes) blue lines, in-creased resolution has only a minor effect on the casewith no damping. More noticeable differences appearfor the case with moderate damping. During the post-linear phase, the magnetic energy decreases more slowlyin HiRes simulation and the streaming velocity decreasesat a faster rate. This effect is due to the higher level ofwave intensity stored in high- k modes during the post-linear phase.As found for the post-linear stage, numerical resolu-tion has also an effect on the saturated state of thesimulations. We expect, and indeed find, that higherresolution simulations achieve a higher level of small-scale waves. In Figure 13 we present the comparisonof the wave spectrum at t Ω = 3 × between Fidand HiRes simulations for two values of ν in = 0 . . ν in there are somedifference around kR L, = 1. But, more importantly,the HiRes simulations exhibit higher level of wave in-tensity at kR L, (cid:29) 1, without regard to the limit oflinear growth (which is equal to 0 in gray-shaded areas).The differences in the wave spectrum between Fid andHiRes simulations also lead to differences in the finalstate of the CR distribution function. In particular,higher amplitudes at large k increase the scattering rateat the resonances close to µ = 0, considering the reso-nance condition k = (1 /µ ) m p Ω /p . Figure 14 presentsthe distribution function at t Ω = 3 × of the simula-tions for Fid runs (upper panels) and HiRes runs (lowerpanels). We show cases with both weak (left) and mod-erate (right) damping, as in Figure 13. For both, HiRessimulations show a better level of isotropization of δf w at p/p ∈ [0 . , k modes in HiRes simulations. No noticeable dif-ference is observed at p > p , corresponding to quitesimilar spectra at kR L, < . Streaming velocity We now consider evolution of the CR bulk or stream-ing velocity, a key parameter since this characterizes thenet CR flux. As mentioned in Section 2.4, in the astro-physical literature the most commonly-adopted assump-tion is that the wave amplitude and streaming velocityare consistent with a state in which the linear growthrate balances the damping rate (Kulsrud & Cesarsky1971). From Equations 8 and 24, this would imply thatif the CRs remain isotropic in a frame moving at V st osmic ray streaming with ion-neutral damping Figure 11. Final spectrum (top panels), distribution function in the simulation frame δf ( p, cos θ ) /f (middle panels), anddistribution function in the wave frame δf w ( p w , cos θ w ) /f (bottom panels) in three simulations with different ν in . Left, middle,and right columns have ν in / Γ max , = 0 . 1, 0 . . 5, respectively. The gray shaded regions in the upper panels delimit theregions where waves cannot grow during the linear phase, according to Equation 24. Thick dashed lines in middle and bottompanels are the contours in momentum space resonant with the limiting k (i.e. the inner border of gray-shaded areas in the upperpanels). The gray dot-dashed lines illustrate resonant contours for a few different wavenumbers. relative to the gas, this velocity will decline until Equa-tion 27 is satisfied.Conceptually, the asymptotic streaming velocity ob-tained through balancing damping and growth relates tothe situation where the instability is driven, i.e., througha CR pressure gradient. In contrast, our numerical setupwith periodic boundary conditions corresponds moreto a ‘transient’ situation in which instability dies out.This affects the astrophysical relevance of the streamingspeed derived in our simulations. Nevertheless, it is in-teresting to follow the measured time-evolution of CRstreaming from the simulations.In the bottom panel of Figure 15 we show the evolu-tion of measured CR streaming velocity v rel (defined inEquation 28), including its dependence on numerical res-olution, for the weak damping case ν in / Γ max , = 0 . 1. Inthe Fid simulation (blue line) the final streaming veloc-ity of CRs is 3 . V A,i . The HiRes simulation isotropizesmore efficiently, showing a decrease of the streaming ve-locity to v rel = 1 . V A,i (but still above full isotropiza- tion v rel = V A,i ). We attribute this difference to bettercapturing of high- k modes in HiRes simulations.If we were to apply Equation 27 to the case with ν in = 0 . max , , the predicted asymptotic streaming ve-locity would be V st /V A , i = 1 . 45. This level is marked bythe dashed horizontal 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 level of waves is present the dis-tribution will continue to isotropize and v rel will con-tinue to decline. In constrast, for the HiRes model with ν in = 0 . max , (see Figure 12, red lines), the late-time v rel remains well above the value 3 . V A,i that would bepredicted by Equation 27, presumably due to the lowerwave amplitudes that reduce the scattering rate (com-pare upper and lower panels of Figure 13, and note alsothe buildup of particles to the right of µ = 0 in the lower-right of Figure 14). Taken together, these results makeclear that the simple approach of setting Γ max , = ν in / v rel instead6 I. Plotnikov, E. Ostriker & X.-N. Bai -8 -7 -6 -5 -4 -3 -2 Figure 12. Effect of numerical grid resolution on the time-evolution of the magnetic wave intensity (top panel) and ofthe CR streaming velocity (bottom panel). Blue lines showthe case with no damping, and red lines show the case with ν in / Γ max , = 0 . 5. Solid lines correspond to fiducial resolutionand dashed lines correspond to high-resolution simulations(see Table 2). depends on the evolution of ( δB/B ) . This is illus-trated for different values of ν in in Figure 12, with morecomplete isotropization in cases of lower ν in . In the realastrophysical case, an additional factor affecting evolu-tion would be the history of (anisotropic) energy sources.Finally, we recall (as pointed out in Section 2.4) thatbalancing the linear growth time and damping time ofwaves does not quantitatively take into account the ad-ditional time required for QLD to isotropize the dis-tribution, suggesting that the criterion Γ max , > ν in / Figure 13. Dependence on spatial resolution of the satu-rated wave spectrum (forward right-handed modes) for mod-els with ν in / Γ max , = 0 . ν in / Γ max , = 0 . DISCUSSIONPerhaps one of the most intriguing results of thepresent work is the evidence of high- k growth duringthe post-linear phase, presented in section 4.3. In Fig-ures 8 and 9 we have shown that these modes do notrequire the driving by CRs to be amplified. Presently,we do not have a satisfactory explanation of this effect,but hypothesize that it is due to some form of mode cou-pling or wave steepening into rotational discontinuities(e.g., Cohen & Kulsrud 1974). One could argue thatthis effect is due to our choice of parameters that leadto relatively high wave amplitude at saturation, up to δB/B (cid:39) . 05. This is likely considerably larger thanis present in the general ISM. However, we also per-formed simulations with n CR /n i as small as 5 × − and V D /V A,i = 2, observing the same effect of fast high- k mode growth during the post-linear phase.We also found that large-wavelength modes with kR L, (cid:28) k ) ∝ k κ − ,and that the IN damping removes the low- k end of the osmic ray streaming with ion-neutral damping Figure 14. Dependence on spatial resolution of the distri-bution function δf w ( p w , cos θ w ) /f , at t Ω = 3 × . Toppanel corresponds to Fid simulations and bottom panel cor-responds to HiRes simulations. Left and right show weakand moderate damping cases. wave spectrum with k < k min , where Γ( k < k min ) <ν in / 2. This can be seen in Figures 7, 9, 11. This implies,assuming only resonant wave-particle interaction, thathigh-momentum CRs with p > p max are not isotropizedand continue to freely stream. Here, p max is deduced byusing Equation 24 with Equation 7 using the resonancecondition kR L, = p / ( pµ ):( pµ ) max p (cid:39) (cid:20) . n CR n i Ω i ν in (cid:18) V D V A,i − (cid:19)(cid:21) , (31)where we fixed κ = 1 . pµ ) max are given in Ta-ble 3, adopting V D /V A,i = 2 and V D = 0 . c as twoextreme cases. Assuming p = m p c for CRs with energy E (cid:39) GeV , and we draw some general conclusions: • DMG : strongly suppressed CRSI; no CRisotropization at energies larger than GeV.Streaming could become extremely large. • MG : marginal instability and isotropy at GeV butnot beyond 100 GeV. The surface layers of molec-ular clouds may therefore be subject to CRSI ifthe initial anisotropy of CRs is not too small. • CNM : the conditions for triggering the CRSI areonly satisfied for highly super-Alfv´enic stream- -8 -7 -6 -5 -4 -3 -2 Figure 15. Dependence on spatial resolution of the time-evolving wave energy (panel a) and streaming velocity (panelb), for ν in / Γ max , = 1 / 10. Blue and red line correspondto Fid and HiRes simulation, respectively. The horizon-tal dashed line delimits the final streaming velocity as ex-pected by balancing the fastest growth rate with dampingrate: Γ max , = ν in / ing. The CR distribution could become quiteanisotropic in this phase. • WNM : most favourable environment for CRSI.Even with tiny anisotropy ( V D /V A,i of order unity)the instability is possible for GeV CRs. CRs with E (cid:29) GeV could also be isotropized if the initialdrift velocity is much larger than V A,i .In our own Milky Way and similar galaxies, the WNMis the dominant component of the ISM by mass; theabove estimate affirms that CRSI is astrophysicallyquite important. We note, however, that in the WNMother wave damping mechanisms can compete with INdamping that could prevent isotropization of CRs withenergies ≥ TeV (see, e.g., Fig.1 and corresponding textin Brahimi et al. 2020), but our values of E max are ingood agreement with similar calculation by Xu et al.(2016) (their Sect. 6.4 and Fig.15).The estimates above are in broad agreement withthose of Kulsrud & Cesarsky (1971) – being based onthe same argument – but are updated here with rep-8 I. Plotnikov, E. Ostriker & X.-N. Bai Table 3. Maximum momentum of cosmic-rays able to self-confine for different phases of the ISM.Phase ( pµ ) max p , V D = 2 V A,i ( pµ ) max p , V D = 0 . c WNM 12.2 345CNM 0.34 22.1MG 1.5 41.7DMG (cid:28) < resentative parameters for different neutral-dominatedmedia.The high-frequency limit was adopted in the presentstudy, i.e., ω A = kV A,i (cid:29) max[ ν in , ν ni ]. In this regime,the neutral fluid and ionized fluid are decoupled, and theion-neutral wave damping rate does not depend on thewavelength, Γ d = ν in / 2. We expect this approximationto hold for any kR L, ≥ . − 1, in general, and to be sat-isfied for CRs at GeV energies and below, which makeup most of the CR energy density and are responsiblefor most of the ionization. This justifies the one-fluidapproach, adopted in the present study. However, inreality very high energy CRs with low resonant frequen-cies are present in the ISM as well. If the whole CR en-ergy spectrum were fully represented, the low- k part ofthe wave spectrum (modes with kR L, (cid:28) max [ ν in , ν ni ])would be damped at slower rate than ν in / 2. This couldhave an interesting effect on the amplification and sur-vival of long-wavelength modes. However, we defer in-vestigation of k -dependent damping to future study.Similar to Paper I, many of our simulations show par-ticle accumulation at µ = 0. Naively, this would beexpected if waves were only present where CRSI is un-damped (at intermediate k near k max = (cid:112) − /κR − L, ),since there would be no large- k waves that are able toscatter small- µ particles subject to the resonance condi-tion kµ = Ω m/p . However, in practice nonlinear MHDeffects populate the large- k region regardless of whetherCRSI is damped or active, and therefore the µ = 0 cross-ing is possible. Nevertheless, in simulations the high- k regime is subject to numerical dissipation, and we findthat higher than standard resolution is required to limitparticle buildup. At the higher values of ν in , particlesstill build up near µ = 0 even with higher resolution (seeFigure 14).As we have previously emphasized, the problem of re-lating the streaming rate and wave amplitude to macro-scopic ambient properties clearly merits further study.The traditional “detailed balance” approach of equatingthe CRSI growth rate to the damping rate has recentlybeen adopted in galaxy formation simulations and otherstudies as a procedure for setting the scattering ratecoefficient, leading to a diffusion coefficient that variesproportional to ν in (see, e.g Hopkins et al. 2021b). To as- sess and quantitatively improve this kind of presciption,however, studies similar to the present one but allow-ing for macroscopic driving (via an imposed CR flux orenergy gradient) will be needed. CONCLUSIONMotivated by the recent progress in coupled fluid-kinetic (MHD-PIC) numerical techniques (Paper I), inthis work we studied the influence of ambipolar diffusion(ion-neutral damping of Alfv´en waves) on the CRSI. Weadopted parameters for a reference model that in theabsence ion-neutral damping lead to full isotropizationof the CR distribution function in the wave frame af-ter the instability saturates. This state corresponds toan asymptoptic state in which the CRs stream at Alfv´enspeed relative to the background fluid. We then exploreddifferent values of ion-neutral damping rate to study theinfluence on the outcome of the CRSI.Our main conclusions are summarized below: • The predicted exponential growth rate of the in-stability including IN damping, Γ( k ), is well-reproduced by the MHD-PIC technique. Thus thelinear theory is in good agreement with analyticalexpectations. • Evolution in a post-linear phase is crucial forisotropization of low- and moderate-energy CRs.During this phase, growth of high- k (wavelength λ (cid:28) R L, ) modes occurs. This growth is notdriven by CR anisotropy (the waves are outsideof the unstable range) but rather appears to be aresult of an MHD wave cascade or wave steepeninginto rotational discontinuities. The high- k wavesare important to scattering when µ is close to 0.The development of the high- k spectrum is bestcaptured in high-resolution simulations. • Systematic comparison between the reference case(no damping) and simulations with IN dampingreveals that the width of the wave spectrum de-creases with increasing ν in . The absence of low- k (compared to R − L, ) waves prevents isotropizationof high-energy CRs. There is no wave growth at allif ν in > max , , consistent with analytic theory. • With stronger IN damping, the maximum am-plitude of waves systematically decreases. When ν in = 0 . max , , the peak wave amplitude is anorder of magnitude below the no-damping case.Lower wave amplitudes reduce the rate of QLDand hence slows isotropization.This work is the second study in a series exploringthe physics of CRSI by means of MHD-PIC simula-tions. This numerical technique has proved to be valu-able because it captures the relevant microphysical CR-gyroscales and allows long-term simulations at low CR osmic ray streaming with ion-neutral damping Amano, T. 2018, Journal of Computational Physics, 366,366Amato, E., & Blasi, P. 2009, MNRAS, 392, 1591—. 2018, Advances in Space Research, 62, 2731Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015,ApJ, 809, 55Bai, X.-N., Ostriker, E. C., Plotnikov, I., & Stone, J. M.2019, ApJ, 876, 60 [Paper I]