Cosmic Ray Transport, Energy Loss, and Influence in the Multiphase Interstellar Medium
DDraft version December 15, 2020
Typeset using L A TEX twocolumn style in AASTeX62
Cosmic Ray Transport, Energy Loss, and Influence in the Multiphase Interstellar Medium
Chad Bustard
1, 2 and Ellen G. Zweibel
1, 3 Department of Physics, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706, USA Kavli Institute for Theoretical Physics, University of California - Santa Barbara, Kohn Hall, Santa Barbara, CA 93107, USA Department of Astronomy, University of Wisconsin - Madison, 475 North Charter Street, Madison, WI 53706, USA
Submitted to ApJABSTRACTThe bulk propagation speed of GeV-energy cosmic rays is limited by frequent scattering off hydro-magnetic waves. Most galaxy evolution simulations that account for this confinement assume the gasis fully ionized and cosmic rays are well-coupled to Alfv´en waves; however, multiphase density inho-mogeneities, frequently under-resolved in galaxy evolution simulations, induce cosmic ray collisionsand ionization-dependent or “plasma-based” transport driven by cosmic ray decoupling and elevatedstreaming speeds in partially neutral gas. How do cosmic rays navigate and influence such a medium,and can we constrain this transport with observations? In this paper, we simulate cosmic ray fronts im-pinging upon idealized, partially neutral clouds and lognormally-distributed clumps, with and withoutplasma-based transport. With these high-resolution simulations, we identify cloud interfaces as crucialregions where cosmic ray fronts can develop a stair-step pressure gradient sufficient to collisionlesslygenerate waves, overcome ion-neutral damping, and exert a force on the cloud. We find that the accel-eration of cold clouds is hindered by only a factor of a few when plasma-based transport is included,with additional dependencies on magnetic field strength and cloud dimensionality. We also probe howcosmic rays sample the background gas and quantify collisional losses. Hadronic gamma-ray emissionmaps are qualitatively different when plasma-based transport is included, but the overall luminosityvaries by only a small factor, as the short cosmic ray residence times in cold clouds are offset by thehigher densities that cosmic rays sample.
Keywords:
Cosmic-rays, Magnetic fields, Gamma-rays, Interstellar medium, Galaxy evolution INTRODUCTIONHighly energetic, non-thermal cosmic rays, despiterepresenting only a microscopic fraction of particles bynumber density ( n CR /n gas ≈ − ), are roughly in en-ergy equipartition with thermal and magnetic energy inthe Milky Way interstellar medium (ISM) (Boulares &Cox 1990). Because cosmic rays do not suffer from ra-diative losses (though hadronic collisions, Coulomb col-lisions, and “collisionless” scattering off magnetic per-turbations can be large energy sinks), and they don’tlose as much energy to adiabatic expansion comparedto thermal gas ( P ∝ ρ / for a relativistic gas instead of P ∝ ρ / for a non-relativistic gas), they fundamentally Corresponding author: Chad [email protected] alter ISM dynamics and are believed to be an importantcomponent of feedback, which regulates star formationin galaxies.For GeV cosmic rays, which make up the peak of thecosmic ray energy spectrum and therefore contain mostof the energy and momentum of the cosmic ray popu-lation, the dominant transport mode is self-confinementvia the streaming instability (Wentzel 1968; Kulsrud &Pearce 1969; Wentzel 1969). In this picture, cosmic rayswith a bulk drift speed greater than the Alfv´en speed canexcite Alfv´en waves through gyroresonance and pitchangle scatter off the waves until the cosmic rays obtainisotropy in the wave frame. Scattering transfers energyfrom the cosmic ray population to the waves, and sub-sequent wave damping heats the background gas. When a r X i v : . [ a s t r o - ph . H E ] D ec Bustard et al. the scattering mean free path is short , cosmic rays arewell-coupled to waves, and the bulk population advectsat the Alfv´en speed down the cosmic ray pressure gra-dient directed along the local magnetic field. Extrin-sic turbulence can alternatively confine cosmic rays, inwhich case cosmic rays scatter off the turbulent cascade.The resulting transport is field-aligned diffusion, but un-like the self-confinement model, there is no net transferof energy to the thermal gas (Zweibel 2017).To-date, most galaxy evolution simulations that in-clude a relativistic cosmic ray fluid assume some com-bination of diffusive and streaming transport in addi-tion to cosmic ray advection with the non-relativisticthermal gas. Given the canonical picture that 10%of supernova energy is converted to cosmic ray energythrough first-order Fermi acceleration, simulations withdiffusive or streaming transport included generally findcolder, smoother, more extended outflows than ther-mally driven winds, and in many cases, the non-thermalsupport is necessary to launch an outflow (e.g. Ipavich(1975); Breitschwerdt et al. (1991); Everett et al. (2008);Uhlig et al. (2012); Hanasz et al. (2013); Salem & Bryan(2014); Ruszkowski et al. (2017); Mao & Ostriker (2018);Buck et al. (2020); Dashyan & Dubois (2020); Bustardet al. (2020); Hopkins et al. (2020a)). Theory and ob-servations alike also suggest the existence of a long-lived cosmic ray reservoir in the circumgalactic medium(CGM), where they may further influence the gas cy-cle in and out of galaxies (Salem & Bryan 2014; Salemet al. 2016; Girichidis et al. 2018; Butsky & Quinn 2018;Kempski & Quataert 2019; Blasi & Amato 2019; Heintzet al. 2020; Ji et al. 2020), but the extent to whichcosmic rays navigate through the disk-halo interface,carry mass-loaded winds, and influence the CGM aretransport-dependent. Given the myriad possibilities forcosmic rays to narrow gaps between observations andsimulations of galaxy evolution, more work is needed toextend and apply first-principles fluid models of cosmicrays to galaxy evolution simulations and, likewise, toconstrain these models with mock observations.An important but mostly overlooked complication,worthy of more detailed study, is that the waves thatscatter cosmic rays are not always guaranteed to bepresent. In particular, density inhomogeneities, fre-quently under-resolved in simulations of galaxy evolu-tion, can induce alternating regions where cosmic raysare coupled or decoupled from waves. Let’s considera front of self-confined cosmic rays impinging upon acloud. As the cosmic rays stream down their pressure In fully ionized Milky Way ISM conditions, this is a goodapproximation, since the mean free path is ≈ . gradient and encounter a drop in Alfv´en speed in thecloud, their pressure build-up at the cloud interface ex-erts a force on the cloud (Wiener et al. 2017, 2019). Up-stream of this “bottleneck”, however, the steady-statesolution is a flat cosmic ray pressure profile unable toexcite hydromagnetic waves and unable to exert a forcealong the direction of the magnetic field (Skilling 1971).Cosmic rays similarly become decoupled when the con-fining waves are heavily damped, forcing steady-statecosmic ray streaming speeds to be much higher than thegas Alfv´en speed in order for the wave excitation rate( ∝ v st ) to be comparable to the damping rate. In par-tially neutral clouds, where damping due to ion-neutralfriction can be strong, this effect is believed to be espe-cially significant (Kulsrud & Pearce 1969).As it’s been shown that cosmic ray bottlenecks cantheoretically accelerate warm, fully ionized clouds ingalaxy halos up to 100s of km/s, it’s of great interest tounderstand the conditions in the ISM rather than CGMwhere cosmic rays are similarly well-coupled and capableof doing work through their pressure gradients. Giventhe added complication that cloud interiors may be par-tially neutral, in which decoupling or fast, “plasma-based transport” serves to flatten cosmic ray pressuregradients, this problem is not trivial and may have im-portant consequences for feedback.Farber et al. (2018) ran magnetohydrodynamic(MHD) ISM “patch” simulations of cosmic ray drivenwinds with a locally adaptive diffusion coefficientboosted in cold, predominantly neutral gas to mockup the effects of ion-neutral damping. This leads todifferent gas properties, a broader spatial distributionof cosmic rays, and higher wind speeds as cosmic rayspreferentially push on the hot gas and escape more eas-ily into the halo. These pioneering “two- κ ” models, inwhich the diffusion coefficient shifts between a low andhigh value depending on gas temperature, underscorethe possible implications of decoupling, but the transi-tion from fully ionized to partially neutral transport ishighly nonlinear with further dependencies on magneticfield strength and cosmic ray pressure gradient (Everett& Zweibel 2011). An idealized study focused on teasingout some of the complicated interplay between bottle-necks, variable transport, and collisional losses in themultiphase ISM is still needed.This study is further motivated by recent full-galaxysimulations from the FIRE collaboration (Chan et al.2019; Hopkins et al. 2020b) that attempt to constraincosmic ray transport with gamma-ray emission andgrammage estimates. Intriguingly, the best-fit simula-tions necessitate very fast cosmic ray transport speeds,especially in diffuse gas, to lower gamma-ray emission ow Cosmic Rays Navigate the Multiphase ISM §
2, we providemore background on the nonlinear cosmic ray transportand energy loss mechanisms we study in this work. Wepresent our assumptions and analytic expectations forcosmic ray transport in § §
4. In §
5, we present and compare 1D and 2Dsimulations of cosmic ray fronts impinging upon rep-resentative ISM clouds. In §
6, we present 2D and 3Dsimulations of cosmic ray transport through a layer oflognormally distributed multiphase clumps in a mockISM. For each “obstacle course”, we quantify cosmicray energy loss (both collisional and collisionless), themomentum imparted to the gas, and probability dis-tribution functions of cosmic ray energy and collisionalenergy loss (with and without taking into account neu-tral gas ). In §
7, we discuss the results and limitationsof our work, and we conclude in § §
6. Readers who are interested in cos-mic ray influence on galaxy evolution should also read § BACKGROUND2.1.
Cosmic Ray Coupling and Decoupling
To first order, cosmic rays comprise a second, rela-tivistic fluid that advects with the nonrelativistic ther-mal gas. We additionally know from the near isotropy of observed cosmic rays and their long confinement timesin the galaxy, estimated from spallation measurements,that cosmic rays are constantly being scattered (see re-cent reviews by e.g. Zweibel (2017); Amato & Blasi(2018); Becker Tjus & Merten (2020)). Their result-ing transport is well-described by a random walk alongtangled magnetic field lines with an energy-dependentdiffusion coefficient. For cosmic rays at energies above afew hundred GeV, scattering off waves generated by anexternal turbulent cascade likely dominates the trans-port. Cosmic rays are coupled to the waves throughresonant interactions, and transport is akin to magneticfield aligned diffusion with a coefficient depending onthe amplitude of scattering waves. Unless the left andright propagating Alfv´en waves in the turbulent cascadehave unequal intensities, though, there is no net transferof energy between cosmic rays and the thermal gas.For GeV-energy cosmic rays, however, the confin-ing waves can be generated by the cosmic rays them-selves. Cosmic rays with even a tiny amount of driftanisotropy can resonate with magnetic perturbationsand become self-confined: they exchange energy withAlfv´en waves through gyroresonance, pitch angle scat-ter off those waves, and become locked to the waveframe, therefore advecting at the local ion Alfv´en speed v ionA = B/ √ πρ ion (in addition to advecting with thebackground gas flow). In a steady-state, the conver-sion of cosmic ray energy to wave energy is balancedby wave damping, which heats the background gas at arate ∝ v ionA · ∇ P CR . We refer to this energy transfer ascollisionless because it is mediated by magnetic fields,not direct collisions between particles.In a multiphase ISM punctuated by density irregu-larities, however, cosmic ray coupling can quickly breakdown. For instance, consider a decrease in Alfv´en speedalong a magnetic field line. One can show that, insteady-state, streaming cosmic rays must conform to aconstant cosmic ray pressure profile for an extended re-gion upstream (Skilling 1971). Under such conditions,no waves are excited , and cosmic rays no longer trans-fer energy or momentum to the thermal gas. We will re-fer to these regions as cosmic ray “bottlenecks” (Wieneret al. 2017; Zweibel 2017; Wiener et al. 2019), a subset Cosmic ray pressure anisotropy may provide an alternativeconfinement mechanism by exciting magnetic waves in the up-stream region, where the cosmic ray streaming instability outlinedhere does not act. This could lock cosmic rays back to the ther-mal gas (instead of allowing them to decouple), but appears not totransfer much momentum or energy to the gas (Zweibel 2020). Forthis work, we will only consider cosmic ray confinement throughdrift anisotropy and save an analysis of pressure anisotropy forfuture work.
Bustard et al. of a larger class of decoupled regions that Skilling (1971)referred to as “free-zones.”This bottleneck effect has been studied with 1D and2D numerical simulations under simplified conditions,namely a cosmic ray front impinging on either a 1Dslab or 2D cylindrical cloud. While these simulationspromisingly show that cosmic ray bottlenecks may ac-tually accelerate these cold clouds (Wiener et al. 2017,2019; Br¨uggen & Scannapieco 2020), thereby helping en-train them into hot outflows, they most importantly ex-pose the complex relationship between cosmic rays andthe multiphase ISM. Although cosmic ray transport isdetermined by local conditions on scales of a cosmic raygyroradius, the global environment determines these lo-cal conditions. In the bottleneck scenario, the presenceof a cold cloud creates a traffic jam, decoupling cosmicrays from thermal gas, extending for macroscopic dis-tances upstream of the cloud. One can then imagine acheckerboard of these clouds, each triggering cosmic raybottlenecks. How a steady source of cosmic rays nav-igates this series of traffic jams depends on the globalcloud and magnetic field structure.Another example of “free-zones” can occur in colder,partially neutral gas, where ion-neutral collisions decou-ple ions and neutrals at scales above the cosmic ray gy-roradius and thereby damp the Alfv´en waves that playthe crucial scattering role (Kulsrud & Pearce 1969). In asteady-state, the resonant streaming instability growthrate, which is proportional to the cosmic ray drift veloc-ity, must balance this increased damping rate, effectivelyincreasing the streaming velocity (sometimes by ordersof magnitude). If ion-neutral damping is strong enough,the cosmic ray mean free path may in fact be larger thanthe partially neutral gas region, leading cosmic rays tofree-stream. Even if the cosmic rays can drive wavesat a rate fast enough to couple them to the gas, cos-mic rays only scatter off Alfv´en waves that propagatein the ions, meaning the streaming velocity is rigorouslythe ion
Alfv´en speed, which can be orders of magnitudelarger than the gas Alfv´en speed when ρ ion << ρ . Collisional Losses
In addition to instigating complex, nonlinear trans-port, cold clouds are also targets for cosmic ray hadronicand Coulomb collisions, the former leading to a catas-trophic decay of pions into gamma-rays. These colli- For the remainder of this paper, we will distinguish betweenthese two important effects, reserving the fequently-used term“fast transport” specifically for elevated diffusive transport aris-ing from cosmic ray decoupling. We will refer to the combinationof decoupling and elevated ion Alfv´en speeds when ion fraction f ion < sional losses, depending on environment, can representlarge energy sinks that suppress cosmic ray influencein galaxies (Hopkins et al. 2020c; Crocker et al. 2020).This is most true for starburst galaxies, which have highvolume-filling factors of dense, neutral gas and by all in-ferences from their gamma-ray emission appear to begood proton calorimeters (Lacki et al. 2011; Tang et al.2014; Yoast-Hull et al. 2015, 2016; Wang & Fields 2018;Krumholz et al. 2020), meaning that all cosmic raysproduced from supernovae are consumed by collisions.Observations of L (cid:63) and dwarf galaxies, however, painta different picture where the vast majority of protonsescape hadronic collisions presumably due to a combi-nation of diffusion and advection in supernova-drivenwinds (Lacki et al. 2011; Ackermann et al. 2016; Fu et al.2017; Lopez et al. 2018). While we only have limitedgamma-ray observations of ≈
10 external star-forminggalaxies (e.g. Ajello et al. (2020)), and while it is diffi-cult to exactly separate the diffuse hadronic gamma-rayemission from point-source emission, the resulting scal-ing relations between gamma-ray luminosity and starformation rate (SFR) are already in severe disagreementwith the current iteration of MHD + cosmic ray dwarfand L (cid:63) galaxy simulations – that is, simulations thatinclude cosmic ray diffusive or streaming transport butassume that transport is not affected by ionization.Chan et al. (2019); Hopkins et al. (2020b) vary cos-mic ray transport models in their state-of-the-art full-galaxy simulations and find that gamma-rays will beoverproduced compared to Local Group dwarf galaxyobservations unless the cosmic ray diffusion coefficientis large ( > cm /s) compared to the canonical valueof a few × cm /s inferred from cosmic ray prop-agation models. Testing a variety of proposed first-principles transport models, Hopkins et al. (2020b) con-strain the issue further: propagation appears to be rate-limited by the warm ionized medium (WIM) and innerCGM, where the current paradigm of self-confinementmodels predicts relatively slow transport speeds. Fasttransport in partially neutral gas is not sufficient to de-crease gamma-ray emission, which is instead dominatedby small patches of the WIM where the authors believea runaway increase in self-confinement could be occur-ring. That is, large cosmic ray pressure gradients couldgenerate confining waves that further trap cosmic rays,increase the cosmic ray pressure gradient, and so on.Clearly, the interplay between cosmic ray transport,wave damping, and collisional losses in the multiphaseISM is very rich and has implications for gamma-rayconstraints as well as general ISM and galactic wind dy-namics. The goal of this paper is to shine a light on thisinterplay with idealized but high resolution simulations. ow Cosmic Rays Navigate the Multiphase ISM ANALYTIC EXPECTATIONSWithin the self-confinement picture of cosmic raytransport, the streaming velocity is v st = B · ( ∇ · P CR ) | B · ( ∇ · P CR ) | v ionA (1)where v ionA = B/ √ πρ i = B/ √ πρf ion , where f ion isthe ion fraction by number, and ρ i = m i n i is the iondensity. Here, and in the following shortened derivation(see also Jiang & Oh (2018); Hopkins et al. (2020b)),we’ve been careful to write the Alfv´en speed as the ionAlfv´en speed, v ionA to reflect that cosmic rays resonatewith Alfv´en waves that propagate in the ions. Notethat although Equation (1) is written in terms of thecosmic ray pressure tensor P CR , we will assume isotropicpressure and replace the divergence by a gradient. Asseen from Equation 1, self-confinement only occurs inthe presence of a cosmic ray pressure gradient directedalong the magnetic field. Equation 1 holds, as written,when there is no wave damping, but more generally, thesteady-state velocity is obtained by equating the growthrate of cosmic ray - excited Alfv´en waves Γ CR with thewave damping rate Γ:Γ CR ( γ ) ≈ π α − α − n CR ( > γ ) n i (cid:18) v D v ionA − (cid:19) = Γ (2) n CR ( > γ ) is the number density of cosmic rays withLorentz factor greater than γ , and we’ll assume α = 4 isthe power-law exponent of the cosmic ray distributionfunction in momentum space.Because the growth rate of the streaming instabilityscales linearly with the bulk drift speed of the cosmicrays, a higher drift speed is needed when wave dampingis present. The net drift with respect to the Alfv´en waveframe is v D − v ionA , and we write this as a diffusive flux(Jiang & Oh 2018): F diffuse = κ || ∇ f ≈ ( v D − v ionA ) f (3)Then the (purely parallel) diffusivity is κ || = f ∇ f π α − α − n i v ionA Ω n CR ( > γ ) (4)Following Hopkins et al. (2020b), we write f / ∇ f = l CR , a cosmic ray scale length, and we can formulatethis diffusion coefficient as an effective boost to thestreaming velocity, v effst = v ionA + κ || / ( γ CR l CR ) where γ CR = 4 /
3. Note that in our simulations, we implementthe diffusive flux κ || and model streaming with speed v st = v ionA , not v st = v effst . We then make some sub-stitutions ( e B = 1 / v ionA ) n i m p , e CR = γ L m p n CR c ) to re-write Equation 4 in terms of more typical codequantities to obtain: κ || = 4 cr L Γ e B l CR πv ionA e CR (5) v effst = v ionA + 3 cr L Γ e B πv ionA e CR (6)where r L = c/ Ω is the cosmic ray Larmor radius, andwe assume cosmic rays occupy a single energy bin at 1GeV.The damping rate Γ can have contributions from mul-tiple processes, including nonlinear Landau damping,which we find to have no bearing on our results (see § in = ν in ≈ − f neutral T / ρ − s − (7)where ν in is the collision frequency between ions andneutrals and f neutral = 1 − f ion is the neutral fraction.Obtaining the ion fraction f ion is quite complicatedas ionization depends on the local radiation field andthe population of low energy (primarily 2 - 10 MeV)cosmic rays that may not follow the same transportas the GeV cosmic rays we model here. Including amore complete treatment of ionization will be the sub-ject of future work, but for now, we remain agnosticabout the exact ion fraction for clouds of varying den-sities and temperatures. We prescribe an ion fractionthat, near 10 K, matches the tabulated Sutherland &Dopita (1993) ionization fraction fairly well. Note that,while the Sutherland-Dopita table assumes collisionalionization equilibrium and only extends down to 10 K,we extend the neutral fraction to lower temperatures fol-lowing a hyperbolic tangent function. We then imposea floor on the ion fraction, f minion , ranging from 10 − to 10 − to parameterize our ignorance and explore thesensitivity of our results to ionization level, althoughnot reaching the f ion characteristic of Galactic molec-ular clouds, which can be as low as f ion ∼ − . Forthe relatively diffuse clouds we focus on, with columndensities ∼ − cm − , it’s likely that there aresufficient UV photons to keep carbon ionized, and ionfractions lower than 10 − are unlikely (Hollenbach &Tielens 1999; Neufeld & Wolfire 2017; Silsbee & Ivlev2019). The equation we use is f ion = 1 − (cid:18)
12 (1 − f minion ) (cid:18) tanh (cid:18) a − Tc (cid:19)(cid:19)(cid:19) (8) Bustard et al. Temperature (K) f i o n Ionization Fraction Function
Sutherland-Dopita 1993 f minion = 10 f minion = 10 f minion = 10 Temperature (K)10 E ff e c t i v e S t r e a m i n g V e l o c i t y ( c m / s ) Solid line: v st ( e B / e CR = 0.1)Dashed line: v st ( e B / e CR = 10)Dotted line: v ionA Solid line: v st ( e B / e CR = 0.1)Dashed line: v st ( e B / e CR = 10)Dotted line: v ionA Solid line: v st ( e B / e CR = 0.1)Dashed line: v st ( e B / e CR = 10)Dotted line: v ionA Solid line: v st ( e B / e CR = 0.1)Dashed line: v st ( e B / e CR = 10)Dotted line: v ionA Pressure = 3.2e-12 dyne/ cm f minion = 10 f minion = 10 g / cm ) Figure 1.
Top panel: Prescribed ion fraction as a func-tion of temperature for three minimum ion values, f minion =10 − , − and 10 − . Bottom panel: Ion Alfv´en velocities( v ionA ) and effective streaming velocities ( v effst from Equa-tion 6) for two different f minion values and varying ratios ofmagnetic to cosmic ray energy density. For high cosmic rayenergy density, v st ≈ v ionA when the ion fraction is low, butthese velocities can diverge significantly when the ion frac-tion is higher or cosmic ray energy density is lower. Thisis especially true near 10 K, where ion-neutral damping in-duces a bump in v effst . where T is the temperature in Kelvin, a = 1 . × K,and c = 2 × K.Figure 1 shows the results of Equation 6 for ionizationfraction functions with various f minion . For a large cos-mic ray energy e CR relative to magnetic energy e B , thestreaming velocity tracks pretty closely to the ion Alfv´envelocity as the last term on the RHS of Equation 6 issmall. We also see, however, that the ion Alfv´en velocityand streaming velocity can diverge considerably whenthe minimum ion fraction is higher or when the cosmicray energy density is lower. Clearly the streaming veloc- ity is sensitive to ion fraction function, and this is espe-cially true at cloud interface temperatures near 10 K.Exploring this sensitivity with more realistic treatmentsof ionization will be the subject of future work, but aswe’ll see in our simulations of cosmic ray fronts, overallcosmic ray penetration into clouds is largely insensitiveto whether we include or exclude ion-neutral damping(the RHS of equation 6). This is because the dynam-ics are set by the cosmic ray energy density outside thecloud, where e CR is by-design comparable to the ther-mal and magnetic energy densities in our simulations;this is the case in the average Milky Way ISM and nearcosmic ray sources, and, rather than situations where e CR is small and therefore won’t significantly affect theISM, this is our regime of interest.Of course, when there is no cosmic ray pressure gradi-ent, cosmic rays are completely decoupled from wavesand free-stream, but when a pressure gradient doesexist, cosmic rays generate self-confining waves. Weshow in our simulations that, in our regime of inter-est when e CR is not much less than e B , a steep pressuregradient can develop when the streaming velocity de-creases in a small region at the cloud interface (compare v st ( T = 3 × K) to v st ( T = 10 K), and subsequenttransport follows v ionA fairly close. The necessary pres-sure gradient to overcome ion-neutral damping can evenbe created from an initially uniform cosmic ray pressurewhen collisional losses are strong (see Appendix), butthese conditions are not as easily met when the cosmicray energy density is lower.This is all consistent with the work of Everett &Zweibel (2011), who study the steady-state interactionbetween a population of GeV energy cosmic ray protons,with non-zero pressure gradient, and cold clouds of vary-ing densities and magnetizations. They find that steeperpressure gradients and higher magnetic field strengthsallow cosmic rays to remain coupled to waves deeper intothe cloud. For a shallow pressure gradient, though, thecosmic rays eventually decouple from waves and free-stream through the cloud, leading Everett & Zweibel(2011) to hypothesize that the cosmic ray energy insideand outside clouds will be roughly the same. Whetherthe diffusive flux can compensate for the drop in ad-vective flux, leading to this scenario, depends on thedynamical and thermodynamical effects of the cosmicray pressure gradient, which is accounted for in our sim-ulations but not in Everett & Zweibel (2011), in whichthe cloud properties are held fixed. COMPUTATIONAL METHODSUntil recently, a large suite of simulations probing theinteraction of cosmic rays with cold, partially neutral ow Cosmic Rays Navigate the Multiphase ISM t ∝ (∆ x ) ). For fast cosmic ray transport (hence,fast signal speed) in partially neutral gas, this conditionmakes high-resolution simulations imprudent.We can overcome this limitation using a new cosmicray treatment based on a two-moment method previ-ously used for radiative transfer (Jiang & Oh 2018)which has been implemented in Athena++ (Stone et al.2020). Unlike the regularization method, this methodhas no dependence on smoothing parameter near cosmicray maxima and boasts a linear stable timestep scal-ing (∆ t ∝ ∆ x ) that makes it computationally tractableto simulate very fast transport speeds and resolve thesmall-scale structures that induce them.We refer the reader to Jiang & Oh (2018) for more de-tails (see also a similar method by Thomas & Pfrommer2019) and use this space just to outline our additions ofplasma-based transport. Namely, we implement a flag,which when turned on, calculates the streaming veloc-ity as the ion Alfv´en velocity everywhere (including forthe v ionA · ∇ P CR heating term), given the temperature-dependent ion fraction function of Equation 8. Thisflag also includes the additional boost due to ion-neutralfriction as a diffusive flux term with diffusion coefficientgiven by Equation 5.We also include sink terms − Λ C and − Λ H , corre-sponding to Coulomb and hadronic collisions, respec-tively, in the cosmic ray energy equation. Corresponding+Λ C and +Λ H / e ± pairs, the majority of which heat the gasthrough Coulomb collisions (Pfrommer et al. 2017). Therest of the hadronic energy loss escapes as gamma-raysand neutrinos. We use the equations of Enßlin et al.(2007); Pfrommer et al. (2017) for the Coulomb andhadronic loss terms:Λ coul = − . × − (cid:16) n e cm − (cid:17) (cid:18) e c ergcm − (cid:19) ergcm − s − (9)Λ hadr = − . × − (cid:16) n cm − (cid:17) (cid:18) e c ergcm − (cid:19) ergcm − s − (10) where n is the number density of nucleons, and n e is thenumber density of free electrons. Note that, while the gas evolution equaations containcollisional and collisionless heating terms, we do not in-clude radiative cooling. This is a practical choice thatallows us to prescribe initially static gas distributionsand isolate the already complex interplay with cosmicrays. We discuss the implications of this further in § MHD SIMULATIONS: SINGLE CLOUDWe ran a large number of 1D and 2D simulations ofdifferent cloud sizes and densities embedded in variousenvironments. We will focus here on the interaction be-tween a cosmic ray front and one representative cloudwith maximum density n cold = 10 cm − in a backgrounddensity of n hot = 0.1 cm − . The initial density profiletakes the same form as in Wiener et al. (2017, 2019): ρ ( r ) = ρ hot + 12 ( ρ cold − ρ hot ) (cid:18) − tanh (cid:18) r − r c t c (cid:19)(cid:19) (11)The fiducial cloud radius, r c , is 10 pc and tapers offwith an interface that is fiducially t c = 5 pc thick. Aswe’ll see, the interface thickness and how many cellsresolve the interface are key parameters. The initialsetup is a 2 kpc long box threaded by a 5 µG magneticfield (plasma β ≈
3) and at constant thermal pressure3 . × − dyne cm − , giving a background temper-ature of 10 K. We do not include radiative cooling,so such an unstable temperature poses no issue. At atemperature of 10 K, the background medium is safelyfully ionized, which isolates the onset of super-Alfv´enicstreaming to the cloud region. The cosmic ray energydensity is set to an initially negligible value.Cosmic rays are injected at the left boundary by grad-ually increasing the cosmic ray energy flux and thengradually shutting it off after 30 Myrs. F CR = F CR (1 − e − t/ )(1 − e ( t − / ) (12)where our fiducial value of F CR is 1 . × − erg cm − s − and t is in Myrs. For pure streaming transport,this yields a cosmic ray front with pressure of order the During the write-up of this project, we found an error in ourcode where we calculated the Coulomb loss rate using the totalgas density instead of the density modified by ionization fraction.Coulomb losses, which should become negligible in cold, neutralgas, are then overestimated in our simulations. Since hadroniccollisions dominate everywhere, this is a small correction to thetotal collisional loss rate and one that is, in reality, partially offsetby our neglect of ionization losses that act primarily in dense gas.
Bustard et al. thermal and magnetic pressures, which in 2D and 3D issignificant enough to warp magnetic field lines in the lat-eral direction when cosmic ray bottlenecks form at coldcloud interfaces. This is exaggerated further in simu-lations with lower magnetic field strength (B = 1 µG instead of the fiducial 5 µG ) or with larger cosmic rayinfluxes, as we show in § ≈ . ≈
37 km/s, and as found in (Wiener et al.2017, 2019), the sound wave that outpaces the cosmicray front distorts and begins to push the cloud before thecosmic rays arrive. We also ran some tests with a back-ground temperature of 3 × K, thereby decreasing thesound speed below the Alfv´en speed, and we found ourmain conclusions about cosmic ray - cloud interactionsto be robust. Of course, decreasing the temperaturewhile keeping the same cloud density means the cloudis now less pressurized, so there are some differences inthe resulting cloud morphology.To properly model streaming with the two-momentmethod, there is a maximum speed of light parameter, V m , that needs to be much greater than the maximumpropagation speed. For simulations with v st = v A (as-suming a fully ionized gas), we set V m = 10 km/s,which is more than sufficient since the gas Alfv´en speedis safely lower than 10 km/s in all simulations. For sim-ulations including plasma-based transport in partiallyneutral gas, we increase V m to 10 km/s. In this case, v ionA approaches 10 km/s in our fiducial cloud setup andlocally exceeds that in the multi-cloud simulations of § V m remains safely higher than the fastestpropagation speed, we also cap the streaming speed at10 km/s and κ || at 3 × cm s − . We tested thevalidity of these caps in 1D simulations by increasingthem each by a factor of 10 (and therefore increasing V m by a factor of 10, as well), and we found negligibledifferences.A separate issue is whether cosmic rays with high ef-fective streaming speeds can be treated as a fluid. Wewill revisit this in the Appendix, but for now, we takea detailed look at cosmic ray - cloud interactions withinthe fluid assumption.5.1.
1D Simulations
We begin with 1D simulations, which therefore forcethe cosmic rays to eventually enter the cloud. We placethe fiducial n cloud = 10 cm − cloud at the center ofthe box a distance x = 1 kpc from the left boundary.We also ran simulations with the cloud placed closer tothe source and find qualitatively similar results. If not explicitly stated, the resolution of all 1D simulations is0.5 pc.Figure 2 shows a time-series of the cosmic ray frontapproaching and entering the cloud for a simulation with f minion set to 10 − . The top panel shows the gas densityand cosmic ray energy density, while the bottom panelshows v ionA and κ || . Although the cosmic ray pressurewithin the cloud is initially constant, spatially depen-dent collisional losses sustain a small pressure gradientdirected towards the cloud center, i.e. l CR >
0. Sincethe cosmic ray energy density within the cloud is ordersof magnitude less than the thermal and magnetic en-ergy densities, the cosmic ray pressure gradient is stillvery small, thereby giving a huge effective diffusivity(see Equation 5).Cosmic ray propagation quickly changes when the cos-mic ray front approaches the cloud. The bottom panelof Figure 2 shows that, in the cloud interface, before thetemperature decreases below ≈ K, v ionA decreasesconsiderably. Here, the gas is still mainly ionized, andthe increase in density outpaces the decrease in ion frac-tion. Once the temperature decreases further, the neu-tral fraction and, hence, v ionA increase significantly, butnot before this decrease in propagation speed has starteda cosmic ray traffic-jam. This bottleneck amplifies thepre-existing cosmic ray pressure gradient (leading tosmall l CR ) at the cloud interface, which suppresses thediffusive flux due to ion-neutral damping. The resultingcosmic ray propagation is then dominated by advectionat the ion Alfv´en speed, as cosmic rays are locked to thewave frame. This is seen in the bottom panel of Figure2, where the diffusivity drops from the capped value of3 × cm s − down to ≈ − cm s − .We check the effects of ion-neutral damping further byvarying f minion between 10 − and 1.0 in Figure 3. Higher f minion values correspond to lower advective fluxes andhigher diffusive fluxes, as expected from § f minion . This appears true in lo-calized regions, but the advective flux still dominantlydetermines the transport. Compared to f minion = 10 − and 10 − , which have almost identical density and cos-mic ray energy profiles, for f minion = 10 − , the diffusivityis boosted, especially at the cloud interface region near T = 10 K. This partially suppresses bottleneck forma-tion at the leading cloud edge (the cloud morphology nolonger shows a pronounced bump) and drives a flat cos-mic ray pressure profile in isolated regions near the cloudedges. Cosmic ray pressure gradients still drive wavesthroughout the cloud interior, though, locking cosmicrays to waves and suppressing the diffusivity. Because v ionA is now lower, there is an overall slower cosmic ray ow Cosmic Rays Navigate the Multiphase ISM D e n s i t y ( g c m ) t = 20 Myr 0 pc 50 pc 100 pct = 28 Myr 0 pc 50 pc 100 pct = 30 Myr 50 pc 100 pc 150 pct = 40 Myr E C R ( e r g c m ) v st = v ionA + Ion-Neutral Damping, f minion = 10 D i ff u s i v i t y ( c m s ) t = 20 Myr 0 pc 50 pc 100 pct = 28 Myr 0 pc 50 pc 100 pct = 30 Myr 50 pc 100 pc 150 pct = 40 Myr v i o n A ( c m s ) Figure 2.
Time series of a cosmic ray front impacting an initially stationary, partially neutral cloud of radius 10 pc andinterface width 5 pc. The blue lines show the density, while the black lines show the cosmic ray energy density. The minimumion fraction we consider is f minion = 10 − , resulting in an increase to the ion Alfv´en speed and time-varying diffusivity within thecloud (both shown in the bottom panel). The cloud initially is at 1.0 kpc away from the left boundary, which corresponds to 0pc in these plot coordinates, but the combination of cosmic ray pressure and preceding acoustic wave pushes the cloud to theright over time. Note that the right-most panel is shifted to follow the cloud as it moves ≈
100 pc over 40 Myrs. Even withplasma-based transport included, cosmic rays bottleneck at the leading edge, where they encounter a dip in transport speed.The ensuing pressure gradient drives waves and suppresses the diffusivity due to ion-neutral damping. The resulting cosmic raypressure gradient is steep at the front and back cloud edges.
50 pc 100 pc 150 pc D e n s i t y ( g c m ) f minion = 10
50 pc 100 pc 150 pc f minion = 10
50 pc 100 pc 150 pc f minion = 10
50 pc 100 pc 150 pc f minion = 1.0 10 E C R ( e r g c m ) Figure 3.
Density (blue lines) and cosmic ray energy density (black lines) for the fiducial cosmic ray front impinging on acloud of radius 10 pc and interface width 5 pc (solid lines) and interface width 0.5 pc (dashed lines). The panels correspond todifferent f minion = 10 − , − , − , and 1 . v st = v A ) going from left to right. Bustard et al.
50 pc 100 pc 150 pc C o lli s i o n a l ( e r g c m s ) f minion = 10
50 pc 100 pc 150 pc D e n s i t y ( g c m )
50 pc 100 pc 150 pc f minion = 1.0 C o lli s i o n l e ss ( e r g c m s )
50 pc 100 pc 150 pc T e m p e r a t u r e ( K ) Figure 4.
Fiducial cosmic ray front impinging on a cloudof radius 10 pc and interface width 5 pc with f minion = 10 − (left) and v st = v A (right). The top row shows energy losses(both collisional and collisionless). The bottom row showsdensity and temperature. passage through the cloud and a cosmic ray profile in-termediate between f minion = 10 − and f minion = 0 (fullyionized). Clearly the nonlinear transport is especiallycomplex in this case, but the cosmic ray energy down-stream of the cloud is surprisingly very similar to thatwith lower f minion values.We also varied the cloud interface width from 5 pc to0.5 pc, which is only one cell thick and shown by thedashed lines in Figure 3. We found that larger widthsinduced more severe bottlenecks, and coupled with theincreased gas mass and, hence, increased collisional en-ergy loss, the cosmic ray energy leaving the back side ofthe cloud decreased. Thinner interfaces do increase theimportance of the diffusive flux, since a quick transitionbetween fully and partially ionized gas doesn’t allow forsuch a build-up of cosmic ray pressure and confiningwaves; however, even with a thin interface resolved byonly one cell, we see a clear drop in cosmic ray pressureat each cloud interface, and the advective flux domi-nates the transport. In light of the limited resolutionavailable in most galaxy-scale simulations, it’s maybemost instructive to show the same cosmic-ray-cloud in-teraction at varying resolution, which exhibits the sametrends as increasing or decreasing the interface width(see Appendix). For now, it is just worth noting thatwe don’t reach convergence until we use a resolution of0.5 pc, our fiducial resolution for all 1D simulations.We ran the same simulation assuming the mediumis fully ionized and compare results in Figure 4. Thetop row shows both collisional and collisionless energyloss. As expected, and as seen in Figure 2, the steeppressure gradient that forms at the front of the cloudcollisionlessly transfers cosmic ray energy to waves at a rate ∝ v ionA · ∇ P CR . Wave damping then heats thecloud edge, as seen in the bottom panel of Figure 4, andbroadens the interface, resulting in a thumb-like densitycompression that moves upstream.In the fully ionized case, this heating tracks the cosmicray front, which at this timestamp, has slowly movedto the middle of the cloud. In the partially neutralcase, fast transport speeds quickly flattens the cosmicray pressure in the cloud, so heating is negligible inthe cloud center. When cosmic rays hit the back ofthe cloud, though, they encounter a drop in streamingspeed as the medium transitions from neutral to ionized(high v ionA to low v ionA ). This second bottleneck persistsfor a long time while cosmic ray energy slowly increasesdownstream of the cloud, and the steep pressure gradi-ent there induces collisionless energy loss, just as it doeson the front cloud edge. This heating may have inter-esting implications for ion abundances and kinematics(Wiener et al. 2017); however, we note that the energyloss rate is only of order 10 − erg cm − s − , which fora density of 1 −
10 cm − , is far less than the expectedcooling rate of order 10 − or higher at temperatures ofa few x 10 K. In future work, we will include radiativecooling to address this directly and determine realisticinterface widths, which are set by a balance betweencooling, conduction, and cosmic ray heating. Whethercosmic rays also modify their environment by increas-ing the ionization fraction is an open question similarlybeyond the scope of this paper.5.2.
2D Simulations
To explore the effects of higher dimensions and buildintuition for our main suite of simulations in §
6, we studythe same cosmic ray front interacting with a 2D cloud.In the case we will show, the cloud density, radius, andinterface width are the same as before in the 1D simula-tions, but we move the cloud to a closer distance of 100pc from the boundary. This is more akin to a cosmicray front impacting dense clouds within an ISM scaleheight of a few hundred pc. To make a fair comparisonbetween 1D and 2D, we re-ran our 1D simulations withthe cloud at a distance of 100 pc instead of 1 kpc. The2D simulation resolution is ≈ µ G (bot-tom row) and 5 µ G (top row). The cloud propertiesare the same as in § ow Cosmic Rays Navigate the Multiphase ISM B = 5 G f minion = 1.0
10 pc B = 5 G f minion = 10 B = 1 G B = 1 G 10 E C R ( e r g c m ) Figure 5.
Snapshots of cosmic ray energy density for asingle-cloud setup in 2D. The cloud density initially peaksat 10 cm − , and white, yellow contours represent densitiesof 1, 10 cm − . Varying the magnetic field strength from 5 µ G(top row) to 1 µ G (bottom row), we see distinct differencesin magnetic field topology created by the pressure wave thatprecedes the cosmic ray front. In the fully ionized case (left),cosmic ray propagation traffic-jams at the front cloud edge,casting a shadow in E CR downstream, while allowing forplasma-based transport through the cloud (right) leads to aslightly less disrupted cloud, greater E CR downstream, anda second bottleneck at the back edge seen most clearly in thelower right panel. enter the cloud, and field lines are allowed to evolve. Forvarying magnetic field strengths (top vs bottom rows)and varying ion fraction (left vs right columns), we seedifferences in how cosmic rays penetrate and pressurizethe cloud and in the cloud morphology (shown by whiteand yellow contours of density 1 and 10 cm − ).In the B = 5 µ G cases, magnetic field line warping isnegligible due to strong magnetic tension, and cosmicrays propagate straight across the cloud. Assuming fullionization, v st decreases sharply in the cloud, inducinga bottleneck on the front edge and a cosmic ray shadowbehind the cloud. This shadow fills in very slowly sincethe collisional loss time in the cloud is significant com-pared to the transport time. The cloud morphology isvery similar to what is seen in Wiener et al. (2019);Br¨uggen & Scannapieco (2020), where the top and bot-tom cloud edges with the lowest gas column get pushedand form tails behind the cloud, followed eventually byfull-cloud acceleration to the right. When plasma-basedtransport is included, the formation of these tails is notas dramatic since the cosmic ray pressure profile be-comes flatter in the cloud and therefore doesn’t push onthe gas as much. The cosmic ray pressure behind thecloud is clearly larger, as well, since the cosmic ray trans-port time is shorter. A second bottleneck, as seen in the 1D simulations, is faintly visible on the back side of thecloud where the transport speeds drops precipitously tothe gas Alfv´en speed, v A instead of v ionA .With a magnetic field strength of B = 1 µ G, corre-sponding to a plasma beta of 75 instead of 3, the mag-netic field topology plays a more significant role. Beforethe cosmic ray front even reaches the cloud, the preced-ing acoustic wave disturbs the magnetic field and warpsit around the cloud. When cosmic rays approach thecloud, they must follow the stretched field lines that aredraping around the cloud and therefore take longer toenter. For field lines that do penetrate the cloud, re-gardless of transport, the slower Alfv´en speed outsidethe cloud leads to a build-up of cosmic ray pressure anda severe bottleneck on the front edge.5.2.1.
Cloud Acceleration: 1D vs 2D
In each simulation, the cloud is eventually acceleratedbut to different extents depending on transport treat-ment, magnetic field strength, and dimensionality. Fig-ures 6 and 7 compare profiles down the cloud midplaneof cosmic ray energy density and gas density after t = 60Myrs. Solid lines show the 1D simulations, dashed linesshow the 2D simulations, and in two cases, we ran pc-resolution 3D simulations without plasma-based trans-port, shown by the dot-dashed lines.Let’s first analyze the 1D results without plasma-based transport (the solid lines in Figure 6), which arethe most straightforward to interpret. Cloud acceler-ation depends on the magnitude and duration of thecosmic ray pressure gradient. Keeping in mind that thecosmic ray pulse only lasts for a finite 30 Myrs, at whichpoint the upstream pressure will fade at a rate propor-tional to the Alfv´en speed, the cloud gets a push onlyuntil the cosmic ray pressures upstream and downstreamequilibrate. For a higher Alfv´en speed (higher magneticfield strength), cosmic rays build up a short-lasting andrelatively low pressure upstream, and since they travelmore quickly through the cloud, they lose less energyto collisions and emerge with a non-negligible pressuredownstream. So the cosmic ray pressure gradient issmall and also short-lived, leading to the least cloudacceleration for the B = 10 µ G simulation. As seen inFigure 6, the distance the cloud travels monotonicallyincreases for decreasing field strength.In 2D, the magnetic field topology is a complicatingfactor, as we see the low and high magnetic field casesboth give very little cloud acceleration. For high mag-netic field strength, field line warping is negligible, andthe setup is effectively 1D. Again, the pressure gradientis relatively small, and the pressures equilibrate quicklyto give very little cloud acceleration. On the contrary,2
Bustard et al.
Figure 6.
Profiles at t = 60 Myrs down the cloud axis showing the evolution of density and cosmic ray energy density for thefiducial cloud in 1D (solid lines), 2D (dashed lines), and for two simulations, 3D (dot-dashed lines). Each panel is for a differentmagnetic field strength, and each simulation is without plasma-based transport ( f minion = 1 . − µ G where the cosmic ray front is most effective at pushing the cloud. for low magnetic field strength, the field lines warp somuch that it allows cosmic rays to move around thecloud, taking pressure from the upstream, while cosmicrays instead converge and build up behind the cloud.The cosmic ray pressure drop across the cloud is small,which again leads to relatively little cloud acceleration.Intermediate magnetic field strengths B ≈ − µ Ggive a sweet spot for cloud acceleration. In these cases,field lines don’t warp as much, so the upstream pressureis higher, while the cosmic ray transit time across thecloud is intermediate, so the downstream pressure is low.This creates a sizeable and fairly long-lasting pressuregradient; however, note that the 3D simulations showflat cosmic ray pressure profiles after 60 Myrs and cloudsthat have better maintained their shape and moved lessdistance. While the differences between 2D and 3D aresmall compared to the differences between 1D and 2D,clearly the extra dimensionality plays a role, motivatinga larger sample of 3D simulations that we will run in thefuture. 5.2.2.
Effects of Plasma-Based Transport
Figure 7 shows the same 1D vs 2D comparison as be-fore, but we now turn on plasma-based transport with f minion = 10 − . Again, for the 1D simulations, the clouddistance increases monotonically for decreasing mag-netic field strength, as the same physical picture of pres-sure equilibration applies. Note, however, that for agiven magnetic field strength, the fully ionized case al-ways gives slightly greater cloud acceleration than whenpartially neutral gas is accounted for. This difference ismost dramatic when the magnetic field strength is highand plasma-based transport through the cloud leads toan especially quick pressure equilibration; for low mag-netic field strengths, the bottleneck at the cloud inter-face is severe and results in a pressure gradient thatpushes the cloud to nearly the same extent regardless of f minion . While the cosmic ray profiles are markedly differ-ent inside the clouds (steep for f minion = 1 . f minion = 10 − ), the total pressure drop from upstreamto downstream is comparably large in both cases, andit is the total cosmic ray pressure difference across thecloud, regardless of the smaller-scale pressure profile dif-ferences, that determines the force.In 2D, efficient cloud acceleration again disfavors highmagnetic field strengths, but clearly plasma-based trans-port has suppressed cloud acceleration in the intermedi- ow Cosmic Rays Navigate the Multiphase ISM Figure 7.
Same as Figure 6 but with f minion = 10 − . Time (Myrs) D M o m e n t u m ( g / s ) f minion = 1.0 B = 1 GB = 2 GB = 3 GB = 4 GB = 5 GB = 10 G 0 10 20 30 40 50 60
Time (Myrs) C l o u d V e l o c i t y ( k m / s ) f minion = 1.0 Time (Myrs) M c l / M c l f minion = 1.0 Figure 8.
For the fiducial 2D cloud, placed 100 pc from the left boundary, these panels show the total 2D momentum, volume-averaged velocity, and change in cloud mass M cl /M cl over time, defining the cloud as all gas with a temperature T < × K.For each simulation, we vary the magnetic field strength between 1 and 10 µ G, and f minion = 1 . ≈ µ G, in which case the field line warping is small, butthe Alfv´en crossing time between the source and the cloud, which partially determines the duration of the cosmic ray pressuregradient, is longer than the 10 µ G simulation. Lower magnetic field strengths yield more field line warping around the cloud,which inhibits acceleration. In all cases, the cloud mass decreases over time, which is to be expected since we don’t includeradiative cooling. Cloud acceleration is decreased by a factor ≈ ate field regime as well, as the density peaks shift to atmost 200 pc after 60 Myrs. In the fully ionized case,intermediate field strengths led to a sweet spot withonly small field line warping and moderate transportspeed across the cloud; combined this created a largepressure difference. Plasma-based transport, however,equilibrates the pressure on a shorter timescale, break-ing down the sweet spot. Low magnetic field strengthsallow for a greater build-up of cosmic ray pressure up- stream, but cosmic ray flow around the cloud once againmakes the pressure drop across the cloud almost negli-gible.Consistent with the morphological differences in Fig-ure 5, plasma-based transport keeps the cloud fairly in-tact, since there is no pressure gradient within the cloud,whereas the f minion = 1 . Bustard et al.
Time (Myrs) D M o m e n t u m ( g / s ) f minion = 10 B = 1 GB = 2 GB = 3 GB = 4 GB = 5 GB = 10 G 0 10 20 30 40 50 60
Time (Myrs) C l o u d V e l o c i t y ( k m / s ) f minion = 10 Time (Myrs) M c l / M c l f minion = 10 Figure 9.
Same as Figure 8 but with f minion = 10 − . Cloud acceleration is decreased by a factor ≈ based transport simulations resulted in high cloud ve-locities, it’s intriguing to wonder whether plasma-basedtransport and a larger cosmic ray flux could lead to sig-nificant acceleration of molecular gas while keeping it in-tact instead of stretched apart. This will be addressedwith future simulations.Our main conclusion is that cosmic ray accelerationof partially neutral clouds is less efficient than acceler-ation of fully ionized clouds, but the differences in im-parted momentum are small (less than a factor of 2).This is quantified in Figures 8 and 9, which show thecloud momentum, velocity, and change in mass, defin-ing the cloud to be all gas with temperature T <
2x 10 K. As expected, the momentum and velocitieswith plasma-based transport included are lower by afactor ≈ t cc ∼ ξ / r cl /v is then shorter, where ξ is the cloud overdensityrelative to the background, r cl is the cloud radius, and v is the velocity of the ambient medium relative to thecloud (Klein et al. 1994). It’s worth noting again thatthere is no radiative cooling in our simulations, whichmay suppress the ablation of gas at the cloud boundaryand even facilitate cloud growth as the surrounding gasmixes and condenses onto the cloud (Gronke & Oh 2018,2020; Kanjilal et al. 2020; Banda-Barrag´an et al. 2020).5.2.3. Comparison to Previous Studies
While we’ve only simulated cosmic ray to thermalpressure ratios of at most a few, therefore acceleratingclouds to only a few km/s, we can already gleam someof the implications for galactic wind driving and com-pare to the broader literature of cloud “wind-tunnel” simulations. Without considering the effects of neutralparticles, we have extended the work of Wiener et al.(2017, 2019) with a faster, more accurate method of cos-mic ray transport (Jiang & Oh 2018) in the Athena++code. This speed-up allowed us to, amongst otherthings, further explore the parameter space of magneticfield strengths and elucidate the implications of field linewarping for cloud acceleration.Our simulations without plasma-based transport aresimilar to another recent simulation suite (Br¨uggen &Scannapieco 2020) but with some differences. Br¨uggen& Scannapieco (2020) simulate a cosmic ray front withratio of cosmic ray pressure to thermal pressure rang-ing from a few to a few tens, impinging upon a warm,fully ionized halo cloud with density 10 − g cm − andradius 100 pc, much larger and lower density than ourrepresentative ISM cloud. Unlike our simulations, whichdefine the cosmic ray flux at the boundary for a finitetime, they maintain a constant cosmic ray energy den-sity at the boundary for the entirety of their simula-tion. While we include cosmic ray heating and colli-sional losses (which would be negligible for a diffuse halocloud), they do not include heating but do include ra-diative cooling.These differences in implementation and setup makeit hard to draw fair comparisons. One differing result,however, is that we find a clear dependence of cloudacceleration on magnetic field strength, while they findonly a small dependence (see e.g. Figure 8 of Br¨uggen& Scannapieco (2020)). This is likely because 1) theyhold the boundary cosmic ray energy density fixed, whileour upstream cosmic ray pressure depends on the Alfv´enspeed, and 2) they implement a switch that constrainscosmic rays to accelerate only the cells with densitygreater than 10% above the ambient density. If im-plemented in our simulations, these differences wouldchange our results, as the amplitude and duration ofthe upstream cosmic ray pressure gradient would haveno magnetic field dependence, and the acoustic pulse ow Cosmic Rays Navigate the Multiphase ISM A MULTIPHASE COSMIC-RAY OBSTACLECOURSEWe now run a suite of simulations of cosmic ray prop-agation through multi-clump “obstacle courses.” Themean density at the bottom boundary in each simulationis ¯ n ≈ cm − and initially falls off exponentially with ascale height H = 0.25 kpc. The magnetic field is con-stant in the y-direction (vertical direction) with a fidu-cial value of 5 µG , and the pressure is constant through- out the simulation box at 3 . × − dyne cm − . Ineach simulation, the domain extends to the y = 4 kpctop boundary, which is sufficiently far away from theclumpy ISM to alleviate worries about boundary effects.On top of this tapered density profile, we impose a log-normal distribution of density perturbations (Br¨uggen2013) with 20 Fourier components in each of the x-, y-,and (if in 3D) z-directions. As with the mean density,the perturbation amplitude also decreases exponentiallywith scale height 0.25 kpc in order for the density vari-ance to attenuate above the disk into a smooth halo.The density profile is then: n = ¯ ne − y/H e αf ( x,y,z ) (13) f ( x, y, z ) = e − y/H (1 − e − y ) (cid:88) k x ,k y ,k z A i sin (2 πk x x/L + φ i ) sin (2 πk y y/L + θ i ) sin (2 πk z z/L + η i )(14)where x, y, z are in units of kpc, H = 0.25 kpc, and k x , k y , k z ∈ [10 , A i and φ i , θ i , η i are randomly chosenfrom a uniform distribution between [0, 1] and [0, 2 π ], re-spectively. The parameter L controls the characteristicsize of perturbations, while α sets the density contrast.The extra factor of (1 − e − y ) in the perturbation equa-tion excludes the clumps away from the bottom bound-ary. We found in test simulations that if a dense clumpis on the boundary, the outflow boundary we imposewill copy the gas density into the ghost cells and there-fore draw in a large amount of artificial mass. Since,as we’ll see, cosmic rays push on clumps differently de-pending on the transport model, this results in differenttotal masses in the simulation box, which muddies ouranalysis of gamma-ray emission, gas momentum, etc.Our clump setup varies, with (L, α ) between (5, 1.5)and (2, 0.5), giving different density probability distri-bution functions (PDFs) shown in Figure 10. For the gaspressure we initialize, at the mean density ∼ cm − , thetemperature is ∼ K. For L = 5, this gives a porousISM dominated by large ( ≈ −
100 pc) over/under-densities. The simulation domain extends to ± ± ≈ −
50 pc) more typical of real ISM clouds, thesedomain widths are sufficient to present cosmic ray inter-actions with clouds of varying densities; however, theynecessitate higher resolution to more properly resolvethe cloud interfaces, which we’ve found to be important( § Bustard et al.
Figure 10.
Density probability distribution function (PDF)for each combination of (L, α ) we present. pc until y = 2 kpc, at which point we make use of staticmesh refinement to decrease resolution to (16, 8) pc andsave computation time.It’s important to note that these simulations do not in-clude cooling or gravity and are not meant to substitutefor more realistic, stratified ISM simulations, thoughthey do replicate the propagation of cosmic rays throughlayers of a fixed thickness and mean column density.Having a constant magnetic field strength out to manyscale heights is also not realistic, so we focus our analysisto the first kpc above the midplane.6.1. Cosmic Ray Propagation and Energy Loss
Figures 11 and 12 show snapshot at t = 40 Myrs ofour B = 5 µ G and B = 1 µ G simulations with meandensity = 1.0 cm − , (L, α ) = (5, 1.5), and f minion = 1 . f minion = 10 − (right). The plots are orientedsuch that cosmic ray flux is input at the lower boundary,and magnetic field lines are vertical. The contours showdensities of n = 1 cm − and n = 10 cm − . The under-dense, inter-cloud region extends down to densities of n ≈ − cm − (see Figure 10). Such a configuration, ifentirely ionized, gives a range of Alfv´en velocities pre-dominantly between 10 and 10 cm/s, but when neutralparticles are accounted for, the lowest Alfv´en velocity is ≈ cm/s, with ion Alfv´en velocities in cloud interiorsreaching beyond 10 cm/s.This difference in propagation speed is clearly re-flected in the top row of Figure 11, which shows thatcosmic rays with plasma-based transport have propa-gated further into the domain and have reached the in-ner “halo”, where the average density has dropped, andcosmic rays from here will free-stream out of the simula-tion domain at high Alfv´en velocities since v A ∝ ρ − / for constant B. In the fully ionized case, the cosmic rays that have managed to propagate furthest are thosethat have avoided cold clouds and squeaked through thecracks. As in the 2D, single cloud simulations, for an ini-tial plasma β = 3, magnetic field line warping plays onlya limited role, meaning most cosmic rays are eventuallyforced into the clouds (and additionally, accelerate theclouds). When accounting for plasma-based transport,cosmic rays leak into the dense cores and quickly reap-pear on the other side. As in Figure 5, there are sharpsteps in cosmic ray energy at cloud edges, both at thefront and the back edges in the f minion = 10 − case, aseach interface induces a bottleneck.For a magnetic field strength of 1 µ G (Figure 12),magnetic field line warping dominantly determines howcosmic rays navigate the medium. Now, cosmic rayswith either transport model reside primarily in the inter-cloud regions. The lower magnetic field strength alsoleads to a greater build-up of cosmic ray pressure thatexerts a force on the medium. As in the single cloudsimulations of § v ionA . This is consistent with previous work that found,using a Monte Carlo method, that cosmic rays in M82-like starburst galaxies would generally sample the meandensity (Boettcher et al. 2013). At temperatures belowT ≈ K, as cosmic rays encounter clumps, the PDFsand concentrations start to diverge. With plasma-based ow Cosmic Rays Navigate the Multiphase ISM y ( kp c ) f minion = 1.0 f minion = 10 y ( kp c ) y ( kp c ) y ( kp c ) y ( kp c ) E c r ( g c m s ) v i o n A ( c m / s ) T ( K ) C o lli s i o n l e ss ( e r g c m s ) C o lli s i o n a l ( e r g c m s ) Figure 11.
Snapshots at t = 40 Myrs for the ISM setup with B = 5 µ G, mean density = 1.0 cm − , and (L, α ) = (5, 1.5). Left panels:
Fully ionized assumption ( f minion = 1 . Right panels:
Plasma-based transport included ( f minion = 10 − ). The toppanel shows cosmic ray energy density and density contours of 1 and 10 cm − , which shows clear variations in how cosmicrays preferentially penetrate or flow around cold clouds. The second row shows v ionA , followed by gas temperature, collisionlessenergy loss rate ( | v ionA · ∇ P CR | ), and collisional ( ∝ e CR n gas ) loss rate, again with dashed contours showing densities of 1and 10 cm − . Collisionless heating at the cloud interfaces is apparent in both transport cases, while collisional energy lossis preferentially higher in dense gas when plasma-based transport is included. A video version of this figure can be found athttps://bustardchad.wixsite.com/mysite/visualizations Bustard et al. y ( kp c ) f minion = 1.0 f minion = 10 y ( kp c ) y ( kp c ) y ( kp c ) y ( kp c ) E c r ( g c m s ) v i o n A ( c m / s ) T ( K ) C o lli s i o n l e ss ( e r g c m s ) C o lli s i o n a l ( e r g c m s ) Figure 12.
Same as Figure 11 but with B = 1 µ G instead of 5 µ G. Magnetic field lines now warp around the cold clumps,squeezing cosmic rays through the gaps in the ISM instead of funneling them through clouds. This effect dominates overdifferences in transport, as the left and right columns (varying f minion ) now show smaller differences. Also note that, comparedto the B = 5 µ G simulations shown in Figure 11, there is now a larger build-up of cosmic ray pressure that more effectivelyaccelerates the gas column. This is especially true for the hot gas, which gets pushed out by the cosmic rays sweeping throughthe under-dense channels. A video version of this figure can be found at https://bustardchad.wixsite.com/mysite/visualizations ow Cosmic Rays Navigate the Multiphase ISM Temperature (K) E C R P D F L = 5, =1.5, B = 5 GL = 2, =1.5, B = 5 GL = 5, =1.5, B = 1 GL = 2, =1.5, B = 1 G Temperature (K) E C R C o n c e n t r a t i o n Temperature (K) G a mm a - R a y E m i ss i o n P D F Temperature (K) G a mm a - R a y C o n c e n t r a t i o n Figure 13.
The left column shows the probability distribution function (PDF) of cosmic ray energy density E CR (top) andgamma-ray emission ∼ E CR n gas (bottom) in 32 temperature bins ranging from T = 10 − × K. The right column showsthe concentrations, obtained by dividing the PDFs in the left column by the temperature PDF, which gives a comparison ofhow concentrated the quantities are compared to the fraction of the ISM at each temperature bin. If the fraction of energyor emission equals the volume fraction at that temperature bin, the concentration will be 1.0. Clearly there are differencesdepending on transport model, with f minion = 1 . f minion = 10 − shown by dashed lines, as well asvarying ISM parameters and magnetic field strengths (different colors). Cosmic rays in all cases preferentially sample the meandensity and temperature ≈ × , where the gas is fully ionized, Alfv´en speeds are low, and cosmic ray residence times arelong. Plasma-based transport leads to preferentially more occupancy in dense, cold clouds but less occupancy at transitiontemperatures T ∼ × K. This trade-off leads to roughly the same total gamma-ray luminosity, regardless of transport( § transport (dashed lines), there is a steep drop in cosmicray occupancy near temperatures T ≈ × K, where v ionA rapidly increases. At similar temperatures, cosmicrays with f minion = 1 . f minion = 1 . f minion = 10 − curves drop more gradually, as cosmic rays uniformly fill cold cloud inte-riors.Lowering the magnetic field strength gives the mostdramatic change. Cosmic rays streaming along thenow-warped magnetic field lines preferentially squeezethrough the under-dense channels between clouds. Thecosmic ray and gamma-ray concentrations are both sup-pressed at low temperatures. Cosmic rays also do notoccupy the high temperature space: this low-density gasgets pushed into the halo, forced out by the cosmic rayfront that sweeps between clouds. In these low mag-netic field cases, then, the great majority of cosmic rayslie within a narrow temperature range near the mean.0 Bustard et al.
The bottom row of Figure 13 shows the PDF ofhadronic gamma-ray emission ∝ e CR n gas . The trendsfrom the top row are clearly imprinted. For both trans-port models, gamma-ray emission is very concentratednear the mean density when B = 1 µ G, but with B= 5 µ G, a substantial fraction of the gamma-ray emis-sion comes from cold clouds at temperatures
T < K. The low-temperature bins that dominate emissiondepend on transport model. At intermediate densities n ≈ −
10 cm − (temperatures around 8 × K), col-lisional energy losses are decreased considerably whenplasma-based transport is included, with more energyloss occurring at the higher densities that fast-movingcosmic rays can leak into. In the fully ionized case, thehighest densities are inaccessible to cosmic rays, as theylose energy instead in cloud interface layers ( T ≈ × K) before penetrating into cloud cores. As we’ll see, thistrade-off still leads to very similar total collisional lossrates ( § v ionA that induces another steep cosmic ray pressure gradient.6.2. Total Cosmic Ray Calorimetry
Figure 14 compiles the total collisional and collision-less energy loss from each of the mock ISM simulationswe ran. Each symbol represents a different transportmodel, including a few simulations where we includedan additional, constant cosmic ray diffusion coefficient;this boosts cosmic ray escape most drastically in thediffuse medium and cloud interfaces where the residencetimes are otherwise very long. Each color is a differentISM configuration, varying L, α , and B. Black dashedlines show contours of constant total cosmic ray energyloss.What one can see most clearly is that, although col-lisionless energy loss can vary by orders of magnitude,collisional energy loss only varies within a factor of a fewregardless of transport model. All models without ad-ditional diffusion are clustered within the outlined box,meaning that, although the collisions occur in very dif-ferent places (see Figure 13), the total collisional energyloss is about the same.Zooming in on the clustered region of Figure 14, wedo see some variation depending on ISM configurationand transport model, but the changes are within a fac- tor of a few. For instance, while a lower magnetic fieldstrength on average leads to lower transport speeds andlonger cosmic ray residence times in the ISM, the net col-lisional energy loss is comparable to that with a highermagnetic field strength. With weak magnetic fields, thecollisionless energy loss ∝ v ionA · ∇ P CR decreases since v ionA ∝ B also decreases; the decrease is only by a factorof ≈
2, however, since slower transport speeds lead to abuild-up in cosmic ray pressure and higher ∇ P CR .Varying α , which determines the density contrast, alsomakes a difference. With higher (lower) density con-trast, the collisional energy loss ∝ n gas increases (de-creases) accordingly. Varying L has almost no conse-quence for collisional loss, however. The slight increasein collisions for L = 2 compared to L = 5 can be at-tributed to the slightly larger density range for L = 2(see Figure 10), which follows the same explanation asvarying α .We also ran three of our clump simulations in 3D, twoof which had f minion = 1 . f minion = 10 − .Each of these simulations had B = 5 µ G, which as wesaw in the single cloud simulations of § β , the differences between 2D and 3D may be-come more important, which we will explore in futurework.It’s insightful to also disentangle the differences in col-lisionless energy loss. If one were to draw a line diag-onally across the plot from lower left to upper right,that would show equal collisional and collisionless en-ergy loss; we see that in all simulations run, collisonallosses exceed collisionless losses, but this gap narrowswhen the density contrast (controlled by α ) decreasesand when plasma-based transport is included. Compar-ing e.g. L = 2, α = 1 . f minion ,we see that plasma-based transport decreases collisionalenergy loss but boosts collisionless energy loss, endingup with a similar total energy loss. The increase in col-lisionless energy loss can be explained by the secondbottleneck that occurs at each back cloud edge (see e.g.Figure 4), therefore giving two interfaces that induce en-ergy loss instead of just one. So plasma-based transportnot only fails to decrease collisional losses substantiallybut can also offset this reduction by increasing collision-less energy loss – the net result is a very similar cosmicray calorimetry.We find larger differences in collisionless and colli-sional energy loss when we include an additional dif-fusion term. Under the assumption of full ionization, anadditional diffusion coefficient of 3 × cm s − sig- ow Cosmic Rays Navigate the Multiphase ISM Figure 14.
Accumulated collisional and collisionless cosmic ray energy loss within one scale height (250 pc) for each of themock ISM simulations. Different symbols represent changes in transport model, while color represents the ISM setup varyingL, α , and B. The black dashed lines, both in the large figure and smaller zoom-in, show contours of constant total energy loss.Despite some significant differences in collisionless energy loss, especially when an additional, constant diffusion term is included,the collisional energy loss (and hence, diffuse gamma-ray emission) only varies within a factor of a few for most models. Changesin ISM setup and 2D vs 3D impart only small variations in total energy loss that are explained in the text. nificantly flattens the cosmic ray pressure gradient andleads to a decrease in collisionless energy loss by a factorof tens. Interestingly, with f minion = 10 − , the same dropin energy loss can be achieved with a smaller diffusioncoefficient of 3 × cm s − . Without this additionaldiffusion, we saw that cosmic rays still bottleneck at thecloud edges, which severely decreases the diffusive fluxthrough the clouds. This bottleneck is effectively wipedout when even a small constant diffusivity is added, as itsmooths the cosmic ray pressure gradient and preventsthe generation of confining waves; this maintains a largediffusivity through the partially neutral cloud interior,and transport no longer simply follows the ion Alfvenspeed. The outcome is that cosmic rays do not sensethe cloud and effectively free-stream through it. Whilethis localized faster transport decreases the collisionlessenergy loss by orders of magnitude, the decrease in col- lisional loss is only a factor of a few since collisions pri-marily come from cosmic rays with long residence timesin the warm ionized gas (see e.g. Figure 13).This finding is consistent with previous arguments inChan et al. (2019) and Hopkins et al. (2020b), who com-pared full-galaxy simulations with varying diffusive andstreaming transport with gamma-ray data. They foundthat, especially for dwarf galaxies, diffuse gamma-rayemission is far over-produced (by a factor of ten or more)compared to observations of Local Group dwarf galax-ies unless the cosmic ray escape time is significantly re-duced; despite strong outflows generated in these dwarfgalaxies that would naively lead to fast advective escape,an additional cosmic ray diffusion coefficient of 3 × cm s − is needed to lower emission to reasonable levels.More specifically, Hopkins et al. (2020b) found that asignificant fraction of gamma-ray emission was coming2 Bustard et al. from the diffuse medium, where long cosmic ray resi-dence times compensated for the low gas density to en-hance collisional losses ∝ e CR n gas . The simulations thatbest fit observed gamma-ray and cosmic ray grammagedata, then, had to allow for fast cosmic ray transport inthe diffuse medium, while fast transport in cold, densegas (where one would naively think most collisions oc-cur) was not sufficient. Our study seems to corrobo-rate this result, as we find collisional energy loss (hence,gamma-ray luminosity) to be fairly degenerate as wevary ISM setups and transport models. The largest de-crease in collisions comes when an extra diffusion termis present – an additional simulation with κ || = 3 × cm s − (not shown in Figure 14) continues the trendtowards lower collisional energy loss.6.3. Influence on Clouds
Farber et al. (2018) probed the implications of cosmicray decoupling for galactic wind driving and star forma-tion feedback. They compared three transport models:cosmic ray advection, field-aligned diffusion with a sin-gle diffusion coefficient (one- κ model), and field-aligneddiffusion with a higher diffusion coefficient in cold gas(two- κ model) to mock the effects of ion-neutral damp-ing. Comparing the one- κ and two- κ models, they foundthat boosting the diffusivity in the cold phase led to afaster, hotter outflow since the cold phase received lessimpulse from cosmic ray fronts. Changing f minion , we findthat this trend partially occurs. Figure 15 shows the mo-mentum in the full simulation box as a function of timefor our (L, α ) = (5, 1.5) ISM setup. The top panel showsthe momentum regardless of gas temperature, while thebottom two panels split the momentum into “cold” gas(T < × K) and “hot” gas (T > × K). Whenplasma-based transport is included, slightly less momen-tum is imparted to both the cold and hot gas, but onlyon the level of tens of a percent. As argued in § Time (Myrs) D M o m e n t u m ( g / s ) Temperature < 8 x 10 K L = 5, = 1.5, B = 5 GL = 2, = 1.5, B = 5 GL = 5, = 1.5, B = 1 GL = 2, = 1.5, B = 1 G
Time (Myrs) D M o m e n t u m ( g / s ) Temperature > 2 x 10 K Figure 15.
Momentum measured in 2D clump simulationswith varying L and α and varying transport, with f minion = 1 . f minion = 10 − (dashed lines). The top panelshows the momentum just in the cold gas (T < × K),and the bottom panel shows the momentum in the hot gas(T > × K). Note the change in y-axis reflecting thatmost of the momentum, in either transport model, is in thecold phase, but lower magnetic field strengths do increasethe hot gas momentum quite substantially as cosmic rayssqueeze out the gas between dense clouds. With a lowerfield strength, the cold momentum increases somewhat, aswell, due to an increased cosmic ray pressure gradient as wesaw in § minishes, and the momentum slightly drops. The coldgas momentum rises a bit more slowly; cosmic ray bot-tlenecks take time to build up (on the order of the Alfvencrossing time between the cloud and the source), and thecosmic ray pressure releases only after a long time sincethe cold clouds move relatively slowly. Pressure gradi-ents on the cold gas then build to greater amplitudesand longer durations than on the hot gas, leading to alarger but more gradual acceleration of the clouds.Rather than plasma-based transport leading tochanges in cloud acceleration, the largest changes in mo- ow Cosmic Rays Navigate the Multiphase ISM f minion = 1 . f minion = 10 − models show an increase in impulse. The cold gas mo-mentum increases because cosmic ray bottlenecks aresevere in either case at the cloud edges. The hot gasmomentum increases because cosmic rays squeeze be-tween the gaps in clouds, sweeping out the hot gas infront of them (see Figure 12).Interestingly, changes in momentum also correlatewith the characteristic size of clumps controlled by our L parameter. More momentum is imparted especiallyto the cold gas when L is smaller; this makes sense asthe limit L → § DISCUSSION AND LIMITATIONSThis work focuses on the interplay between cosmicrays and cold clouds, taking a critical look at the rolesof bottlenecks induced upstream and elevated transportspeeds in partially neutral clouds. We focus on the self-confinement model of cosmic ray streaming, with advec-tive flux modified by varying ion Alfv´en speed and diffu-sive flux resulting solely from ion-neutral damping. Wehave not, until this point in the paper, explored the roleof other damping mechanisms or instabilities that sup-press or enhance confinement. While we mainly leavethis to future work, we did re-run a subset of our 1Dcloud simulations including nonlinear Landau damping(NLLD), which results from thermal particles taking en-ergy from interacting Alfv´en waves. Following Hopkinset al. (2020b), the damping rate can be written asΓ
NLLD = (cid:20) ( γ CR − π / (cid:18) c s v ionA r L l CR (cid:19) (cid:18) e CR e B (cid:19)(cid:21) / (15)This is proportional to the cosmic ray pressure gradi-ent, suggesting that waves may be efficiently damped inthe cloud interface region where a steep pressure gradi-ent forms, but the diffusive flux is inversely proportional: κ || ∝ ( ∇ P CR ) − / (Loewenstein et al. 1991). Because ofthis, our simulations show only a small additional diffu-sion coefficient < cm /s in the cloud interface, andthe resulting cosmic ray and cloud evolution is almostexactly the same whether we account for NLLD or not.Another limitation, which we plan to explore in futurework, is that these simulations do not include cooling orconduction. These play important roles in cloud sur-vival during acceleration (Br¨uggen & Scannapieco 2016; Gronke & Oh 2018), and they also set the thermody-namic state of the gas and the width of intermediatetemperature transition regions, or cloud envelopes. Notethat a typical Fields length in 10 K gas is less than0.1 pc, much smaller than our fiducial cloud setup witha 5 pc interface. Heating from cosmic ray streaming,while it noticeably heats the interface in our simulations,is not likely to set the interface width in real, coolingclouds. From our simulations, we find collisionless lossrates (hence, heating rates) less than 10 − erg cm − s − , subdominant compared to the expected radiativecooling and conductive heating rates if those were takeninto account (see also Figure 9 of Everett & Zweibel(2011)).Models and observations of molecular clouds exposedto interstellar radiation fields, however, do show a grad-ual transition between hot and cold phases, with warm,ionized envelopes surrounding the cold cores (Goldsmith& Li 2005). Indeed, a growing literature on cosmic raypenetration into molecular clouds has come to similarconclusions about the trapping of cosmic rays by eitherextrinsic or self-generated turbulence in interface regions(Morlino & Gabici 2015; Schlickeiser et al. 2016; Dogielet al. 2018; Ivlev et al. 2018; Inoue 2019), especially nearcosmic ray sources where it appears that cosmic rayscan be confined even in the presence of neutral particles(Nava et al. 2016; Brahimi et al. 2020).Constraining these theoretical results with observa-tions is important and already underway to truly diag-nose cosmic ray transport in dense, cold clouds. Silsbee& Ivlev (2019) considered whether the low-energy (lessthan GeV) cosmic ray spectrum from the Voyager probewould be better fit by diffusive or free-streaming prop-agation of cosmic rays penetrating through a large col-umn density into the Local Bubble. While confinementwill be stronger at these energies than for GeV cosmicrays, they find that diffusive propagation is a better fitto the spectrum than free-streaming, which would addi-tionally require an unreasonably high column density tofit the Voyager data.Fujita et al. (2020) model free-streaming and diffu-sive propagation (due to preexisting MHD turbulence)in molecular clouds, obtaining profiles of the ionizationrate, 6.4 keV line flux, and gamma-ray emission. Whilethere are only a few observed supernova remnants forwhich all three of these diagnostics are available, thequalitative comparison in their Discussion section findsthat both diffusive and free-streaming propagation arelikely realized.Joubaud et al. (2020) studied a diffuse cloud near theOrion-Eridanus superbubble that was observed to havea 34% lower cosmic ray flux compared to the local flux4 Bustard et al. estimate. As the authors note, the estimates of Everett& Zweibel (2011) suggest only a slight drop in cosmicray flux within cold clouds and would require a muchhigher than observed magnetic field strength to accountfor such a 30% drop. Our simulations assume cloudparameters closer to the Eridu estimates ( n ≈ − and B ≈ µG ) and find that the cosmic ray pressuredrops by a factor ≈ CONCLUSIONSMany MHD + cosmic ray simulations, aimed at mod-eling galaxy evolution, now take into account the self-confinement of GeV-energy cosmic rays through fre-quent scattering off cosmic ray - generated Alfv´en waves.While this kinetic process occurs at scales of order thecosmic ray gyroradius, many orders of magnitude belowthe current resolution limit of galaxy evolution simu-lations, fluid theories of this streaming transport havebeen developed and, with the advent of various numeri-cal techniques, implemented in MHD codes for use in agalaxy-scale context. With these tools now available,our understanding of cosmic ray influence in the in-terstellar, circumgalactic, and intracluster medium hasrapidly progressed; however, as of this writing, two im-portant trends are evident in most published simula-tions: 1) it’s assumed that GeV-energy cosmic ray trans-port proceeds as if the background thermal gas is every-where fully ionized. In reality, the waves the cosmic raysresonate with propagate only in the ions, and streamingis therefore at the ion Alfv´en speed v ionA = v A / √ f ion .Also, ion-neutral friction can damp the confining waves,thereby decoupling cosmic rays and greatly increas-ing their diffusive flux; 2) partially neutral ISM cloudsand other density irregularities that induce decouplingand nonlinear cosmic transport are frequently under-resolved. If we imagine the multiphase ISM as an obsta-cle course of coupled vs decoupled regimes, further com-plicated by hadronic and Coulomb collisions betweencosmic rays and ambient gas, then how do cosmic raysnavigate this medium, and can we constrain this trans- port with observations? The answers have implicationsfor ISM dynamics, interpretation of observations, andthe acceleration of cold interstellar gas in galactic windsand fountains.In this paper, we presented a suite of high-resolution,idealized simulations of cosmic ray fronts navigating themultiphase ISM. We ran each simulation twice: first,assuming the medium is fully ionized and, second, ac-counting for changes in ionization that induce “plasma-based transport”. We began with 1D simulations of anenergetically significant cosmic ray front hitting a sin-gle, partially neutral cloud with radius and density typ-ical of ISM conditions. We then explored how cosmicray transport, cosmic ray energy loss, and cloud accel-eration and morphology varied in multiple dimensionsand varying magnetic field strength. We then ran thesame cosmic ray fronts through mock ISM setups: slabsof lognormally-distributed multi-phase clumps on top ofan exponentially decaying mean density. We assessedhow cosmic rays sampled the ISM and how plasma-basedtransport affected cloud acceleration and collisional andcollisionless energy loss.Our main findings were: • Cloud interfaces can play a crucial role in settingthe cosmic ray propagation through the cloud in-terior. The drop in Alfv´en speed while the gasis still fully ionized leads to a bottleneck in theinterface; this reinforces the cosmic ray pressuregradient that generates confining waves and lockscosmic rays to the wave frame. The effects of ion-neutral damping are thereby suppressed, and theadvective flux rather than diffusive flux generallydominates transport ( v st ≈ v ionA ). • The stair-step structure that develops at cloudinterfaces facilitates collisionless energy transferfrom the cosmic rays to waves ∝ v ionA ∇ P CR . Un-like for a fully ionized cloud where a single cos-mic ray gradient forms (Wiener et al. 2017, 2019),there’s now a bottleneck at both the front and back cloud edges. In our mock ISM simulations,this extra energy loss partially offsets a reductionin collisional energy loss, keeping the total cos-mic ray calorimetry similar to when full ionizationis assumed. The additional gas heating may alsohave implications for ion abundances and kinemat-ics at cloud interfaces, though we estimate that therole of heating would be secondary to cooling andconduction. We plan to explore this further bymore self-consistently modeling interface proper-ties rather than leaving the interface width as afree parameter. ow Cosmic Rays Navigate the Multiphase ISM • Cosmic rays sample the ISM differently when vary-ing ionization is accounted for. Plasma-basedtransport allows cosmic rays to leak into densecloud interiors while, assuming full ionization, cos-mic rays are confined to thin cloud boundary lay-ers. These differences are imprinted on gamma-rayemission maps. • Total collisional losses are very similar for the twotransport models. With plasma-based transport,short cosmic ray residence times in clouds are off-set by the denser gas they access. So while gamma-ray emission maps with the two transport modelslook very different, the gamma-ray luminosity isdominantly set by the cosmic ray transport in thewarm ionized medium, which is the same in eachtransport model. This result is consistent withHopkins et al. (2020b). • When ionization is accounted for, cloud acceler-ation is decreased since plasma-based transportin cloud interiors flattens the cosmic ray pressureprofile. The change is not as large as one mightexpect, however, as our simulations exhibit a fac-tor < κ ” models depending on gas temperature, this methodwill more fully capture the generation of bottlenecks andsteep pressure gradients at cloud interfaces, as well asplasma-based transport in cloud interiors that is robust,at least in these simulations, to changes in ionizationfraction. More detailed modeling of ionization from lo-cal sources and lower energy cosmic rays could alter thispicture and should be a next step.ACKNOWLEDGMENTSThe authors are greatly indebted to Yan-Fei Jiangfor providing us with his two-moment method of cos-mic ray transport in Athena++. We also thank JoshWiener, Mateusz Ruszkowski, Ryan Farber, Peng Oh,and Max Gronke for stimulating discussions. This re-search was supported in part by the National ScienceFoundation under Grant No. NSF PHY-1748958 andby the Gordon and Betty Moore Foundation throughGrant No. GBMF7392. EGZ is happy to acknowl-edge support from NSF Grant AST 2007323. This workmade extensive use of the yt (Turk et al. 2011) softwarepackage, a toolkit for analyzing and visualizing quantita-tive data. Computations were performed using the Ex-treme Science and Engineering Discovery Environment(XSEDE), which is supported by National Science Foun-dation grant number ACI-1548562 (Towns et al. 2014).Specifically, our computing resources stemmed from al-location TG-AST190019 on the Stampede2 supercom-puter.APPENDIX.1. Effect of Resolution
In this section, we carry out a small resolution study, primarily in 1D, to test convergence and inform larger-scalesimulations of galaxy evolution, which may under-resolve cosmic ray interactions with multiphase structures. Figure16 shows our fiducial 1D cloud of radius 10 pc and interface width 5 pc simulated at resolutions of 2, 1, 0.5, and 0.25pc. As noted in Section 5.1, not resolving the cloud interface is akin to having a physically smaller interface width.Lower resolution leads to a smaller total drop in cosmic ray energy through the cloud, less heating at the front edge,and different cloud morphologies (note the disappearance of the upstream “thumb” caused by heating). Convergenceis achieved at a resolution of 0.5 pc, which resolves the 5 pc - thick interface by a sufficient number of cells. Similarly,our simulations with a 10 pc interface width reached convergence at 1 pc.6
Bustard et al.
50 pc 100 pc 150 pc D e n s i t y ( g c m ) E C R ( e r g c m ) v st = v ionA + Ion-Neutral Damping, f minion = 10 Figure 16.
Density and cosmic ray energy density for the fiducial cosmic ray front impinging on a cloud of radius 10 pc andinterface width 5 pc. Each panel corresponds to a different resolution. When the interface region is well-resolved ( ≈
10 cells perinterface), the results are well-converged. With poorer resolution, the cosmic ray front doesn’t respond to the drop in Alfv´enspeed at the interface, therefore barreling through the cloud and appearing on the other side with a larger cosmic ray energy.The spike in density on the front cloud edge, caused by cosmic ray heating, is also not as noticeable at low resolution.
Because the total pressure drop determines the force on the cloud, we expect that under-resolved simulations willunderestimate cloud acceleration. To test this, we ran our 2D cloud simulations (with f minion = 10 − ) at resolutionsof 4, 2, and 1 pc. The resulting cloud momenta and velocities are shown in Figure 17. As expected, deterioratingresolution decreases cloud acceleration and, hence, the perceived effectiveness of cosmic ray-driven feedback.While this pc or sub-pc resolution requirement is initially daunting, it’s worth noting that, if transport is dominatedby the advective rather than diffusive flux, such problems become more computationally tractable. Our choice to set ahigh V m = 10 cm/s, advective cap of 10 cm/s, and diffusive cap of 3 × cm /s in the first place was a pre-emptivemeasure anticipating that large diffusive fluxes, in particular, would need to be accommodated. As we found, that israrely the case, and transport speeds are generally limited by v ionA . We found only small differences, then, when V m and the diffusive and advective caps are each decreased by a factor of 5-10, which should provide a glimmer of hopefor convergence of future large-scale simulations with plasma-based transport. Time (Myrs) D M o m e n t u m ( g / s ) f minion = 10 B = 5 GB = 1 G 0 10 20 30 40 50 60
Time (Myrs) C l o u d V e l o c i t y ( k m / s ) f minion = 10 Figure 17.
Cloud momentum (left) and velocity (right) at varying resolutions of 1 pc (solid lines), 2 pc (dashed lines), and 4pc (dot-dashed lines) for the 2D cloud simulation of § µ G (black lines) and B = 1 µ G (red lines). Deterioratingresolution leads to less cloud acceleration, since cosmic ray pressure gradients are smoothed at the unresolved cloud boundary.Even at the fiducial 1 pc resolution, we don’t appear to be reaching convergence, as expected from Figure 16.
Finally, we test how cosmic rays sample the ISM and transfer energy in our clumpy simulations with lower resolutions.We again focus on plasma-based transport with f minion = 10 − , and we test resolutions of 32, 8, 4, and 2 pc in our (L, α ) = (2, 1.5) setups. The resulting cosmic ray energy density PDFs and gamma-ray luminosity histograms are shownin Figure 18. Lowering the resolution smooths cloud interfaces, allowing cosmic rays to leak more easily into cold, ow Cosmic Rays Navigate the Multiphase ISM Temperature (K) E C R P D F L = 2, =1.5, B = 5 GL = 2, =1.5, B = 1 G Temperature (K) G a mm a - R a y L u m i n o s i t y ( e r g / s ) L = 2, =1.5, B = 5 GL = 2, =1.5, B = 1 G
Figure 18.
Probability distribution function (see also Figure 13) of cosmic ray energy density (left) and histogram of gamma-rayluminosity (right) at varying resolutions of 2 pc (solid lines), 4 pc (dashed lines), 8 pc (dot-dashed lines), and 32 pc (dotted lines).For either magnetic field strength tested, lowering resolution leads to higher cosmic ray occupancy in cold gas, especially atinterface temperatures T ≈ K. This leads to slightly more gamma-ray emission coming from cold clumps. This is especiallynoticeable in the B = 1 µ G simulations, though the majority of emission still comes from the T ≈ × K band. Note thechange in x-axis limits in the right panel to focus on the lower temperature region where most gamma-ray emission comes from. dense gas. This trend is evident in both the B = 5 µ G and B = 1 µ G simulations. This leads to more gamma-rayemission coming from the dense, cold gas, especially in the interface temperature range near 8 × K.The total gamma-ray luminosities over all temperature bins are relatively unchanged when the resolution is varied,though, and even appear to decrease in the B = 5 µ G case with 32 pc resolution. This is consistent, also, with ourfindings of § α parameter. In Figure 14, one can seethat lower α , akin to lower resolution that would smooth out density perturbations, actually leads to a decrease incollisional losses. Poor resolution, then, appears not to be the dominant factor that grossly overestimates gamma-rayemission in large-scale simulations, but more high-resolution simulations of cosmic rays in different environments (e.g.superbubbles vs mean ISM) and for a range of mean gas densities are needed..2. Validity of a Fluid Cosmic Ray Model
A possible limitation of this work concerns the validity of a fluid model for cosmic rays. This assumption restsupon the cosmic ray mean free path λ mfp being shorter than other length scales of interest, namely the cosmic rayscale length l CR and is, in fact, baked into the diffusivity we’ve implemented: Assuming a short mean free path andbalancing scattering against acceleration down the pressure gradient leads to a scattering rate ν ∼ c / ( l CR v s ), where v s is the streaming velocity at which wave damping balances wave growth. This implies that the ratio of mean freepath to scale length λ mfp /l CR ∼ v s /c <
1. So even at cloud interfaces where the diffusivity is initially very large andthe scale length becomes very short as the cosmic ray front approaches, λ mfp < l CR by design.In Figure 19, we plot the mean free path λ mfp = κ/c versus l cr = P CR / ∇ P CR at a time of 28 Myrs when the cosmicray front is just hitting the cloud we presented in Section 5.1. At every snapshot we output to data, the mean freepath is below the cosmic ray scale length, but clearly, as in this snapshot, there are locations upstream of the cosmicray front where the mean free path is greater than the scale length a few cells away from it..3. Cosmic Ray Exclusion via Collisional Losses
One of the main findings of our work is that the interface between hot and cold gas plays a crucial role in cosmicray propagation and dynamical effects on the cloud. Our result is broadly consistent with a growing literature oncosmic ray penetration into molecular clouds, the history of which is quite rich (see e.g. Skilling & Strong (1976) andCesarsky & Volk (1978) for seminal papers) and outlined in Everett & Zweibel (2011). Specifically, recent works havenoted that cosmic ray penetration into molecular clouds is highly nonlinear and self-modulated by the excitation ofwaves in the cloud interface (see the Discussion section). While we have focused on a cosmic ray front moving witha preferred direction toward a cloud, it may be that molecular clouds sit in a “sea” of cosmic rays – an initially flatcosmic ray energy profile. In this case, with no pre-existing pressure gradient, the cosmic ray profile is modulated bycollisional loss rates in the dense molecular cloud. The resulting pressure difference sucks fresh cosmic rays into the8
Bustard et al. x (kpc) m f p ( p c ) v st = v ionA + Ion-Neutral Damping t = 28 Myr l c r ( p c ) f minion = 10 f minion = 10 Figure 19.
The cosmic ray mean free path λ mfp = κ || /c vs the cosmic ray scale length l CR for the 1D simulations of § cloud, and how much those exterior cosmic rays can pressurize the cloud interior depends on the transport throughthe cloud interface.We model this scenario for two different clouds each of radius 50 pc, with varying interface widths of 25, 5, and 0.5pc, and with densities of 10 and 100 cm − . Initially, the cosmic ray energy density is uniform everywhere, but it dropswithin the cloud as collisions decrease the cosmic ray pressure. Figure 20 shows our results at snapshots when thecosmic ray energy density has “bottomed out”, i.e. when a steady-state has been reached between cosmic ray inflowthrough the interface and energy loss in the interior. The cloud morphology changes somewhat due to the pressureimbalance, but the cosmic ray pressure we start with is at least a few times smaller than the thermal pressure, so theeffect is small. For the fiducial cloud density of 10 cm − (top row) with minimum ion fractions of 10 − and 10 − , thecosmic ray energy stays fairly level. As expected, the interface is the biggest bottleneck for inflowing cosmic rays, soa thicker interface leads to a larger drop in cosmic ray energy within the cloud.For a denser cloud (100 cm − ), which is the fiducial value in Everett & Zweibel (2011), the energy drop is moresubstantial. There are also more apparent differences when we start with a lower vs higher cosmic ray energy. Thelower cosmic ray energy leads to more decoupling and higher diffusive flux. This is especially true when there is a thininterface. Without a bottleneck, the pressure gradient doesn’t steepen and generate confining waves, so the diffusiveflux remains large.This study addresses the assertion in Everett & Zweibel (2011) that cosmic ray energy will be roughly uniform insideand outside cold clouds and shows that this very much depends on interface properties. It also shows that the extentto which cosmic rays penetrate clouds and the relative importance of diffusive and advective fluxes depends somewhaton cosmic ray content already in the cloud. This is different from the single cosmic ray front simulations we havefocused on in this paper, where the presence of a steep front very much negates the diffusive flux and locks cosmic raysto Alfv´en waves. In the now uniform pressure setup, the collisional loss time is not necessarily short enough to cause asteep pressure gradient, especially when the interface is narrow, and the diffusive flux can be significant. Future workwill address more realistic scenarios with multiple cosmic ray bursts from different directions, as well as cosmic rayproduction within cold clouds. REFERENCES Ackermann, M., Albert, A., Atwood, W. B., et al. 2016,A&A, 586, A71Ajello, M., Di Mauro, M., Paliya, V. S., & Garrappa, S.2020, ApJ, 894, 88Amato, E., & Blasi, P. 2018, Advances in Space Research,62, 2731 Banda-Barrag´an, W., Br¨uggen, M., Heesen, V., et al. 2020,arXiv e-prints, arXiv:2011.05240Becker Tjus, J., & Merten, L. 2020, PhR, 872, 1Blasi, P., & Amato, E. 2019, PhRvL, 122, 051101Boettcher, E., Zweibel, E. G., Yoast-Hull, T. M., &Gallagher, J. S., I. 2013, ApJ, 779, 12 ow Cosmic Rays Navigate the Multiphase ISM -100 pc 0 pc 100 pc10 D e n s i t y ( g / c m ) f minion = 10 E c ( e r g c m ) -100 pc 0 pc 100 pc10 D e n s i t y ( g / c m ) f minion = 10 E c ( e r g c m ) -100 pc 0 pc 100 pc10 D e n s i t y ( g / c m ) f minion = 10 E c ( e r g c m ) Figure 20.