Cosmic Ray Transport in the Ionized and Neutral ISM: MHD-PIC Simulations and Effective Fluid Treatments
DDraft version February 25, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Cosmic Ray Transport in the Ionized and Neutral ISM: MHD-PIC Simulations andEffective Fluid Treatments
Christopher J. Bambic, Xue-Ning Bai,
2, 3 and Eve C. Ostriker Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA; [email protected] Institute for Advanced Study, Tsinghua University, Beijing 100084, People’s Republic of China; [email protected] Department of Astronomy, Tsinghua University, Beijing 100084, People’s Republic of China Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA; [email protected] (Received XXXX; Revised XXXX; Accepted XXXX)
Submitted to ApJABSTRACTCosmic rays (CRs) have critical impacts in the multiphase interstellar medium (ISM), driving dy-namical motions in low-density plasma and modifying the ionization state, temperature, and chemicalcomposition of higher-density atomic and molecular gas. We present a study of CR propagation be-tween the ionized ISM and a neutral cloud. Using one-dimensional magnetohydrodynamic particle-in-cell simulations which include ion-neutral drag to damp Alfv´en waves in the cloud, we self-consistentlyevolve the kinetic physics of CRs and fluid dynamics of the multiphase gas. By introducing the cloud inour periodic domain, our simulations break translational symmetry and allow the emergence of spatialstructure in the CR distribution function. A negative spatial gradient forms across the fully-ionizedISM region while a positive gradient forms across the neutral cloud. We connect our results with CRhydrodynamics formulations by computing the wave-particle scattering rates as predicted by quasilin-ear, fluid, and Fokker-Planck theory. For momenta where the mean free path is short relative to thebox size, we find excellent agreement among all scattering rates. By exploring different cloud sizes andion-neutral collision rates, we show that our results are robust. Our work provides a first-principlesverification of CR hydrodynamics when particles stream down their pressure gradient, and opens apathway toward comprehensive calibrations of transport coefficients from self-generated Alfv´en wavescattering with CRs.
Keywords: cosmic rays – instabilities – magnetohydrodynamics (MHD) – ISM: clouds INTRODUCTIONThe gas in the interstellar medium (ISM) spans a hugerange of temperatures and densities, with the coldestatomic and molecular phases taking the form of cloudsembedded within a diffuse medium (Dickey & Lockman1990; Wolfire et al. 1995; Draine 2011). Pervading boththe diffuse ISM and dense clouds are relativistic cos-mic rays (CRs), with a (broken) power-law distributionextending over more than ten orders of magnitude in en-ergy (Grenier et al. 2015). With the total energy densityin CRs comparable to the thermal, kinetic, and mag-netic energy densities of the thermal ISM plasma, CRsmay have important dynamical consequences, especiallyin extraplanar regions where they may help drive galac-tic winds (Zweibel 2017; Recchia 2021). While CRs ofenergy ∼ GeV dominate dynamical effects, the lower-energy portion of the distribution is the most importantsource of ionization in regions shielded from UV radi-ation by dust (Padovani et al. 2020). This ionizationdrives both chemistry and heating. Cold atomic and molecular clouds have low ionizationfraction ( x i (cid:46) − ; Draine 2011), and short-wavelengthAlfv´en waves propagating in these clouds are subjectto strong damping through ion-neutral collisions. Be-cause CR transport is governed by scattering off of thesewaves (see below), CR mean free paths become quitelong within clouds. Furthermore, since the envelopesof clouds represent transition regions to the more ion-ized warm, diffuse gas, the propagation of CRs throughclouds may depend strongly on the structure of thisboundary layer (Padovani et al. 2009; Everett & Zweibel2011; Ivlev et al. 2018; Silsbee & Ivlev 2019).Alfv´en waves in the ISM are subject to growth via thegyroresonant cosmic ray streaming instability (CRSI;Lerche 1967; Kulsrud & Pearce 1969). If the CR driftvelocity v D (the mean velocity of the CR distributionin the rest frame of the thermal gas) exceeds the Alfv´envelocity v A , Alfv´en waves grow rapidly, feeding off thefree energy from the momentum space anisotropy ofthe drifting distribution. Magnetic perturbations exert a r X i v : . [ a s t r o - ph . GA ] F e b C. J. Bambic, X. Bai, & E. C. Ostriker
Lorentz forces on particles, pitch angle scattering thosein gyroresonance with Alfv´en waves. The resulting dif-fusion in momentum space tends to isotropize the CRdistribution in the frame of the waves. The CRSI satu-rates once the momentum distribution function becomesisotropic in the wave frame, implying a drift velocity v D ∼ v A .A nearly isotropic CR distribution lends itself to afluid description. While charged particles genericallyobey the 6D Vlasov-Maxwell system of equations, self-generated Alfv´en waves limit the anisotropy of the CRdistribution function, ensuring that the first two mo-ments (energy density and flux) describe the evolution.In essence, the near isotropy of the drifting distributionprovides a closure for a hierarchy of moment equations,transforming the CR pressure tensor into a scalar andmaking kinetic particles behave as a continuous fluid.The fluid approach to CR transport, CR hydro-dynamics (McKenzie & Voelk 1982; Zweibel 2017),is particularly useful on astrophysically macroscopicscales. CRSI-generated waves grow on the CR gyroscale( ∼ − − − pc) while the overall distribution func-tion varies on much larger scales, ∼ pc, a sepa-ration of up to 12 orders of magnitude in length scale.CR hydrodynamics overcomes this scale separation bywrapping the microphysics of Alfv´en wave-particle in-teractions into momentum-dependent “diffusion coef-ficients” which are themselves functions of the wave-particle scattering rate (Shalchi 2009; Pfrommer et al.2017). Unfortunately, computing the diffusion coeffi-cients is nontrivial, and has remained a subject of intensedebate with relevance for the solar wind (Jokipii 1971),molecular clouds (Padovani et al. 2020), galaxy forma-tion/evolution (Buck et al. 2020; Hopkins et al. 2021a,b),galactic winds and outflows (Farber et al. 2018; Zweibel2017; Bustard et al. 2020), galactic halos (Kempski &Quataert 2020; Ji et al. 2020), and even the intraclustermedium (ICM; Guo & Oh 2008; Ruszkowski et al. 2017).Low-ionization regions further complicate transport.In weakly-ionized atomic and molecular clouds, colli-sions between ions and neutral particles impose a dragforce on the ions, damping the short-wavelength Alfv´enwaves that are gyroresonant with CRs at GeV and lowerenergies (Kulsrud & Pearce 1969). Without significantwave-particle scattering, CR transport is ballistic, withparticles freely streaming through clouds. In this situa-tion, the CR mean free path to wave-particle scatteringcan approach or exceed the scale of clouds (Abdo et al.2010). Thus, the cloud-ISM interface represents a tran-sition in the physics of CR transport.The non-trivial large-scale CR dynamics induced bythis interface has been the subject of a number of previ-ous works. Skilling & Strong (1976) proposed that scat-tering of CRs by their self-generated Alfv´en waves coulddivert their trajectories, “excluding” them from enter-ing molecular clouds with the large fluxes inferred fromthe diffusive medium. This conclusion is countered by Cesarsky & Volk (1978) who argued that diffusion dueto small-scale magnetic irregularities and self-generatedwaves is ineffective at excluding CRs from GMCs, exceptfor particles with energies (cid:46)
50 MeV.Everett & Zweibel (2011) revisited this system usingfluid theory with an imposed spatial diffusion coeffi-cient. Using the steady state CR hydrodynamics equa-tions, they found that variation of the diffusion coef-ficient across the interface could decrease the numberdensity within clouds by an up to an order of magni-tude. Ivlev et al. (2018) argued that such modulationcan be explained as a ratio of the CR fluxes due to dif-fusion and free streaming. Silsbee & Ivlev (2019) pro-vided a means of connecting these numerical and semi-analytic models to observations of the CR ionization ratein GMCs (Indriolo et al. 2007), finding that variationsin the transport parameters yielded significant changesto the ionization rate vs. column density relation. Re-cently, Fujita et al. (2020) argued that X-ray and gammaray observations could be used to constrain the trans-port regime of CRs in GMCs. In galaxy-scale fluid sim-ulations, Semenov et al. (2020) found that changing CRdiffusivity around GMCs acts to regulate star formationand the structure of galactic disks. Simple analytic andfluid models consistently point to observational conse-quences imposed by variations in CR transport physics.Fluid models with a constant spatial diffusion coeffi-cient are often used to examine the role of CRs on galaxyproperties (Booth et al. 2013; Salem & Bryan 2014; Pak-mor et al. 2016; Pfrommer et al. 2017; Wiener et al. 2017;Chan et al. 2019; Hopkins et al. 2021b). While thesemodels can yield useful constraints and comparisons toobservations (Hopkins et al. 2021a), a full treatment ofCRs around GMCs requires calculating the diffusion co-efficient from first principles.In this paper, we take the first steps toward a ki-netic study of CR propagation through the multi-phase ISM. Using magnetohydrodynamic particle-in-cell(MHD-PIC) simulations developed within the
Athena code framework (Bai et al. 2015), we study a drift-ing population of CRs within a toy model of a mag-netized ISM-cloud system: a one-dimensional box withdistinct cloud and ISM regions separated by a sharpboundary. Our simulations self-consistently evolve thegyroresonant CRSI and subsequent CR-wave interac-tions in both regions while including the effects of ion-neutral damping of Alfv´en waves in the cloud. Thus,we extend the simulations of Bai et al. (2019) (here-after BOPS19), which studied the gyroresonant CRSIin a uniform plasma without damping, to a model sys-tem relevant for the multiphase ISM. In separate work,we have also used MHD-PIC to investigate the effect ofion-neutral drag on development of the CRSI in linear,post-linear, and saturated stages, considering a uniformmedium but different levels of the drag coefficient (Plot-nikov et al 2021, in prep.).
HD-PIC Simulations of Molecular Clouds THEORETICAL PRELIMINARIESThe fluid-like behavior of cosmic ray (CR) particlesas described by CR hydrodynamics is a consequence ofwave-particle scattering from Alfv´en waves excited bythe gyroresonant CRSI. In this section, we discuss keyresults of the CRSI and how CR hydrodynamics is con-structed from quasilinear theory. The resulting momentequations will motivate the results in this work, namelythe emergence of a spatial gradient in the CR energy density as well as the computation of the wave-particlescattering rate.2.1.
Wave-Particle Scattering Rates
Throughout this paper, we work in spherical coordi-nates in momentum space, ( p, θ, φ ), where p is the CRmomentum, θ is the pitch angle, and φ is the gyrophase.We define the pitch angle of a particle relative to themagnetic field direction through its cosine, µ ≡ cos( θ ).With this definition, particles moving forward, parallelto the magnetic field correspond to µ = 1, particles mov-ing backward, parallel to the field have µ = −
1, andparticles gyrating about the magnetic field completelyperpendicular to the field direction have µ = 0. Thislatter case corresponds to a 90 ◦ pitch angle.A simple analysis of a particle propagating along abackground magnetic field B = B ˆ x perturbed by anAlfv´en wave with frequency ω shows that the particlewill be in resonance with the wave when the particlemomentum p = γm v satisfies ω − kvµ ± Ω = 0 . (1)Here, k represents the wavenumber along the back-ground field, and the ± symbol represents the differencein resonance condition for left (+) and right ( − ) handedwaves. The magnitude of the particle velocity v is de-noted by v , and Ω is the relativistic CR gyrofrequency,related to the non-relativistic cyclotron frequency Ω c byΩ = Ω c γ = qB γmc , (2)where q and m are the CR proton charge and ion massrespectively, and the Lorentz factor γ is given by, γ = 1 (cid:114) − (cid:16) | v | c (cid:17) = (cid:115) (cid:18) | p | mc (cid:19) . (3)Where appropriate, we differentiate between momentumor wave components parallel ( (cid:107) ) and perpendicular ( ⊥ )to the equilibrium field. For an Alfv´enic perturbation δ B ⊥ such that the amplitude satisfies | δ B ⊥ | (cid:28) | B | , aparticle in gyroresonance with the wave will experiencean average parallel Lorentz force, leading to a change inpitch angle, ∆ θ ≈ π (cid:18) | δ B ⊥ || B | (cid:19) cos φ, (4)over one cyclotron orbit (Kulsrud 2005). Within thequasilinear approximation, uncorrelated Alfv´en wavepackets each scatter a resonant particle by ∆ θ , leadingto random walk diffusion of the pitch angle at a rate, ν scat = N (cid:88) i =1 (cid:104) ∆ θ (cid:105) i t = N (∆ θ ) πN Ω − = π (cid:18) δB ⊥ B (cid:19) Ω , (5) C. J. Bambic, X. Bai, & E. C. Ostriker where we sum over the number of scattering events N .Working out the full quasilinear theory refines thisscattering rate. On slow timescales, a particle distri-bution f ( t, x, p ) in a background magnetic field underthe influence of a spectrum of weak electromagneticfluctuations undergoes quasilinear diffusion (Kennel &Engelmann 1966). In the frame comoving with Alfv´enwaves, the wave electric field vanishes since the mag-netic field is stationary in this frame. Because magneticfields can do no work on the particles, diffusion in mo-mentum space must occur through pitch angle scatter-ing (Jokipii 1966). Thus, particles obey the quasilineardiffusion equation, ∂f w ∂t + v w µ w ∂f w ∂x = ∂∂µ w (cid:20) − µ w ν QL ( µ w ) ∂f w ∂µ w (cid:21) , (6)which describes spatial advection along field lines as wellas pitch angle diffusion through momentum space. Here,we introduced the quasilinear scattering rate ν QL as wellas the subscript “ w ” for quantities measured in the waveframe. Note that this equation is only valid in the waveframe since any other frame would require correctionsto the RHS due to wave electric fields.The quasilinear scattering rate is given by ν QL ( p w , µ w ) = π Ω k res I ± ( k res ) , (7)where the Alfv´en wave power spectrum I ± ( k ) is normal-ized through, (cid:90) I ± ( k ) dk = (cid:28) δB ⊥ B (cid:29) x , (8)and the notation (cid:104) ... (cid:105) x indicates a spatial average.Above, we introduced the resonant wavenumber, k res = m Ω c p w µ w , (9)which corresponds to the wavenumber of an Alfv´en wavein resonance with a particle with wave-frame momen-tum p w and pitch angle µ w . At µ w = 0, the resonantwavenumber is formally infinite, and no wave can scat-ter the particle. This situation is referred to as the 90 ◦ pitch angle problem (Felice & Kulsrud 2001) and canbe alleviated by resonance broadening, mirror scatter-ing (Holcomb & Spitkovsky 2019), or nonlinear wave-particle interactions (BOPS19).Wave-particle scattering relaxes the CR distributionto an isotropic distribution “self-confined” to co-movewith Alfv´en waves. Thus, Equation 6 is a specific in-stance of a Fokker-Planck equation for the pitch angleevolution of the distribution (Shalchi 2009), ∂f w ∂t + v w µ w ∂f w ∂x = ∂∂µ w (cid:20) D µµ ∂f w ∂µ w (cid:21) , (10)where the Fokker-Planck (FP) diffusion coefficient isgiven by Mertsch (2020) as D µµ = (cid:10) (∆ µ ) (cid:11) t . (11) The angle brackets here represent an ensemble average(see Equations 50 and 51 for our numerical implemen-tation of this calculation).By tracking individual particles, we can compute thechange in pitch angle ∆ µ ( t ) for each particle and averageover the distribution function to compute D µµ , which isrelated to the quasilinear scattering rate through D µµ ( p w , µ w ) = 1 − µ w ν QL ( p w , µ w ) . (12)Thus, the two rates presented here (Equations 7 and 11)provide independent means of computing the wave-particle scattering rate. BOPS19 previously showed thatthe distribution function evolves consistent with Equa-tion 6 for particles away from µ = 0. Here, we studyindividual particle orbits and compare their ensembleaveraged scattering rates (via Equation 11) to those pre-dicted from quasilinear theory (via Equation 7). Sincethe Fokker-Planck rate captures nonlinear effects lost inthe quasilinear approximation, comparison of these scat-tering rates provides insight into the strength of nonlin-ear versus quasilinear scattering.2.2. Fluid Theory
Fluid theories reduce the dimensionality of statisticalsystems by averaging over moments of the underlyingmomentum space distribution. This procedure gener-ates a hierarchy of moment equations, the truncationof which requires a closure. In our case, the closure willcome from an assumption about the near-isotropy of theCR distribution.We define the particle energy as E ( p ) = γmc . TheCR energy density for a given momentum p w is then E CR ,w ≡ (cid:90) dµ w E ( p w ) f w ( p w ) , (13)and the corresponding parallel CR energy flux in thewave frame is F CR ,w ≡ (cid:90) dµ w E ( p w ) v (cid:107) ,w ( p w ) f w ( p w ) , (14)where the parallel velocity v (cid:107) ,w = v w µ w = p w µ w / ( mγ w ).Note that formally, these quantities are energy densityand flux per momentum density . In Figures 1-3, we plotthe total energy density and flux, i.e. the result of ap-plying an additional integration 2 π (cid:82) dp w p w .Multiplying the Fokker-Planck equation (10) by theenergy and energy flux per particle ( E ( p ) and E ( p ) v (cid:107) ( p )respectively) and integrating over µ w and gyrophase φ yields ∂∂t E CR ,w + ∂∂x F CR ,w = 0 , (15) ∂∂t F CR ,w + ∂∂x (cid:90) dµ w v w E µ w f w = − (cid:90) dµ w E v w D µµ ∂f w ∂µ w , (16) HD-PIC Simulations of Molecular Clouds f w is close to isotropic, we can assume that anyanisotropy in the distribution function does not con-tribute substantially to the integral over µ w for the sec-ond term on the left-hand side of Equation 16. We there-fore approximate the expression inside the gradient as (cid:90) dµ w v w µ w E f w ≈ v w E ( p w ) (cid:90) dµ w f w = v w E CR ,w . (17)This assumption of approximate isotropy of the distri-bution relies upon v D /c (cid:28)
1, which is certainly true forrelativistic CRs. Thus, CR hydrodynamics may still beapplied for a distribution with rapid streaming as longas v A (cid:28) v D (cid:28) c .We now introduce an effective, “fluid” scattering rate ν eff , such that the pitch angle dependence of the dis-tribution function and wave-particle scattering rate isentirely absorbed into ν eff , ∂∂t F CR ,w + v w ∂∂x E CR ,w = − ν eff F CR ,w . (18)Here, ν eff = 1 F CR , w (cid:90) dµ w E v w D µµ ∂f w ∂µ w , (19)and for quasilinear theory substitution of Equation 12results in ν eff , QL = 1 F CR , w (cid:90) dµ w E v w (cid:18) − µ w (cid:19) ν QL ∂f w ∂µ w . (20)For a flat spectrum with k res I ± ( k res ) = constant,Equation 20 reduces to ν eff , QL = ν QL . We treat the fluxmoment equation (18) as a function of momentum, com-puting the effective scattering rate for each momentumbin. This procedure allows us to compare the quasilin-ear, Fokker-Planck, and fluid scattering rates as a func-tion of particle momentum. METHODSWe wish to study the spatial and temporal evolution ofthe CR distribution function and compare this evolutionto predictions from fluid theory. This evolution arisesnaturally at the interface of two transport regimes: anambient ISM where particle transport is diffusive anda cloud where Alfv´en waves are damped through ion-neutral collisions and CRs stream freely.3.1.
The MHD-PIC Method
The magnetohydrodynamic particle-in-cell (MHD-PIC) method is a plasma model which evolves a kinetic species (CRs) under the influence of force fields calcu-lated using the equations of MHD (Bai et al. 2015). Forthe gyroresonant CRSI, the thermal plasma, referred toas the “gas” with the subscript “g” is described by theequations of ideal MHD with source terms: ∂ρ∂t + ∇ · ( ρ v g ) = 0 , (21) ∂ ( ρ v g ) ∂t + ∇ · ( ρ v g v g − BB + P Tot )= − (cid:18) qn CR E + J CR c × B (cid:19) − ν IN ρ v ⊥ g , (22) ∂ E Tot ∂t + ∇ · [( E Tot + P Tot ) v g − ( B · v g ) B ] = − J CR · E . (23)Here, ρ is the gas density, I is the unit tensor, the to-tal pressure P Tot = P g + B /
2, and the electric field E = − v g × B /c . The total energy is given by E Tot = P g γ ad − ρv g + 12 | B | , (24)where γ ad is the adiabatic index. We work in units wherethe magnetic permeability is unity such that 4 π = 1.Ion-neutral damping is implemented as a simple ex-ponential attenuation factor, v ⊥ g → v ⊥ g e − Γ IN ∆ t , (25)where ∆ t is the simulation time step and Γ IN is theion-neutral damping rate. The details of this damp-ing including numerical properties and the effect on thegrowth and saturation of the gyroresonant CRSI is dis-cussed in detail by Plotnikov et al. (2021 in prep.).Particles are pushed by the Lorentz force, d p j dt = (cid:16) qmc (cid:17) j ( c E + v j × B ) , (26)where “j” refers to the j th particle, p j = γ j v j , and thecharge to mass ratio q/mc ≡ δf -method (Dimits & Lee 1993; Parker & Lee 1993; Hu &Krommes 1994; Denton & Kotschenreuther 1995; Kunzet al. 2014b; Bai et al. 2019). In the δf -method, thedistribution function f ( t, x, p ) is split into a uniformstatic background f ( p ) and the particles evolved in thesimulation, represented by a perturbation on this back-ground, δf ( t, x, p ). The PIC method itself pushes par-ticles in phase space in order to compute a weightingfunction for the jth particle, w j : w j = δf ( t, x j ( t ) , p j ( t )) f ( t, x j ( t ) , p j ( t )) = 1 − f ( x j ( t ) , p j ( t )) f (0 , x j (0) , p j (0)) . (27) C. J. Bambic, X. Bai, & E. C. Ostriker
We can then straightforwardly compute the CR numberand current densities in Equations 22 and 23, n CR = n CR , + (cid:90) δf ( t, x, p ) d p (cid:39) n CR , + N p (cid:88) j =1 w j S ( x − x j ) , (28) J CR ( t, x ) = J CR , + q CR (cid:90) v δf ( t, x, p ) d p (cid:39) J CR , + q CR N p (cid:88) j =1 w j v j S ( x − x j ) , (29)where the summations are over the total number ofparticles N p , the subscript “0” represents moments ofthe background distribution f ( p ), and S ( x − x j ) is theshaping function for interpolating a point particle at x -coordinate x j to the grid. We use a triangular-shapedcloud (TSC; Birdsall & Langdon 1985) for S ( x − x j ).The simulations are run in a one-dimensional box withperiodic boundaries and an equilibrium magnetic field B = B ˆ x . Particles are limited by the numerical speedof light C = 300 v A , which reduces the separation be-tween wave velocities ( v A ) and particle velocities ( ∼ C ).Any analysis involving the speed of light c (computationof energy density or flux), uses this numerical speed oflight c = C . Similarly, Lorentz transformations use thisnumerical speed of light. Since the error in these trans-formations scales as O ( v A / C ), our choice of C ensures v A / C (cid:28) Athena
MHD code (Stone et al.2008), with constrained transport to ensure ∇ · B = 0(Evans & Hawley 1988). Time integration is donethrough the corner transport upwind method (Gardiner& Stone 2005, 2008) and we use the Roe Riemann solver(Roe 1981) with third-order reconstruction for spatialintegration. The details of the MHD-PIC method canbe found in Bai et al. (2015) and BOPS19.3.2. Modifications to Previous Simulations
BOPS19 studied a suite of magnetohydrodynamicparticle-in-cell (MHD-PIC) simulations of the gyroreso-nant CR streaming instability (CRSI) which varied theequilibrium number density ( n CR , ) and drift velocity( v D ) to study the growth of the instability and quasi-linear evolution of the particle distribution. Their sim-ulations were performed in a frame co-moving with theCR drift velocity such that the momentum space distri-bution of the particles was described by an isotropic κ distribution (Summers & Thorne 1991), f κ ( p ) = n CR , ( πκp ) / Γ( κ + 1)Γ( κ − ) (cid:34) κ (cid:18) pp (cid:19) (cid:35) − ( κ +1) , (30)where p = 300 v A is the peak momentum of the dis-tribution and κ = 1.25. Waves in the BOPS19 simula-tions were able to grow unimpeded by physical dampingmechanisms, succumbing only to numerical diffusion.The cloud problem investigated in this paper requiresthree modifications to the BOPS19 simulations:1. Rest Frame rather than Drift Frame: For the cloudto be stationary within the simulation (lab) frame,the simulation must be performed in the rest frameof the thermal plasma rather than in a frame drift-ing with the CRs. Because particles are relativis-tic, transformations of the κ distribution to thisnew frame require a Lorentz transformation.2. Large Box Size: Particles are allowed to traverse ∼
10 mean free paths to scattering with waveswithin the ambient ISM region (hereafter “ISM”).The cloud is made comparable in size to the ISMto capture large-scale evolution.3. Ion-neutral Damping of Waves: Alfv´en waves areattenuated in the cloud region while being allowedto grow unimpeded in the ISM.Related to (1), we have verified that we can recoverthe BOPS19 results, including wave growth rates, powerspectra, and distribution functions for simulations withno wave damping independent of frame. Related to(2), based on calculations of mean free paths in runsof the BOPS19 problem with varying box size (see Ap-pendix A), we choose the minimum ISM length, L ISM = 10 v A Ω − c , (31)which corresponds to ≈ n CR , / n i = 10 − . For higher CR numberdensities, more mean free paths are present in the ISM.Table 1 summarizes parameters for the 4 simulationspresented in this work. For all simulations in this pa-per, we use a constant initial drift velocity v D /v A = 10to increase the wave growth rate, saturation amplitude,and initial CR flux relative to the fiducial model ofBOPS19. The Fiducial simulation, where the ISMlength L ISM and the cloud length L Cloud are equal and n CR , /n i = 10 − , is our primary focus. NCR2 studies thecase where the waves grow more rapidly and the scatter-ing rate is twice that in the
Fiducial run. By extendingthe cloud region in the simulation
Long Cloud , the CRenergy flux decreases more slowly and becomes a sub-dominant term in the fluid equation. Finally, by extend-ing the ISM region relative to the cloud, the
Long ISM simulation studies the situation where the CR distri-bution is isotropized more rapidly. All simulations are
HD-PIC Simulations of Molecular Clouds Table 1.
Summary of SimulationsSimulation n CR , /n i v D /v A ν IN L ISM /L Cloud L (Ω c ) ( v A Ω − c )Fiducial 10 −
10 5.08 × − × NCR2 2 × −
10 10.16 × − × Long Cloud 10 −
10 5.08 × − × Long ISM 10 −
10 5.08 × − × computed in a box at least 20 times longer than thatstudied in BOPS19, with each grid cell spanning 10 v A Ω − c , where we work in units of v A = Ω c = 1. Forthe Long Cloud and
Long ISM simulations, the box is 50times longer. Thus, all simulations presented are com-puted with only 16 particles per cell per type (where 8different particle types are used which span the full rangeof momentum, − ≤ log ( p d ) ≤
2) until t = 10 v A Ω − c .3.2.1. Rest Frame vs. Drift Frame
BOPS19 did their calculations in the initial drift frameof the CRs, implying that their background thermalplasma initially had a velocity − v D ˆ x . We shall insteadwork in the frame where the initial background fluid(both ISM and cloud) are stationary, so that the CRdistribution is isotropic in a frame moving with veloc-ity v D ˆ x . We shall refer to the initial rest frame of thebackground plasma (the simulation frame) as either the“rest” or “lab” frame, and the frame where the initialCR distribution is isotropic as the “drift” frame.The 4-momentum of a particle p µ is defined as p µ = γ ( v ) ( c, v x , v y , v z ) , (32)where γ is the Lorentz factor (Equation 3). The mag-nitude of a 4-vector, p µ p µ , is invariant under a Lorentzboost along the x -direction, Λ µν : Λ µν [ v D ] = γ ( v D ) − v D c γ ( v D ) 0 0 − v D c γ ( v D ) γ ( v D ) 0 00 0 1 00 0 0 1 , (33) p µ drift = Λ µν [ v D ] p ν rest . (34)The inverse of Λ µν , needed to boost from drift frame tolab frame, is trivially computed with the substitution v D → − v D .In initializing the particles, we employ the κ distribu-tion, but because this applies in the drift frame we mustboost to the lab frame. In the lab frame, the CR dis-tribution is not initially isotropic; the anisotropy resultspurely from the drift velocity of the distribution. Thetransformation of the κ distribution f κ ( p drift ) is then f κ ( p drift ) = f κ ( | Λ µν [ v D ] p ν rest | ) = f ( p rest ) , (35) where the underlying distribution is unchanged asLorentz transformations preserve phase space volume.Here, the magnitude symbol refers to magnitude of the3-momentum, i.e. (cid:112) p + p + p . Because of conserva-tion of particle number, the Lorentz boost transformsthe isotropic κ distribution into a drifting prolate distri-bution unstable to the gyroresonant CRSI.In the isotropic frame, the distribution has no net ve-locity, so there is no net current; however, in the restframe, the net drift of the distribution creates significantcurrent in the x -direction. This current is computed byboosting the 4-current density J µ CR = ( q CR n CR , , J CR , )from the drift frame of the CRs to the gas rest frame.3.2.2. CRSI Growth Rate & Ion-Neutral Damping
The atomic and molecular ISM has ionization frac-tions ranging from ∼ − in the warm atomic gas to ∼ − in the cold atomic gas to (cid:46) − in the moleculargas (Draine 2011). Many neutral atoms and moleculesare present, serving as targets for the few ions carryingAlfv´en waves. These neutrals provide a collisional dragforce on the ions, which removes their momentum anddamps the waves.In this work, we parameterize the ion-neutral damp-ing rate based on the peak growth rate of the gyrores-onant CRSI in the absence of damping. The growthrate as a function of k can be worked out from thefull Vlasov-Maxwell system with a drifting populationof CRs (Zweibel 2003; Amato & Blasi 2009). Using thenotation introduced in BOPS19, the growth rate isΓ CR ( k ) = 12 n CR , n i Ω c (cid:18) v D v A − (cid:19) Q ( k ) , (36)where for a κ distribution, Q ( k ) = √ πκ / Γ( κ + 1)Γ( κ − ) 1 s [1 + 1 / ( κs )] κ , (37)and s = kp / ( m Ω c ). The growth rate is maximizedat s = (cid:112) − /κ . For κ = 1.25, p = 300 v A , and v D = 10 v A used throughout this work, the peak growthrate is Γ CR ( k peak ) ≈ n CR , /n i ) Ω c .In low-ionization regions, the damping rate of Alfv´enwaves depends in general on the wave frequency ω . For C. J. Bambic, X. Bai, & E. C. Ostriker trans-relativistic CRs and typical conditions in atomicand molecular gas, ω = kv A,i (cid:29) ν IN , where v A,i is theAlfv´en speed considering just the ions , and ν IN is theion-neutral collision rate (see Plotnikov et al. 2021, inprep. for detailed discussion of ISM conditions). In thislimit, the damping rate of waves is simplyΓ IN = ν IN . (38)We define the critical collision frequency as the ratewhere the peak wave growth equals damping, ν crit = 2Γ CR ( k peak ) = 5 . (cid:18) n CR , n i (cid:19) Ω c . (39)Throughout this work, we set the ion-neutral dampingrate (in the cloud) to be this critical rate. We performedan exploration of various damping rates, concluding thatany damping rate above the critical rate (even 10 timesthis rate) produces the same particle evolution as the ν IN = ν crit case. This is perhaps not surprising, sinceany damping rate above ν crit will completely suppresswave growth, resulting in ballistic propagation throughthe cloud. Plotnikov et al (2021, in prep.) considersimulations with a range of ν IN /ν crit = 0 . − EMERGENCE OF SPATIAL STRUCTUREWe are studying a time-dependent problem whereina distribution of cosmic rays (CRs) drives the growthof Alfv´en waves through the gyroresonant CRSI. Simul-taneously, wave-particle scattering slowly removes theanisotropy of the drifting CR distribution, decreasingthe drift velocity v D . Because CRs diffuse through theISM region and free stream through the cloud, spatialstructure inexorably emerges in the distribution func-tion. A negative spatial gradient in the CR energy den-sity spans the ISM while a positive gradient spans thecloud. The minimum in the CR energy density is lo-cated at x = 10 v A Ω − c in the Fiducial simulation,which we define as the “leading ISM-cloud interface.”In this section, we explore the origin of this energydensity gradient, which appears universally in our cloudsimulations. First, we introduce the diagnostics used tostudy the CRs and waves.The Alfv´en wave energy density E A is defined, E A = 12 ρ | δ v ⊥ | + | δ B ⊥ | | δ B ⊥ | , (40)where v ⊥ = v y ˆ y + v z ˆ z and B ⊥ = B y ˆ y + B z ˆ z are theperpendicular components of the velocity and magneticfields respectively. Using the δf -method formalism, we For convenience, in general we omit the “i” subscript on theAlfv´en speed, but here we include it to emphasize that the rele-vant Alfv´en speed is v A,i = B/ρ / i . can define the CR number and energy densities and theparallel energy flux at the i th position x i , n CR ( t, x i ) = (cid:90) d p [ f ( p ) + δf ( t, x i , p )] , (41) E CR ( t, x i ) = (cid:90) d p E ( p ) [ f ( p ) + δf ( t, x i , p )] , (42) F CR ( t, x i ) = (cid:90) d p v (cid:107) ( p ) E ( p ) [ f ( p ) + δf ( t, x i , p )] . (43)We can convert between energy densities and fluxes indifferent frames by constructing the energy flux 4-vectorin the rest frame of the CRs, F µ rest = ( E CR c, F CR , , , (44)and boosting into the frame of the waves, F µ wave = Λ µν [ v A ] F ν rest . (45)This procedure can be used to find the background fluxand energy density, E CR , = 3 . × n CR , ρ v A , (46) F CR , = 3 . × n CR , v D ρ v A , (47)where we work in units of the background gas density ρ and Alfv´en velocity v A .The temporal evolution of the Alfv´en wave energydensity and CR flux as measured in the wave frame isshown in Figures 1a and 1b respectively. While wavesgrow at the rate predicted by Equation 36, the satura-tion amplitudes both in the ISM and cloud can differfrom the Fiducial run by up to a factor of 4. Thisis not surprising since the
Long Cloud simulation has2.5 × the number of particles as the Fiducial run andthus 2.5 × the CR momentum. Particle scattering trans-fers this momentum to Alfv´en waves (Kulsrud 2005; Baiet al. 2019); however, strongly-suppressed waves in thecloud are unable to grow sufficiently to scatter and re-ceive momentum at a significant rate. Nearly all mo-mentum lost by the CRs must be transferred to wavesin the ISM. Similarly, the saturation amplitude in the Long ISM simulation is lower than that of the
Fiducial run since the ISM length is increased by a factor of 4while the total CR momentum to be deposited increasesby only a factor of 2.5.4.1.
Wave and Cosmic Ray Co-evolution
For all simulations, the average wave-frame energyflux decreases with time; however, no simulation reaches
HD-PIC Simulations of Molecular Clouds )10 W a v e E ne r g y D en s i t y ( A / v A ) FiducialNCR2Long CloudLong ISM CR ( k peak ) )0246810 CR E ne r g y F l u x i n W a v e F r a m e ( F w a v e C R / C R , v A ) FiducialNCR2Long CloudLong ISM
Figure 1.
Left : Alfv´en wave energy densities E A in ISM (solid) and cloud (dotted) regions respectively for all simulations. Inthe ISM region, all simulations show growth rates consistent with the theoretical expectation, 2Γ CR ( k peak ) (black dashed line).The saturation amplitudes vary substantially depending on the ISM to cloud length ratio and n CR , /n i . Right : Flux of CRs, asmeasured in the wave frame, in units E CR , v A . For all simulations, the flux declines in time as the CR distribution is isotropizedby particle-wave scattering. a fully isotropic state within the simulation time. Be-cause waves grow most rapidly in the NCR2 simula-tion and grow to an amplitude larger than all but the
Long Cloud run, scattering is frequent and the distri-bution is more rapidly evolved toward isotropy. The
Long Cloud simulation isotropizes slowly since at anygiven time, the majority of particles are located in thecloud region and experience no scattering. In Figure 2,we show the profiles of the Alfv´en wave amplitude to-gether with the CR number density, energy density, andflux, at time 5 × Ω − c . At this stage, Figure 1 showsthat the wave energy is saturated and the CR flux hassignificantly declined from its initial value. Figure 2 con-veys one of the main results of this paper: the presenceof a cloud region allows for the formation of a spatial gra-dient in the CR number and energy densities. The max-imum is on the upstream side of the ISM (near x = 0)while the minimum is where the streaming CRs enterthe cloud from the ISM (near x = 10 v A Ω − c for thefiducial model).The gradients in the profiles of E CR and n CR can beunderstood within the context of the CR fluid flux equa-tion, Equation 18. In the cloud, Alfv´en wave growth issuppressed and the scattering rate remains small. Thus,the RHS of Equation 18 can be ignored. Since the en- ergy flux is always decreasing in time (Figure 1b), a posi-tive spatial energy density gradient must form across thecloud to balance the temporal drop in flux.In the ISM, a negative gradient forms from particlediffusion as CRs traverse the region. Since the simula-tion is periodic, the CR energy density (pressure) gra-dients must balance one another on the upstream sideof the ISM (which is also the downstream side of thecloud) located at the domain boundary. The temporaldecrease in flux and negative spatial energy density gra-dient must work together to balance the negative scat-tering rate term on the right-hand side of Equation 18.As we shall show, the magnitude of the negative num-ber/energy density gradient is consistent with fluid the-ory for the measured CR flux based on the predictedeffective scattering rate.Since our analysis relies on constructing an ISM wherewave-particle scattering is sufficient for particles to prop-agate diffusively, we verify that multiple mean free pathsare present in the ISM upon saturation. The mean freepath to scattering λ mfp for a particle with momentumcoordinates ( p , µ ) can be estimated as λ mfp ( x ) = (cid:90) dµ | µ | vν scat ∼ π p Ω c (cid:18) δB ⊥ ( x ) B (cid:19) − , (48)0 C. J. Bambic, X. Bai, & E. C. Ostriker A / v A × + × n CR, w / n CR, 0 × + × w / CR, 0
25 50 75 100 125 150 175x ( v A )6.006.05 F CR, w / CR, 0 v A = 0.26 L n CR, 0 n i vv A ( BB ) × + × w / v A F CR, w / v A Time = 5.00e+04 ISM Cloud
Figure 2.
Spatial distribution of wave and particle energy and flux for
Fiducial simulation, shown at time t = 5 × v A Ω − c .All quantities are smoothed using convolution with a smoothing length of 0.01 L . The mean free path for the ISM is shown(light blue bar), and the cloud region is shaded. At this time, the energy density and number density of the CRs form a negativegradient across the ISM region and a positive gradient across the cloud. The gradient develops in the ISM due to spatial diffusionas CRs are scattered by waves (top row). The gradient in the cloud is a consequence of the periodic boundary conditions whichmust balance the CR pressure across the ISM-cloud interface. The profile is consistent with predictions from CR fluid theoryin a periodic domain. where we have used the approximate scattering ratefrom Equation 5. We set p = p above since we aremost interested in the peak of the distribution wherethe majority of particles are present. In the top panelof Figure 2, we show the mean free path in the ISM, λ = 0 . L , based on wave amplitudes averaged overthis region. The mean free path in the cloud is ≈ L ;particles free stream with no significant wave-particlecollisions.At the time t = 5 × v A Ω − c in the Fiducial simula-tion, particles traverse approximately 2 mean free pathsin the ISM. This mean free path is not particularly short,perhaps calling into question the validity of a diffusionapproximation for particle transport. Based on trackingparticle trajectories, we find that indeed even thoughCRs undergo many small angle scattering events, theyexperience few direction reversals over the simulationtime. Yet, gradient structures are nonetheless able todevelop in the ISM and cloud. Thus, while particle propagation may not be strictly diffusive in real space,exhibiting a random walk in x , CRs in our simulationsevolve diffusively in momentum space. In this way, even2 mean free paths in our simulations can roughly capturethe diffusion process.We can gain further insight into the formation of thegradients by studying the space-time evolution of theCR moments (Figure 3). Before t = 4 × Ω − c , thesimulations exhibit transient behavior, with fluctuationsin the energy density and flux on scales of the cloud.Fluctuations in the energy density are asymmetric aboutthe ISM-cloud interface. Furthermore, the cosmic raypressure gradient develops slightly earlier in the ISMregion, indicating that the cloud gradient is a responseto the pressure gradient in the diffusive transport region.For a brief interval around t = 3 × Ω − c , a peakin the number/ energy density forms on the upstreamside of the cloud. This peak forms immediately after thewave amplitudes exceed E A /ρ v A ∼ − , i.e. when the HD-PIC Simulations of Molecular Clouds Figure 3.
Space-time evolution of wave energy density E A , perturbed particle number δn CR , w /n CR , and energy δ E CR , w / E CR , densities, and flux ∆ F CR ,w . The fluctuation in the CR energy flux is defined ∆ F CR ,w ( t ) ≡ F CR ,w ( t ) − (cid:104) F CR ,w ( t ) (cid:105) x , where theangle brackets denote an average over the simulation box at a given time t . The growth of Alfv´en waves leads to a reduction inthe mean free path in the ISM region, resulting in modulations in CR energy flux and densities. After a short transient from ∼ − × Ω − c , the distribution adjusts to form an approximately spatially-constant flux that declines over time. In thissteady state configuration, a negative (positive) energy density gradient forms in the ISM (cloud), consistent with fluid theory. scattering rate becomes shorter than the current simu-lation time. Thus, this peak is a transient effect whichforms just as particles begin scattering back into thecloud. The peak is soon erased as forward-propagatingparticles crossing the ISM are scattered away from theleading interface, unable to replenish the deficit left byparticles scattered forward into the cloud.4.2. Origins of the Gradient: Kinetic Perspective
The advantage of treating CRs as kinetic particles inthe PIC method is that we are able to obtain the fulldistribution function at any given time or location in thesimulation. We choose to output slices in the distribu-tion function in position bins averaged over a width of1.25 × v A Ω − c , corresponding to 16 position bins inthe Fiducial and
NCR2 simulations and 40 bins in the
Long Cloud and
Long ISM runs. Figure 4 shows the space-time evolution of momen-tum bins of the wave-frame distribution function. Wechoose to study three momentum bins, corresponding tothe peak momenta of the n CR ,w , E CR ,w , and F CR ,w inte-grals, i.e. log ( p w /p ) = 0.19, 0.45, and 0.55 respectively.For ease of comparison, we study the quantity δf w /f κ ,where δf w = f w − f κ . This perturbation is slowly erasedas the distribution function relaxes toward an isotropic,approximately κ distribution in the wave frame.Initially, δf w /f κ increases linearly with µ w , antisym-metric about µ w = 0, i.e. the distribution has anet flux (Row 1). As waves grow and scattering be-gins (Row 2), the distribution evolves differently atthe leading and trailing interfaces, respectively down-stream and upstream from the ISM. At the leadinginterface (central two columns), particles with posi-tive pitch angle arrived by crossing through the ISM.These forward-propagating particles were scattered to-2 C. J. Bambic, X. Bai, & E. C. Ostriker -0.100.000.10 f w f Time = 0.00e+00 x [0.0,12.5] x [87.5,100.0] x [100.0,112.5] x [187.5,200.0]-0.100.000.10 f w f Time = 2.00e+04 -0.100.000.10 f w f Time = 3.00e+04 -0.050.000.05 f w f Time = 5.00e+04 w -0.030.000.03 f w f Time = 7.50e+04 w w w ISM Cloud
Figure 4.
Slices of the wave frame perturbed distribution function δf w = f w − f κ at four locations near the ISM-cloud interfaces.Three momentum bins are shown, corresponding to the peak of the number density (light blue), energy density (orange), andenergy flux (black) integrals, at log ( p w /p ) = 0.19, 0.45, and 0.55 respectively. As time increases (top to bottom), the waveframe distribution moves from an anisotropic drifting state toward an isotropic, µ -independent state except near µ = 0. ward smaller pitch angle, leading to the development ofa shelf-like structure near µ w = 0 . µ w = 1 (Row 3). A similar phenomenaoccurs at the trailing interface. Backward-propagatingparticles are scattered toward more negative pitch an-gle, increasing the deficit near µ w = − . µ w = −
1. The deficit near µ w = − . ◦ barrier.As time progresses, a common story emerges in theISM. We focus on forward-propagating particles sincethey are the dominant population. Particles enteringfrom the trailing cloud-ISM interface ( x = 0) are scat-tered as they cross the ISM. Upon arriving at the lead-ing cloud interface ( x = 100 in units 10 v A Ω − c ), theyare suddenly able to free stream into the cloud, rapidlyvacating the interface region. This sudden transitionleads to a net deficit just upstream of the leading ISM-cloud interface, and a steep gradient forms near x = 100(Row 4). Particles built-up near µ w (cid:38) ◦ barrier as they crossthe ISM, providing fewer particles to replenish the deficitleft by the now free-streaming CRs entering the cloud;the gradient is enhanced. These effects are compoundedby the continuing decrease in the CR flux, providingfewer and fewer particles to fill in the deficit with time. Thus, we are left with an excess at x = 0 and a deficitat x = 100. This process is the essence of the spatialdiffusion which yields a gradient in the CR number andenergy densities.The positive gradient in the cloud is formed by a differ-ent process. We see the emergence of a peak in the dis-tribution function of forward-propagating particles near µ w = 0 as particles cross the cloud. Since low amplitudewaves in the cloud cannot scatter particles at any sub-stantial rate, the distributions in the cloud are purelyinherited from the ISM. The peak apparent in the right-most column of Figure 4 can be explained by consid-ering the crossing time of particles through the cloud.CRs with pitch angles of µ w = 0.1 have 10% of the par-allel velocity of those with µ w = 1 for fixed momentum.Thus, particles near µ w = 0 entered the cloud at an ear-lier time in the simulation, when the overall CR flux washigher. CRs near | µ w | = 1 entered the cloud recentlyby comparison, after the CR flux dropped substantially.Thus, we are seeing an overlap of wave-frame distribu-tions at different times, and this overlap effect becomesincreasingly clear near the trailing interface ( x = 0).By studying the distribution functions, we find a ki-netic explanation for the same behavior predicted byfluid theory. A negative spatial gradient across the ISMforms in response to wave-particle scattering. The pos- HD-PIC Simulations of Molecular Clouds = 0.26 L × + × w / CR, 0
Time = 5.00e+04 Fiducial = 0.10 L × + × w / CR, 0
Time = 2.50e+04 NCR2 = 0.06 L × + × w / CR, 0
Time = 5.00e+04 Long Cloud = 0.14 L v A 1c )3.0 × + × w / CR, 0
Time = 5.00e+04 Long ISM
Figure 5.
Energy density profiles and mean free paths for all simulations. Times are selected such that the gradient ismaximized in the simulation. In all simulations, the peak of the energy density perturbation is near the upstream side of theISM (trailing end of cloud) x = 0 while the minimum appears downstream from the ISM (leading end of cloud). The NCR2 and
Long Cloud gradients in the ISM are the largest, implying that larger wave amplitudes (scattering rates) yield more substantialgradients, consistent with fluid theory. itive gradient across the cloud is a response to the dropin CR energy flux, since particles further in the cloud forthe same momentum and pitch angle entered the cloudat an earlier time when the overall CR flux was higher.Together, these processes work to create a pressure bal-ance of the CR “fluid” across the contact discontinuitybetween the ISM and cloud.4.3.
Comparing Spatial Structure among Simulations
Before using our simulations to compute the wave-particle scattering rates, we turn our attention to thestructure and magnitude of the spatial gradient in CRenergy density, which provides an insight into the mag-nitudes of the wave-particle scattering rates.Figure 5 displays the energy density spatial structureand mean free paths in the ISM for all simulations attimes when the gradient is near maximal. The
Fiducial simulation contains the fewest mean free paths acrossthe ISM ( ∼ Long ISM simulation containsnearly 6 mean free paths.The overall magnitude of the fluctuation in CR energydensity is connected to the amplitudes of Alfv´en waves.This is evidenced by the fact that the simulations withthe largest amplitude waves (
NCR2 and
Long Cloud )have energy density fluctuations more than a factor of 3 larger than the
Fiducial run, despite possessing an ISMregion of equal size. Similarly, the fluctuations in the
Fiducial and
Long ISM runs are approximately equal,consistent with the similarity in Alfv´en wave energy den-sity between these simulations (Figure 1).Comparing the steepness of gradients in the ISM andcloud for a given simulation provides insight into therelative importance of each term in the fluid equation(18). For the
Fiducial and
NCR2 runs, the gradientsmust be equal in magnitude in the ISM and cloud, be-cause they are of equal length, L ISM = L Cloud . In thecloud region, the energy density gradient must be equalin magnitude (and opposite in sign) to the time deriva-tive in energy flux because the scattering rate is negligi-ble. Therefore, the two terms on the LHS of Equation 18must be equal in magnitude for these two models. Fol-lowing similar reasoning, one would expect the flux timederivative term to be ∼ / Long Cloud model, and the flux time derivativeterm to be 4 times the pressure gradient term for the
Long ISM model. Thus, modulating the ISM-to-cloudlength ratio, L ISM /L Cloud explores the relative impor-tance of each term and provides a test of robustness forthe fluid theory when we compute the fluid scatteringrate in Section 5.2.4
C. J. Bambic, X. Bai, & E. C. Ostriker k I ( k ) Critical
Time = 5.0 ×10 Forward RightForward LeftBackward RightBackward Left
NCR2
Time = 2.5 ×10 k x k I ( k ) Long Cloud
Time = 5.0 ×10 k x Long ISM
Time = 5.0 ×10 Figure 6.
ISM Alfv´en wave power spectra at times when the CR energy density gradient is near maximal. For plotting purposes,these spectra are smoothed on a length scale of 100 k L , where k L v A / Ω c = 2 π / L = π and π × − for the Critical / NCR2 and
Long Cloud / Long ISM simulations respectively. Note that this smoothing length is 5 × larger than that used for computingthe effective quasilinear scattering rate. Forward propagating waves contain the majority of power, shared equally between leftand right polarizations. The self-generated wave spectrum from our κ distribution of CRs has a slope I ( k ) ∼ k − at scales justbelow the peak scale near kv A / Ω c = v A / p ≈ × − . So far we have focused on the energy density struc-ture in the ISM, which always has a negative gradient(as expected for a diffusive region with positive flux);however, a different process operates in the cloud. Forall of our models, the energy density shows a consis-tent positive gradient across the cloud. This increase isa consequence of the continual decrease in energy fluxwith time in all the simulations. As the CRs approachisotropy, the rate of change of the flux decreases (Fig-ure 1b), and the energy density gradient in the cloud isslowly erased (Figure 3). In the present simulations, thetemporal decrease in flux is a consequence of periodicboundary conditions. In the real ISM, we would insteadexpect an approximate steady state for the flux to be reached. A state with a temporally constant flux andnegligible scattering within a dense, neutral cloud (dueto strong local damping) would not have a spatial gra-dient in the energy density within the cloud. Rather, adownward energy density “ramp” in the diffusive ISMregion would be followed by an energy plateau withinthe cloud. The emergence of the CR density gradientsin the cloud, while entirely consistent with fluid theory,is an artifact of the periodic boundary condition in thepresent simulations. SCATTERING RATE COMPARISONIn this section, we present a numerical verification ofthe moment equation which underpins CR hydrodynam-
HD-PIC Simulations of Molecular Clouds ( p w )( c ) Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 QuasilinearFluidWave
NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 p w / p ( p w )( c ) Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 p w / p Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Figure 7.
Effective scattering rates computed from quasilinear theory (Equation 20; orange) and fluid theory (Equation 49;dark blue) as a function of momentum for all simulations. The naive estimate based on the average Alfv´en wave energy density(Equation 5; light blue) is over-plotted for comparison. Curves are smoothed in momentum for plotting purposes. We show thetime evolution of the scattering rates by varying the transparency of the curves with time; the most transparent lines are fromthe earliest measurement time and the darkest lines are from the latest measurement time within the range displayed. ics. We compare the wave-particle scattering rates aspredicted by quasilinear and Fokker-Planck theory toshow that these methods agree well with CR hydrody-namics when the mean free path to scattering is suffi-ciently short. Since the scattering rate encodes spatialdiffusion, this section serves as first-principles confirma-tion of the quasilinear calculation for the parallel CRdiffusion coefficient in fluid theory. Deviations from thequasilinear prediction point to nonlinear effects, partic-ularly near the µ = 0 barrier (see Sections 6.2 and 6.3).5.1. Effective Quasilinear Scattering Rate
We compute both the effective quasilinear and fluidscattering rates using the full distribution function mea-sured in the wave frame f w . Equation 20 provides the procedure for computing the effective scattering ratefrom quasilinear theory. Following BOPS19, we de-compose our Alfv´en waves into left/right-handed andforward/backward-propagating modes, 4 power spectrain all. We compute these power spectra based on thewaves in the ISM only (zero-padding the cloud) andnormalize the spectra according to Equation 8. Thus,the integrated power spectra together equal the quantity( δB ⊥ /B ) averaged over the ISM. Power spectra arethen smoothed over a scale of 20 k L , where k L = 2 π/L is the wavenumber corresponding to the box size L . Wefind that less smoothing of the power spectra in k -spaceyields too much noise in the scattering rate which causesthe integral in Equation 20 to be poorly behaved. Thechosen smoothing does not change the shape of the spec-6 C. J. Bambic, X. Bai, & E. C. Ostriker trum. Smoothed Alfv´en wave spectra for all simulationsare shown in Figure 6. This procedure allows us to com-pute the quasilinear scattering rate through Equation 7.To obtain the effective scattering rate from ν QL and f w , we smooth f w and average the distribution over theISM, omitting the spatial bins nearest the ISM-cloudboundaries. This omission eliminates particles trappednear µ = 0 (see Figure 4 bottom right panels) whichcan bias the µ -derivative near µ = 0; however, includingthese bins simply increases the noise in the scatteringrate integral. We compute the µ -derivative of f w usinga centered finite difference and perform the integrals inEquation 20 over pitch angle and gyrophase only, leavingscattering rates as a function of momentum. The sameprocedure holds for the energy flux in the denominatorof ν eff according to Equation 20. The results are shownas the orange curves in Figure 7.5.2. Fluid Scattering Rate
The fluid scattering rate is computed from the mo-ments of f w . First, the energy-weighted moments of thedistribution function are calculated according to Equa-tions 42 and 43. The integrations are performed overpitch angle and gyrophase alone such that the resultingexpressions represent infinitesimal moments as a func-tion of momentum. These moments are then insertedinto the fluid equation, Equation 18. The time deriva-tive of the energy flux is computed by measuring theflux at 200 output times over the simulation time of10 Ω − c , corresponding to an interval of 500 Ω − c be-tween outputs. Then, a second order accurate centeredfinite difference is used to compute the time derivative.The energy density gradient is computed by performinga linear fit to the energy density as a function of position x in the ISM, comprising 8 spatial bins for the wave spec-trum output in the Fiducial and
NCR2 simulations, and20 spatial bins in the
Long Cloud and
Long ISM runs.The slope of this fit is the energy density gradient. Wehandle the energy flux term on the RHS of Equation 18by averaging the energy flux over the ISM. RearrangingEquation 18 then yields the following form for the fluidscattering rate, ν fluid ( p w ) = − ∂F CR ∂t ( p w ) + v w ∂ E CR ∂x ( p w ) F CR ( p w ) , (49)where ν eff has been replaced by ν fluid to eliminate ambi-guity with the effective quasilinear scattering rate. Thisfluid scattering rate as a function of momentum for allsimulations is shown in Figure 7. Note that each term inEquation 49 can be computed from the CR energy den-sity gradients (Figure 5) and CR flux. The time rate-of-change in the flux must be balanced by the energy den-sity gradient in the cloud region, where ν eff ≈
0. Thusthe cloud gradient provides the time derivative termwhile the ISM gradient provides the position derivative.In general, we find good agreement between the ef-fective scattering rates computed via fluid and quasi- linear theory. The scattering rates do not differ bymore than 50% over more than an order of magnitudein momentum, 1 < p w /p <
20. For much of this range(1 < p w /p <
10) the deviation is within 25%.Except in the case of very strong scattering (simula-tion
NCR2 ), the fluid and effective quasilinear scatteringrates are lower than the naive estimate for the scatter-ing rate (Equation 5; light blue line in Figure 7) by afactor of order unity. The naive estimate, which relieson a strict random walk in pitch angle, is most accuratein the regime of strong scattering when particles are lessaffected by the µ = 0 barrier (see Section 6.2).Given that the quasilinear scattering rate ν QL ( p w , µ w )and certainly the distribution function suffer from par-ticle discreteness noise while the fluid scattering rateincludes nonlinear effects (see Sections 6.3), these minordeviations are within the range of accuracy we expectto achieve from our PIC simulations. In addition, wenote that the integral in Equation 20 is poorly behavedand runs into a “sign problem” for lower particle mo-menta, yielding a negative effective scattering rate—anunphysical result. These scattering rates are omitted inFigure 7. The agreement between the scattering ratesindicates that, for particles with momenta near the peakof the distribution where the mean free path is shorterrelative to the tails of the distribution, the quasilinearprediction is borne out by the fluid behavior of the CRs.5.3. Fokker-Planck Coefficient
The quasilinear diffusion equation (6) is a Fokker-Planck equation with the quasilinear scattering rateweighted by (1 − µ ) / D µµ . Because of this correspondence, weexpect particles to diffuse in pitch angle at a rate givenby Equation 11. In this subsection, we present our pro-cedure for tracking individual particles and computingtheir scattering rates. This method provides a means ofcomparing the quasilinear scattering rate (Equation 7)with the Fokker-Planck rate. In this way, we test thevalidity of quasilinear theory in a discrete sense, ratherthan the statistical ensemble captured by fluid theory.We track a total of 3200 particles in each simula-tion. These particles are sampled equally from 200spatial bins (16 particles per spatial bin). The rangeof momenta in the drifting (initially isotropic) frame, − < log ( p d /p ) <
2, is divided into 8 bins. We samplean equal number of particles from each of these 8 bins,drawn from the initial κ distribution. Thus, the overallsampled distribution is not a κ distribution, but rather8 distinct momentum ranges within which the particlesobey a κ distribution (400 particles per momentum bin).In pitch angle, we sample from the initially flat driftframe distribution, tracking forward and backward par-ticles of the same | µ d | . For instance, if we track particle i with phase space coordinates ( p d,i , µ d,i , x i ), we also tracka particle with coordinates ( p d,i , − µ d,i , x i ). This proce-dure ensures an equal number of forward and backward- HD-PIC Simulations of Molecular Clouds Figure 8.
Left: Tracked particle orbits from
Fiducial simulation evolving in space (top), pitch angle (middle), and scatteringrate (bottom). Right: Fokker-Planck scattering rate at time t = 5 × Ω − c as a function of wave-frame momentum p w andpitch angle µ w . Note that the top left plot is a space-time plot, where position is on the x -axis, consistent with the simulation.Particles are selected such that they appear near the trailing cloud-ISM interface ( x = 0) at time t = 5 × Ω − c . Orbits evolvewithin the periodic domain, scattering in the ISM region and free streaming through the cloud (grey region). All scatteringrates shown are restricted to the ISM only, as cloud scattering is negligible. propagating particles in the drift frame, which whenboosted to the wave frame, yields more particles at pos-itive than negative pitch angle in the wave frame.Particles’ full phase space coordinates (position andmomentum) are output every ∆ t F P = 50Ω − c , implyingthat particles with p ≈ p undergo 50/ γ ≈
35 gyrationsabout the magnetic field between outputs. If a particle i is measured at a time t , we compute the scattering rateaccording to Equation 11 as ν F P ( t, p w,i , µ w,i ) = ( µ w,i ( t + ∆ t F P ) − µ w,i ( t )) t F P . (50)Figure 8 displays particle positions, pitch angles, andscattering rates as a function of time for a selection of6 particles in the Fiducial simulation which undergo a90 ◦ pitch angle crossing. We select particles such thatthey arrive near the peak of the CR energy density gradi- ent at t = 5 × Ω − c from the bin 0 < log ( p w /p ) < × − Ω c , althoughrates can fluctuate up to 10 − Ω c for the particlestracked. A better sense of the full scattering rate distri-bution is shown in the right panel of Figure 8. Clearly,3200 particles is a sparse sampling of the full distribu-tion function (this is seen most clearly by juxtaposingthis figure with Figure 9 of BOPS19). Yet, we still seeclear patterns emerging, namely an increase in scatter-ing rate near the peak of the distribution around p anda decrease in scattering rate with increasing momentum.Unfortunately, the sparseness does not allow us to seeconstant scattering rates along resonant lines k res = con-stant (dotted black lines); however, we see an approxi-8 C. J. Bambic, X. Bai, & E. C. Ostriker ( p w )( c ) Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 Fiducial
Time = 4.5 - 5.0 ×10 QuasilinearFokker-PlanckWave
NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 NCR2
Time = 2.0 - 2.5 ×10 p w / p ( p w )( c ) Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 Long Cloud
Time = 4.5 - 5.0 ×10 p w / p Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Long ISM
Time = 4.5 - 5.0 ×10 Figure 9.
Ensemble averaged scattering rates computed from Fokker-Planck theory (Equations 50 and 51; dark red) andquasilinear theory (Equations 7 and 52; orange) as a function of momentum for all simulations. These scattering rates agreeremarkably well across nearly 3 orders of magnitude in momentum, with exquisite concordance at the peak of the CR momentumdistribution. Similarly, these scattering rates agree well with the naive estimate (Equation 5; light blue) at this peak. mately constant scattering rate as a function of µ w forfixed momentum. The exception to this pattern is near µ w = 0, where the scattering rate drops to ∼ − Ω c .With the scattering rate distribution shown in theright panel of Figure 8, we are able to compute a pitchangle-averaged scattering rate which we can compare tothat from quasilinear theory (Equation 7). We definethis ensemble average, (cid:104) ν F P (cid:105) ( p w ) = (cid:82) (cid:82) ν F P ( p w , µ w , x ) f w ( p w , µ w , x ) dµ w dx (cid:82) (cid:82) f w ( p w , µ w , x ) dµ w dx . (51)Here, the integral in µ w is over the full range of pitchangle, − x is over the ISMregion. Note that as was done with the effective quasilin-ear scattering rate, we omit the spatial bins nearest the cloud-ISM boundaries since particles retain their distri-bution from the free-streaming cloud in these regions.Because the distribution is sparse, we set all momen-tum and pitch angle bins in f w and ν F P to 0 where noparticle is present in that element of phase space.Comparing this Fokker-Planck rate to quasilinear the-ory requires a different weighting of the quasilinear scat-tering rate ν QL . Transforming to pitch angle coordinatesintroduces a geometric factor of (1 − µ w ) / (cid:104) ν QL (cid:105) ( p w ) = (cid:82) (cid:16) − µ w (cid:17) ν QL ( p w , µ w ) f w ( p w , µ w ) dµ w (cid:82) f w ( p w , µ w ) dµ w , (52) HD-PIC Simulations of Molecular Clouds ν QL is given by Equation 7. Note that because the quasi-linear scattering rate is a function of the Alfv´en wavepower spectrum as measured over the entire ISM, thescattering rate is already averaged over x . In the aboveintegral, we average f w over x , again omitting spatialbins nearest the cloud-ISM boundaries.The ensemble averaged Fokker-Planck scattering rateis compared to the suitably weighted quasilinear scatter-ing rate in Figure 9. Across a wide range of momenta,nearly three orders of magnitude, these scattering ratesagree remarkably well, with near perfect agreement atthe peak of the CR momentum distribution. In addition,these scattering rates agree well with the naive scatter-ing rate estimate from Equation 5. Some deviations arepresent, most notably at momenta below p w /p = 0.5and above p w /p = 20, where the Fokker-Planck scat-tering rate is higher than the quasilinear prediction. Weaddress the question of the validity of quasilinear theoryas well as the robustness of fluid models of CR transportin the next section. DISCUSSIONBy breaking translational symmetry with spatially-dependent ion-neutral damping, our simulations haveenabled an exploration of fluid behavior in a collisionlessCR population under the sole influence of self-inducedwave-particle interactions. This result is one of a grow-ing number of studies in which collisionless wave-particleinteractions act to replace the particle-particle collisionsunderpinning transport in MHD fluids. Whether in col-lisionless shocks (Spitkovsky 2008), kinetic turbulence(Howes et al. 2008; Meyrand et al. 2019), or MRI-unstable shear flows (Kunz et al. 2016), wave-particleinteractions regulate transport, heating, and fluid-scalestructure. Yet, while fluid behavior might emerge, theprecise value of transport coefficients may deviate sub-stantially from collisional or weakly collisional predic-tions (Spitzer 1962; Braginskii 1965) or take on a func-tional form ill-suited to fluid models (Arzamasskiy et al.2021, in prep.; Kunz et al. 2014a).Using our simple toy model for a sharp boundary be-tween a mostly-neutral cloud and well-ionized plasma,we directly demonstrate spatiotemporal behavior consis-tent with that expected from the energy flux equation ofCR hydrodynamics. The effective scattering rates in thefluid theory are consistent with the fully kinetic quasi-linear prediction and supported by studies of individualparticle motions. Thus, for the restrictive case of field-parallel CR diffusion subject only to scattering fromsmall-scale, self-generated Alfv´en waves for our chosenparameters, the quasilinear prediction for the diffusioncoefficient is likely accurate. In this section, we addressthe validity of the quasilinear approximation, the roleof nonlinear effects in transport, and how fluid theoryremains so robust in our simple system. 6.1.
Validity of Quasilinear Theory
Quasilinear theory, when applied to the gyroresonantCRSI, treats the growth of Alfv´en waves based on astatic, drifting distribution function f ( p ), i.e. through linear theory. The growing waves induce scatteringamong the particles, modifying the distribution func-tion. In general, the evolution of the distribution func-tion is a complex, fully nonlinear problem computed nu-merically in our simulations; however, if particle andwave dynamics occur on disparate timescales, we canuse a scale separation technique to evolve the distri-bution function. Quasi -linear theory assumes that theoverall distribution function evolves slowly compared tothe dynamical timescale for the waves. In essence, itposits a separation of timescales, kv A (cid:29) | ∂ ln f /∂t | .This relation is usually well satisfied for small waveamplitudes, as the ratio between the two scales is( v A /c )( δB/B ) − (cid:29) δB/B ∼ − , c = 300 v A .A critical assumption underlying quasilinear theory isthe “random phase approximation”, i.e. fields are delta-correlated in k -space and time (Kulsrud 2005). Whenthis assumption is satisfied, particles experience Gaus-sian white noise forcing and freely diffuse in momentumspace as a central limit effect. Quasilinear theory ap-plies as long as particles encounter sufficiently many un-correlated wave packets such that they undergo chaoticdiffusion in momentum space (see Besse et al. (2011)for further discussion). In our simulations, this condi-tion is met when the mean free path is short relative tothe ISM scale. As long as multiple mean free paths fitwithin the ISM, particles experience multiple scatteringevents while traversing the ISM and chaotic diffusionensues. Nonlinear wave-wave interactions may distortthe waves and lead to complications which violate therandom phase approximation (see Section 6.3).Away from the peak of the CR momentum distribu-tion, the mean free path to scattering is longer thanthe ISM scale. The quasilinear approximation is inap-plicable for these particles, a fact best exemplified inthe flattened structure in the fluid scattering rate near p w /p = 10 in the Fiducial simulation (Figure 7). Thistime-dependent feature moves toward higher momentumas the simulation is run longer, indicating that the struc-ture is a consequence of high momentum particles lack-ing the time to respond to wave growth through manyencounters with waves. Similarly, the deviations at lowmomentum between the fluid and quasilinear rates inFigure 7 as well as the Fokker-Planck and quasilinearrates in Figure 9 point toward sources of scattering notcaptured by a quasilinear treatment.Despite these deviations for particles away from thedistribution peak, (1) our tracked particles undergochaotic diffusion in pitch angle, (2) the ensemble aver-aged Fokker-Planck and weighted quasilinear scatteringrates agree well near the peak, and (3) our Alfv´en wave0
C. J. Bambic, X. Bai, & E. C. Ostriker amplitudes remain small. Thus, we argue that quasilin-ear theory is a valid approximation for the growth andsaturation of the gyroresonant CRSI in our simulationsfor CRs near the peak of the distribution.6.2.
Implications for CR Hydrodynamics
Transport in CR fluid theory is encoded in the spatialdiffusion coefficient, the inverse of which (times c ) cor-responds to the fluid scattering rate that we measure.For CR fluid theory to accurately describe the system’sdynamics, this diffusion coefficient should incorporateboth quasilinear and nonlinear scattering mechanisms.The agreement between our quasilinear predictions andfluid scattering rates requires efficiently crossing the 90 ◦ pitch angle ( µ = 0) barrier.Overcoming this barrier goes beyond quasilinear the-ory and requires nonlinear effects (see Section 6.3). Pa-rameters chosen in our simulations are far from realis-tic. In particular, the ratio of n CR , /n i ∼ − is highlyexaggerated compared to the ratio of CR to thermalparticle density in the Galaxy, ∼ − (although in low-ionization regions n CR /n i is much higher). While thisexaggeration is necessary to make our simulations com-putationally feasible, the large wave amplitudes unreal-istically enhance nonlinear effects. A more quantitativeunderstanding of the potential dependence of the scat-tering rate on wave amplitude is needed to confirm thatthe quasilinear rate is applicable in CR fluid treatmentsfor realistic ISM environments.In fact, we can already identify deviations of the effec-tive scattering rate (Equation 20, based on the quasilin-ear theory) from the fluid scattering rate (Equation 49,which incorporates nonlinear effects) in Figure 7, for CRmomenta p (cid:46) p . We have noted that calculation of theeffective scattering rate from Equation 20 suffers fromnumerical noise, but the trend is already evident. Thiseffective scattering rate is weighted by a gradient in thepitch angle distribution, which is maximized near µ = 0(Figure 4) where quasilinear scattering rates vanish.Since this effective rate would be formally equal to thefluid scattering rate if quasilinear theory fully describedwave-particle interactions, deviations between the effec-tive and fluid rates point to the key role nonlinear effectsplay in determining the total fluid scattering rate, andsubsequently, diffusion. We can conclude that nonlin-ear effects start to dominate the overall fluid scatteringrate for particles with p (cid:46) p in our simulations. Eventhough realistic wave amplitudes would be lower thanthose in our simulations, nonlinear effects may also beimportant in realistic ISM environments, as has alreadybeen suggested from some recent simulations of galaxyformation (e.g., Hopkins et al. 2020).We also comment that when studying the Fokker-Planck scattering rate, weighting in pitch angle is uni-form; the role of nonlinear effects in overcoming µ = 0 isnot as significantly manifested. Consequently, there isbetter agreement between the quasilinear and Fokker- Planck rates in Figure 9. Despite this agreement, thedeviation between the two rates becomes more signifi-cant for particles with p (cid:46) . p (see Section 6.3), indi-cating that nonlinear effects dominate over a wide rangeof pitch angles.Further research is necessary to elucidate the signifi-cance and nature of these nonlinear effects in the hopeof developing transport equations faithful to all sourcesof wave-particle scattering. A first step can come fromMHD-PIC simulations which yield numerical estimatesof spatial diffusion; however, analytic or semi-analytictechniques may be necessary to extend numerical insightinto the regimes relevant for the Galaxy.6.3. Nonlinear Sources of Scattering
We know some nonlinearity must be present for par-ticles to overcome the µ = 0 barrier. Here, we discusssome of the most widely considered mechanisms thatmay contribute to the enhancement of scattering overthe quasilinear prediction in our simulations.Our wave spectrum is determined by the initial CRdistribution, not a scale-by-scale transfer of energy aswould be expected within an MHD turbulent cascade.While in general, wave-wave interactions are weak atlow amplitude, we do observe, as reported in Plotnikovet al. (2021, in prep.), that the spectrum evolves toinclude high- k modes, eventually achieving a power-lawspectrum in intensity. The origin of this cascade is yetto be understood; however, the effect is present evenfor wave amplitudes which remain well within the linearregime.One important consequence of this spectral evolutionis generation of abrupt features analogous to rotationaldiscontinuities in transverse magnetic fields. Such fea-tures, reported in Plotnikov et al. (2021, in prep.) andpresent in BOPS19 as well as this work, might be a con-sequence of nonlinear wave steepening into rotationaldiscontinuities (Cohen & Kulsrud 1974). Particles en-countering such abrupt features effectively see a suddenchange of field direction and hence a sudden change in µ relative to the perturbed field. This scattering mecha-nism is generally insignificant given our small wave am-plitudes, but it can become significant when µ is closeto zero. A reflection can in principle be achieved by | µ | (cid:46) δB/B . This is identified in BOPS19 as the domi-nant mechanism for overcoming the µ = 0 barrier, and islikely also the mechanism responsible in our simulations.More generally, wave-wave interactions induced bycascades couple modes in k -space, yielding correlatedspectral and temporal structure in the waves, relevanton the small spatial (high- k ) scales of low-momentumgyroresonant particles. These correlations break the as-sumption of delta-correlated fields (the random-phaseapproximation) underlying quasilinear theory (Kulsrud2005), and particles no longer experience Gaussian whitenoise forcing. The distribution function and wave spec-trum co-evolve on similar timescales, and wave-particle HD-PIC Simulations of Molecular Clouds µ = 0 (Dupree 1966; Weinstock1969; V¨olk 1973; Achterberg 1981). Since diffusion co-efficients are time integrals over correlation functions(Shalchi 2009), nonlinear effects such as wave-wave in-teractions and spatially-localized rotational discontinu-ities modify the scattering rates measured in Sections 5.2and 5.3 away from the quasilinear prediction.As is pointed out by Holcomb & Spitkovsky (2019),the singularity in the resonant wavenumber in Equa-tion 9 is strictly artificial and can be removed by relax-ing the magnetostatic approximation ( ω/k ∼ µ = 0.Mirror scattering may play some role in crossing µ = 0(Felice & Kulsrud 2001; Holcomb & Spitkovsky 2019);however, this mechanism is likely subdominant in ourone-dimensional simulations. Gradients in the magneticfield can form magnetic mirrors which adiabatically scat-ter particles via the mirror force. This effect is non-resonant and does not break conservation of magneticmoment, in contrast to the aforementioned rotationaldiscontinuities. As variation in the field strength in oursimulations is due only to wave motions, we never formlarge spatial gradients in the magnetic field and thusnever generate significant mirror-like structures. Mirrorscattering would likely gain greater significance in mul-tiple dimensions and requires future studies of the CRSIwhich go beyond 1D. In addition, mirror structure maynaturally be present in background MHD turbulence,associated with transit time damping from fast magne-tosonic modes (e.g. Schlickeiser & Miller (1998), Yan &Lazarian (2004)).6.4. Future Directions
While we have extended the BOPS19 MHD-PIC sim-ulations to a more realistic environment, we remainfar from the conditions relevant to the real ISM. Per-haps most crucially, our problem was studied with pe-riodic boundary conditions. Thus, rather than achievea steady state energy flux of CRs, we were forced toaddress a time-dependent problem in which the energyflux is continuously decreasing, leading to the formationof a non-physical positive spatial gradient in the cloudregion. While this structure is consistent with CR hy-drodynamics, it would not develop for a (quasi) steadystate in which wave growth is balanced by damping. Modifying the boundary condition to a constant fluxof CRs entering the simulation domain is a key next steptowards greater realism for the problem we have studied.This boundary condition would establish a constant CRenergy density gradient and therefore a constant scat-tering rate for the ISM region. If waves are stronglydamped within the cloud, we would expect constant en-ergy density in this region. Such a system could providea laboratory for measuring the spatial diffusion coeffi-cient as a function of n CR , /n i , CR pressure gradient, ν IN /ν crit , and L ISM /L Cloud .Even with this modification, the problem of CR trans-port in and around GMCs remains a challenge. Threedimensional turbulent magnetic fields, CR energy lossesdue to H impact ionization (McCall et al. 1998), andhadronic losses from CR impacts which produce gammarays through pion decay in the GeV (Yang et al. 2014;Tibaldo et al. 2015) and TeV (Aharonian et al. 2006;HESS Collaboration et al. 2016; H. E. S. S. Collabora-tion et al. 2018) bands may all contribute to controllingCR transport, energy densities, and ionization rates inGMCs. Thus, a full treatment of this problem requiresnot only the plasma physics of CR transport, but GMCchemistry and CR energy loss mechanisms as well.In the absence of three dimensional effects, additionalexternally-driven MHD turbulence, or CR energy losses,and for numerically expedient but unrealistic parameterregimes, our work remains but a first step in understand-ing CR transport in the multiphase ISM. Yet, by takinga first-principles approach for CR diffusion and compar-ing to fluid treatments, our work opens a path towardfuture studies of CR transport coefficients in realisticenvironments. In this way, studies of CRs in systemsfrom GMCs to galaxy clusters can explore astrophysi-cal macroscales while remaining firmly grounded in theplasma physics of the CR microscales. CONCLUSIONWe have presented the first self-consistent kineticsimulations of CR transport across an inhomogeneousdomain, modeling an embedded neutral cloud withinthe ionized ISM. By breaking translational symmetrythrough ion-neutral damping of Alfv´en waves in thecloud region, our simulations enable us to see aspectsof fluid behavior in the spatial structure of the CR dis-tribution function. In particular, we show that the sim-ulation results are consistent with the predictions of CRhydrodynamics, in which an energy density gradient andtime-dependent energy flux work together to balancewave-particle scattering throughout the ISM. In the ISMregion, the gradient in the energy density is in the op-posite direction to the net CR flux. We can understandthe ISM energy gradient as the consequence of diffusivepropagation imposed by wave-particle scattering.In the cloud region, where there is negligible scatter-ing, the gradient in the energy density is in the same di-rection as the net CR flux, since the decrease in time of2
C. J. Bambic, X. Bai, & E. C. Ostriker the flux must be directly balanced by a spatial gradientof the pressure. This behavior is a consequence of ourperiodic boundary conditions which ensure that energyflux is a constantly decreasing function of time. Thecloud contains a superposition of propagating CR dis-tributions, unchanged after they first entered the cloud.Structure in the ISM energy density allows us to com-pute a wave-particle scattering rate based on the spatio-temporal evolution of the CR moments—a fluid ap-proach. We compare this rate to the quasilinear pre-diction and Fokker-Planck theory based on particle tra-jectories. Suitably weighted, all of these scattering ratesagree near the peak of the CR distribution where mul-tiple CR mean free paths fit within the ISM region.The agreement we find among scattering rates servesas a first-principles verification of CR hydrodynamics. Adiffusion coefficient computed from quasilinear theory isan accurate description of field-parallel transport due toself-generated wave-particle pitch angle scattering. Wenote that for the parameters studied in this work, non-linear effects are exaggerated, which enables crossing ofthe µ = 0 barrier at affordable spatial resolution.For more realistic environments with lower wave am-plitudes, broadening of gyroresonances would be re-duced, but other nonlinear wave-particle interactions,both resonant and non-resonant, may become impor-tant. Despite these uncertainties, the evidence from ourwork suggests that CR hydrodynamics is a valid model for CR transport in the ionized ISM, with wave-particlescattering naturally leading to fluid behavior of the col-lisionless distribution of CRs. Our work thus opens apathway toward first-principles calibration of CR fluidtransport coefficients in the multiphase ISM.ACKNOWLEDGMENTSWe are grateful to Illya Plotnikov for his contribu-tions to numerical tools for this project. CJB is thank-ful to Ellen Zweibel, Christoph Pfrommer, Matt Kunz,Anatoly Spitkovsky, Cole Holcomb, Bruce Draine, PhilHopkins, and Russell Kulsrud for advice and encour-agement. CJB is supported by the NSF Graduate Re-search Fellowship and the Churchill Foundation of theUnited States. XNB acknowledges support by NSFCgrant 11873033. The work of ECO was supported bygrant 510940 from the Simons Foundation. This workbegan at the Multiscale Phenomena in Plasma Astro-physics program at KITP in Santa Barbara, CA. Com-putational resources for our simulations were providedby the Princeton Institute for Computational Scienceand Engineering (PICSciE) and the Office of Informa-tion Technology’s High Performance Computing Centerat Princeton University.
Software:
Athena (Stone et al. 2008; Stone & Gar-diner 2009), yt (Turk et al. 2011)REFERENCES Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ,710, 133, doi: 10.1088/0004-637X/710/1/133Achterberg, A. 1981, A&A, 98, 161Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R.,et al. 2006, Nature, 439, 695, doi: 10.1038/nature04467Amato, E., & Blasi, P. 2009, MNRAS, 392, 1591,doi: 10.1111/j.1365-2966.2008.14200.xBai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015,ApJ, 809, 55, doi: 10.1088/0004-637X/809/1/55Bai, X.-N., Ostriker, E. C., Plotnikov, I., & Stone, J. M.2019, ApJ, 876, 60, doi: 10.3847/1538-4357/ab1648Besse, N., Elskens, Y., Escande, D. F., & Bertrand , P.2011, Plasma Physics and Controlled Fusion, 53, 025012,doi: 10.1088/0741-3335/53/2/025012Birdsall, C. K., & Langdon, A. B. 1985, Plasma PhysicsVia Computer Simulation (USA: McGraw-Hill, Inc.)Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y.2013, ApJL, 777, L16, doi: 10.1088/2041-8205/777/1/L16Boris, J. P. 1970, Proceedings of 4th Conference onNumerical Simulation of Plasmas, 3Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205 Buck, T., Pfrommer, C., Pakmor, R., Grand, R. J. J., &Springel, V. 2020, MNRAS, 497, 1712,doi: 10.1093/mnras/staa1960Bustard, C., Zweibel, E. G., D’Onghia, E., Gallagher, J. S.,I., & Farber, R. 2020, ApJ, 893, 29,doi: 10.3847/1538-4357/ab7fa3Cesarsky, C. J., & Volk, H. J. 1978, A&A, 70, 367Chan, T. K., Kereˇs, D., Hopkins, P. F., et al. 2019,MNRAS, 488, 3716, doi: 10.1093/mnras/stz1895Cohen, R. H., & Kulsrud, R. M. 1974, Physics of Fluids,17, 2215, doi: 10.1063/1.1694695Denton, R. E., & Kotschenreuther, M. 1995, Journal ofComputational Physics, 119, 283,doi: 10.1006/jcph.1995.1136Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215,doi: 10.1146/annurev.aa.28.090190.001243Dimits, A. M., & Lee, W. W. 1993, Journal ofComputational Physics, 107, 309,doi: 10.1006/jcph.1993.1146Draine, B. T. 2011, Physics of the Interstellar andIntergalactic Medium (Princeton University Press)
HD-PIC Simulations of Molecular Clouds Dupree, T. H. 1966, Physics of Fluids, 9, 1773,doi: 10.1063/1.1761932Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659,doi: 10.1086/166684Everett, J. E., & Zweibel, E. G. 2011, ApJ, 739, 60,doi: 10.1088/0004-637X/739/2/60Farber, R., Ruszkowski, M., Yang, H. Y. K., & Zweibel,E. G. 2018, ApJ, 856, 112,doi: 10.3847/1538-4357/aab26dFelice, G. M., & Kulsrud, R. M. 2001, ApJ, 553, 198,doi: 10.1086/320651Fujita, Y., Nobukawa, K. K., & Sano, H. 2020, arXive-prints, arXiv:2009.13524.https://arxiv.org/abs/2009.13524Gardiner, T. A., & Stone, J. M. 2005, Journal ofComputational Physics, 205, 509,doi: 10.1016/j.jcp.2004.11.016—. 2008, Journal of Computational Physics, 227, 4123,doi: 10.1016/j.jcp.2007.12.017Grenier, I. A., Black, J. H., & Strong, A. W. 2015, ARA&A,53, 199, doi: 10.1146/annurev-astro-082214-122457Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251,doi: 10.1111/j.1365-2966.2007.12692.xH. E. S. S. Collaboration, Abdalla, H., Abramowski, A.,et al. 2018, A&A, 612, A1,doi: 10.1051/0004-6361/201732098HESS Collaboration, Abramowski, A., Aharonian, F., et al.2016, Nature, 531, 476, doi: 10.1038/nature17147Holcomb, C., & Spitkovsky, A. 2019, ApJ, 882, 3,doi: 10.3847/1538-4357/ab328aHopkins, P. F., Chan, T. K., Squire, J., et al. 2020, arXive-prints, arXiv:2004.02897.https://arxiv.org/abs/2004.02897—. 2021a, MNRAS, 501, 3663,doi: 10.1093/mnras/staa3692Hopkins, P. F., Squire, J., Chan, T. K., et al. 2021b,MNRAS, 501, 4184, doi: 10.1093/mnras/staa3691Howes, G. G., Cowley, S. C., Dorland, W., et al. 2008,Journal of Geophysical Research (Space Physics), 113,A05103, doi: 10.1029/2007JA012665Hu, G., & Krommes, J. A. 1994, Physics of Plasmas, 1, 863,doi: 10.1063/1.870745Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007,ApJ, 671, 1736, doi: 10.1086/523036Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O., et al. 2018,ApJ, 855, 23, doi: 10.3847/1538-4357/aaadb9Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, MNRAS,496, 4221, doi: 10.1093/mnras/staa1849Jokipii, J. R. 1966, ApJ, 146, 480, doi: 10.1086/148912 —. 1971, Reviews of Geophysics and Space Physics, 9, 27,doi: 10.1029/RG009i001p00027Kempski, P., & Quataert, E. 2020, MNRAS, 493, 1801,doi: 10.1093/mnras/staa385Kennel, C. F., & Engelmann, F. 1966, Physics of Fluids, 9,2377, doi: 10.1063/1.1761629Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445,doi: 10.1086/149981Kulsrud, R. M. 2005, Plasma physics for astrophysics(Princeton University Press)Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014a,Physical Review Letters, 112, 205003,doi: 10.1103/PhysRevLett.112.205003Kunz, M. W., Stone, J. M., & Bai, X.-N. 2014b, Journal ofComputational Physics, 259, 154,doi: 10.1016/j.jcp.2013.11.035Kunz, M. W., Stone, J. M., & Quataert, E. 2016, PhRvL,117, 235101, doi: 10.1103/PhysRevLett.117.235101Lerche, I. 1967, ApJ, 147, 689, doi: 10.1086/149045McCall, B. J., Geballe, T. R., Hinkle, K. H., & Oka, T.1998, Science, 279, 1910,doi: 10.1126/science.279.5358.1910McKenzie, J. F., & Voelk, H. J. 1982, A&A, 116, 191Mertsch, P. 2020, Ap&SS, 365, 135,doi: 10.1007/s10509-020-03832-3Meyrand, R., Kanekar, A., Dorland, W., & Schekochihin,A. A. 2019, Proceedings of the National Academy ofScience, 116, 1185, doi: 10.1073/pnas.1813913116Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A,501, 619, doi: 10.1051/0004-6361/200911794Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, SSRv,216, 29, doi: 10.1007/s11214-020-00654-1Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V.2016, ApJL, 824, L30, doi: 10.3847/2041-8205/824/2/L30Parker, S. E., & Lee, W. W. 1993, Physics of Fluids B, 5,77, doi: 10.1063/1.860870Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., &Springel, V. 2017, MNRAS, 465, 4500,doi: 10.1093/mnras/stw2941Recchia, S. 2021, arXiv e-prints, arXiv:2101.02052.https://arxiv.org/abs/2101.02052Roe, P. L. 1981, Journal of Computational Physics, 43, 357,doi: 10.1016/0021-9991(81)90128-5Ruszkowski, M., Yang, H.-Y. K., & Reynolds, C. S. 2017,ApJ, 844, 13, doi: 10.3847/1538-4357/aa79f8Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312,doi: 10.1093/mnras/stt2121Schlickeiser, R. 1989, ApJ, 336, 243, doi: 10.1086/167009Schlickeiser, R., & Miller, J. A. 1998, ApJ, 492, 352,doi: 10.1086/305023 C. J. Bambic, X. Bai, & E. C. Ostriker
Semenov, V. A., Kravtsov, A. V., & Caprioli, D. 2020,arXiv e-prints, arXiv:2012.01427.https://arxiv.org/abs/2012.01427Shalchi, A. 2009, Nonlinear Cosmic Ray Diffusion Theories,Vol. 362, doi: 10.1007/978-3-642-00309-7Silsbee, K., & Ivlev, A. V. 2019, ApJ, 879, 14,doi: 10.3847/1538-4357/ab22b4Skilling, J., & Strong, A. W. 1976, A&A, 53, 253Spitkovsky, A. 2008, ApJL, 673, L39, doi: 10.1086/527374Spitzer, L. 1962, Physics of Fully Ionized GasesStone, J. M., & Gardiner, T. 2009, NewA, 14, 139,doi: 10.1016/j.newast.2008.06.003Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., &Simon, J. B. 2008, ApJS, 178, 137, doi: 10.1086/588755Summers, D., & Thorne, R. M. 1991, Physics of Fluids B,3, 1835, doi: 10.1063/1.859653Tibaldo, L., Digel, S. W., Casandjian, J. M., et al. 2015,ApJ, 807, 161, doi: 10.1088/0004-637X/807/2/161Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS,192, 9, doi: 10.1088/0067-0049/192/1/9 V¨olk, H. J. 1973, Ap&SS, 25, 471, doi: 10.1007/BF00649186Weinstock, J. 1969, Physics of Fluids, 12, 1045,doi: 10.1063/1.2163666Wentzel, D. G. 1971, ApJ, 163, 503, doi: 10.1086/150794Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467,906, doi: 10.1093/mnras/stx127Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens,A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152,doi: 10.1086/175510Yan, H., & Lazarian, A. 2004, ApJ, 614, 757,doi: 10.1086/423733Yang, R.-z., de O˜na Wilhelmi, E., & Aharonian, F. 2014,A&A, 566, A142, doi: 10.1051/0004-6361/201321044Zweibel, E. G. 2003, ApJ, 587, 625, doi: 10.1086/368256—. 2017, Physics of Plasmas, 24, 055402,doi: 10.1063/1.4984017
HD-PIC Simulations of Molecular Clouds L /10 v A c Figure 10.