Cosmic-Ray Acceleration of Galactic Outflows in Multiphase Gas

We investigate the dynamical interaction between cosmic rays (CRs) and the multiphase interstellar medium (ISM) using numerical magnetohydrodynamic (MHD) simulations with a two-moment CR solver and TIGRESS simulations of star-forming galactic disks. We previously studied transport of CRs within TIGRESS outputs using a"post-processing"approach, and we now assess the effects of the MHD backreaction to CR pressure. We confirm our previous conclusion that there are three quite different regimes of CR transport in multiphase ISM gas, while also finding that simulations with"live MHD"predict a smoother CR pressure distribution. The CR pressure near the midplane is comparable to other pressure components in the gas, but the scale height of CRs is far larger. Next, with a goal of understanding the role of CRs in driving galactic outflows, we conduct a set of controlled simulations of the extraplanar region above $z=500$ pc, with imposed boundary conditions flowing from the midplane into this region. We explore a range of thermal and kinematic properties for the injected thermal gas, encompassing both hot, fast-moving outflows, and cooler, slower-moving outflows. The boundary conditions for CR energy density and flux are scaled from the supernova rate in the underlying TIGRESS model. Our simulations reveal that CRs efficiently accelerate extra-planar material if the latter is mostly warm/warm-hot gas, in which CRs stream at the Alfv\'en speed and the effective sound speed increases as density decreases. In contrast, CRs have very little effect on fast, hot outflows where the Alfv\'en speed is small, even when the injected CR momentum flux exceeds the injected MHD momentum flux.

1. INTRODUCTION Cosmic rays (CRs) are charged particles moving with relativistic speeds.While their origins are believed to be localized to supernova remnants where they are accelerated in shocks (e.g., Bell 2004;Morlino & Caprioli 2012;Blasi 2013), CRs spread throughout the interstellar medium (ISM) thanks to their quasi-collisionless nature.Direct observations of CRs in the solar system indicate that their spectrum extends over more than ten orders of magnitude in kinetic energy, from at least 10 6 eV to 10 20 eV.For the protons, which comprise most of the CR energy, the spectrum is well approximated by a broken power law that peaks at kinetic energies near 1 GeV (e.g., Grenier et al. 2015;Cummings et al. 2016).
In the Milky Way, the energy density of CRs in the local ISM is comparable to the thermal, turbulent, and magnetic energy densities (e.g., Boulares & Cox 1990;Beck 2001), suggesting that CRs can significantly contribute to the dynamics of gas in the ISM.
Explaining the multiphase nature of galactic outflows is a longstanding theoretical issue (see review by Heckman & Thompson 2017).While the hot component is interpreted as supernova-heated gas that accelerates under its own thermal pressure gradients out of the host galaxy, the presence of a high-velocity cooler component is still puzzling.One possible explanation is that the warm/cold outflowing gas consists of ISM material that is accelerated via ram pressure from the surrounding hot wind.A caveat for this mechanism is that the hot wind can rapidly destroy the embedded cool gas due to a combination of shocks and hydrodynamical instabilities (e.g., Cooper et al. 2009;Scannapieco & Brüggen 2015;Schneider & Robertson 2017;Zhang et al. 2017).However, idealized simulations of cool clouds travelling thorough a hotter medium have shown that, under appropriate cloud and environmental conditions, cool gas can survive its journey by mixing with the surrounding hot gas: for sufficiently high gas density, mixing between the two gas phases reduces the cooling time of the hot gas, triggering its condensation (e.g., Armillotta et al. 2016;Gronke & Oh 2020;Sparre et al. 2020;Banda-Barragán et al. 2021).Schneider et al. (2020) explored this possibility by means of high-resolution simulations of starburst dwarf galaxies, finding that mixing between cool and hot gas and subsequent cooling is an effective way of transferring momentum from the hot to the cool phase of the wind, accelerating the latter to very high velocity.At the same time as more massive cool clouds are accelerated by mixing, lower mass clouds can be destroyed and add to the mass flux of the hot wind (see e.g.Fielding & Bryan 2022, and references therein).
In addition to gaining momentum from the hot wind, cool gas can in principle be accelerated by interactions with CRs.The dynamical role of CRs in driving galactic winds has been investigated in both analytic models (Dorfi & Breitschwerdt 2012;Recchia et al. 2016;Mao & Ostriker 2018;Quataert et al. 2021a,b) and numerical simulations of isolated galaxies or cosmological zoom-ins (e.g., Uhlig et al. 2012;Salem & Bryan 2014;Pakmor et al. 2016;Ruszkowski et al. 2017;Chan et al. 2021;Girichidis et al. 2022;Peschken et al. 2023;Thomas et al. 2023) and portions of ISM (e.g., Girichidis et al. 2016;Simpson et al. 2016;Farber et al. 2018;Girichidis et al. 2018;Rathjen et al. 2021).All these studies found that CR pressure gradients can drive galactic outflows; however, the efficiency of this process is strongly dependent on the CR transport prescription adopted in the model because this affects the CR pressure gradient and therefore the momentum transfer to thermal gas.
Modeling CR transport on galactic scale is a challenging task due to the complex physical processes that couple CRs to the thermal gas, which are not yet fully understood (see review by Amato & Blasi 2018).The interaction between CRs and thermal gas is primarily mediated by magnetic fields.CRs stream along magnetic field lines, while scattering off small-scale (of order the CR gyroradius) magnetic fluctuations.This scattering process couples CRs with the gas and sets the effective diffusive propagation speed of the CR distribution.Current understanding suggests that the origin of waves is energy dependent (Blasi et al. 2012;Zweibel 2017;Evoli et al. 2018), with CRs at moderate (GeV) and low (sub-GeV) energy scattered by Alfvén waves excited by the CRs themselves via streaming instability (the self-confinement scenario; e.g., Kulsrud & Pearce 1969;Wentzel 1974), while high-energy (ultra-GeV) CRs are scattered by extrinsic turbulent fluctuations cascading down to small scales (the extrinsic turbulence scenario; e.g., Chandran 2000;Yan & Lazarian 2002).Nevertheless, to date these ideas have not yet been fully tested with realistic physical models of the underlying multiphase and magnetized ISM.
In most studies of ISM dynamics, the CR mean free path is much smaller than the spatial scales of interest, and it is reasonable to treat CRs as a fluid.Since the CR-wave interaction is not resolved on macroscopic scales, this interaction is treated in the transport equations via a scattering (or diffusion) coefficient.Most commonly, a constant diffusion coefficient which ignores the multiphase structure of the underlying gas is adopted, with a value that is motivated by observational constraints.Considering the vastly varying conditions within the multiphase ISM, however, a constant diffusion coefficient is unrealistic.To address this issue, we recently developed a more detailed physical prescription for the transport of the CR fluid, in which the scattering coefficient varies with the properties of the ambient gas, with a functional form motivated by the theory of self-confinement (Armillotta et al. 2021, hereafter Paper I).Our model focuses on GeV CRs as they contain most of the energy and momentum of the CR population and are therefore more relevant for the gas dynamics.
According to the self-confinement paradigm, a CR distribution with a bulk drift speed greater than the Alfvèn speed can excite Alfvèn waves through gyro-resonance, and individual CRs scatter off these waves as the CR distribution drifts in the direction of decreasing CR density.The processes of wave excitation through the streaming instability and isotropization via scattering off the waves so produced have been validated directly via MHD-PIC simulations (e.g., Bai et al. 2019).In principle, scattering by resonant Alfvén waves can prevent the CR drift (or "streaming") speed from exceeding the local Alfvén speed if wave amplitudes are sufficiently large.However, wave amplitudes and therefore scattering rates are reduced by wave damping, and this damping varies with the local properties of the ISM (e.g., Kulsrud 2005;Plotnikov et al. 2021;Bambic et al. 2021).
In line with the self-confinement picture, our model treats CR fluid transport as a combination of advection by the thermal gas, streaming along the magnetic field at the local ion Alfvèn speed, and diffusion relative to the wave frame.The CR scattering coefficient is determined by the local balance between wave excitation and damping mediated by local gas properties (considering both ion-neutral damping and non-linear Landau damping, Kulsrud & Cesarsky 1971;Kulsrud 2005).Paper I describes the incorporation of this prescription in the algorithm for CR transport implemented by Jiang & Oh (2018) in the Athena++ MHD code (Stone et al. 2020), as well as our model for ionization of warm and cold gas as set by CRs in the tens of MeV regime.In Paper I and in Armillotta et al. (2022, hereafter Paper II) we present the application of our model to computing the propagation of CRs in the TIGRESS MHD simulations modeling kpc-sized portions of star-forming galactic disks for a range of conditions (Kim & Ostriker 2017;Kim et al. 2020a).The advantage of the TIGRESS simulations is that star formation and feedback (including both supernovae and photoelectric heating) are modeled in a self-consistent manner, thus providing a realistic representation of the multiphase star-forming ISM.Our work demonstrates that the transport of CRs is quite different in different thermal phases of the gas, with the CR scattering coefficient varying over more than four orders of magnitude depending on the properties of the underlying gas (e.g., density and ionization fraction).This challenges the common assumption of uniform scattering, and highlights the importance of an accurate representation of the multiphase ISM in CR transport modeling.
One limitation of our previous studies is that they took a postprocessing approach instead of selfconsistently computing the simultaneous evolution of thermal gas, magnetic fields, and CRs.This prevented us from examining the effects of the MHD backreaction to the CR pressure.In the present work, we move beyond the postprocessing approach and investigate the impact of CRs on multiphase ISM dynamics, with a particular focus on the role of CRs in accelerating galactic outflows.We pursue our investigation in two steps.First, starting from the TIGRESS solar-neighborhood snapshots postprocessed for CR transport in Paper II, we perform new simulations in which MHD and CRs are evolved together.These simulations allow us to analyse how the transport and distribution of CRs in realistic multiphase gas are affected by "live" MHD, thereby testing the conclusions of our previous studies.
In the second part of this paper, we perform MHD simulations of galactic outflows including CRs.For the MHD variables, the initial conditions are obtained from the original TIGRESS solar-neighborhood simulation outputs, with a domain extracted from the extra-planar region (|z| > 500 pc).During the simulation, thermal gas and CRs are injected as boundary conditions at the bottom of the simulation box.We explore different thermal and kinematic conditions of the injected gas, cov-ering the typical properties of hot and fast-moving outflows, as well as cool and slow-moving outflows, both of which are expected to be driven by supernovae into the extra-planar region (Kim & Ostriker 2018;Vijayan et al. 2020).Our goal is to understand differences in the ways CR pressure gradients act to accelerate gas within different phases; the extra-planar regions are where this is expected to occur, since the CR scale height is larger than that of the gas.
The paper is organized as follows.In Section 2, we briefly summarize our numerical algorithms.In Section 3, we present the results of the simulations with self-consistent MHD and CRs that use the TIGRESS snapshots postprocessed in Paper II as initial conditions.In Section 4, we describe the setup of the galactic outflow simulations; outcomes of these models for hot wind and warm winds are presented in Section 5 and Section 6, respectively.In Section 7, we discuss our work in relation to observational findings and other recent theoretical works.Finally, in Section 8, we summarize our main results.

MHD EQUATIONS WITH CR TRANSPORT
In this work, we use the MHD code Athena++ (Stone et al. 2020) including the two-moment algorithm for CR transport developed by Jiang & Oh (2018).The full set of ideal MHD equations including CR transport is (2) Here, ρ is the gas density, v is the gas velocity, B is the magnetic field, e = (1/2)ρv 2 + P t /(γ − 1) + B 2 /(8π) is the gas energy density with P t the gas thermal pressure and γ = 5/3 the gas adiabatic index, e c is the CR energy density, F c is the CR energy flux, and P ↔ c is the CR pressure tensor.We assume approximately isotropic pressure, so that P ↔ c ≡ P c I ↔ , with P c = (γ c −1) e c = e c /3, where γ c = 4/3 is the adiabatic index of the relativistic fluid, and I ↔ is the identity tensor.In Equation 6, the speed v m represents the maximum velocity for CRs propagation.In principle, v m is equal to the speed of light c for relativistic CRs.However, here we adopt the "reduced speed of light" approach (see e.g.Skinner & Ostriker 2013, for the corresponding two-moment radiation implementation) with v m = 10 4 km s −1 ≪ c, as it has been demonstrated that the result is not sensitive to the exact value of v m as long as v m is much larger than any other speed in the simulation (Jiang & Oh 2018).This enables larger numerical timesteps based on the CFL condition for this set of hyperbolic equations.
In Equation 2-3, Φ is the "external" gravitational potential from the old stellar disk and dark matter halo (Kim & Ostriker 2017, Equation 6).Self-gravity is not included in the simulations performed for this paper.ρL = n H (n H Λ(T ) − Γ) is the net cooling function, where n H is the hydrogen number density.The cooling coefficient Λ(T ) is computed using the fitting formula in Koyama & Inutsuka (2002) for T < 2 × 10 4 , and the tabulated values in Sutherland & Dopita (1993) with solar metallicity for T > 2 × 10 4 .For warm and cold gas we apply heating to represent the photoelectric effect on grains.The heating rate Γ scales with the instantaneous far-ultraviolet luminosity from star particles.For reference solar-neighborhood values, we adopt a heating rate of Γ 0 = 2 × 10 −26 erg s −1 , and a mean FUV intensity of 4πJ FUV,0 = 2.7 × 10 −3 erg s −1 cm −2 (see Kim et al. 2020a for further details).
In Equation 2 and Equation 6, the term σ represents the rate of momentum density exchanged between the CR population and the thermal gas.In the energy equations, Equation 3 and Equation represents the direct CR pressure work done on or by the gas, while represents the rate of energy transferred to the gas via damping of gyro-resonant Alfvén waves.In the above, v s is the CR streaming velocity, defined to have the same magnitude as the local Alfvén speed in the ions v A,i ≡ B/ √ 4πρ i , oriented along the local magnetic field and pointing down the CR pressure gradient.Here, ρ i is the ion mass density, where the fractional ionization is calculated assuming collisional ionization equilibrium for T > 2 × 10 4 and assuming CR ionization for T < 2 × 10 4 (see Section 2.2.5 of Paper I).The diagonal tensor σ ↔ tot is the wave-particle interaction coefficient, defined to allow for both scattering and streaming along the direction parallel to the magnetic field, and only scattering in the directions perpendicular to the magnetic field, For the relativistic case, σ ∥ = ν ∥ /c 2 , where ν ∥ is the scattering rate parallel to the magnetic field direction due to Alfvén waves that are resonant with the CR gyro-motion.In the simulation, we compute σ ∥ in a self-consistent manner (see below) and set σ ⊥ = 10 σ ∥ .The latter can be understood as scattering by unresolved fluctuations of the mean magnetic field.In Paper I, we explored the transport of CRs in the absence of perpendicular scattering (σ ⊥ ≫ σ ∥ ) as well as the case σ ⊥ = 10 σ ∥ , and did not found any substantial difference in the CR distribution.
The scattering coefficient parallel to the magnetic field direction σ ∥ is derived based on the predictions of the self-confinement picture and assuming that, in steady state, the excitation of Alfvén waves by streaming CRs is balanced by some form of wave damping (Kulsrud & Pearce 1969;Kulsrud & Cesarsky 1971).We consider nonlinear Landau damping (NLL) and ion-neutral damping (IN).The former occurs when thermal ions have a resonance with the beat wave formed by the interaction of two resonant Alfvèn waves, while the latter arises from friction between ions and neutrals in partially ionized gas, where neutrals are not tied to magnetic fields.We refer to Section 2.2.3 in Paper I and Section 2.2.2 in Paper II for a detailed description of the derivation of σ ∥ .Here, we only report the final solution.The scattering coefficient σ ∥ reduces to in well ionized, low-density gas where nonlinear Landau dominates, and in primarily neutral, denser gas where ion-neutral damping dominates.In Equation 10, Ω 0 = e|B|/(m p c) is the cyclotron frequency for e the electron charge and m p the proton mass, v t,i is the ion thermal velocity (which we set equal to the gas sound speed c s ), m i is the ion mass, n i is the ion number density, v p = 1 − (m p c2 /E) 2 is the CR proton velocity, where m p is the proton mass and E ≡ E k + m p c 2 is the total relativistic energy, with E k the kinetic energy.For CRs with E k ≃ 1 GeV, which are the focus of this work, v p ≈ c.The quantity n 1 , which has units of a density, depends on the CR distribution function F (p) in momentum space as where p 1 = m p Ω 0 /k is the resonant momentum for wavenumber k.In Appendix A1 of Paper I, we demonstrate that, assuming F (p) ∝ p −4.7 (e.g., Aguilar et al. 2014Aguilar et al. , 2015)), Finally, in Equation 5-6, the terms Λ coll n H e c and Λ coll n H F c /v 2 p represent, respectively, the rates of CR energy density and CR flux energy lost due to collisional interactions with the ambient gas.GeV CRs interact with the ambient gas through ionization of the neutral atomic/molecular gas and hadronic collisions leading to decays of pions into γ-rays (e.g., Padovani et al. 2020).In the above, Λ coll is the CR collisional coefficient.In Paper I, we estimate that Λ coll is equal to 4 × 10 −16 cm 3 s −1 for CRs with E k = 1 GeV (see also Padovani et al. 2020).

FROM POSTPROCESSED TO SELF-CONSISTENT SIMULATIONS
With the goal of incorporating the MHD backreaction to the CR force, in this paper we move from the postprocessing approach employed in Paper I and Paper II to a self-consistent approach computing the MHD together with the CRs.Among the three TIGRESS models investigated in Paper II, here we focus on the model representative of the solar neighborhood environment only.
In the postprocessing simulations, the MHD quantities are frozen in time, while the energy and flux density of CRs are evolved through space and time given the background distribution of thermal gas, magnetic field, and star cluster particles extracted from the TIGRESS simulation snapshots.During the post-processing simulations, CR energy is injected near star cluster particles to model effects from supernova events.For each star cluster particle, the rate of injected energy is computed as where ϵ c is the fraction of supernova energy that goes into production of CRs, assumed to be equal to 0.1 (e.g.Morlino & Caprioli 2012;Ackermann et al. 2014), E SN = 10 51 erg is the energy released by an individual supernova event, and ṄSN is the number of supernovae per unit time determined from the Starburst99 code based on the age and mass of the star cluster.We note that the current postprocessing simulation approach differs from that in Paper I in that we omit the term ) in the energy equation (Equation 5); we explain the reason for this modification in Section 3.1.
Figure 1 displays the distribution on the grid of several MHD and CR quantities in one sample postprocessed TIGRESS snapshot.The simulation has box size L x = L y = 1024 pc and L z = 7168 pc, with a uniform spatial resolution ∆x = 8 pc, sufficient to achieve robust convergence of several ISM and outflow properties, as well as convergence of CR properties (Kim & Ostriker 2017;Armillotta et al. 2021).The first four columns of Figure 1 show slices at y = 0 (upper panels) and z = 0 (lower panels) of hydrogen number density n H , gas temperature T , magnitude of gas velocity v, and magnitude of ion Alfvén speed v A,i .This particular snapshot is representative of an outflowing phase of the original TIGRESS simulation, in which part of the gas heated and accelerated by supernova blast waves breaks out of the galactic plane, generating large-scale outflows in the coronal region.These outflows are composed of warm-cold fountain gas reaching several kpc from the midplane, and hot winds that escape from the disk (Kim & Ostriker 2018;Vijayan et al. 2020;Kim et al. 2020a).As a consequence, most of the computational volume is occupied by hot (T > 10 6 K) and rarefied gas, with a decrease in the hot-gas volume filling factor near the midplane, mostly composed of warm/cold (T ≲ 10 4 K) gas.We note that the gas velocity v exceeds the ion Alfvén speed v A,i in the hot phase of the gas, while v A,i exceeds v in the warm phase (especially in the neutral gas where x i ≪ 1).
The three rightmost panels of Figure 1 display some outputs of the CR transport algorithm: scattering coefficient σ ∥ , CR pressure 1 P c /k B , with k B the Boltzmann constant, and magnitude of CR flux in the vertical direction, F c,z .The scattering coefficient distribution closely follows the distribution of the background MHD quantities.In particular, σ ∥ is relatively high (above 10 −28 cm −2 s) in hot, highly ionized regions and quite low (below 10 −29 cm −2 s) in cooler, neutral regions.Intermediate-density regions at the interface between neutral and ionized gas are characterized by the highest values of σ ∥ (∼ 10 −27 cm −2 s).Similarly to the scattering coefficient distribution, the distribution of CR pressure also reflects the gas distribution: CRs accumulate in high-density regions, where the relatively-low gas velocities (v < 50 km s −1 ) do not foster their removal, while CRs in regions with hot and fast-moving winds (v ≫ 100 km s −1 ) are rapidly advected away from the mid-plane.One can note that the CR-flux streamlines mostly align with the velocity streamlines in regions with hot winds, meaning that CRs coupled to the hot gas escape the disk through these "chimneys".
Starting from the CR-postprocessed TIGRESS simulation snapshots (as in Figure 1), we now perform new simulations where the distributions of thermal gas, magnetic field, and CRs are evolved together (see Section 2).While these new simulations properly compute the MHD together with the CRs, they are not fully self-consistent in that they do not include self-gravity to follow new star formation, and they do not include injection of energy and momentum in the thermal gas from radiation and supernova feedback.Given that new star formation and feedback to create hot gas are not included in the current simulations, we run only for a timescale (a few Myr) shorter than the time for the hot gas to be advected out of the domain.We turn off CR energy injection for these "MHD relaxation" restart simulations.Our primary goal with these simulations is to allow for the MHD backreaction on the CRs, and to test how this modifies the CR distribution and resulting transport characteristics.
Figure 2 displays how the MHD and CR quantities plotted in Figure 1 have evolved in 3 Myr.A visual comparison clearly shows that the distribution of CR pressure becomes much smoother once the MHD backreaction is included.In the next section, we investigate what leads to the redistribution of CR pressure.

CR pressure redistribution
In Figure 3, we zoom in on a filament of warm (T ≃ 10 4 K) and dense gas located above the disk region to show time evolution of the spatial distributions of gas temperature T , hydrogen number density n H , CR pressure P c , gas velocity magnitude |v|, and ion Alfvén speed |v A,i |, as well as the orientation of the velocity and magnetic field.At t = 0, CRs are trapped in the dense gas.In the postprocessing simulation, the expansion of the hot and fast-moving gas into the warm gas leads to the transport of CRs towards the dense filament (see velocity streamlines in the hot gas).Once in the filament -where CR streaming dominates over advection (v A,i > v) -CRs are unable to escape due to the magnetic field lines mostly parallel to the filament edge.This CR trapping results in a significant CR pressure gradient at the filament edges.At t > 0, the backreaction of the CR pressure on the gas rearranges the velocity and magnetic field topology, allowing CRs to propagate out of the dense gas.The velocity field direction and the evolution of the temperature and density distributions clearly indicate that the warm and dense gas is now expanding into the hot gas as a consequence of the CR pressure gradients initially present at the edge of the filament.
A more quantitative analysis of the evolution of the velocity and magnetic field topology at interfaces between cold/warm and hot gas is presented in Figure 4. Here, the left and right panels display the probability distributions of B • ∇P c /|∇P c | and v • ∇P c /|∇P c |, respectively.These are, respectively, the cosine of the angle between the magnetic field direction and the CR pressure gradient direction, and the cosine of the angle between the velocity direction and the CR pressure gradient direction.In order to track the interface regions, the distributions are weighted by the normalized density gradient |∇ρ|/ρ.At t = 0, the distribution of B • ∇P c /|∇P c | in peaks at values near zero, meaning that the CR pressure gradient is mostly perpendicular to the magnetic field direction.As discussed above, this prevents CRs from streaming out of the dense gas in the postprocessing simulations.Once the backreaction of the CR pressure on the gas is included, however, the angle between the CR pressure gradient direction and the magnetic field direction decreases, as indicated by the slight shift of the probability distribution towards larger absolute values of B • ∇P c /|∇P c |.
The CR backreaction is even more evident in the right panel of Figure 4, showing how the distribution of the angles between the velocity direction and the CR pressure gradient direction varies with time.At t = 0, the distribution of v • ∇P c /|∇P c | peaks at values near 1.As shown in Figure 3, CR pressure gradients are initially aligned with the velocity streamlines at interfaces, and pointing in the same direction.Due to the excess of CR pressure in the cold/warm dense gas, when MHD becomes "live", CRs exert a back-reaction force that leads the dense gas to expand into the hot gas.As a result, the distribution of v • ∇P c /|∇P c | in Figure 4 peaks at values near −1 at t = 0.5 Myr, meaning that the velocity field direction becomes mostly anti-parallel to the CR pressure gradient (see also panels at t = 0.5 Myr in Figure 3).We conclude that the rearrangement of velocity and magnetic field topology enable CRs initially trapped in the dense gas to propagate out of it, leading to a more uniform CR pressure distribution (see Fig- Figure 5 quantitatively shows the redistribution of CR pressure in the different phases of the gas once the backreaction of CRs is included.The left and right plots display the volume-weighted probability distributions of CR pressure vs. gas density at t = 0 and t = 3 Myr, respectively.Overall, P c increases in the low-density hot gas: at t = 0, the distribution of CR pressure spans across more than two orders of magnitude, from P c /k B ∼ 10 4 K cm −3 to as low as P c /k B ∼ 10 1.5 K cm −3 ; by contrast, at t = 3 Myr, the  distribution is much narrower and P c /k B consistently remains above 10 2.5 K cm −3 .While increasing in the hot and rarefied gas, P c decreases in the cooler high-density gas: for n H ≳ 10 −2 cm −3 , P c decreases by a factor 1.5 from t = 0 to t = 3 Myr.Finally, we analyse how the mean vertical profiles of CR and MHD pressures vary in "live MHD" simulations.We first apply post-processing and then run with "live MHD" for ten TIGRESS solar neighborhood snapshots taken between t = 200 to 550 Myr, when the mean surface density is equal to Σ gas = 9.5 M ⊙ pc −2 and the mean star formation rate per unit area is In order to study mean trends, all averages are taken over a given time in the 10 simulations initiated from these 10 TIGRESS snapshots.
The left and right plots of Figure 6 display the horizontally and temporally averaged profiles of CR pressure P c , thermal pressure P t , vertical kinetic pressure P k,z , and and vertical magnetic stress P m,z , at t = 0 and t = 3 Myr, respectively.The vertical kinetic pressure and magnetic stress are defined as These are the ẑ ẑ components of the Reynolds and Maxwell stress tensors, respectively, which appear in the conservation-law form of the momentum equation.
With shearing-periodic boundary conditions, contributions to the mean vertical force per unit volume are obtained from the vertical gradient of horizontal averages over the domain of these terms (e.g., Piontek & Ostriker 2007).P m,z can in principle be negative if the vertical magnetic field dominates, but in practice, the horizontal field dominates in the midplane region where magnetic stresses are dynamically important.For z ≲ 500 pc, the horizontal average of P m,z is positive and is smaller than the horizontal average of the magnetic pressure )/8π by a factor 1.3 − 1.5.At large |z|, the magnetic stress is dynamically unimportant compared to other terms, but magnetic fields still mediate the interaction between CRs and the thermal gas.
Near the midplane, both the CR pressure and the MHD pressures slightly decrease (by 20% and 30% at the midplane, respectively) going from t = 0 to t = 3 Myr.From a detailed examination of the simulations, we find that the decrease of MHD pressures is mostly in the hot gas, as the latter does not gain new energy and momentum in the absence of supernova feedback.By contrast, the decrease of CR pressure is mostly in the warm/cold and dense gas, due to CRs propagating out of it (see above).At higher |z|, P k,z increases in time (from t = 0 to t = 3 Myr) as the vertical CR pressure-gradient forces accelerate the gas away from the disk.We shall explore the dynamical effect of CRs in the extra-planar region in the next sections.
By comparing the left panel of Figure 6 with the right panel of Figure 14 in Paper I or Figure 4 in Paper II, one can notice that while in the present work the midplane CR pressure is in equipartition with the thermal and kinetic pressures, in Paper I and Paper II the midplane CR pressure was a factor of ∼ 2 higher than the other pressures.This is because in the postprocessing simulations carried out for this paper, the adiabatic work term on the RHS of Equation 5 In the present post-processing simulations, we omit this adiabatic work term because it can introduce spurious sources of CR energy at interfaces between cold/warm and hot gas when there is no CR backreaction on the gas.To understand this effect, first note that assuming negligible time-dependent term and collisional loss term in Equation 6, σ Second, at hot/cool interfaces, in post-processing simulations the velocity field and the CR pressure gradient are generally in the same direction (since CRs are flowing from hot gas into cooler gas; see Figure 3 and Fig ) near these interfaces is positive, and overall the adiabatic work done at these interfaces, where |∇P c | is large, can be considerable.
In Paper I and Paper II, we limited the contribution from adiabatic work by setting −v• σ ↔ tot •(F c −4/3ve c ) = 0 in cells with e c > e k .Even with this limit, the total adiabatic work (integrated over the whole simulation domain) was comparable to the total injected CR energy (see Table 2 of Paper II).As noted above, however, as soon as we turn on CR backreaction in "live MHD" simulations, the velocity vectors reorient and the adiabatic work becomes initially negative, and then within a few Myr the velocity field becomes preferentially perpendicular to ∇P c (Figure 4).In "live MHD" simulations, CRs lose energy overall by doing work on the gas.By setting the adiabatic work term to zero to avoid unphysical enhancement of CR energy in postprocessing simulations, for a 10% energy injection rate from SNe to CRs, the CR midplane pressure is comparable to the other MHD pressures, similar to the observed situation.

Role of Streaming, Diffusive and Advective
Transport In this section, we verify that the main findings of our previous work regarding transport still hold in the "live MHD" simulations.Those are 1) the scattering coefficient varies over more than four orders of magnitude depending on the properties of the background gas, with the maximum value reached at intermediate densities, and 2) CR advection dominates in the high-velocity, low-density hot gas, while CR diffusion and streaming are more important in higher-density, cooler gas.We emphasize that the results for scattering coefficients we report on here are valid for CRs with E kin = 1 GeV only.In fact, σ ∥ depends on E kin (and p) through the n 1 term (see Equation 12): n 1 decreases with increasing E kin (see Appendix A.1 in Paper I), which implies the values of σ ∥ as produced by the streaming instability would be lower for CRs with E kin > 1 GeV.
In Figure 7, we show the variation of σ ∥ with density n H (left) and temperature T (right) at t = 0 and t = 3 Myr The analysis is based on all the simulations performed in this section (same set of simulations used for Figure 6); we show median values as well as 16th and 84th percentiles.The profiles predicted by the postprocessing simulations and by the "live MHD" simulations are overall similar, thus confirming our find- ing from Paper I. In the left panel, σ ∥ is fairly flat at low densities up to n H ∼ 10 −2 cm −3 , where the gas is well ionized and nonlinear Landau damping dominates, and then rapidly decreases at higher densities, where the gas is mostly neutral and ion-neutral damping becomes stronger than nonlinear Landau damping (see Paper I for a detailed analysis of the dependence of σ ∥ on n H ).
In the right panel, at low temperatures (T < 10 4 K) where gas is mainly neutral and ion-neutral damping dominates, σ ∥ < 10 −31 cm −2 s and increases mildly with T (as n H decreases). Near T ≈ 10 4 K, where gas becomes mostly ionized (Sutherland & Dopita 1993), σ ∥ abruptly increases by more than four orders of magnitude to a few ×10 −28 cm −2 s; in this high temperature regime where nonlinear Landau damping dominates, σ ∥ slowly decreases with T .Despite the overall similar profiles, σ ∥ is smaller at t = 3 Myr than at t = 0.This difference can be attributed to the smaller CR pressure gradients in the "live MHD" simulations.In the low-density, hightemperature regime, where 10), σ ∥ decreases by a factor ∼ 2 under "live" MHD.In the high-density, low-temperature regime, where Equation 11), σ ∥ decreases by almost one order of magnitude under "live" MHD.As σ IN,∥ decreases with density faster than σ NLL,∥ , the regime where ion-neutral damping dominates extends down to lower densities at t = 3 Myr compared to t = 0.This results in slightly different turnovers of the scattering coefficient profiles: this is at n H ≃ 1 − 2 × 10 −2 cm −3 at t = 0, and n H ≃ 5 − 6 × 10 −3 cm −3 at t = 3 Myr.
In Figure 8, we show the relative contributions of advection, streaming, and diffusion to the overall CR transport.If we define the vertical flux as F c,z ≡ v c,eff,z (4/3)e c , each of these contributions is characterized by a speed, where the first two are the vertical components of the gas flow speed and ion Alfvén speed.The diffusive flux F d is obtained from Equation 6 in the limit of negligible time-dependent term and collisional losses by subtracting the advective and streaming fluxes from the total CR flux: F d ≡ F c − 4/3(v + v s )e c .Using Equation 8, the diffusive velocity becomes The left and right panels of Figure 8 show the median profiles of the vertical components of the advection speed v, ion Alfvén speed v A,i , and diffusive speed v d as a function of gas density and temperature, respectively.The velocity profiles are consistent at t = 0 and t = 3 Myr: the advection speed largely dominates in hot, low-density, well ionized gas, while the diffusive speed largely dominates in cold, high-density, mostly neutral gas.Streaming is the most important transport mechanism in warm and intermediate-density regions.Since the overall transport speed is the largest of the three, and this maximum speed is smallest for the streaming-dominated regime, the primary overall limit on CR transport is from the ion Alfvén speed in gas at n H ∼ 0.01 − 0.1 cm −3 with T ∼ 10 4 K.We emphasize that our conclusion regarding the overall limit on transport set by v A,i holds only for CRs with E kin ∼ 1 GeV.At higher CR energy where the scattering rate is expected to decrease, the diffusion velocity (Equation 16) would increase, pushing the blue curves in Figure 8 upward and squeezing the regime where streaming dominates transport.For CRs of sufficiently high energy, diffusion would dominate in all warm and cold gas.
We point out that even though the transport speed in the colder, higher-density gas is comparable to or larger than the transport speed in the hotter, lower-density gas, the CR residence time is significantly greater in the dense gas, as indicated by the higher mean CR pressure there (see Figure 5).In fact, the propagation of CRs out of the higher-density, primarily neutral gas (which comprises most of the ISM mass) is limited by the low transport speed in the surrounding intermediate-density streaming-dominated gas.Given the low-scattering (but high-column) midplane and high-scattering (but low column) exterior illustrated by the slice depicted in Figure 2, realistic galaxies have physical elements that connect to both traditional "leaky box" and diffusive pictures.
As the distribution of CR pressure becomes smoother from t = 0 to t = 3 Myr, the diffusive speed decreases in the low-density, high-temperature regime, where v d ∝ √ ∇P c /e c .By contrast, the diffusive speed increases in the high-density, low-temperature regime, where v d ∝ 1/e c , as the energy density decreases in the dense gas (see above).However, these changes do not significantly affect CR transport overall, since the main transport mechanism at low density and high temperature is advection, and since the CR transport is already extremely rapid (producing quite uniform P c ) in dense gas.

SIMULATIONS OF GALACTIC OUTFLOWS WITH COSMIC RAYS
In the second part of this work, we investigate under what physical conditions CRs are able to accelerate galactic outflows.To this end, we use the same code described in Section 2 to perform controlled simulations of galactic outflows with imposed inflow boundary conditions.The initial conditions of these simulations are taken from the original TIGRESS solar neighborhood simulation outputs (Kim et al. 2020a), with a domain extracted from the upper extra-planar region (z > 500 pc).
We select two different TIGRESS snapshots -hereafter called M1 (Model 1) and M2 (Model 2) -from which to extract the initial conditions.M1 is representative of an outflow-dominated period, while M2 is representative of a quiescent period.The total gas surface density and star formation rate (averaged over 40 Myr) for these two snapshots are Σ gas = 9.6 M ⊙ pc −2 and Σ SFR = 1.1 × 10 −2 M ⊙ kpc −2 yr −1 for M1, and Σ gas = 9.9 M ⊙ pc −2 and Σ SFR = 4 × 10 −3 M ⊙ kpc −2 yr −1 for M2.For the two snapshots, the main gas and magnetic field properties, measured at z = 500 pc and averaged in the horizontal directions, are listed in Table 1.Here, we define three different gas phases based on temperature: warm (T < 2 × 10 4 K), warm-hot (2 × 10 4 K < T < 5 × 10 5 K), and hot (T > 5 × 10 5 K) phase.For each phase, the horizontally-averaged quantities are computed as qph (z) = x,y q(x, y, z)Θ ph (T )∆x∆y with ∆x = ∆y = 8 pc the grid resolution, and Θ ph (T ) the top-hat function that returns 1 for gas at temperatures within the temperature range of each phase (ph = warm, warm-hot, or hot) or 0 otherwise.In columns 3-4-5, the mass, MHD momentum, and MHD energy fluxes of the thermal gas and magnetic fields along the vertical direction are respectively defined as In column 6, f A , the area fraction occupied by a given phase, is defined as In columns 7-8-9-10-11, the "typical" quantities q are defined as the horizontally-averaged quantities divided by the fractional area occupied by each gas phase, We note that vz averages over all gas, while vz,out averages just over outflowing (v z > 0) gas.
From the properties listed in the table, it is clear that M1 and M2 differ in terms of thermal and dynamical gas properties.For M1, the hot gas occupies most of the volume in the extra-planar region (f A (z = 500 pc) = 0.74), and is characterized by very high velocities (ṽ z,out (z = 500 pc) = 350 km s −1 ).By contrast, in M2 the warm component fills up a considerable part of the volume (f A (z = 500 pc) = 0.59), and the hot gas is characterized by a relatively low velocity (ṽ z,out (z = 500 pc) = 100 km s −1 ).While the temperature of hot gas is more than a factor of two higher in M1 than in M2, the densities within the hot gas are similar for the two models.
The simulation domain has size L x = L y = 1024 pc in the horizontal directions and extends from z = 0 to z = L z = 3584 pc in the vertical direction.For z ≥ 500 pc, the initial conditions for the thermal gas are extracted from the TIGRESS snapshot; for z < 500 pc, the hydro variables are initialized to spatially uniform values.These are the same values assigned to the hydro variables in the ghost zones at the bottom of the z-axis throughout the entire simulation (see below and Table 2).The initial conditions for the magnetic field are entirely extracted from the TIGRESS snapshot.The CR pressure is initialized to be uniform in the horizontal directions and decreasing in the vertical direction as with z i = 500 pc.Here, the CR scale height H c = 1.5 kpc is determined through a linear fit of lnP c versus z for z ≥ 500 pc in Figure 6.For z < 500 pc, P c is initialized to the normalization factor P c,0 .The CR energy flux is initialized to 0. All the simulations employ periodic boundary conditions in the horizontal directions and open boundary conditions at the top of the vertical direction, while the boundary conditions at the bottom of the vertical direction are set to user-defined values.In doing so, we control the injection of gas and CRs in the simulation box.Table 2 lists the model names and the inflow boundary conditions.Our CR simulation model nomenclature is such that the first two letters refer to the TIGRESS snapshot from which the initial conditions are extracted.Depending on the thermal phase of the injected gas, the set of simulations is divided into two main groups, HW and WW, standing for hot wind and warm wind, respectively.The third letters denote the different CR pressure adopted for the HW models (see column 6) or the different inflow density adopted for the WW models (see column 3).We describe the inflow boundary conditions employed in the HW and WW models in Section 5 and Section 6, respectively.
For convenience in comparing the MHD and CR fluxes, Table 2 includes for each model the mass flux F M,0 , the MHD momentum flux F p,MHD,0 , and the CR momentum flux F p,c,0 injected in the boundary conditions (note that F p,c,0 ≡ P c,0 ).

COSMIC RAYS AND HOT WINDS
We start with the analysis of the HW models where hot gas is injected at the bottom of the simulation box.In the HW models, the inflow boundary conditions of the MHD variables are set to their initial mean values at z = 500 pc in the hot outflowing gas (except for the horizontal components of the gas velocity and magnetic field, which are set to 0).This means that the injected gas is hot and fast-moving in the M1 -HW models (T 0 = 10 7 K, v z,0 = 350 km s −1 ), and mildly hot and relatively slow in the M2 -HW models (T 0 = 5 × 10 6 K, v z,0 = 110 km s −1 ).
For both the M1 -HW and the M2 -HW cases, we run three different simulations, characterized by different conditions for CRs (Table 2).The NoCR models are control models that do not include CR physics, while the LPCR and HPCR models include CR physics, but differ in the inflow boundary conditions for the CR pressure, P c,0 , and the vertical CR energy flux, F c,z,0 (the horizontal components of the CR energy flux are set to 0).
In the LPCR models (where LPCR stands for low CR pressure), F c,z,0 is set to the total CR energy flux injected from supernova events in the associated TIGRESS snapshot: F c,z,0 ≡ F c,inj = Nsp sp=1 Ėc,sp /(L x L y ), with the sum taken over all the star cluster particles N sp , and Ėc,sp the rate of CR energy injected from a given star cluster particle (Equation 13).Since the gas advection speed in hot gas dominates over the other CR propagation speeds (see v z,0 vs. v Ai,z,0 in Table 2), we compute2 P c,0 as F c,z,0 /(4v z,0 ).This choice assumes that 1) once injected, CRs do not experience additional gains/losses of energy during their propagation, and 2) CRs propagate at the same velocity of the hot gas.The first assumption is generally valid in the region at z ≲ 500 pc, where streaming and collisional losses are smaller than the injected energy density.However, in the disk region, which contains mostly cold/warm gas, the transport of CRs is slower than would be the case if all the gas were hot and high velocity, so the second assumption provides a lower limit to the CR pressure at the disk-halo interface (z ∼ 500 pc).
In the HPCR models (where HPCR stands for high CR pressure), P c,0 is set to the mean value of CR pressure (in the total gas) measured at z = 500 pc in the corresponding TIGRESS snapshot postprocessed with CRs, while F c,z,0 = 4P c,0 v z,0 .This is an upper limit to the actual CR pressure in the hot gas at z = 500 pc because hot gas tends to have somewhat lower CR pressure than warm gas (and therefore, lower than the average) even in simulations with time-dependent MHD and CR physics, where the distribution of CR pressure is smoother across different gas phases (see Figure 2).
We note that the injected MHD momentum flux exceeds the injected CR momentum flux in all models except for M2 -HW -HPCR.In particular, in M1 -HW -LPCR, the injected MHD momentum flux is more than one order of magnitude larger than the injected CR momentum flux (see Table 2).2. The initial conditions for the MHD variables are extracted from the M1 and M2 TIGRESS snapshots, respectively.From top to bottom, the panels show the slices through y = 0 of hydrogen number density nH, gas temperature T , gas vertical velocity vz, vertical magnetic field Bz, and CR pressure Pc.
Figure 10 displays the distribution on the grid of hydrogen number density n H , gas temperature T , vertical velocity v z , and CR pressure P c in three different models, M1 -HW -LPCR, M1 -HW -HPCR, and M2 -HW -LPCR, at t = 70 Myr.We note that compared to the M1 simulations, the M2 simulations are characterized by a slower injected wind (see Table 2 and Figure 9).
Despite their difference in CR pressure (the CR pressure imposed in the boundary conditions differs by more than one order of magnitude; see Table 2), the two M1 -HW models exhibit overall similar gas distribution: most of the volume is occupied by hot (T > 10 6 K), rarefied and fast-moving gas, which surrounds higherdensity cooler clouds/filaments with both positive and negative velocities.In model M2 -HW -LPCR, the gas distribution is quite different, with most of the volume occupied by gas with temperatures between 10 4 − 10 6 K ("warm and warm-hot").Within this gas, the cooler and denser component flows inward, while the hotter and more rarefied component flows outward.We emphasize that only hot gas is injected in these models.Therefore, the warm/warm-hot gas present in the simulation domain is either warm/warm-hot gas which was already present at t = 0 and has not escaped the box, or injected gas that has cooled down.We estimate the net amount of hot gas that has cooled down throughout the simulation by adding the mass of warm/warm-hot gas that has left the simulation box to the difference between the final (t = 100 Myr; end of the simulation) and the initial (t = 0) mass of warm/warm-hot gas within the simulation domain.We find that the amount of cooled-down hot gas is about 30% of the final mass of warm/warmhot gas in the M1 models, and less than 20% in the M2 model, implying that most of the warm/warm-hot gas present in the simulation domain was already present at t = 0.The M1 and M2 models also differ in terms of CR pressure distribution.Unlike M2 -HW -LPCR (which has hot gas inflowing at 110 km s −1 ), neither M1 -HW -LPCR nor M1 -HW -HPCR (which both have hot gas inflowing at 350 km s −1 ) present a clear decrease of CR pressure with increasing z.Rather, both of the highvelocity HW models exhibit an upward trend in CR pressure with increasing z.Since the large-scale vertical CR pressure gradient in the M1 -HW models points away from rather than toward the disk midplane, we expect that CRs will be ineffective in accelerating gas outward.
Figure 11 shows, for all models, the time evolution of the outward (v z > 0) mass flux F M,out (defined as in Equation 18, but just considering cells with v z > 0) for each thermal phase (separated into three columns) taken at z = 3 kpc.We show the M1 -HW and M2 -HW models in the top and bottom panels, respectively.For all phases, the evolution of the outward mass flux is similar in the three M1 -HW models, regardless of the presence or absence of CRs.This implies that the contribution of CRs to the acceleration of the ambient gas is negligible in these models.Interestingly, the model with higher CR pressure (M1 -HW -HPCR) exhibits the lowest outward mass flux at later times, especially in the hottest component of the gas.This makes the time-averaged (average over t = 20 − 100 Myr) mass flux (dotted line) lower in the M1 -HW model.By contrast, the mass flux evolution is significantly different among the three M2 -HW models, which have a factor ∼ 3 lower hot gas injection velocity than the M1 -HW models (110 km s −1 vs. 350 km s −1 ).The model without CRs (M2 -HW -NoCR) is characterized by the lowest mass flux, indicating that CRs are dynamically important in these models.Figure 11 shows that, at earlier times, for the M2 models, F M,out is about one order of magnitude higher in the model with higher CR pressure (M2 -HW -HPCR) than in the model with lower CR pressure (M2 -HW -LPCR) for all thermal phases.However, while F M,out remains roughly constant with time in M2 -HW -LPCR, it decreases with time in M2 -HW -HPCR, becoming comparable to M2 -HW -LPCR towards the end of the simulation.
In the following sections, we quantify to what extent the dynamical impact of CRs on the gas varies across different models, exploring the reasons behind the different trends observed in Figure 11.

Momentum exchange between gas and cosmic rays
In this section, we measure the exchange of vertical momentum flux between gas and CRs in different models.We start from the momentum equation in the vertical direction (Equation 2), considering a shearingperiodic box and taking horizontal and temporal aver- ages: with ⟨q⟩ the average over time of q(z; t), the horizontal average of a quantity q.In Equation 23, we assume that the time-dependent term and the loss term of Equation 6are on average negligible, so that ⟨σ 23, F p,MHD (z) ≡ ⟨P k,z + P t + P m,z ⟩ is the contribution to the momentum flux from the MHD pressures (see Equation 19), while F p,c (z) ≡ ⟨P c ⟩ is the contribution to the momentum flux from the CR pressure.We may integrate Equation 23 from an initial height z i to an arbitrary height z, and express the momentum equation in terms of momentum flux differences along the vertical direction (see Kim & Ostriker 2018;Vijayan et al. 2020;Armillotta et al. 2022): Here, ⟨ ṗz ⟩(z) is the volume-integrated rate of change in z-momentum normalized to the area of the horizontal plane, and the momentum flux and weight differences are defined as for W(z) the weight of gas above z, with L z /2 the top of the simulation box.Note that −∆ z W(z) > 0 is the weight of gas in the range [z i , z].
In steady state, ⟨ ṗz ⟩(z) ≈ 0, and Equation 24 reduces to: this says that in steady state, the momentum flux differences ∆ z F p,MHD (z) and ∆ z F p,c (z) together must balance the weight of gas in a given range of z.We can rearrange Equation 28 so that the LHS represents the "net" momentum flux difference of the thermal gas that arises from its interaction with CRs: We note that in the absence of CRs (∆ z F p,c (z) = 0), since the weight −∆ z W(z) > 0 (for z > z i ), the momentum flux of thermal gas must always decrease as the flow moves outward (∆ z F p,MHD (z) < 0) to compensate for losses in climbing out of the gravitational potential.
In the presence of CRs, ∆ z (F p,MHD (z) − W(z)) may be either positive or negative depending on whether CRs give momentum to (∆ z F p,c (z) < 0) or receive momentum from (∆ z F p,c (z) > 0) the thermal gas.
Figure 12 shows the time-averaged vertical profiles of the LHS and RHS terms of Equation 29 as measured in the three hot wind models displayed in Figure 10, M1 -HW -LPCR, M1 -HW -HPCR, and M2 -HW -LPCR.The momentum flux and weight differences are calculated from an initial height z i = 500 pc (Equation 26), and the temporal average is over t = 20 − 100 Myr.In Figure 12, we also plot ⟨ ṗz ⟩ to show that this term is negligible and the system is in quasisteady state in the three displayed models.We omit M2 -HW -HPCR from this analysis as this model does not exhibit a steady behaviour, that is |⟨ ṗz ⟩| ≫ 0.
In both of the M1 simulations, ∆ z F p,c is positive and increases towards larger z, meaning that CRs gain momentum from the gas (∆ z (F p,MHD (z) − W(z)) < 0) while moving outward.Over the height range we consider here, CRs gain only a few % of the initial (i.e., at z = 500pc) MHD momentum flux in the M1 model with lower CR pressure (M1 -HW -LPCR), and about 20% in the M1 model with higher CR pressure (M1 -HW -HPCR); for M1 -HW -HPCR this amounts to a 25% increase in the CR pressure.As a consequence, the MHD momentum flux decreases faster in the latter model, thus explaining why the outward mass flux is lower here than in M1 -HW -LPCR (see Figure 11).In the M2 simulation, the exchange of momentum flux between gas and CRs goes in the opposite direction: ∆ z F p,c is negative and decreases towards larger z as a consequence of CRs transferring momentum to the gas (∆ z (F p,MHD (z) − W(z)) > 0).The change in F p,c from z = 500 pc to z = 3 kpc is about 20% of the z = 500 pc MHD momentum flux and about 30% of the z = 500 pc CR momentum flux.
From this analysis, we conclude that the dynamical impact of CRs on hot outflows varies depending on the conditions of the thermal gas: in the M2 -HW -LPCR model, characterized by a relatively slow hot wind, CRs efficiently transfer momentum flux to the gas.By contrast, in the M1 -HW models, where the velocities of the wind are much higher, momentum is transferred from the thermal gas to the CRs, causing the gas to decelerate.In the next section, we investigate how the different gas properties in the two groups of models control the different contribution of CRs to the overall gas dynamics.

Vertical profiles
In Figure 13, we present the time-averaged vertical profiles of the HW simulations, showing MHD pressure, CR pressure, fractional area occupied by individual thermal phases, outward gas velocity, and gas density in the simulations shown in Figure 12.The top panels of Figure 13 plot the time-averaged vertical profiles of MHD pressure P MHD and CR pressure P c in the total gas.In M1 -HW -LPCR, P c is much lower than P MHD , thus Temporally averaged vertical profiles of CR pressure Pc and MHD pressure PMHD (first row ), fractional areas for each thermal phase f A,ph (second row ), typical velocities ṽout,ph (third row ), and hydrogen number densities n H,out,ph (forth row ) of the outflowing gas (average is over t = 20 − 100 Myr).Different columns refer to different hot wind models: M1 -HW -LPCR (left column), M1 -HW -HPCR (middle column), and M2 -HW -LPCR (right column).In the top panels, the profiles of CR and MHD pressures are indicated with orange and purple lines, respectively.In the central and bottom panels, different colors represent different gas phases: cyan for warm/warm-hot, and crimson for hot.In the second-row panels, different line styles are used to separate outflowing (vz > 0) from inflowing (vz < 0) gas: dotted for outflowing, dashed for inflowing, solid for total.In the bottom-row panels, solid and dotted lines represent the typical (Equation 22) hydrogen densities ñH,out,ph and the horizontally-averaged (Equation 17) hydrogen densities, respectively.
explaining the negligible contribution of CRs to the gas dynamics in this model (see Figure 11 and Figure 12).In M1 -HW -HPCR, the two pressures are comparable.The MHD pressure decreases outward, while the CR pressure increases outward (in agreement with Figure 12, where F p,c = ⟨P c ⟩).With a CR pressure gradient pointing upward along the z axis, CRs exert a negative vertical force (d⟨P c ⟩/dz > 0 in Equation 23) that deprives the surrounding gas of momentum (see analysis in Section 5.1), causing its deceleration (see third row).Finally, in M2 -HW -LPCR, both P MHD and P c decrease with z, with the CR pressure profile flatter than the MHD pressure profile.In this model, CRs exert an overall positive vertical force that increases the momentum of the gas.In the range z = [0, 1.5 kpc], the hot gas accelerates, and beyond this, the velocity decreases slightly.At larger z, most of the volume is filled with warm and warm-hot temperature gas (see also Figure 10).
The different CR distribution in M1 -HW and M2 -HW -LPCR depends on the different thermal and dynamical properties of the gas.In M1 -HW, the fastmoving hot gas helps to push the lower-velocity, warmer gas present in the simulation domain towards higher altitudes.The interaction with a faster-moving fluid shreds the dense cloudlets during their motion, triggering mixing between the two fluids (see Figure 10).This results in a small increase in the mean velocity of warm/warmhot gas with z, although this remains < 100 km s −1 .As shown in the second row of Figure 13, the fractional area occupied by the warm/warm-hot gas increases from about 0.01% to 20% within z = 3 kpc.Since the volume remains predominantly occupied by the hot gas in both M1 -HW models (with f A,w/w−h < 0.2), the CR transport is primarily controlled by hot gas advection.The decrease in the hot gas velocity at higher z in the M1 -HW -HPCR model leads to an increase in CR pressure with z (in the advection-dominated limit, ), and at the same time the positive CR pressure gradient and interaction with warm gas decelerates the hot gas.Because CR transport is slower in warm gas than hot gas, an additional factor contributing to CR pressure buildup at large z is the increase of f A,w/w−h with height.
In M2 -HW -LPCR, the relatively slow hot wind is less able to accelerate the warm and warm-hot gas.Also, the fraction of warm and dense gas initially present in the simulation domain is much higher than in M1. Figure 13 shows that for z ∼ 0.5 − 1 kpc, ∼ 30% of the volume is occupied by inflowing (v z < 0) warm/warmhot temperature gas.CRs stall in the high-density region at low altitudes slightly increasing their pressure.Similarly to the M1 models, f A,w/w−h increases with increasing z in M2 -HW -LPCR.In this model, though, warm/warm-hot gas occupies most of the volume for z > 1 kpc, and therefore it is this material that pri-marily governs the CR transport.This warm/warm-hot gas is mainly outflowing, with an average velocity increasing with z (see the right panel in the third row of Figure 13).This implies that the advection of CRs out of the simulations domain increases with z, thus explaining the decrease of P c with z.We note that the large fraction of warm gas present in M2 -HW -LPCR makes this model akin to the warm wind models discussed in the next section.We therefore refer to Section 6 for a detailed discussion of the role of CRs in accelerating outflows in the warm gas.
The middle panel in the third row of Figure 13 shows that the typical (see Equation 22) outflow velocity of the hot gas ṽout,hot decreases with z in M1 -HW -HPCR, where CRs act to decelerate the wind.Because of this deceleration, despite the identical initial gas conditions, the typical outflow velocities are lower in M1 -HW -HPCR than in M1 -HW -LPCR, where the CR pressure is negligible compared to the other pressures.By contrast, in M2 -HW -LPCR, ṽout,hot increases with z for z < 1.5 kpc, while flattening at higher altitudes.The typical outflow velocity of the warm/warm-hot gas always increases with z, even in the M1 -HW models where CRs decelerate the hot gas.This happens both because the warm/warm-hot gas gains momentum from the hotter and faster-moving gas, and because lowervelocity cloudlets increasingly "drop out" of the flow with height (see also Vijayan et al. 2020;Armillotta et al. 2022).
Finally, the bottom panels of Figure 13 show that, in the M2 -HW -LPCR model, the warm/warm-hot outflowing gas exhibits lower typical densities compared to the M1 -HW models.Part of the reason for the lower densities in the M2 model is that, with a hot gas mass injection rate four times smaller than in the M1 -HW models (see Table 2), the warm/warm-hot gas can expand in area, dropping its density.Furthermore, unlike the M1 -HW models, in the M2 model the horizontally-averaged density of the warm/warmhot gas (dashed line) decreases outward with z due to the inefficiency of the lower-velocity hot gas in sweeping up the warm cloudlets.With the lower-velocity gas dropping out, the higher-velocity warm/warm-hot gas expands in area, further reducing its density.

COSMIC-RAY DRIVEN WARM WINDS
In this section, we analyse the simulations in which slow-moving warm (T = 10 4 K) gas is injected at the bottom of the simulation box (referred to as WW models in Table 2).In the WW models, the injected gas is warm (T 0 = 10 4 K) and slow-moving (v z,0 ∼ 10 − 20 km s −1 ).The inflow boundary conditions for the vertical components of gas velocity and magnetic field are set to their initial mean values in the warm outflowing gas at z = 500 pc.As in the HW-LCPR models, we fix F c,z,0 to the total CR energy flux injected in the associated postprocessed TIGRESS snapshot.P c,0 is then com- puted as F c,z,0 /[4(v z,0 + v Ai,z,0 )], assuming that both streaming and advection contribute to the propagation of CRs in the warm gas (diffusion is negligible at the gas densities explored here, except for the high-density HD model).We note that for the WW models, P c,0 is comparable to the average value of CR pressure measured in the postprocessed snapshot at z = 500 pc.Hence, there is no need for a second "HPCR" model defined using the postprocessed CR pressure for P c,0 , as in the HW models.
For the WW models, we explore the effects of three different densities for the injected gas: n H,0 = 10 −1 , 10 −2 , 10 −3 cm −3 for models HD, ID, and LD, respectively, which stand for high, intermediate, and low density.The HD models have density closest to that of the warm fountain gas in the TIGRESS simulation (see Table 1).For each set of M1 and M2 simulations, we employ the same value of P c,0 .This is computed for v Ai,z,0 from the ID model, rather than for the value of v Ai,z,0 of the corresponding model (in that case, P c,0 would vary by less than a factor 1.5).In the M2 -WW models, the injected CR pressure is almost a factor of two lower compared to the M1 -WW models.We note that the injected CR momentum flux is larger than the MHD momentum flux for all models except for M2 -WW -HD.We find that our conclusions regarding the effect of varying the injected density in WW models are independent of the initial conditions employed in the simulation (M1 vs M2).Hence, in the following, we mostly focus on the analysis of the M1 simulations.
Figure 14 displays the distribution on the grid of hydrogen number density n H , gas temperature T , vertical velocity v z , and CR pressure P c , in the three M1 -WW models at t = 100 Myr.Unlike the HW models, P c now decreases outward along the vertical direction in all models, suggesting a non negligible contribution of CRs to counteracting gravity and accelerating gas in this extraplanar region.A comparison between the three panels shows that the large-scale vertical gradient of CR pressure becomes smaller going from the model with higher injected gas density (M1 -WW -HD) to the model with lower injected density (M1 -WW -LD).At z > 1 kpc, the gas density and velocity distributions are qualitatively similar in the three models: lower-density gas flows outward, while higher-density gas mostly flows inward.The outflowing gas is clearly accelerated by CR pressure forces as thermal pressure forces are negligible in the warm gas populating most of the simulation box.A difference that emerges from a visual comparison is that the fraction of gas at T ≳ 10 5 K (color-coded in yellow and red in the second row of Figure 14) increases going from M1 -WW -HD to M1 -WW -ID to M1 -WW -LD.At the time the snapshots are taken (t = 100 Myr), the fast-moving hot gas initially present in the simulation domain has already escaped the box.As we shall discuss in Section 6.4, gas at T > 10 4 K currently present in the simulation box is low-density gas heated up by CRs.
In the top panels of Figure 15, we plot the time evolution of the outward (v z > 0) and inward (v z < 0) mass flux F M (Equation 18 (Equation 19), and MHD energy flux F E,MHD (Equation 20) taken at z = 3 kpc in the three models.Unlike our approach in Figure 11, here we do not separate the fluxes of the three different thermal phases since the fraction of hot gas is negligible (in terms of both volume and mass) in the WW models.We do, however, separate outward and inward MHD fluxes based on the sign of v z .In the three models shown, F M , F p,MHD , and F E,MHD decrease after ∼ 10 Myr, as the fast-moving hot gas initially present in the simulation domain escapes the box.Overall, the fluxes at z = 3 kpc are comparable in the three models despite two orders of magnitude difference in the injection density.In all models, the outward fluxes are on average higher than the inward fluxes (see also Section 6.2).We note that, at z = 3 kpc, F M in the M1 -WW models is lower than in the M1 -HW models (see Figure 11; the ratio is ∼ 6 − 10 for time-averaged mass fluxes in HW vs. WW models).This is true even for M1 -WW -HD, which has an order of magnitude higher mass injection rate than the M1 -HW models; the dense inflowing gas in M1 -WW -HD simply builds up at low z (see Figure 14).
The bottom panels of Figure 15 display the time evolution of the CR momentum flux F p,c ≡ ⟨P c ⟩, and CR energy flux F E,c ≡ ⟨F c (z)⟩ at z = 3 kpc.The CR fluxes are separated into outward and inward fluxes based on the sign of the vertical effective CR propagation velocity v c,eff,z ≡ F c /(4P c ).The energy flux of CRs is completely dominated by the outflowing term.Figure 14 also shows that pressure is fairly uniform at any given height.The regions with lower gas density generally correspond to the locations which have v c,eff,z > 0; this is the majority of the area, such that the momentum flux of the outflowing CRs is larger than the momentum flux of the inflowing CRs at z = 3 kpc.In all models, the CR fluxes at z = 3 kpc remain more than one order of magnitude larger than the corresponding MHD fluxes, although they are lower than their injection values (see Table 2).In all panels, the profiles are divided by the MHD momentum flux Fp, MHD at zi = 0.5 kpc.The time average is taken starting from t = 30 Myr to t = 150 Myr.On average, the system is in quasi-steady state (⟨ ṗ⟩ is much smaller than the other terms).

Transfer of momentum to the gas
Similarly to what we have done in Section 5.1 for the HW models, in this section we quantify the exchange of momentum flux between gas and CRs outward along the vertical direction in different models.Figure 16 displays the time-averaged vertical profiles of the "net" MHD momentum flux difference ∆ z (F p,MHD (z)−W(z)), and CR momentum flux difference ∆ z F p,c (z) (see Equation 29) in the three M1 -WW models.The differences are calculated from an initial height z i = 500 pc, and the temporal averages are over t = 50−150 Myr.We discard the first 50 Myr from the analysis to avoid transients affected by the initial conditions, where a fast-moving hot wind is present.In all panels, ⟨ ṗz ⟩ is negligible, meaning that the systems are in statistical steady state.
The main evidence of Figure 16 is that unlike the situation shown for hot winds in Figure 13, ∆ z F p,c is always negative and decreasing towards larger z in the WW models, confirming that CRs effectively transfer momentum to the gas as the latter flows towards higher altitudes.Over the height range we consider here, the gas gains more than 200% of the initial (i.e.z = 500 pc) MHD momentum flux (corresponding to about 50% of the z = 500 pc CR momentum flux) in M1 -WW -HD and M1 -WW -ID, while the gas gains about 100% of the initial MHD momentum flux (corresponding to less than 20% of the z = 500 pc CR momentum flux) in M1 -WW -LD.In the next section, we use additional information from vertical profiles to analyse the main differences, in terms of CR and MHD properties, among the three models.

Vertical profiles
In Figure 17, we analyse the time-averaged vertical profiles of the outward and inward mass and MHD momentum fluxes from the three M1 -WW models.Between ∼ 0.5 − 2.5 kpc, the inward mass and momentum fluxes in all models are fairly similar; this reflects population of inflowing gas present in the initial conditions.In the model with higher gas density (M1 -WW -HD), there is more than an order of magnitude decrease in |F M (v z > 0)| within the first ∼ 0.5 kpc, due to rapid dropout of gas (which is injected with a velocity of only 15 km s −1 ).Beyond this point, |F M (v z > 0)| decreases more slowly with z, since CRs effectively transfer momentum to the gas (see Figure 16).In M1 -WW -ID, |F M (v z > 0)| decreases by a factor of a few within z ∼ 1 kpc, after which it remains flat.In M1 -WW -LD, the outward mass flux |F M (v z > 0)| initially increases at small z, and then plateaus at a value slightly lower than in M1 -WW -ID.The vertical profiles of the MHD momentum fluxes are qualitatively similar to the mass flux profiles.After a steep decrease at z < 0.5 kpc in model M1 -WW -HD and a steep increase in model M1 -WW -LD, the outward MHD momentum flux F p, MHD (v z > 0) remains nearly constant with z, at similar levels for all models.
To understand the warm wind models in more detail, in the left panel of Figure 18, we plot the temporally and horizontally averaged vertical profiles of CR pressures.The CR pressure profile is less steep in M1 -WW -LD than in the other two models, consistent with the smaller fractional exchange in momentum flux between gas and CRs discussed in Section 6.1 (20% vs. 50% of the CR flux at z = 0.5 kpc, where F p,c ≡ ⟨P c ⟩). Models M1 -WW -HD and M1 -WW -ID exhibit similar CR pressure profiles, except for the region at z ≲ 300 pc, where the profile in model M1 -WW -HD is relatively flat due to the effective CR diffusion at high density (n H > 0.1 cm −3 at z ≲ 300 pc; see middle panel of Figure 18).At higher z, CR diffusion becomes negligible as the gas density rapidly drops.
In the middle panels of Figure 18, we display the temporally averaged typical density of the outflowing gas ñH,out as a function of z.In all models, the value of ñH,out at z = 0 is larger than the density of the injected gas n H,0 .Due to the inflow boundary conditions employed at the bottom of the vertical axis, gas flowing inward cannot leave the simulation box, but rather accumulates near the boundary, increasing its density.In M1 -WW -HD, ñH,out has to drop by more than one order of magnitude, to ∼ 10 −2 cm −3 , before CR pressure gradient forces begin to have an effect.This is confirmed by the right panel of Figure 18, showing that the outward vertical velocity ṽz,out initially decreases until z ∼ 0.3 kpc.At higher z, CR forces effectively accelerate the gas, as demonstrated by the vertical velocity increasing with z.In the lower-density models, much less gas builds up, and the flow accelerates outward starting near z = 0.For z ≳ 0.5 kpc, the profiles of CR pressure and gas density in M1 -WW -HD and M1 -WW -ID become roughly similar.This results in similar gas acceleration, as shown in the right panel of Figure 18, where the offset in gas velocity is due to the different velocity near the inflow (at z ≲ 500 pc).The mean outward velocities at z = 3 kpc are ≃ 38 and 45 km s −1 in M1 -WW -HD and M1 -WW -ID, respectively.The outward velocities reach v z,out ≳ 100 km s −1 (well above the averages in Figure 18) in low-density (n H ≲ 10 −3 cm −3 ), high-temperature (T ≫ 10 5 K) regions (see Figure 14).
In M1 -WW -LD, the injected gas has sufficiently low density to be very efficiently accelerated by CR forces.This is clearly visible in the right panel of Figure 18, showing that the outward velocity rapidly increases with z in the region at z ≲ 1 kpc.At high altitudes, the outward gas velocity becomes relatively flat, meaning that CR forces merely balance gravity.Due to the effective acceleration, the outward mass flux does not drop with z as in the other models (see Figure 17), leading to a smaller fractional variation of ñH,out .As we shall see in Section 6.3, streaming dominates the transport of CRs at low altitudes in the WW models.In this regime, . With a smaller fractional drop in gas density, the drop in CR pressure is smaller in the LD model.
In summary, the transfer of momentum flux from the CR population to the gas is overall greater in the models with high/intermediate density at the base of the extra- planar region.For these models, the larger fractional drop in gas density leads to a larger drop in CR pressure, and, as a result, a larger momentum flux transfer.However, the differences in the gas flows at high altitudes (z ∼ 2 − 3 kpc) are in fact not so large among the models: the density and velocity of the outflowing gas are within a factor 2 of each other.Outward mass, momentum, and energy fluxes for the three models are therefore also similar (see Figure 15 and Figure 17).This is true despite two orders of magnitude difference in the injection density.The implication is that mass loading of CR-driven galactic winds adjusts to match the carrying capacity allowed by the momentum flux available for transfer from the CR distribution (see Mao & Ostriker 2018 and further analysis in Section 6.3).

Different cosmic-ray pressures at the wind base
So far, we have analysed the dynamical effect of CRs in models with different density of the injected gas.Here, we compare the results from models M1 -WW -ID and M2 -WW -ID.These two models have the same injected gas density, and nearly the same injected velocity and momentum flux.The injected CR pressure P c,0 and CR flux F c,z,0 are however lower by a factor two in M2 -WW -ID (see Table 2).
The mean vertical profiles of CR pressure, outflowing gas density, and outward gas velocity, are displayed in Figure 19.In M2 -WW -ID, CR forces are ineffective in accelerating the gas near the injection region.As a consequence, the outward vertical velocity decreases until z ∼ 100 pc, and the gas accumulates near the boundary; the resulting gas density is a factor of 2 higher than in model M1 -WW -ID.Away from the injection region, CRs are overall able to accelerate the gas.Even though the CR pressure gradients in the two models are comparable for z ≲ 1 kpc, the gas density drops faster in M2 -WW -ID due to the smaller initial outward velocity.The gas acceleration produced by M2 -WW -ID is initially higher than the acceleration produced by M1 -WW -ID, while becoming comparable for z > 1 kpc.Beyond ∼ 1 kpc, the velocity profiles are similar, although the density and therefore the mass flux remains higher in M1 -WW -ID than in M2 -WW -ID.As we shall see in the next section, the extent to which CRs accelerate the ambient gas depends on the so-called CR sound speed, C c = dP c /dρ.

Comparison with analytic predictions
In this section, we compare the results from our warm wind simulations with the predictions of the onedimensional (1D) steady-state model developed by Mao & Ostriker (2018) to describe properties of warm galactic winds driven by CRs.We focus on model M1 -WW -ID since the density of the injected warm gas in this model is most consistent with the typical gas density in SN-driven fountain flows at the base of the extra-planar region for solar neighborhood conditions (n H ∼ 0.01 − 0.1 cm −3 , Kim & Ostriker 2018).Mao & Ostriker (2018) treat the gas as a single-phase isothermal medium, and assume that velocity, magnetic field and pressure gradient streamlines follow the largescale gravitational potential gradient.The former condition is approximately satisfied in our WW simulations, where most of the gas has a temperature of 10 4−5 K.The latter condition is also a good approximation of the mean trend in our WW simulations, where gravity is in the vertical direction, as are the mean pressure gradients and the mean magnetic field.Although our simulations have outflows and inflows simultaneously, both are primarily in the vertical direction.The Mao & Ostriker (2018) model includes CR advection and streaming only, under the assumption that CR diffusion is negligible in the extra-planar region.This assumption is verified in our simulations.
From the CR energy conservation law, Mao & Ostriker (2018) derive the steady-state expression for the   Mao & Ostriker (2018).The left panel shows the vertical profile of CR pressure Pc from the simulation M1 -WW -ID (red line) compared with the profile predicted by the 1D model (black line).The middle panel shows the gas velocity variation along the z axis dzvz produced by the simulation M1 -WW -ID (red line) in comparison to the 1D model prediction (black line).The third panel shows the vertical profiles of the root mean square gas velocity ⟨v 2 z ⟩ 1/2 (purple line) and the effective sound speed C 2 eff (light green line).The 1D model predicts that, in the Cartesian geometry of the present simulations, C 2 eff has to be larger than v 2 z in order to have a steady-state outflow.
CR pressure along each streamline: Here, the "0" subscript denotes values at the streamline footpoint, γ c = 4/3 is the CR adiabatic index, and A is the wind cross-sectional area.In the local Cartesian box of our models, A does not vary with z, i.e.A(z)/A 0 = 1 for the flow.3Since in our simulations gas and CR motions are mostly along the z-axis, i.e. ⟨v⟩ ≈ ⟨v z ⟩ and ⟨v A,i ⟩ ≈ ⟨v A,i,z ⟩, we can estimate how the mean CR pressure ⟨P c ⟩ should vary along the vertical direction according to Equation 30.In the left panel of Figure 20, we compare the time-averaged vertical profile of CR pressure in M1 -WW -ID with the analytic prediction, and find an excellent agreement for z > 500 pc. Figure 21 shows that the mean ion Alfvèn speed is always larger than the mean gas advection velocity, implying that streaming dominates over advection in the warm extra-planar gas (see also Figure 8 and Section 3.2).We note that ⟨v z ⟩ in Figure 21 is the average between the mean inflowing and outflowing velocities, and is therefore slightly lower than ⟨v z,out ⟩ shown in Figure 18.Following Mao & Ostriker (2018) (see their eqs.26 and 27 for more general 1D geometry, and see also Breitschwerdt et al. 1991;Everett et al. 2008), the steadystate momentum equation along vertical streamlines in a local cartesian box is: where the third term within the brackets is the CR contribution to the squared sound speed, Equation 32 can be easily demonstrated writing v and v A,i as a function of ρ in Equation 30.From the mass conservation law and B z = const.(for ∇ • B = 0), we know that in steady state, the gas velocity and ion Alfvén speed along vertical streamlines can be written as and respectively.In Equation 34, ρ i /ρ i,0 can be approximated to ρ/ρ 0 since gas is mostly ionized (x i ∼ 1) at the typical temperatures reached for extraplanar gas in our simulations (T ≳ 10 4 K).In the streaming-dominated regime (as in Figure 21), P c ∝ ρ 2/3 from Equation 32and Equation 34, so that the CR sound speed increases with decreasing ρ, i.e. dP c /dρ ≈ 2/3 P c /ρ ∝ ρ −1/3 .By contrast, in the advection-dominated regime, the CR sound speed decreases with decreasing ρ, i.e. dP c /dρ ≈ 4/3 P c /ρ ∝ ρ 1/3 .When advection dominates, as is true for fast, hot gas in our HW models, the scaling behavior of the effective sound speed does not encourage wind acceleration.
From Equation 31, the ordinary differential equation describing a steady-state vertical wind is given by where C eff is the effective sound speed, defined as In the middle panel of Figure 20, we compare the average outward velocity variation along the z axis ⟨d z v out,z ⟩ in M1 -WW -ID with the predictions of the analytic model (Equation 35), finding an overall good agreement.The gradient ⟨d z v out,z ⟩ is mostly positive, as the outward velocity increases with z (see Figure 18).From Equation 35, for an accelerating wind with d z v z > 0, it must be true that C 2 eff > v 2 z (d z Φ is positive by definition).The third panel of Figure 20 shows that this condition is satisfied in our simulation.Of the two terms on the RHS of Equation 36, the CR term C 2 c = dP c /dρ dominates over the thermal term c 2 s (c s ∼ 10 km s −1 for T ∼ 10 4 K).Hence, in the streaming-dominated regime of M1 -WW -HD, C 2 eff increases with z as ρ decreases.The condition C 2 eff > v 2 z is satisfied in all WW models (at z ≳ 500 pc), except for M1 -WW -LD.Here, the outward velocity exceeds the effective sound speed at z ≳ 1 kpc.In fact, the wind is not accelerating in this region for this model (see right panel of Figure 18).We note also that for the subsonic regime (C z ), it is straightforward to show that the acceleration obeys v z /v z,0 = ∆Φ/(3C 2 c,0 ) + 1 3 for ∆Φ the potential difference.Even though both density and CR pressure are higher in M1 -WW -ID than in M2 -WW -ID, the velocity and C c are similar in the two models at z ∼ 1 kpc, which explains the similarity of the acceleration profiles in Figure 19.Finally, we note from Equation 35 that there is no steady-state solution for a wind passing through a sonic point (C 2 eff = v 2 z ) in a local cartesian box.This would be possible only if the cross-sectional area A increases as the flow moves outward.With the opening of the streamlines, a term C 2 eff d z ln A is subtracted from the gradient of the gravitational potential so that the numerator of Equation 35 passes through zero and winds go through a critical point to reach an asymptotic velocity (Mao & Ostriker 2018).

Cosmic-ray heating
In addition to their direct dynamical effects, energy loss from the CR bulk flow can also be an important source of heating for the ambient gas (e.g., Jacob & Pfrommer 2017;Kempski & Quataert 2020). 4In the self-confinement regime, damping of streaming-excited Alfvén waves mediates transfer of energy from the CR population to the gas at a rate v s • σ ↔ tot • (F c − 4/3ve c ) (see RHS of Equation 3).If we subtract the sum of the work equation (obtained by the dot product of v and Equation 2) and the magnetic energy equation (obtained by the dot product of B and Equation 4) from the total energy equation (Equation 3), we can see that such energy transfer leads to an increase of the thermal energy of the gas: with e t = P t /(γ − 1) the thermal energy.We note that the last term is always positive as streaming always drains energy from CRs (see definition of v s in Equation 7); with σ ↔ tot • (F c − 4/3ve c ) → −∇P c in steady state, the contribution to the heating rate is |v A,i • ∇P c |.
To evaluate the extent to which CRs affect the thermal state of the gas, we quantify the three source terms on the RHS of Equation 37 in our fiducial model M1 -WW -ID. Figure 22 shows the magnitude of the average rates of radiative cooling n 2 H Λ, photoelectic heating n H Γ, and CR heating v s • σ ↔ tot • (F c − 4/3ve c ) as a function of n H and T .For n H ≳ 3 − 4 × 10 −3 cm −3 and T ≲ 10 4 K, most of the gas is in thermal equilibrium with photoelectric heating balanced by cooling, n H Γ ≃ n 2 H Λ, while CR heating is negligible.At higher temperature and lower density, the rate of photoelectric heating quickly decreases, while the rate of CR heating only slowly decreases with density, and increases at the highest temperatures.In the density range 10 −4 < n H < 10 −3 cm −3 and temperature range 2×10 4 < T < 10 5 K, most of the gas is in thermal equilibrium with CR heating balanced by cooling, H Λ. At higher temperature and lower density, gas becomes thermally unstable with the rate of CR heating overall larger than the rate of radiative cooling.CRs heat the gas up to temperature of order 10 6 K.In Paper I and Paper II, we found that for the TI-GRESS model representing a solar-neighborhood environment, the average CR pressure near the midplane exceeded the other ISM pressures by a factor of 2−3.This can be compared with the local Milky Way, where the estimated pressure of CRs with E kin ≥ 1 GeV is closer to equipartition with the other pressures.We speculated that this disagreement might be due to the absence of CR backreaction on the gas in the purely-postprocessing simulations of Paper I. Our expectation was that, with the rearrangement of the velocity and magnetic field topology in simulations with time-dependent MHD and CR physics, CRs trapped in the ISM would be able to propagate out of the disk, leading to a CR pressure decrease near the midplane.In Section 3 of the present paper, we show that rearrangement of the velocity and magnetic field topology indeed enables CRs trapped in the dense gas to propagate out of it.The decrease of pressure in the warm gas is compensated by the increase of pressure in the hot gas, leading to an overall more uniform pressure between hot and warm/cold gas (Figure 2, Figure 5).There is therefore only a marginal decrease in the mean midplane CR pressure with "live MHD" compared to pure post-processing CR transport (see Figure 6).This suggests that the overall CR confinement in the disk is not controlled by local magnetic geometry, but by the combination of advection in hot gas and Alfvénic streaming in warm gas, as limited by scattering in low density, ionized gas (Figure 8).We find the scattering rate does not significantly change when "live MHD" is turned on (Figure 7).
Based on further tests, we have found that the higherthan-expected CR pressure reported from the postprocessing simulations of Paper I and Paper II was due to spurious adiabatic work exerted from the thermal gas on the CR fluid (see Section 3.1).The effect is particularly strong at interfaces between cold/warm and hot gas, where CR pressure gradients are high and oriented in the same directions as the velocity vectors, which results in significant gains of CR energy density (≈ v • ∇P c ).We find that with the rearrangement of the velocity field topology in self-consistent "live MHD" simulations, the adiabatic work diminishes in magnitude (CR pressure gradients are smaller) and becomes mostly negative (at interfaces the velocity vectors are mostly in the direction of decreasing CR pressure).We conclude that since adiabatic work can significantly boost the effective CR energy injection in postprocessing simulations, it is best to zero out this term, and we do so for the postprocssing simulations of Section 3. We find that, in steady state, the average pressure of CRs near the midplane is in equipartition with the other pressures, and in good agreement with the observed value.
10 5 10 6 T (K) In the ISM, approximate equipartition is attained by balancing energy gains from star formation feedback and energy losses due to dissipative processes (Ostriker et al. 2010;Ostriker & Shetty 2011;Kim & Ostriker 2015a;Ostriker & Kim 2022): balancing far-UV photoelectric heating and cooling for thermal pressure, balancing momentum flux injection from supernovae with turbulent dissipation for kinetic pressure, and applying turbulent driving in combination with shear to sustain the pressure in the magnetic field.Likewise, the pressure of CRs is determined by the balance between CR energy injection from supernovae and propagation of CRs out of the disk (Armillotta et al. 2021(Armillotta et al. , 2022)).Other energy sink terms play a minor role in setting the midplane CR pressure: in the post-processing runs of Section 3, the average rates of collisional and streaming losses relative to the rate of CR energy injection are 0.1 and 0.2, respectively. 5 In the case of negligible losses, the midplane CR pressure P c (0) can be written as a function of the vertical flux of CR energy above the supernova input layer F c,inj (see also Armillotta et al. 2021Armillotta et al. , 2022)): an effective CR scale height (measured in the simulation through a linear fit of lnP c vs. z within 1.5 kpc), and κ eff ≡ σ −1 eff an effective diffusion coefficient.We note that κ eff may be understood as a measurement of the efficiency of CR propagation, including not only CR diffusion, but also advection and streaming.For H c,eff ≃ 1.3 kpc as measured in our post-processing simulations, the midplane pressure-flux relation is satisfied for κ eff ≃ 3.6 × 10 28 cm 2 s −1 .In the selfconsistent simulations, where the distribution of CR pressure is smoother, H c,eff slightly increases, and so does κ eff .After t = 3 Myr, H c,eff ≃ 1.5 kpc, while κ eff ≃ 5.6 × 10 28 cm 2 s −1 .
We emphasize that the value of κ eff is an order of magnitude larger than the value of σ −1 ∥ in the low-density ionized gas that fills space outside of the midplane (see Figure 2 and Figure 7), showing that advection and streaming (rather than pure diffusion) are essential to CR transport.With F c ≡ (4/3)v c,eff e c = 4v c,eff P c , the effective transport speed is v c,eff = (1/4)κ eff /H c,eff , which for the post-processing or self-consistent simulations is v c,eff = 22 or 30 km s −1 , respectively.This range is consistent with the (vertical) ion Alfvén speed in gas at T ≲ 10 4 K, as shown in Figure 8. Indeed, we conclude in Section 3.2 that the primary overall limit on transport of GeV CRs in the multiphase ISM is the ion Alfvén speed in warm gas at n H ∼ 0.01 − 0.1cm −3 , since this regime is where the maximum over the diffusion, advection, and streaming speeds has the smallest value.
In addition to the midplane CR pressure, another observable quantity that can be predicted by our simulations is the CR grammage, which is a measure of the column of gas traversed by CRs during their propagation.
In our simulations, the grammage can be computed as with µ H = 1.4,and Ėinj the total CR energy injected per unit time (see Section 3.2.2 and Equation 27 in Paper I).
In good agreement with observations (X ∼ 10 g cm −2 for GeV energy, Amato & Blasi 2018), we find that the average CR grammage is X ≃ 10.8 g cm −2 in the postprocessing simulations, and X ≃ 7.2 g cm −2 in the simulations with MHD at t = 3 Myr.Provided that the CR scale height is larger than that of the gas, it is straightforward to show that Equation 39 is equivalent to X ≈ (3/8)Σ gas v p /v c,eff for Σ gas the total gas surface density.With an average Σ gas = 9.5 M ⊙ pc −2 in the TIGRESS solar neighborhood simulation, the resulting estimated values of the grammage are in good agreement with computed values when adopting v c,eff listed above.
Based on the success of our simulations in reproducing key observed properties in the local ISM, we conclude that the self-confinement paradigm is able to explain transport of GeV CRs.Provided that the underlying MHD simulation has a realistic representation of the multiphase, star-forming ISM, it does not appear necessary to include sources of small-scale MHD waves other than self-excitation by CRs, or more generally to adopt an ad hoc scattering coefficient.When combined with realistic flow and Alfvén velocities in hot and warm ionized gas, scattering set by the balance between growth from CR streaming instability and nonlinear Landau damping produces realistic transport.We emphasize, however, that the quantitative results presented here apply only to CRs with energy E kin ∼ 1 GeV.Because the scattering coefficient decreases with increasing CR particle momentum (through n 1 , see Equation 12), we expect the CR confinement time within the disk and, as a consequence, the CR grammage (and measured B/C ratio) to decrease with increasing momentum, in agreement with the observations (e.g., Ptuskin et al. 2009).It will be very interesting to extend the current models to higher CR energy in order to quantify the scaling behavior of the grammage, and to identify when self-excited waves fail to produce results consistent with observations.
To date, the only other numerical studies that have tested CR transport with a variable scattering model based on self-confinement are those of Hopkins et al. (2021Hopkins et al. ( , 2022) ) employing cosmological zoom-in FIRE simulations (Hopkins et al. 2018a) with CRs, adopting the same CR injection efficiency ϵ c = 0.1 that we do.In contrast to our conclusions, they suggest that the standard self-confinement model is incompatible with observations.In their simulations, CRs are excessively confined, with grammage, residence time, and CR pressure up to an order of magnitude higher than observations.They argue that the failure is the consequence of selfconfinement instability: if P c rises, σ ∥ increases, which slows down CR diffusion, further increasing P c and σ ∥ in a runaway process that stops when the effective propagation speed decreases to the local ion Alfvèn speed.While we agree with their argument that streaming dominates diffusion in low-density regions, and we find that streaming is indeed the main limit in some regimes, this is not the only factor controlling CR transport.In our simulations, advection is the dominant transport mechanism in hot, rarefied gas at T ≳ 2 × 10 4 K and n H ≲ 10 −2 cm −3 , while waves are strongly damped by ion-neutral collisions so that diffusion is dominant at n H ≳ 0.5 cm −3 and T ≲ 5 × 10 3 K (see Figure 8).Thus, in supernova remnants where CRs are injected, and in superbubbles and chimneys that channel hot gas from the disk to the halo, gas velocities are so high that CRs are quickly advected away from the midplane.Furthermore, the high diffusion in the predominantly-neutral gas prevents CRs from being trapped there.As a result, the typical CR pressure we find is compatible with ISM observations in the Milky Way.
One possible reason for the difference between our results and those of Hopkins et al. (2021Hopkins et al. ( , 2022) ) is that the mass resolution adopted in cosmological zoom-in simulations (∆M ∼ 10 3 − 10 4 M ⊙ ) is insufficient to properly resolve the hot phase of the ISM.In our localdisk simulations, the mass resolution at the typical densities (n H ∼ 10 −4 − 10 −3 cm −3 ) of the hot ISM is ∆M = ρ∆x 3 ∼ 10 −3 − 10 −2 M ⊙ , far lower than the value 7000 M ⊙ employed by Hopkins et al. (2021) in their Milky Way-like model.Since the Sedov stage of SNR evolution ends when remnant mass is ∼ 1000 M ⊙ (Kim & Ostriker 2015b), it is not possible to capture the Sedov stage of evolution of supernova remnants with mass elements 7000 M ⊙ , and therefore supernova feedback in FIRE is realized via momentum injection (Hopkins et al. 2018b,a) without resolving the hot ISM.Interestingly, Hopkins et al. (2021) remark that the densities of their superbubble regions are ∼ 0.01 cm −3 (one or two orders of magnitude denser than in our simulations), and that these regions often have very low effective diffusion.Though they do not report the advection velocities within superbubbles, it is possible that streaming limits CR transport there, unlike the case for our higher resolution models.Hopkins et al. (2022) explore a variety of alternatives to standard self-confinement, but our conclusion is that at least for GeV CRs, it appears that the standard model is satisfactory.Crucially, this requires that the underlying MHD simulation of the multiphase ISM is sufficiently accurate, which demands much higher resolution than is possible in any present cosmological simulations (or zoom simulations of massive galaxies) to follow the creation and propagation of hot gas.Potentially, high resolution MHD + CR simulations, using TIGRESS or similar frameworks, can be employed to calibrate effective CR scattering rates for use in lower resolution cosmological simulations.

Cosmic-ray driven outflows
As reviewed in Section 1, the last decade has seen an increasing number of studies, both analytic and numerical, exploring the role of CRs in accelerating largescale galactic outflows.However, the predictions of these studies are highly sensitive to the physics of CR transport included in the model.
A key distinction is between the case when CR transport is dominated by diffusion, and when it is dominated by streaming.Quataert et al. (2021a,b) have used a combination of analytic estimates and time-dependent numerical simulations to study the properties of galactic winds driven in each of these limits, considering spherically symmetric geometry as would be appropriate for a wind emerging from the center of a galaxy.Their diffusive solutions primarily focus cases with κ ≡ 1/σ ∥ > r 0 c i , where r 0 is the base radius of the wind and c i is the isothermal gas sound speed, and they argue that when this inequality holds, diffusive CR transport yields mass-loss rates and wind powers higher than streaming transport.When κ ≡ 1/σ ∥ < r 0 c i , they instead find much reduced wind fluxes for the diffusive case.For r 0 of a few kpc in a galactic disk, and c i ∼ 10 km s −1 for warm ionized gas, r 0 c i ∼ 10 28 cm 2 s −1 .This can be compared to the results shown in Figure 7 for σ ∥ , based on our self-consistent computation of the scattering coefficient assuming that waves are excited by the streaming instability.There, for typical densities n H ≃ 10 −3 − 10 −2 cm −3 (see Figure 18 from our warm wind simulations), σ −1 ∥ ∼ few × 10 −27 cm 2 s −1 .This, together with our previous results in Armillotta et al. (2022) showing that σ −1 ∥ is even smaller in inner-galaxy environments, suggests that the large-diffusion regime κ > r 0 c i would not be applicable.Although quantitative results would presumably vary somewhat with environment, Figure 8 shows that in fact only relatively dense and cool gas is diffusion dominated, with hot gas advection-dominated and warm, intermediatedensity gas streaming-dominated.Note that these statements refer to the ∼ GeV protons that comprise the bulk of the CR energy density; high energy CRs have lower scattering rates and would therefore be diffusion dominated down to lower density (i.e. the blue curves in Figure 8 would shift up), but they could not power significant winds.
In our warm wind simulations, diffusion is minimal, and the primary limit on CR transport is initially the Alfvén speed, with advection increasing as gas accelerates outward.We show in Section 6.3 that the flow profiles from the simulations (e.g.profiles of CR pressure and acceleration rate) are in very good agreement with the predictions of the 1D analytic model of Mao & Ostriker (2018), which includes CR advection and streaming only, and assumes an isothermal flow.Mao & Ostriker (2018) proposed (see their Eq.46) that the mass-loss rate per unit area in a CR-driven warm wind is expected to scale as for P c,0 the CR pressure at the base of the wind, V H the characteristic halo velocity, and u 0 the velocity of warm gas in the fountain flow in the region where the wind originates.The limited spatial domain in our present simulations precludes testing whether this predicted scaling holds, since the CR-driven wind does not pass through a critical point.The acceleration of warm gas by CRs within our domain is, however, comparable to that achieved by transfer of momentum from hot gas to warm cloudlets when a powerful hot wind is present (compare velocities profiles in Figure 13, Figure 18, Figure 19).Also, our warm wind simulations show that the fluxes carried by the CR-accelerated wind are very insensitive to the density of the gas flowing into the extraplanar region: Figure 17 shows only a factor of a few variation in the wind fluxes among models where the inflow density varies by two orders of magnitude.This suggests that the mass flux of a CR-driven warm wind may be set by the velocity u 0 of warm gas in the fountain region, consistent with Equation 40, rather than the mean density ρ 0 of the fountain.An alternative way of writing the predicted CR-driven wind mass is While several recent studies have focused on the role of CR pressure in launching warm outflows, Everett et al. (2008) presented a semi-analytic 1D model of hot outflows driven by both thermal and CR pressure, aimed at explaining the diffuse soft X-ray emission observed at the Galactic Center.Everett et al. (2008) concluded that CR pressure is crucial to driving a hot wind from the Milky Way's Center, and that cases with CR pressure comparable to thermal gas pressure produce the best fit to the observations.Our simulation that is most similar to those of Everett et al. (2008) is M1 -HW -HPCR, in which most of the volume is occupied by hot, fast-moving gas, and CR and gas pressures are comparable (see Figure 13).However, our conclusion is that CRs do not contribute to the hot wind acceleration, but rather slow down the gas.
We identify two main differences between our hot wind simulations and the models of Everett et al. (2008).First, they treat the gas as a single-phase medium, while gas is multiphase in our simulation.Even though subdominant in terms of volume, the increasing fraction of warm, slow-moving gas cloudlets at high altitude will help to slow down CRs and in doing so help to create a CR gradient pointing upward.Furthermore, while CR transport is advection-dominated in our simulation (since advection dominates in the volume-filling hot gas; see Figure 8), it is streaming-dominated in Everett et al.'s best model (see their Fig. 8 for z < 3 kpc).This is because the (imposed) magnetic field at the wind base in their model is two orders of magnitude higher than the average value for the hot wind of our simulation (see Table 1).Although we reach different conclusions regarding the role of CRs in hot winds, we agree with Everett et al. that streaming transport is crucial to accelerating gas.Indeed, in Section 6, we demonstrate that CRs can efficiently accelerate flows of gas in the warm medium, where CR transport is streamingdominated.
Several other recent studies (e.g., Girichidis et al. 2016;Simpson et al. 2016;Farber et al. 2018;Girichidis et al. 2018;Rathjen et al. 2021) have investigated CR driving of outflows using simulations of galactic disk patches, similar in geometry and numerical resolution to our models.Similar to our conclusions, they find that CRs are able to drive outflows of warm gas that reach moderate velocity (v ≲ 100 km s −1 ) within a few kpc, and which are also somewhat denser than hot winds driven by SNe.Of course, densities in both hot SNdriven winds and warm CR-driven winds vary (in space and time) depending on the input energy and momentum fluxes driving them, and hot winds also help in accelerating warm embedded clouds.Comparing our M1 -HW models with our M1 -WW models, which have comparable hot-gas and CR pressure injected (see Table 2), the density of warm gas driven by CRs is a factor of a few higher in the M1 -WW models than the density of hot gas in the M1 -HW models (compare center panel of Figure 18 with bottom panels of Figure 13).However, the density and mass flux of warm gas driven by interaction with the hot gas in the M1 -HW models is higher than the density and mass flux of warm gas driven by CRs in the M1 -WW models (mass fluxes of warm gas for the HW and WW cases are shown in top-left panels of Figure 11 and Figure 15, respectively).We point out that this difference in density and, most importantly, in pressure (given that the temperature is ∼ 10 4 K in both cases) could be used to distinguish between warm outflows generated by interaction with hot gas and those driven by CRs in observational studies.
An important difference between our warm wind simulations and the above disk-patch simulations is that the latter do not include CR streaming, only diffusion.We find that streaming dominates CR transport in our WW simulations (and in warm ionized gas generally), and, in doing so, it controls the dynamical and thermal evolution of the outflow.The absence of gas heating via CR streaming might explain why, in the above simulations, most of the outflowing gas has a temperature of 10 4 K (see e.g., Simpson et al. 2016;Farber et al. 2018;Girichidis et al. 2018), while, in our WW simulations, a significant fraction of gas is heated up to temperature T > 10 5 K.
In work to date, CR streaming has been included only in larger-scale simulations of isolated galaxies or cosmological zoom-ins only (e.g., Ruszkowski et al. 2017;Chan et al. 2019;Hopkins et al. 2020Hopkins et al. , 2021;;Thomas et al. 2023), which have much lower resolution than our simulations.While most of these works acknowledge the importance of CR streaming in accelerating galactic outflows, they do not discuss the impact of CR-streaming heating on the thermal state of the gas.Ji et al. (2020) analyse the simulation of a Milky Way-like galaxy performed by Hopkins et al. (2020), finding that heating from CRs via streaming losses is much weaker than radiative cooling.Unlike us, they do not compare heating and cooling rates in different gas phases, but rather they compute their average values as a function of the distance from the galactic plane.Similarly, in Section 6.4, we find that on average radiative cooling dominates over heating from CRs (see Figure 22).However, gas cooling rates become much lower than CR heating rates in the low-density gas.
The right panel of Figure 22 shows that CR heating exceeds cooling gas when T ≳ 10 5 K, which is consistent with the prediction for the onset of thermal instability in Kempski & Quataert (2020, see their Eq. 47).Their prediction is that under conditions where the CR pressure is large compared to the gas pressure (true in the extraplanar region -see e.g. Figure 6 or the F p panels of Figure 15), the onset of thermal instability is where dΛ/dT ∼ 0. For our cooling function (Fig. 1 of Kim & Ostriker 2017), there is a maximum at T ∼ 10 5 K.We note however that thermal instability only develops in selected lower-density regions within our warm winds (see Figure 14), rather than throughout the flow.It is unknown whether on larger scales this could ultimately lead to a thermally-driven wind, as has been suggested by Modak et al. (2023); within our simulation domain, the CR pressure remains quite large compared to the thermal gas pressure.

CONCLUSIONS AND PROSPECTS
The primary focus of this paper is on the dynamical interaction of CRs with multiphase gas in extraplanar regions, in order to understand when and how CRs are able to drive galactic winds.Our study makes use of the TIGRESS MHD simulations of portions of star-forming galactic disks (Kim & Ostriker 2017;Kim et al. 2020a) for initial conditions of the gas and magnetic fields, and also to set CR injection rates based on SN rates in the simulations.
As a prelude to our investigation of wind driving, we follow up on Paper I and Paper II to analyse how the MHD backreaction modifies the transport and distribution of CRs, thereby testing the conclusions of our previous studies.In the postprocessing CR simulations of Paper I and Paper II (with frozen MHD variables), the prevailing orientations of the velocity and magnetic field lines cause CRs to be confined within the warm+cold, dense gas.However, once "live MHD" is turned on, the excess of CR pressure in these regions forces dense gas to expand into the hot gas.The resulting rearrangement of the velocity and magnetic field topology enables CRs to propagate out of the dense gas.The net result is a more uniform CR distribution between the hot and warm+cold gas regions when MHD is "live" (compare Figure 1, Figure 2).Even though the "live MHD" simulations produce a smoother CR distribution compared to the postprocessing simulations, our other previous conclusions regarding CR transport remain valid.In particular, we confirm that the scattering coefficient varies over more than four orders of magnitude depending on the gas properties (e.g., gas density and ionization), and that the transport of CRs is different in different thermal phases of the gas (Figure 7, Figure 8).In ionized gas at T > 10 4 K, the scattering coefficient is σ ∥ ∼ 10 −28 − 10 −27 cm −2 s, whereas in mostly-neutral gas at T < 10 4 K, σ ∥ < 10 −31 cm −2 s.As a result of the high scattering rate, CR transport is primarily by advection in the high-velocity, low-density hot gas, and by streaming at the Alfvén speed in intermediate-density warm ionized gas.Given the extremely low scattering rate in neutral warm and cold gas (which is the majority of mass in the ISM), diffusion dominates CR transport in these regions.The CR pressure near the midplane is comparable to the kinetic, magnetic, and thermal pressure in the gas, but the scale height of CRs is far larger than that of the warm+cold gas (Figure 6).
In our studies of CR-driven winds, for the initial conditions we extract a domain from the extra-planar region (z > 500 pc) of the TIGRESS solar neighborhood MHD simulation outputs, considering both an epoch of strong and weaker outflow (see Table 1 for conditions of each case).We conduct simulations (with duration ∼ 100 Myr) in which we inject both gas and CRs with prescribed boundary conditions at the bottom of the simulation box.In the hot wind models (entries 1-6 in Table 2), the inflow boundary conditions for the MHD variables are set to their initial mean values from TI-GRESS in the hot outflowing gas at z = 500 pc, while exploring different conditions of CR pressure.In the warm wind models (entries 7-12 in Table 2), the injected gas is warm (T = 10 4 K) and slow-moving; the gas velocity and magnetic field inflow boundary conditions are set to the mean values within the warm outflow from TIGRESS at z = 500 pc, while exploring different values of injected gas density and CR pressure.
We find that CRs in fast-moving, predominantly-hot winds do not aid in the outflow acceleration, but instead slow down the gas.In both of our M1 -HW models with fast-moving hot gas outflows, the CR pressure evolves to increase outward because the advection velocity decreases outward, so that momentum is transferred from gas to CRs (see Figure 12).When a powerful hot wind is present, its interaction with embedded denser/cooler cloudlets shreds and accelerates them, effectively transferring momentum to create large outward mass fluxes of warm gas (Figure 11).
By contrast, our warm wind models demonstrate that CRs can effectively transfer momentum to the surrounding gas as the latter moves towards higher altitudes, with the MHD momentum flux increasing by ∼ 200% and warm gas accelerating to ∼ 50 − 60 km s −1 within 3 kpc (Figure 16, Figure 18).The reason that CR-driven acceleration is effective for the warm wind case is that streaming at the Alfvén speed controls CR transport for the conditions of warm gas in the extra-planar region (Figure 21).When streaming dominates, the CR pressure decreases as P c ∝ ρ 2/3 , and therefore the effective sound speed increases upward as C 2 eff ≈ dP c /dρ ∝ ρ −1/3 ).This condition ensures that the CR sound speed remains larger than the gas vertical velocity, which is an essential requirement for the gas to be accelerated (see Equation 35).In our WW simulations, mass, momentum, and energy fluxes of the gas flowing out of the box are nearly independent of the injected gas density, while the fluxes increase with the level of CR pressure at the base of the extra-planar region.We also find that CRs affect the thermal state of the outflow: in the low-density (n H < 10 −4 cm −3 ), higher-temperature (T > 10 5 K) regime, gas becomes thermally unstable with the rate of heating due to CR streaming losses larger than the rate of radiative cooling.
In conclusion, our findings indicate that under realistic scattering rates that are self-consistently determined, CRs can effectively drive flows of warm gas away from the galactic disk.This process is likely to be important in locations where there has not been a recent superbubble breakout powering a hot wind.Our results suggest that CRs remain a promising candidate to help remove gas from galaxies and limit star formation over cosmic timescales.
It is important to emphasize, however, that while the simplified simulations presented in this study provide valuable insights into the conditions under which CRs are able to accelerate the surrounding gas, there is still much work to be done.In particular, to fully understand the initiation of multiphase outflows and the relative importance of the CR pressure force vs. transfer of momentum from hot gas in accelerating cooler and denser outflows, it is necessary to self-consistently include CR energy injection together with thermal/kinetic energy injection as part of the supernova feedback process.Incorporation of CRs within local star-forming ISM simulations such as TIGRESS offers a practical way to explore parameter space with sufficient numerical resolution to follow all of the ISM phases and their interaction with CRs.For comprehensive understanding of the complete outflow formation and evolution process, including acceleration to velocities exceeding the escape speed as streamlines open up, it will also be necessary to conduct simulations of entire galactic disks.For simulations of dwarf galaxies, it will be possible to employ sufficiently high resolution to directly follow star formation and feedback processes including CR physics; for simulations of massive galaxies at lower resolution, a wind launching model may be calibrated based on resolved local-disk simulations (as in Kim et al. 2020b) but including CRs.

Figure 1 .Figure 2 .
Figure 1.Sample TIGRESS snapshot post-processed with the algorithm for CR transport (Jiang & Oh 2018; Armillotta et al. 2021).The upper (lower) row of panels shows x-z (x-y) slices through the center of the simulation box, where x, y, and z are the local radial, azimuthal, and vertical directions.From left to right, columns show hydrogen number density nH, gas temperature T , gas speed |v|, ion Alfvén speed |vA,i|, scattering coefficient σ ∥ , CR pressure Pc, and magnitude of the vertical CR flux |Fc,z|.The arrows overlaid on the gas velocity, Alfvén speed, and vertical CR flux slices indicate the projected directions of the gas velocity, Alfvén speed, and CR flux, respectively, in each slice.The CR pressure is divided by the Boltzmann constant kB = 1.38 × 10 −16 erg K −1 .
ure 5).For t ≥ 1 Myr, the distribution of v • ∇P c /|∇P c | peaks near values ≳ 0, as CR pressure gradients are smoothed out at interfaces (see panels at t = 1 Myr in Figure 3).Both the distribution of B • ∇P c /|∇P c | and v • ∇P c /|∇P c | achieve convergence for t ≥ 1 Myr.

Figure 3 .
Figure3.Zoom-in on a filament of warm gas located just above the disk in the simulation shown in Figure2.The snapshots are taken at t = 0 Myr (top panels), t = 0.5 Myr (middle panels), and t = 1 Myr (bottom panels).Each row shows, from left to right, slices of gas temperature T , hydrogen number density nH, CR pressure Pc, magnitude of gas velocity v, and magnitude of ion Alfvén speed vA,i.The arrows overlaid on the hydrogen density and CR pressure slices indicate the directions of the magnetic field and gas velocity, respectively.

Figure 5 .
Figure 5. Two-dimensional distribution of CR pressure Pc vs. gas density nH at t = 0 (left panel ) and t = 3 Myr (right panel ).The distributions are computed as a two-dimensional probability density function showing the volume of gas within each logarithmic bin, normalized by the bin area.

Figure 6 .Figure 7 .
Figure6.Horizontally-averaged vertical profiles of CR pressure Pc (purple), thermal pressure Pt (cyan), vertical kinetic pressure P k,z (yellow), vertical magnetic stress Pm,z (orange) at t = 0 (left panel ) and t = 3 Myr (right panel ).The average is computed over 10 snapshots output from TIGRESS and evolved with time-dependent MHD and CR physics in this work.The shaded areas cover the 16th and 84th percentiles of fluctuations.

Figure 8 .
Figure 8. Medians of the three vertical components of the effective CR propagation speed -advection speed v (orange), ion Alfvén speed vA,i (green), and diffusive speed v d (blue) -as a function of hydrogen density nH (left panel ) and gas temperature T (right panel ) at t = 0 (dashed lines) and t = 3 Myr (solid lines).

)Figure 9 .
Figure 9.Initial conditions of two sample galactic outflow simulations: M1 -HW -LPCR (left panels) and M2 -HW -LPCR (right panels) -see Table2.The initial conditions for the MHD variables are extracted from the M1 and M2 TIGRESS snapshots, respectively.From top to bottom, the panels show the slices through y = 0 of hydrogen number density nH, gas temperature T , gas vertical velocity vz, vertical magnetic field Bz, and CR pressure Pc.

Figure 10 .
Figure 10.Snapshots from three hot wind simulations: M1 -HW -LPCR (fast wind and low CR pressure; left column), M1 -HW -HPCR (fast wind and high CR pressure; middle column), M2 -HW -LPCR (relatively slow wind and low CR pressure; right column).From top to bottom, the panels show the slices through y = 0 of hydrogen number density nH, gas temperature T , gas vertical velocity vz, and CR pressure Pc divided by the injected CR pressure Pc,0.The snapshots are taken at t = 70 Myr.

Figure 11 .
Figure 11.Time evolution of outward mass fluxes FM, out for the M1 higher velocity hot wind models (top panels) and the M2 lower velocity hot wind models (bottom panels) simulations; note that upper and lower panels have different scales.The fluxes are taken at z = 3 kpc.Each column represents different phases: warm (left column), warm-hot (middle column), and hot (right column).In each panel, different colors refer to different models: green for models without CRs (NoCR), coral for models with low Pc (LPCR), and blue for models with high Pc (HPCR).The dashed gray lines indicate the mass flux injected at the bottom of the simulation box (z = 0).The dotted lines indicate the time-averaged fluxes, with the averages taken over t = 20 − 100 Myr.

Figure 12 .
Figure 12.Temporally averaged vertical profiles of the "net" MHD momentum flux difference ∆(Fp, MHD − W) (purple lines) and of the CR momentum flux difference ∆Fp, c (orange lines) in the total gas for hot wind models (average is over t = 20 − 100 Myr).The difference is calculated between an arbitrary height z and the initial height zi = 0.5 kpc.Different panels refer to different models: M1 -HW -LPCR (left panel ), M1 -HW -HPCR (middle panel ), and M2 -HW -LPCR (right panel ).In each panel, the gray line shows the momentum change rate per unit area ⟨ ṗ⟩.The vertical profiles are divided by the MHD momentum flux Fp,MHD measured at zi = 0.5 kpc.
Figure13.Temporally averaged vertical profiles of CR pressure Pc and MHD pressure PMHD (first row ), fractional areas for each thermal phase f A,ph (second row ), typical velocities ṽout,ph (third row ), and hydrogen number densities n H,out,ph (forth row ) of the outflowing gas (average is over t = 20 − 100 Myr).Different columns refer to different hot wind models: M1 -HW -LPCR (left column), M1 -HW -HPCR (middle column), and M2 -HW -LPCR (right column).In the top panels, the profiles of CR and MHD pressures are indicated with orange and purple lines, respectively.In the central and bottom panels, different colors represent different gas phases: cyan for warm/warm-hot, and crimson for hot.In the second-row panels, different line styles are used to separate outflowing (vz > 0) from inflowing (vz < 0) gas: dotted for outflowing, dashed for inflowing, solid for total.In the bottom-row panels, solid and dotted lines represent the typical (Equation22) hydrogen densities ñH,out,ph and the horizontally-averaged (Equation17) hydrogen densities, respectively.

Figure 15 .
Figure14displays the distribution on the grid of hydrogen number density n H , gas temperature T , vertical velocity v z , and CR pressure P c , in the three M1 -WW models at t = 100 Myr.Unlike the HW models, P c now decreases outward along the vertical direction in all models, suggesting a non negligible contribution of CRs to counteracting gravity and accelerating gas in this extraplanar region.A comparison between the three panels shows that the large-scale vertical gradient of CR pressure becomes smaller going from the model with higher injected gas density (M1 -WW -HD) to the model with lower injected density (M1 -WW -LD).At z > 1 kpc, the gas density and velocity distributions are qualitatively similar in the three models: lower-density gas flows outward, while higher-density gas mostly flows inward.The outflowing gas is clearly accelerated by CR pressure forces as thermal pressure forces are negligible in the warm gas populating most of the simulation box.A difference that emerges from a visual comparison is that the fraction of gas at T ≳ 10 5 K (color-coded in yellow and red in the second row of Figure14) increases going from M1 -WW -HD to M1 -WW -ID to M1 -WW -LD.At the time the snapshots are taken (t = 100 Myr), the fast-moving hot gas initially present in the simulation domain has already escaped the box.As we shall discuss in Section 6.4, gas at T > 10 4 K currently present in the simulation box is low-density gas heated up by CRs.In the top panels of Figure15, we plot the time evolution of the outward (v z > 0) and inward (v z < 0) mass flux F M (Equation18), MHD momentum flux F p,MHD

Figure 16 .
Figure 16.Temporally averaged vertical profiles of momentum flux difference for the warm wind simulations M1 -WW -HD (left panel ), M1 -WW -ID (middle panel ), and M1 -WW -LD (right panel ); average is over t = 50 − 150 Myr.The purple line shows the "net" MHD momentum flux difference ∆(Fp, MHD − W) taking into account the flux loss due to gas climbing out of the potential, the orange line shows the CR momentum flux difference ∆Fp, c, while the gray line shows the momentum change rate per unit area ⟨ ṗ⟩/A.The difference is calculated between an arbitrary height z and the initial height zi = 0.5 kpc.In all panels, the profiles are divided by the MHD momentum flux Fp, MHD at zi = 0.5 kpc.The time average is taken starting from t = 30 Myr to t = 150 Myr.On average, the system is in quasi-steady state (⟨ ṗ⟩ is much smaller than the other terms).

Figure 17 .
Figure 17.averaged mass flux FM (upper panels) and MHD momentum flux Fp, MHD (bottom panels) as a function of z for the warm wind simulations M1 -WW -HD (left panels), M1 -WW -ID (middle panels), and M1 -WW -LD (right panels); average is over t = 50 − 150 Myr.The green and red lines represent the outward and inward mass fluxes, respectively.

Figure 18 .
Figure 18.Temporally averaged vertical profiles of CR pressure Pc (left panel ), typical hydrogen number density ñH,out (middle panel ), and typical vertical velocity ṽz,out (right panel ) of the outflowing gas.Different colors represent different warm wind models: green for M1 -WW -HD, coral for M1 -WW -ID, blue for M1 -WW -LD.

Figure 19 .
Figure 19.Temporally averaged vertical profiles of CR pressure Pc (left panel ), typical hydrogen number density ñH,out (middle panel ), and typical vertical velocity ṽz,out (right panel ) of the outflowing gas.Different colors represent different warm wind models: coral for M1 -WW -ID (with higher CR pressure and flux), magenta for M2 -WW -LD (with lower CR pressure and flux).

Figure 20 .
Figure20.Comparisons of warm wind simulation results with the predictions of the one-dimensional (1D) steady-state model for CR-driven winds byMao & Ostriker (2018).The left panel shows the vertical profile of CR pressure Pc from the simulation M1 -WW -ID (red line) compared with the profile predicted by the 1D model (black line).The middle panel shows the gas velocity variation along the z axis dzvz produced by the simulation M1 -WW -ID (red line) in comparison to the 1D model prediction (black line).The third panel shows the vertical profiles of the root mean square gas velocity ⟨v 2 z ⟩ 1/2 (purple line) and the effective sound speed C 2 eff (light green line).The 1D model predicts that, in the Cartesian geometry of the present simulations, C 2 eff has to be larger than v 2 z in order to have a steady-state outflow.

Figure 21 .
Figure 21.Horizontally and temporally averaged profiles of the vertical components of the gas velocity vz (orange line) and ion Alfvèn speed vA,i,z (gray line) in M1 -WW -ID.
Cosmic-ray pressure in the disk

Figure 22 .
Figure 22.Medians of the rates of radiative cooling n 2 H |Λ|, photoelectric heating nHΓ, and CR streaming heating vs • σ ↔ tot • (Fc − 4/3vec) as a function of hydrogen density nH (left panel ) and gas temperature T (right panel ) in M1 -WW -ID.For both panels, the shaded areas cover the 16th and percentiles of fluctuations.

5
) with Σ SFR the SFR surface density, m ⋆ the total mass of new stars per supernova (95.5M ⊙ in Kim & Ostriker 2017, from a Kroupa IMF), H c,eff = ⟨|d ln P c /dz|⟩ −1 These numbers are obtained by integrating the energy sink terms over the simulation domain within |z| = H c,eff , with H c,eff an effective CR scale height, and taking averages over multiple snapshots.
In Equation 11, x i = n i /n H is the ion fraction, m n is the neutral mass, n n the neutral number density (see Section 2.2.5 of Paper I for the derivation of x i , n n , m n , and m i ), ⟨σv⟩ in ∼ 3 × 10 −9 cm 3 s −1 is the rate coefficient for collisions between H and H + (Draine 2011, Table2.1).

Table 2 .
Model list and inflow boundary conditions.