# The Launching of Cosmic Ray Driven Outflows

Xiaoshan Huang (黄小珊)1 Shane W. Davis1
1Department of Astronomy, University of Virginia, Charlottesville, VA 22904,
USA E-mail<EMAIL_ADDRESS>0000-0003-2868-489X 0000-0001-7488-4468

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Cosmic rays (CRs) are thought to be an important feedback mechanism in star-
forming galaxies. They can provide an important source of pressure support and
possibly drive outflows. We perform multidimensional CR-magnetohydrodynamic
simulations including transport by streaming and diffusion to investigate wind
launching from an initially hydrostatic atmosphere by CRs. We estimate a
characteristic Eddington limit on the CR flux for which the CR force exceeds
gravity and compare it to simulated systems. Scaling our results to conditions
in star-forming galaxies, we find that CRs are likely to contribute to driving
outflows for a broad range of star formation environments. We quantify the
momentum and energy transfer between CRs and gas, along with the associated
mass outflow rates under different assumptions about the relative importance
of streaming and diffusion for transport. In simulations with streaming, we
observe the growth and saturation of the CR acoustic instability, but the CRs
and gas remain well coupled, with CR momentum transferred efficiently to the
gas even when this instability is present. Higher CR fluxes transferr more
energy to the gas and drive stronger outflows. When streaming is present, most
of the transferred energy takes the form of Alfvén wave heating of the gas,
raising its pressure and internal energy, with a lower fractional contribution
to the kinetic energy of the outflow. We also consider runs with radiative
cooling, which modifies gas temperature and pressure profiles but does not
seem to have a large impact on the mass outflow for super-Eddington CR fluxes.

###### keywords:

galaxies: ISM — CR-magnetohydrodynamics — ISM: jets and outflows —method:
numerical simulation

††pubyear: 2021

##  1 Introduction

The general picture that stars are formed from collapsing gas is widely agreed
upon. However, it has long been a puzzle that the observed star formation rate
(SFR) is usually below what is predicted by the simple assumption of
gravitational collapse of gas clouds on dynamical timescales (Kennicutt Jr,
1998; Krumholz & Tan, 2007; Kennicutt Jr & Evans, 2012). Cosmological
simulations resolving galaxy formation indicate that the energy and momentum
feedback from star formation is a key component for explaining the inefficient
conversion from gas to stars (Somerville & Davé, 2015). Such feedback can
happen via numerous channels. Galactic outflow, which is observed in wide
range of galaxies (Veilleux et al., 2005; Spilker et al., 2020; Veilleux et
al., 2020), is thought to be an important feedback mechanism for removing gas
from galaxies, shaping the circum-galactic/inter-galactic medium, and
modifying galaxy evolution.

The origin of galactic outflow is not completely understood. There are several
promising mechanisms for launching outflows, but their relative importance in
different types of galaxies is still a matter of debate (Zhang, 2018).
Although there is evidence AGN may be relevant in some cases, the outflows
correlate with star formation activity (Martin, 2005; Veilleux et al., 2005),
implying that mechanisms tied to star-formation operate. Hot supernova-driven
thermal winds are almost certainly present in some galaxies (Chevalier &
Clegg, 1985). Entrainment of colder gas in these hot winds (Strickland &
Heckman, 2009; Scannapieco & Brüggen, 2015; Fielding et al., 2018; Gronke &
Oh, 2020), and radiation pressure on dust (Krumholz & Thompson, 2012; Zhang et
al., 2018; Huang et al., 2020; Kim et al., 2021) are promising mechanisms for
explaining the molecular and weakly ionised gas that is observed. We refer to
the latter mechanism as non-thermal feedback since relativistic particles
(photons) are doing the driving rather than hot gas.

Another non-thermal feedback candidate is cosmic ray (CR) pressure. Being
roughly in energy equipartition in our local universe (Boulares & Cox, 1990;
Ferriere, 2001), CRs are thought to be an attractive source to drive outflows
in Milky Way-like systems (Ipavich, 1975; Breitschwerdt et al., 1991a, 1993;
Everett et al., 2008), as well as in star-bursting systems (Socrates et al.,
2008). Meanwhile, numerical works also suggest that CR can drive outflow with
significant mass-loading and effectively regulate the star formation rate
(Jubelgas et al., 2008; Uhlig et al., 2012; Booth et al., 2013; Hanasz et al.,
2013; Salem & Bryan, 2014).

Recent semi-analytical and numerical simulations further elaborate CR’s
feedback role in terms of driving outflow, modifying galaxy evolution and
shaping circum-galactic environment at different scales. For example,
cosmological and galaxy formation simulations suggest that when compared to
simulations only including thermal feedback, adding CRs changes the star
formation history, especially in massive galaxies (Ruszkowski et al., 2017;
Butsky & Quinn, 2018; Hopkins et al., 2020). CRs usually enhance outflow,
constrain wind morphology and thicken the galaxy disk (Bustard et al., 2020;
Buck et al., 2020). Interestingly, both analytical work and numerical studies
suggest that the CR driven wind is likely distinct from a typical momentum or
energy driven wind, leaving a unique wind profile (Salem & Bryan, 2014; Mao &
Ostriker, 2018; Girichidis et al., 2018). In the circum-galactic medium (CGM)
or intra-galactic medium (IGM), CRs are an important pressure source that
supports gas in addition to thermal pressure and potentially allows a cooler
halo (Butsky & Quinn, 2018; Buck et al., 2020). The collisionless and
collisional loss of CR energy is also an important energy component in the CGM
(Bustard et al., 2020; Crocker et al., 2021).

The cosmological and galaxy evolution simulations including CRs provide
straightforward comparisons with observation, showing that CRs are an
important feedback mechanism. In addition, high resolution local simulations
with resolved multi-phase gas structure are another channel to diagnose CR
physics and its interaction with radiative cooling, conduction, and magnetic
field evolution. For example, CRs are able to drive an initially cold cloud
that embedded in hot background to roughly the observed outflow velocities
without complete disruption by local dynamical instability (Wiener et al.,
2019; Brüggen & Scannapieco, 2020). In the CGM, CRs modify the onset of
thermal instability and subsequent multi-phase gas formation Butsky et al.
(2020); Kempski & Quataert (2020).

To connect the cloud-wind interaction scale and the galaxy formation scale, in
this work, we investigate the CR’s role in launching an outflow from dense
atmosphere against gravity near the disk with high resolution simulations.
Similar studies of radiative feedback found that the wind-launching process
could be impacted by the presence of dynamical instabilities (Krumholz &
Thompson, 2012; Davis et al., 2014). Resolving local dynamics and potential
instability (Begelman & Zweibel, 1994; Heintz & Zweibel, 2018) is an essential
step for bridging the gap between well-resolved simulations of individual
cloud evolution and incorporating CR to galaxy evolution framework employing
sub-grid models.

Despite the importance of CR feedback, our knowledge about the CR transfer
mechanism is limited by the lack of direct extra-galactic observation. The
canonical treatment of CR transfer usually includes two processes: streaming
and diffusion.The latter models CR transfer as classical diffusive process
with effective diffusivity (Quataert et al., 2021; Semenov et al., 2021).
Transport by streaming means that CRs stream along the magnetic field and
excite waves, with subsequent wave damping enabling energy and momentum
exchange with gas. The streaming velocity is usually determined by the balance
between wave growth and wave damping, as well as local ionisation condition,
ranging from Alfvénic to super-Alfvénic streaming (Skilling, 1971; Holguin et
al., 2019; Bai et al., 2019; Bustard & Zweibel, 2020; Zweibel, 2020; Bambic et
al., 2021). Extensive studies show that the transfer mechanism has significant
impact the disk and halo morphology, the outflow properties and the energetics
of CGM and IGM (Wiener et al., 2013; Ruszkowski et al., 2017; Wiener et al.,
2017; Buck et al., 2020; Ramzan et al., 2020).

This work aims to study the efficiency of CR driving from a disk atmosphere as
a function of a prescribed CR flux, under different assumptions about the
transport mechanism.We adopt the two-moment scheme proposed by Jiang & Oh
(2018), which solves CR transfer equations including both streaming and
diffusion self-consistently. We perform a series of adiabatic simulations
primarily with Alfvénic streaming and isotropic diffusion to explore how
relative importance of streaming and diffusion changes CR feedback. We also
briefly consider the role of radiative cooling on the launching CR winds. A
detailed study of radiative cooling in CR driven multi-phase outflow and its
impact on cold gas survivability will be presented in a companion paper.

The plan of this paper is as follows: in Section 2 we describe the equations
solved, simulation set-up and our prescription for a CR Eddington limit. The
simulation results are reported in Section 3. We discuss the primary
implication from our results and relevance to star-forming galaxies in Section
4, and summarise the main conclusions in Section 5.

##  2 Numerical Method

###  2.1 Cosmic Ray Transportation Equations

The equations we solve are ideal magneto hydrodynamics (MHD) equations and
cosmic rays (CR) transportation equations, which is based on the two-moment
scheme proposed by Jiang & Oh (2018):

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0,$
(1a) $\displaystyle\frac{\partial(\rho\mathbf{v})}{\partial
t}+\nabla\cdot(\rho\mathbf{v}\mathbf{v}-\textbf{BB}+\textsf{P*})=$
$\displaystyle\sigma_{\rm c}\cdot[\mathbf{F}_{\rm c}-\mathbf{v}\cdot(E_{\rm
c}\textsf{I}+\textsf{P}_{\rm c})]+\rho g$ , (1b) $\displaystyle\frac{\partial
E}{\partial
t}+\nabla\cdot[(E+P^{*})\mathbf{v}-\textbf{B}(\textbf{B}\cdot\mathbf{v})]=$
$\displaystyle-(\mathbf{v}+\mathbf{v}_{\rm{s}})\cdot(\nabla\cdot\textsf{P}_{\rm
c}),$ (1c) $\displaystyle\frac{\partial\textbf{B}}{\partial
t}-\nabla\times(\mathbf{v}\times\textbf{B})=0,$ (1d)
$\displaystyle\frac{\partial E_{\rm c}}{\partial t}+\nabla\cdot\mathbf{F}_{\rm
c}=$
$\displaystyle-(\mathbf{v}+\mathbf{v}_{\rm{s}})\cdot\sigma_{c}\cdot[\mathbf{F}_{\rm
c}-\mathbf{v}\cdot(E_{\rm c}\textsf{I}+\textsf{P}_{\rm c})]+Q$ (1e)
$\displaystyle\frac{1}{V_{\rm m}^{2}}\frac{\partial\mathbf{F}_{\rm
c}}{\partial t}+\nabla\cdot\textsf{P}_{\rm c}=$
$\displaystyle-\sigma_{c}\cdot[\mathbf{F}_{\rm c}-\mathbf{v}\cdot(E_{\rm
c}\textsf{I}+\textsf{P}_{\rm c})]$ . (1f)

Here $\rho$, $\mathbf{v}$, $E$ are fluid density, velocity, total energy, B is
magnetic field strength. $P^{*}$ is the sum of gas pressure and magnetic
pressure with permeability $\mu=1$.

Equation 1e and Equation 1f are the CR momentum and energy equations. They are
the first and second moment integration of the CR advection-diffusion
equation, which was derived by Skilling (1971).
$\mathbf{v}_{\rm{s}}=-\textbf{sgn}(\textbf{B}\cdot\nabla P_{\rm
c})\mathbf{v}_{A}$ is the CR streaming velocity. It has the magnitude of
Alfvén velocity and directs down the CR pressure gradient. This description of
CR transportation is based on the “self-confinement” picture. In this
framework, when CRs stream through plasma, they excite Alfvén waves. When the
streaming velocity exceeds the Alfvén velocity, the Alfvén wave is amplified
and eventually saturates. As the Alfvén wave grows and perturb the magnetic
field, the irregularities in magnetic field will scatter the CRs and change
their pitch angles, reducing the CR streaming velocity to near the Alfvén
speed. So the CRs stream at local Alfvén speed and down the CR pressure
gradient. The wave damping enables the momentum and energy exchange between
CRs and gas.

$E_{\rm c}$ and $\mathbf{F}_{\rm c}$ are the cosmic ray energy density and
flux. Equation 1e is widely adopted in CR-MHD simulations. Compared to the
one-moment scheme, the two-moment scheme solves the time-dependent Equation 1f
to update $\mathbf{F}_{\rm c}$, instead of adopting steady-state CR flux
description. Solving separate CR energy and momentum equations makes the code
less diffusive because it better handles streaming velocity where the CR
pressure gradient is flat. $V_{\rm m}$ is the maximum CR propagation velocity,
it is constant through the whole domain. This is analogous to adopting the
reduced speed of light approximation in an explicit two-moment method to solve
radiation transfer equation (Skinner & Ostriker, 2013), $V_{\rm m}\ll c$
replaces the speed of light $c$ in Equation 1f to relax the timestep. As long
as $V_{\rm m}$ is significantly larger than the maximum flow or Alfvén
velocity, the impact on the dynamics is limited (Jiang & Oh, 2018).

$\sigma_{c}$ quantifies the CR-gas interaction:

$\sigma_{\rm c}^{-1}=\sigma_{\rm
c}^{{}^{\prime}-1}+\frac{\textbf{B}}{|\textbf{B}\cdot(\nabla\cdot\textsf{P}_{\rm
c})|}\mathbf{v}_{\rm A}\cdot(E_{\rm c}\textsf{I}+\textsf{P}_{\rm c}),$ (2)

where $\sigma_{c}^{\prime-1}$ is the classical CR diffusion coefficient
$\kappa_{\rm diff}$. The second terms is an effective streaming coefficient
$\sigma_{\rm str}^{-1}$. When $d\mathbf{F}_{\rm c}/dt$ terms disappears,
Equation 1e reduces to the classic CR energy transport equation (Breitschwerdt
et al., 1991b). Because the source terms for cosmic rays (Equation 1e,1f) and
the fluid (Equation 1b,1c) are the same except for the opposite signs, the
conserved variables are the total energy $E+E_{\rm c}$ and momentum
$\rho\mathbf{v}+\mathbf{F}_{\rm c}/V_{\rm m}^{2}$.

###  2.2 Eddington flux for Cosmic Rays

A primary goal of this study is to examine the impact of the CR transport
mechanism (specifically streaming versus diffusion) on the efficiency of wind
launching and acceleration. But, the launching is also sensitive to the flux
of CRs through the atmosphere, so we begin by estimating what magnitude of CR
flux we expect to be necessary to drive winds in the limit that either
streaming or diffusion dominates. Motivated by Socrates et al. (2008), we
formulate this in terms of a CR Eddington flux, where the acceleration due to
CR forces just balance the local gravitational acceleration $g$, which we take
to be constant. For a given flux $F_{\rm CR,x}$ at the base of the atmosphere,
we define the CR Eddington flux $F_{\rm Edd}$ as the the flux that satisfies
$-\partial P_{\rm c}/\partial x=-\rho g$, where $P_{\rm c}$ is the CR
pressure. We assume that the CRs are nearly isotropic so that $E_{\rm
c}=3P_{\rm c}$.

First we derive an expression for $F_{\rm edd}$ in the CR streaming limit.
When streaming dominates the CR transport, $v_{\rm s}=v_{\rm A}$ if $E_{\rm
c}$ is monotonic. Because $V_{\rm m}$ is usually significantly larger than any
velocity in the simulation, we further assume that the time dependent term in
Equation (1f) is negligible, so that $F_{{\rm c},x}\approx 4(v_{\rm
s}+v_{x})P_{\rm c}$, where $v_{x}$ is the gas velocity along $x$ direction. In
other words, the CR flux is the sum of advective flux and streaming flux. In
the launching region of the flow, $v_{x}$ is small and $v_{\rm A}\gg v_{x}$,
so we further assume that $F_{{\rm c},x}\approx 4v_{\rm A}P_{\rm c}$. The CR
energy Equation (1e) in 1D geometry becomes:

$\frac{\partial E_{\rm c}}{\partial t}+3v_{\rm A}\frac{\partial P_{\rm
c}}{\partial x}+4P_{\rm c}\frac{\partial v_{\rm A}}{\partial x}=0,$ (3)

In steady state, the time dependent term vanishes, resulting $\partial(P_{\rm
c}v_{\rm A}^{4/3})/\partial x=0$.

With these assumptions, we can approximate

$\frac{\partial P_{\rm c}}{\partial x}\simeq\frac{4P_{\rm c}}{3v_{\rm
A}}\frac{\partial v_{\rm A}}{\partial x}=-\frac{2P_{\rm
c}}{3\rho}\frac{\partial\rho}{\partial x},$ (4)

where the second equality follows from assuming $\mathbf{B}$ is constant.
Defining $H_{\rho}\equiv(-\partial\rm ln\rho/\partial x)^{-1}$, and setting
$\partial P_{\rm c}/\partial x=\rho g$, we can solve for

$F_{\rm Edd,str}=6v_{\rm A}g\rho H_{\rho},$ (5)

which we define as the Eddington flux for CRs in the streaming limit. Note
that unlike the expression for radiation with electron scattering opacity,
this expression is not constant, depending on the background density and
density gradient. This means that for a given flux, the acceleration of the
flow against gravity will generally be most effective at lower densities. In
contrast, to radiation dominated atmospheres, we might expect that streaming
CRs will more generically drive outflows, but possibly only in superficially
low density regions of the atmosphere. This may lead to effective
acceleration, but only for a relatively small amount of gas.

The diffusion limit is simpler, with $F_{\rm CR,x}\approx-\kappa_{\rm
diff}(\partial P_{\rm c}/\partial x)$, where $\kappa_{\rm diff}$ is a
characteristic CR diffusivity. So the Eddington flux in the diffusion limit is

$F_{\rm Edd,dif}\equiv-\kappa_{\rm diff}\frac{\partial P_{\rm c}}{\partial
x}=\kappa_{\rm diff}\rho g.$ (6)

For $\kappa_{\rm diff}$ approximately constant, this is somewhat more like the
radiative Eddington limit, but with an opacity that is inversely proportional
to density, so the result again suggest that lower density gas will be more
readily accelerated.

###  2.3 Simulation Setup

####  2.3.1 Units and Scaling

When cooling is neglected, Equations (1a)-(1f) can be non-dimensionalised and
rescaled via three free parameters, which we will take to be the
characteristic length scale $l_{0}$, temperature $T_{0}$, and density
$\rho_{0}$. The characteristic velocity $v_{0}$ is set to an isothermal sound
speed using $v_{0}\equiv c_{s,0}=\sqrt{k_{b}T_{0}/(0.6m_{p})}$, where $k_{\rm
B}$ is Botzmann’s constant, $m_{p}$ is the proton mass. This leads to
dimensionless time $t_{0}=l_{0}/v_{0}$, energy density
$E_{0}=\rho_{0}v_{0}^{2}$, flux $F_{0}=\rho_{0}v_{0}^{3}$, and acceleration
$a_{0}=v_{0}^{2}/l_{0}$. The choice of $T_{0}$, $\rho_{0}$, and $l_{0}$ can be
rescaled to match varying galactic environments, as discussed in section 4.2.
The exception will be where cooling is discussed in section 3.3, as cooling
depends explicitly on the choice of temperature and density.

In the rest of paper, we report the dimensionless numbers unless otherwise
noted. The fiducial vertical ($x$) and horizontal ($y$ or $z$) sizes are
$L_{x}=400$ and $L_{y}=50$, with $8000\times 1000$ total cells for the
fiducial resolution. The aspect ratio of calculation domain varies between the
simulations to adjust for specific problems as well as to minimize the
computation cost. We discuss the effect of domain aspect ratio in Section
3.4.1. In all simulations, we set the pressure floor $P_{\rm floor}=10^{-7}$
and density floor $\rho_{\rm floor}=10^{-8}$. We set the fiducial $V_{\rm
m}=1000$ and discuss the effect of $V_{\rm m}$ in Section 3.4.2.

####  2.3.2 Initialization and boundary conditions

We envision our simulation as an atmosphere situated above a star forming
environment (e.g. a galactic disk). The gravity $g$ is created by the gas and
stars below the atmosphere and a CR flux is provided by the stellar winds and
supernovae occurring in the (unsimulated) environment below. Although
idealised, this setup allows us to isolate and explore the impact of varying
the CR flux in the launching region of the outflow, which is uncertain due to
open questions about CR production, transport and destruction within the star
forming environment. The downside of this setup is that it does not admit
steady state wind solutions of the form present in spherical symmetry or other
geometry that allows the wind to expand into a larger cross sectional area as
it moves out (Ipavich, 1975; Mao & Ostriker, 2018; Quataert et al., 2021).

We initialise the gas density to an isothermal profile in hydrostatic
equilibrium with static background gravity. The gas density is uniform in the
$y$ direction. In the $x$ direction, the distribution is the maximum of
$\rho_{0}\exp(-x/h)$ and $\rho_{\rm bkgd}=7.988\times 10^{-5}$. We choose a
scale height $h=2.018$. The initial isothermal temperature $T_{\rm
init}=1.667$. Hydrostatic equilibrium for $\rho>\rho_{\rm bkgd}$ requires
$g=-0.826$. The magnetic field is initialised to be uniform with
$B_{x}=B_{0}=2$ and all other components set to zero.

In the streaming limit, we set an isotropic diffusivity with $\kappa_{\rm
diff}=10^{-8}$. In the streaming-diffusion and non-streaming runs, we adopt
isotropic diffusion with the values reported in Table 1. The only exception is
HSE_1F_hd_aniso, where the diffusivity parallel to the field $\kappa_{\rm
diff,\parallel}=100$ and the diffusivity perpendicular to the magnetic field
$\kappa_{\rm diff,\bot}=10^{-8}$.

The boundaries are periodic in the horizontal ($y$ and $z$) directions. For
the $x$ direction the top MHD boundary is outflow, while the bottom MHD
boundary is reflecting for hydro variables and outflow for magnetic field
variables. In order to prescribe a uniform CR flux, we fix the CR flux
$F_{\rm{CR},x}$ in the bottom ghost zones and copy the CR energy density, $y$
(and $z$ for three dimensional runs) component of the CR flux to the ghost
zones from the last active zone. For the top CR boundary condition, we copy
the CR flux and energy from the last active zone. For the diffusion-only
simulations, we implemented an alternate right MHD boundary that prevents any
gas inflow, avoiding unphysical oscillations in the non-streaming simulations.
We confirmed that the streaming runs were insensitive to this choice for the
top MHD boundary condition.

Although we target a constant CR flux, we cannot precisely control the CR flux
entering the domain when streaming dominates due to impacts from the CR
“bottleneck” effect. After an initial transient phase, the boundary settles to
a nearly constant CR flux that is usually slightly below the target $F_{\rm
CR,x}$, with temporal variations of several percent or less of the average
value. Higher target fluxes provide results closer to the target value. For
the fiducial expected $F_{\rm CR,x}=15.0$, the time average between
$5.0<t<25.0$ of CR flux across the boundary in the simulation is
$\overline{F}_{\rm CR,in}=12.2$. We define the average difference

$\overline{\Delta F}_{\rm CR,in}=\frac{\int_{t}\left|<\mathcal{F}_{E_{\rm
CR}}>-\bar{F}_{\rm CR,in}\right|dt}{\int_{t}dt},$

where $<X>=\int_{\rm in}XdA/\int_{\rm in}dA$ is volume average of $X$ at the
bottom boundary, $\mathcal{F}_{E_{\rm CR}}$ is the flux of CR energy returned
by Riemann solver at the bottom boundary. For a target flux $F_{\rm
CR,x}=15.0$, the fluctuation $\overline{\Delta F}_{\rm
CR,in}/\overline{F}_{\rm CR,in}=7.5\%$. For target CR fluxes $F_{\rm
CR,x}=60.0$ and $F_{\rm CR,x}=135.0$, the average CR flux at the boundary is
$\overline{F}_{\rm CR,in}=59.074$, $\overline{F}_{\rm CR,in}=134.451$, and
fluctuation are only $1.1\%$ and $0.5\%$ accordingly.

##  3 Results

Table 1: Summary of Simulation Parameters

Name | $\kappa$ | $F_{\rm CR,x}$ | $B_{0}$ |  $L_{x}\times L_{y}$ b |  $N_{x}\times N_{z}$ b | $V_{\rm m}$
---|---|---|---|---|---|---
HSE_1F_str | $10^{-8}$ | 15.0 | 2.0 | $400\times 50$ | $8000\times 1000$ | 1000
HSE_hF_str | $10^{-8}$ | 7.5 | 2.0 | $400\times 50$ | $8000\times 1000$ | 1000
HSE_4F_str | $10^{-8}$ | 60.0 | 2.0 | $800\times 50$ | $16000\times 1000$ | 1000
HSE_9F_str | $10^{-8}$ | 135.0 | 2.0 | $800\times 50$ | $16000\times 1000$ | 1000
HSE_20F_str | $10^{-8}$ | 300.0 | 2.0 | $800\times 50$ | $16000\times 1000$ | 1000
HSE_1F_ld | $1.0$ | 15.0 | 2.0 | $200\times 5$ | $4000\times 100$ | 1000
HSE_1F_cd | $10.0$ | 15.0 | 2.0 | $200\times 25$ | $4000\times 500$ | 1000
HSE_1F_hd | $100.0$ | 15.0 | 2.0 | $200\times 5$ | $4000\times 100$ | 1000
HSE_1F_hd_aniso | $100.0$ | 15.0 | 2.0 | $200\times 5$ | $4000\times 100$ | 1000
HSE_1F_hd_nsa | $100.0$ | 15.0 | 2.0 | $200\times 25$ | $4000\times 500$ | 1000
HSE_1F_ld_nsa | $1.0$ | 15.0 | 2.0 | $200\times 25$ | $4000\times 500$ | 1000
HSE_1F_str_wb | $10^{-8}$ | 15.0 | 1.0 | $200\times 5$ | $4000\times 100$ | 1000
HSE_1F_str_sb | $10^{-8}$ | 15.0 | 4.0 | $200\times 5$ | $4000\times 100$ | 1000
HSE_1F_str_scc | $10^{-8}$ | 15.0 | 2.0 | $400\times 50$ | $8000\times 1000$ | 1000
HSE_1F_str_wcc | $10^{-8}$ | 15.0 | 2.0 | $400\times 50$ | $8000\times 1000$ | 1000
HSE_9F_str_scc | $10^{-8}$ | 135.0 | 2.0 | $400\times 50$ | $8000\times 1000$ | 1000
HSE_1F_str_VM | $100.0$ | 15.0 | 2.0 | $400\times 50$ | $8000\times 1000$ | 2000
HSE_1F_str_LR | $10^{-8}$ | 15.0 | 2.0 | $400\times 50$ | $4000\times 500$ | 1000
HSE_1F_str_HR | $10^{-8}$ | 15.0 | 2.0 | $400\times 50$ | $16000\times 2000$ | 1000
HSE_1F_str_3D | $10^{-8}$ | 15.0 | 2.0 | $400\times 25\times 25$ | $4000\times 250\times 250$ | 1000

  * a

Streaming is turned off, equivalently setting $\mathbf{v}_{\rm{s}}=0$ in
Equation 1c, Equation 1e, $\mathbf{v}_{\rm A}=0$ in Equation 2.

  * b

The reported dimension corresponds to the largest simulation domain for each
set of parameters. Notice that some of the simulation has lower aspect ratio
$L_{y}/L_{x}$, we discuss the effect of aspect ratio in Section 3.4.1

  * c

The three adiabatic simulations with different cooling strength are scaled to
different physical units.

###  3.1 The Effect of CR Flux on Streaming: from sub-Eddington to super-
Eddington

The five streaming dominated simulations are HSE_hF_str HSE_1F_str,
HSE_4F_str, HSE_9F_str and HSE_20F_str, the diffusivity is set to be
$\kappa_{\rm diff}=10^{-8}$, so that $\sigma_{\rm
c}^{-1}\approx\frac{\textbf{B}}{|\textbf{B}\cdot(\nabla\cdot\textsf{P}_{\rm
c})|}\mathbf{v}_{\rm A}\cdot(E_{\rm c}\textsf{I}+\textsf{P}_{\rm c})$. We
estimate the Eddington flux according to Equation 5. Assuming $v_{\rm s}$ is
calculated using the initial density $\rho_{\rm init}=\rho_{0}$ at the base of
the atmosphere and the initial uniform magnetic field $B_{0}$, the estimated
Eddington flux $F_{\rm edd,str}=20.0$.

We select HSE_1F_str with a target CR flux $F_{\rm CR,x}=15.0$ as our fiducial
run. This is somewhat sub-Eddington at the base of the domain, but becomes
super-Eddington as the density drops.

Figure 1: The density snapshots of HSE_1F_str. From left to right, t=1.5, 7.5,
18.0 and 24.0. The picture shows only half of the simulation domain in the $x$
direction.

Figure 1 shows density snapshots from HSE_1F_str for the bottom half of the
domain. At early time $t=1.5$, while the CR is streaming in the dense
atmosphere, the gas density profile remains close to the initial distribution.
This is shorter than the timescale for CRs to stream out of the gas, which we
approximate as the Alfvén-wave crossing time:

$\displaystyle t_{\rm str}=\int_{0}^{2l_{0}}\frac{dx}{v_{\rm s}(x)}$
$\displaystyle=\frac{2h_{0}}{B_{0}/\sqrt{\rho_{0}}}(1-e^{-\frac{l_{0}}{h_{0}}}),$

which gives $t_{\rm str}\sim 2.0$ for fiducial values. At later times, we see
that a substantial fraction of the gas is lofted up as the atmosphere expands
outward under a combination of CR and gas pressure. Since mass is conserved,
the density drops near the base. Snapshots at $t=18.0$ and $t=24.0$ are
similar, suggesting the gas is reaching a quasi equilibrium density profile,
even though a modest outflow is present at these later times. There is also
evidence of abrupt jumps in the density that show up more clearly in the
horizontally average plots below. We attribute these inhomogeneities to the
development of an acoustic instability, as discussed in section 4.3. Figure 1
shows that this instability does not lead to large amplitude variations on the
horizontal dimensions, at least for the magnetic field strength and geometry
used here.

Figure 2: The $y$ direction averaged density (the first row), CR pressure (the
blue lines in the second row), gas pressure (the red lines in the second row)
and gas velocity (the third row). The line style corresponds to different
time, the thick solid lines in each panel are the initial profile. Other
sampled times are t=1.8 (the thin solid lines), t=5.0 (the dashed lines) and
t=18.0 (the dashed dotted lines).

Figure 2 show horizontally averaged quantities and illustrates the
“bottleneck” effect. When the CR force acts on gas and streaming dominates the
CR flux, the streaming velocity $v_{\rm s}\propto 1/\sqrt{\rho}$, increasing
as the gas density drops. When CRs stream into the atmosphere, CRs tend to
pile up in the dense region where they stream slower, resulting in a CR
pressure gradient from dense gas to diffuse gas, producing the “bottleneck”.
Given the monotonic density profile, $P_{\rm c}$ reaches its maximum at base
of the atmosphere and decreases outward.

Around $t=1.8$, the gas density is similar to the initial profile and nearly
unperturbed. The gas pressure closely follows the CR pressure, and the
combined pressure gradients exceed the gravitational force, launching the
atmosphere upward. By $t=5.0$, the CRs begin to stream out of the dense gas
and drive the expansion, enhancing both the CR and gas pressure in the
initially low density background region ($x\gtrsim 20$). The large enhancement
in gas pressure is driven by CR heating. Gas velocities are quite high for
$x<100l_{0}$, due to the combined acceleration by CR and gas pressure
gradients.

As more gas is pushed towards the upper boundary, the density gradient becomes
shallower and the Alfvén velocity gradient consequently decreases (the $B$
field remains more constant). Therefore, the CR streaming “bottleneck” is less
effective, the CR pressure gradient decreases, providing less heating and
acceleration. The velocity decreases as more mass is swept up in the outflow
and continues to do work against the constant gravitational acceleration. The
combined gas and CR pressure gradients still exceed gravity, with velocity
increasing upward due to the resulting acceleration. At $t=18.0$, the diffuse
gas starts to flow out of the domain, carrying away a modest faction of the
initial mass. After this point and before we end the simulation at $t=25.0$,
all the quantities shown in this plot are changing slowly, with pressure and
density gradients gradually becoming more shallow with time.

Even after horizontal averaging, the effects of the instabilities noted above
can be discerned as the growth of a series of shocks in the gas density,
pressure, and velocity profiles for $45\lesssim x\lesssim 75$ at $t=18.0$. We
identify them as the result of CR acoustic instability and discuss them in
more detail in Section 4.3. The instability, however, saturates at a level
where it does not significantly change the volume averaged CR force. Despite
the instability, the magnetic field is relatively unperturbed, making CRs
effectively stream along the uniform magnetic field. Hence, there is no
analogue of the Rayleigh-Taylor instability that regulates the acceleration of
gas in radiatively driven atmospheres (Krumholz & Thompson, 2012).

The overall dynamical picture is that CRs accelerate gas against gravity even
with this moderate CR flux, but a significant fraction of the energy goes into
heating. According to Equation 1c, CRs heat the gas at the rate of
$\mathbf{v}_{\rm{A}}\cdot\nabla P_{\rm c}$. Since $v_{\rm A}\propto
1/\sqrt{\rho}$, the lower density regions tend to experience greater heating.
Gas is in the low density background can be heated by over two orders of
magnitude with the lowest density gas experiencing the largest heating. These
results suggest radiative cooling could play a significant role and we
consider its impact below.

Figure 3: The density snapshots at $t=10.0$ for four streaming dominated
simulations with different CR flux. From left to right, the injected CR flux
$F_{\rm CR,x}=7.5$(HSE_hF_str), 15.0(HSE_1F_str), 60.0(HSE_4F_str),
135.0(HSE_9F_str). We only show part of the domain in the $x$ direction (see
$L_{x}$ in Table 1).

Figures 3 and 4 show the impact of varying the target CR flux injected at the
base of the atmosphere. We consider simulations with CR fluxes 0.5
(HSE_hF_str), 4 (HSE_4F_str), 9 (HSE_9F_str) and 20 (HSE_20F_str) times the
flux in the fiducial run. Although the overall dynamics are similar,
HSE_4F_str and HSE_9F_str produce more significant gas outflow in shorter
time. In contrast, in the sub-Eddington case HSE_hf_str with half of the
fiducial CR flux, the atmosphere expands outward over a longer timescale, but
by the time we finish other simulations at $t=25$, there is no significant gas
outflow. Figure 3 compares the gas density snapshot at $t=10.0$ of these
simulations. As we increase the CR flux from sub-Eddington to super-Eddington,
the atmosphere is more dispersed and the outflow carries larger mass and
momentum. We see that as the CR flux increases, the shocks driven by acoustic
instability become more prominent, but the amplitude never becomes large
enough to have a significant feedback on the gas acceleration.

Figure 4: The gas outflow from simulations in the streaming limit with
different injecting CR flux: HSE_hF_str (green), HSE_1F_str (red), HSE_4F_str
(red), HSE_9F_str (blue) and HSE_20F_str (purple). In this plot, we report the
results from simulations with the same box size $L_{x}=400$ instead of the
reported $L_{x}$ in Table 1. Upper panel: the mass outflow rate per unit area.
Lower panel: the average outflow velocity. Both defined as Equation 7

Figure 4 shows the volume average of the mass outflow rate per area
$\dot{M}_{\rm out}$ (the upper panel) and averaged outflow velocity $v_{\rm
out}$ (the lower panel) at $x=400l_{0}$ as a function of time. $\dot{M}_{\rm
out}$ and $v_{\rm out}$ are calculated as:

$\dot{M}_{\rm out}=\frac{\int_{t}\int_{\rm
out}\mathcal{F}_{\rho}dydt}{\int_{t}\int_{\rm out}dydt},\quad v_{\rm
out}=\frac{\int_{t}\int_{\rm out}\rho vvdxdydt}{\int_{t}\int_{\rm
out}\mathcal{F}_{\rho}dydt},$ (7)

where $\mathcal{F}_{\rho}$ is the density flux returned by Riemann solver on
the outer boundary, $\rho vv|_{\rm out}$ is the $\rho vv$ calculated at the
outer boundary, $\int_{\rm out}$ denotes integration on the outer boundary. At
early times, there is modest inflow and velocities are negative due to the
fact that the isothermal and constant density background has no pressure
gradient to support the gas against gravity. For the sub-Eddington run, the
combination of CR and gas pressure is never sufficient to overcome gravity and
this inflow persists to the end of the simulation. For the marginal and super-
Eddington cases, the CR and gas pressure forces eventually exceed gravity and
drive outflows. The general trend is that both the velocities and mass outflow
rates increase as the injected CR flux increases.

Figure 5: Force terms in the $x$ direction as a function of time for the
nearly-Eddington (the left panel) HSE_1F_str and the Super-Eddington case
HSE_9F_str (the right panel). The force terms are summed over across the
domain. In both panels, the red dashed line is gas pressure gradient force,
the purple dashed line is CR pressure gradient force, and the gray solid line
is total $x$ direction momentum. The rate is sampled every $\Delta t=0.001$.

The outflow rates in simulations reflect the different force balance in sub-
Eddington systems and super-Eddington systems. Figure 5 shows the $x$
direction forces as a function of time for the fiducial near-Eddington case
(HSE_1F_str) and a super-Eddington case (HSE_9F_str). Here we only explicitly
plot the gas pressure gradient force (the red dashed line) and CR pressure
gradient force (the purple dashed line), the sum with other forces including
gravity, magnetic field force and advected momentum is the grey solid line. A
direct difference between near-Eddington and super-Eddington is the inversion
of CR and gas pressure gradient force. We also exam other simulations with
different CR flux, the super-Eddington runs always have a higher CR force than
gas pressure gradient force, and vice versa for the sub-Eddington and near-
Eddington cases.

Despite the different interplay between gas pressure and CR forces, the CR
force (dashed purple lines) follows similar time evolution: it peaks at the
beginning, then relaxes to a lower value in a few $t_{\rm str}$, and slowly
decreases. This is because the “bottleneck” effect relies on gas density
gradient. However, as gas is accelerated and dispersed, the density gradient
becomes shallower and even vanishes, so the CR pressure gradient forces drops
commensurately at later time. The gas pressure gradient is modified by CR-gas
interactions, and contributes to support against gravity. Although the total
force tends to be relatively constant and close to zero at later times, there
is no steady state wind developed in our simulations since the initial gas is
redistributed and slowly depleted. Nevertheless, our results confirm that the
CR Eddington flux is a good criterion for distinguishing the presence of
significant levels of CR driven outflow.

###  3.2 The Effect of CR Transport Mechanism: Streaming versus Diffusion

In this section, we consider the relative importance of CR diffusion and
streaming. Although all the streaming runs reported above have a small nominal
diffusive component ($\kappa_{\rm diff}=10^{-8}$), the CR flux is completely
dominated by streaming. We first consider three hybrid cases where both
streaming and diffusion are present: low diffusivity (HSE_1F_ld $\kappa_{\rm
diff}=1.0$), comparable diffusivity (HSE_1F_cd, $\kappa_{\rm diff}=10.0$), and
high diffusivity (HSE_1F_hd, $\kappa_{\rm diff}=100.0$) but parameters are
otherwise set to match the fiducial run (HSE_1F_str). We also consider two
diffusion (non-streaming) simulations with high diffusivity (HSE_1F_hd_ns,
$\kappa_{\rm diff}=100.0$) and low diffusivity (HSE_1F_ld_ns, $\kappa_{\rm
diff}=1.0$). These non-streaming simulations have $\mathbf{v}_{\rm{s}}=0$ and
$\mathbf{v}_{\rm{A}}=0$ in Equation 1a-1f and Equation 2, so CRs cannot heat
the gas via wave damping.

Figure 6: The density snapshots at $t=10.0$ for different CR transport models.
From left to right, the first two columns are streaming and diffusion, the
diffusivity is $\kappa_{\rm diff}=10.0$(HSE_1F_cd) and $\kappa_{\rm
diff}=100.0$ (HSE_1F_hd). The fourth and fifth columns are pure diffusion
without streaming. The diffusivity are accordingly $\kappa_{\rm diff}=100.0$
(HSE_1F_hd_ns) and $\kappa_{\rm diff}=1.0$ (HSE_1F_ld_ns). The aspect ratio of
plots are adjusted for the purpose of visual comparison with each other, the Y
labels indicate the actual size of calculation domain.

We first look at the diffusion-only simulations. The CR flux $\mathbf{F}_{\rm
c}\sim-\kappa_{\rm diff}\nabla P_{\rm c}$, so $\nabla P_{\rm
c}\sim-\mathbf{F}_{\rm c}/\kappa_{\rm diff}$. If we ignore the time-dependent
term in Equation (1f), the CR contribution to gas momentum at RHS of Equation
(1b) is $\nabla P_{\rm c}$. So for the a given CR flux, larger $\kappa_{\rm
diff}$ gives the smaller acceleration since CRs diffuse through the gas more
quickly without much interaction. The gas density for HSE_1F_hd_ns and
HSE_1F_ld_ns at $t=10.0$ is shown in the third and fourth column in Figure (6)
respectively. In HSE_1F_ld_ns CRs cannot diffuse through the gas rapidly, so a
strong CR pressure gradient develops at the base of the atmosphere, with the
CR force exceeding gravity and lofting the entire atmosphere. Presumably,
there would be feedback on the injected CR flux in a more realistic set-up but
the key result here is that CRs are capable of collectively launching and
accelerating the gas against gravity, with no instabilities (e.g. Parker or
Rayleigh-Taylor like) present to disrupt the shell.

In contrast, HSE_1F_hd_ns experiences a more gentle acceleration. CRs diffuse
through the gas expanding the atmosphere in a manner more similar to the
streaming runs but with very little outflow. This difference between
$\kappa_{\rm diff}=100.0$ and $\kappa_{\rm diff}=1.0$ dynamics also can be
interpreted in terms of the CR Eddington flux. For the injected CR flux
$F_{\rm CR0,in}\approx 15.0$, the low diffusivity ($\kappa_{\rm diff}=1.0$)
run is super-Eddington at the base ($F_{\rm edd}\approx 0.83$), but the high
diffusivity ($\kappa_{\rm diff}=100.0$) run is substantially sub-Eddington
($F_{\rm edd}\approx 82.62$) so the atmosphere relaxes to steady state with a
larger scale height rather than driving outflow.

We turn on both streaming and diffusion in HSE_1F_ld, HSE_1F_cd and HSE_1F_hd.
The gas momentum source term (RHS of Equation 1b) is proportional to the
interaction coefficient $\sigma_{\rm c}$. The first term in the interaction
coefficient is the classical diffusivity, the second term is $\sigma_{\rm
str}^{-1}$ is an effective “diffusivity” for streaming. We estimate it based
on the initial value at the base of atmosphere: $\sigma_{\rm
str,init}^{-1}\approx 4v_{\rm A}/(\partial\ln P_{\rm c}/\partial
x)\approx-6v_{\rm A}H_{\rho}\approx 24.2$ (see similar derivation in Equation
5). In HSE_1F_ld we set $\kappa_{\rm diff}=\sigma_{\rm c}^{\prime-1}=1.0$, so
that the interaction coefficient is dominated by the effective streaming
coefficient, while in HSE_1F_hd $\kappa_{\rm diff}=100.0$, so the diffusivity
dominates. The isotropic diffusivity in HSE_1F_cd is $\kappa_{\rm diff}=10.0$,
which is about the same order of magnitude as $\sigma_{\rm str,init}^{-1}$.

The first two columns in Figure 6 compare the density snapshots for the
streaming-diffusion simulations: HSE_1F_cd and HSE_1F_hd at $t=10.0$. The gas
distribution in HSE_1F_ld (not shown) is almost identical to the streaming
simulation HSE_1F_str with the same CR flux shown in Figure 1. In HSE_1F_ld,
the small diffusivity does not lead to a large CR pressure gradient as the
pure diffusion case because $\sigma_{\rm c}^{-1}$ is now dominated by the
streaming component. By preventing the development of large gradients, CR
streaming limits the CR force when CR diffusivity is low. At the other end
where the diffusivity is large, the gas density profile in HSE_1F_hd is
similar to the non-streaming case HSE_1F_hd_ns with the same diffusivity,
except for the presence of shocks formed from the acoustic instability. In
HSE_1F_cd, where streaming and diffusion are comparable, the result falls
between these limits, but most similar to the HSE_1F_ld run.

When $\kappa_{\rm diff}$ is small, the interaction coefficient is dominated by
the streaming component, and the dynamics is similar to streaming limit.
Increasing $\kappa_{\rm diff}$ will ensure that
$\sigma_{c}^{-1}\sim\kappa_{\rm diff}$ but also lower the interaction
coefficient, which usually leads to a reduction of the CR pressure gradient,
reducing the CR force on the gas. The upper panel of Figure 7 illustrates the
decrease in the CR pressure gradient as $\kappa_{\rm diff}$ increases in the
streaming-diffusion runs at early time ($t=1.5$), approximately when the CR
force peaks.

Figure 7: $y$ direction averaged CR pressure (the first row), total CR flux
(solid lines in the second row), the streaming flux component (dashed lines in
the second row) for $0<x<50$ at $t=1.5$ for the three streaming-diffusion
simulations: HSE_1F_str_ld (red lines), HSE_1F_str_cd (blue lines) and
HSE_1F_str_hd (green lines).

The bottom panel shows the total CR flux (the solid lines) and the streaming
flux (the dashed line) at the same time. Since gas has not yet received a
significant acceleration, the advection flux ($4v_{\rm x}P_{\rm c}$) is small.
Thus, the difference between total CR flux and streaming flux is dominated by
the diffusive flux. When diffusivity is low, the CR flux is dominated by the
streaming flux ($4v_{\rm A}P_{\rm CR}$), but the diffusive flux is significant
wherever the pressure gradient is present and diffusivity is large. Although
changing $\kappa_{\rm diff}$ modifies the balance of streaming and diffusive
flux, the total flux is fixed by our boundary condition, so $\sigma_{\rm c}$
roughly controls the amount of momentum that CR transfers to gas. The upshot
is that for a given CR flux, the overall momentum transfer from CRs to gas is
generally reduced when streaming operates.

###  3.3 The Impact of Cooling

In this section, we study the effect of radiative cooling by comparing three
streaming-dominated simulations with different cooling times. A more detailed
treatment of cooling and the effects of thermal instability on CR driven
outflows will be presented in a companion paper (Huang, Jiang & Davis, in
preparation). We adopt an optically thin radiative model with cooling rate
$n_{H}^{2}\Lambda(T)$ and a constant supplemental heating rate $n_{H}\Gamma$.
The cooling function $\Lambda(T)$ is taken from Wiener et al. (2019), which is
an approximation to fit the CLOUDY data from Wiersma et al. (2009):

$\Lambda(T)=1.1\times 10^{-21}\times 10^{\Theta(\log(T/10^{5}\rm K))}\rm
erg~{}cm^{3}~{}s^{-1},$ (8)

where

$\Theta(x)=0.4x-3+\frac{5.2}{e^{x+0.08}+e^{-1.5(x+0.08)}}$ (9)

We cut off the radiative cooling for gas with temperature below $10^{4}$K.

Because cooling is proportional to $n_{H}^{2}$, the cooling time significantly
increases as the gas density decreases with height. The initial cooling
timescale at two scale heights is larger than the cooling timescale at the
base by more than two orders of magnitude. We report our simulations in terms
of $t_{\rm cool}/t_{0}$, the ratio of the cooling time to sound crossing time
at the base of the simulation, but this will overestimate the importance of
cooling higher in the simulation where densities are lower. If $t_{\rm
cool}<t_{0}$ we expect cooling to impact the gas pressure gradient, but if
$t_{\rm cool}>t_{0}$, cooling may not be important for the gas dynamics.

In contrast to our non-radiative simulations, introduction of cooling fixes
the density and temperature of the simulation, removing the scaling freedom
that was previously present. HSE_1F_str_wc has $t_{\rm cool}/t_{0}=0.7$
(initially) at the base. The density scaling unit is $\rho_{0}=2.35\times
10^{-27}\rm g~{}cm^{-3}$, the length unit is $l_{0}=1.54\times 10^{20}\rm cm$,
and temperature unit $T_{0}=9.9\times 10^{4}\rm K$. For HSE_1F_str_sc and
HSE_9F_str_sc the density scaling unit is $\rho_{0}=2.35\times 10^{-25}\rm
g~{}cm^{-3}$, the length unit is $l_{0}=1.54\times 10^{20}\rm cm$, and the
temperature unit $T_{0}=9.9\times 10^{3}\rm K$, yielding $t_{\rm
cool}/t_{0}=0.06$ at the base. However, we continue to report results in
dimensionless to facilitate comparison to non-radiative simulations.

Figure 8: The $y$ direction averaged gas density (the first row), gas pressure
(the second row), CR pressure (the third row) for the marginal-Eddington flux:
HSE_1F_str (blue solid lines, no cooling), HSE_1F_str_wc (orange solid lines,
weak cooling) and HSE_1F_str_sc (green solid lines, strong cooling), and the
super-Eddington flux HSE_9F_str (red dashed lines, no cooling) and
HSE_9F_str_sc (purple dashed lines, strong cooling) The snapshots are sampled
at $t=20$.  Figure 9: The mass outflow rate per unit area (the upper panel)
and average velocity (the lower panel) for the marginal-Eddington flux:
HSE_1F_str (blue solid lines, no cooling), HSE_1F_str_wc (orange solid lines,
weak cooling) and HSE_1F_str_sc (green solid lines, strong cooling), and
super-Eddington flux HSE_9F_str (red dashed lines, no cooling) and
HSE_9F_str_sc (purple dashed lines, strong cooling). Comparing the simulations
without cooling, with short cooling timescale and long cooling timescale, the
three marginal-Eddington runs shows similar outflow properties, the two super-
Eddington runs have different mass outflow rate, but similar velocity.

Figure 8 shows snapshots of horizontally ($y$) averaged gas density (top
panel), CR pressure (middle panel) and gas pressure (bottom panel) at
$t^{\prime}=20$. The blue solid line shows the profile of the non-radiative
run HSE_1F_str, and the other two solid lines are the radiative simulations.
The red dashed lines shows the non-radiative super-Eddington flux simulation
HSE_9F_str, and the purple dashed line is the corresponding radiative
simulation: HSE_9F_str_sc.

For simulations with near-Eddington CR flux the orange lines show the profiles
from the radiative run HSE_1F_str_wc. The gas pressure at the base is lower
due to the cooling, but because the cooling is relatively moderate, it has
limited effect on the overall dynamics. The gas density and CR pressure are
similar to non-radiative simulations.

The green lines are the profiles from HSE_1F_str_sc, where the short cooling
time significantly alters the gas profile at the base compared to HSE_1F_str.
Since the internal energy is radiated away, the gas lacks pressure support and
falls toward the base, forming a much sharper density gradient and is less
dispersed compared to the non-radiative and moderate cooling runs. The sharper
density gradient increases the bottleneck effect for CR streaming, leading to
a stronger CR force on the gas.

The horizontally averaged gas density profile is more strongly clumped in both
radiative runs. The shocks created by acoustic instability are enhanced by
cooling. Since the shocks are relatively isothermal, the enhanced density in
the shock increases the cooling rate, lowering gas pressure relative to the
shocks in non-radiative runs. The relatively lower pressure support allows the
gas to condense more effectively.

Despite the differences in the gas density profiles, the outflow properties
from the three near-Eddington runs shows in Figure 9 are similar. The outflow
we observed in the simulations with near-Eddington CR flux are composed
primarily of low density gas, where the cooling time is significantly longer
than dynamical timescales. An interesting contrast is provided the super-
Eddington simulation HSE_9F_str_sc, which also has a short cooling time and
corresponds to the purple dashed line in Figure 8 and Figure 9.

In the lower panel of Figure 8, the gas pressure is significantly reduced by
cooling and the gas pressure gradient is inverted compared to the initial
profile. However, the gas density profile is similar to its non-radiative
counterpart simulation HSE_9F_str (red dashed line) because the CR force
dominates over the gas pressure force. The momentum transfer between CRs and
gas that primarily determines the gas density and velocity distribution in
super-Eddington systems is only indirectly impacted by radiative cooling.
Figure 9 shows that the mass outflow rates are reduced by cooling but the
effect is less than 50% and velocities are relatively unimpacted. The upshot
is that cooling is likely to be important in setting density profiles of near-
Eddington and sub-Eddington atmospheres, but is probably less important in the
super-Eddington cases, where CR forces are dominant.

###  3.4 The Effect of Simulation Parameters

In this section, we briefly discuss a number of factors in our setup and their
impact on the numerical results, including the simulation dimensions and
aspect ratio in Section 3.4.1, along with resolution, dimensionality, and
choice of $V_{\rm max}$ in Section 3.4.2.

####  3.4.1 Aspect Ratio and Domain Size

Some of the simulations listed in Table 1 are performed with different aspect
ratios. While the fiducial box size is $L_{x}\times L_{y}=400\times 50$, we
have some simulations with longer $L_{x}$ or shorter $L_{y}$. First we discuss
the effect of varying $L_{x}$, then the effect of varying $L_{y}$.

We increase $L_{x}$ in streaming dominated runs with larger CR flux, because
the larger acceleration requires large box size to capture the gas - CR
interactions. For the higher CR flux runs, the efficiency of momentum and
energy transfer between CR and gas depends on $L_{x}$. We define the $x$
direction momentum transfer efficiency $\dot{p}_{\rm{\rightarrow
g}}/\dot{p}_{x,\rm{CR}}$ and energy transfer efficiency
$\dot{E}_{\rm{\rightarrow g}}/\dot{E}_{\rm CR}$ as:

$\displaystyle\dot{p}_{x,\rm{\rightarrow
g}}/\dot{p}_{x,\rm{CR}}=\frac{\int\sigma_{\rm
c}\hat{\mathbf{n}}_{x}\cdot[\mathbf{F}_{\rm c}-\mathbf{v}\cdot(E_{\rm
c}\textsf{I}+\textsf{P}_{\rm c})]dV}{\int_{\rm in}\mathcal{F}_{F_{\rm
CR,x}}dA}$ (10) $\displaystyle\dot{E}_{\rm{\rightarrow g}}/\dot{E}_{\rm
CR}=\frac{\int(\mathbf{v}+\mathbf{v}_{\rm{s}})\cdot(\nabla\cdot\textsf{P}_{\rm
c})dV}{\int_{\rm in}\mathcal{F}_{E_{\rm CR}}dA}$ (11)

where $\mathcal{F}_{F_{\rm CR,x}}$ and $\mathcal{F}_{E_{\rm CR}}$ are the flux
of CR momentum in $x$ direction and CR energy returned by Riemann solver. So
the $x$ direction momentum (or energy) transfer efficiency describes what
fraction of total injected $x$ direction CR momentum (or CR energy) is given
to the gas in the domain.

Figure 10: The effect of box size to momentum transfer efficiency. The lines
are time evolution of the ratio between total momentum transferred to gas from
CR and total injected CR momentum for HSE_hF_str (green lines), HSE_1F_str
(red lines), HSE_4F_str (orange lines) and HSE_9F_str (blue lines). The line
style shows the domain width $L_{x}$, where solid lines are for $L_{x}=200$,
dashed lines are for $L_{x}=400$ (fiducial), and dotted lines are for
$L_{x}=800$.

Figure 10 shows the momentum transfer efficiency a function of time for
different choices of box length $L_{x}$. In general, $L_{x}$ has a relatively
insignificant effect on the simulations with sub-Eddington fluxes (red and
green curves). However, the box size largely impacts the result of super-
Eddington cases. For example in HSE_9F_str (blue curves), both momentum and
energy transfer efficiency between CR and gas increases significantly as we go
from shorter ($L_{x}=200$) to longer ($L_{x}=800$) boxes. the larger domain
allows more time for the CR forces to act on the gas, transferring a larger
momentum fraction.

Figure 11: The gas density snapshots from HSE_1F_cd at $t=10.0$. From left to
right: $L_{y}=5,25,50$. The aspect ratio of the third row is adjusted for the
purpose of visual comparison, the $y/l_{0}$ labels denote the actual width.

For our nearly plane-parallel geometry and nearly vertical magnetic field, the
width of the simulation domain $L_{y}$ seems to have negligible results on the
overall dynamics. Figure 11 shows the gas density snapshots at $t=10.0$ for
HSE_1F_cd, in which both streaming and diffusion are operating with
$\kappa_{\rm diff}=10.0$. There is a modest impact on the growth of the
acoustic instability, but the differences have insignificant impact on the
dynamics in the vertical direction.

####  3.4.2 Resolution, Dimensionality and $V_{\rm max}$

Figure 12 shows a comparison of results at $t=10$ from high resolution
(HSE_1F_str_HR), moderate resolution (HSE_1F_str), and low resolution
(HSE_1F_str_LR) runs. We fix the vertical and horizontal dimensions,
increasing the number of zones in both directions by a factor of two for each
increased resolution run. Aside from a modest effect on the development of the
acoustic instability, the runs are insensitive to the resolution adopted here.
The momentum and energy transfer efficiencies between gas and CR are also
insensitive to resolution. Note that the amplitude of the shocks depends on
resolution, with lower amplitude in the HSE_1F_str_LR run, and nearly
disappearing in even lower resolution test runs (not shown). But, the
amplitude seems to be converging between HSE_1F_str_HR and HSE_1F_str.

Figure 12: The $y$ direction averaged gas density (upper panel) and $x$
direction velocity (lower panel) of HSE_1F_str_LR (blue), HSE_1F_str (orange),
HSE_1F_str_HR (green), HSE_1F_str_VM (red). HSE_1F_str_3D (purple) is the 3D
simulation that is performed with the same resolution as HSE_1F_str_LR, in
which we show the profile averaged throught $y$ and $z$ directions.

To study the impact of our two-dimensional assumption, we also perform a
three-dimensional simulation, which is shown as the purple curve in Figure 12.
As with resolution, there is no evidence of significant dependence on
dimensionality, at least for the streaming simulations and magnetic field
geometry employed here.

Finally, we consider the dependence on $V_{m}$, which acts like a reduced
speed of light for the time dependent term in Equation 1f. This approximation
allows us to take larger timesteps, which scale inversely with $V_{\rm m}$.
Previous tests have indicated that this is a reasonable approximation as long
as $V_{\rm m}$ is larger than other physical velocities in the simulation. In
simulations presented here, the maximum physical speed tends to be the Alfvén
speed in the diffuse background medium, which is about $v_{\rm A,max}\approx
224$, while we use $V_{\rm m}=1000$ in all the simulations except
HSE_1F_str_VM, where $V_{\rm m}=2000$. In our dynamical problems, the time
dependent term $(1/V_{\rm m}^{2})\partial\mathbf{F}_{\rm c}/\partial t$ in
Equation 1f is much smaller than the other terms so we expect our choice of
$V_{\rm m}$ should not impact our results. Jiang & Oh (2018) point out that
the main physical effect of using $V_{\rm m}$ is overestimation of the
momentum carried by CRs by $c^{2}/V_{\rm m}^{2}$, but this has little impact
on the CR force, which determines momentum transfer to the gas.

The red curve in Figure 12 shows the density snapshot of HSE_1F_str_VM, also
at $t=10.0$. The gas distribution and velocity are similar to others, with
only modest differences in the development of the acoustic instability. The
momentum and energy transfer efficiency agree with the $V_{\rm m}=1000$
simulations.

##  4 Discussion

###  4.1 Momentum and Energy Transfer

Our simulation results described in section 3 explore a number of factors that
impact the momentum and energy exchange between CRs and gas. To further
quantify this exchange, we define the momentum efficiency $\dot{P}_{\rm
x,togas}/\dot{P}_{\rm x,CR}$ and energy efficiency $\dot{E}_{\rm
togas}/\dot{E}_{\rm CR}$ in Equation 10 and Equation 11 respectively by taking
the ratio of total momentum (energy) transferred from CRs to gas and the total
injected CR momentum (energy). Table 2 tabulates the the time averaged
efficiencies at a point when these simulations roughly reach quasi-steady
state.

Table 2: Energy and Momentum Budgeta

Name | $\dot{P}_{x,\rm{\rightarrow g}}/\dot{P}_{\rm x,CR}$ | $\dot{E}_{\rm{\rightarrow g}}/\dot{E}_{\rm CR}$ | $\dot{E}_{\rm{\rightarrow g,KE}}/\dot{E}_{\rm{\rightarrow g}}$ | $\dot{E}_{\rm{\rightarrow g,H}}/\dot{E}_{\rm{\rightarrow g}}$
---|---|---|---|---
HSE_hF_str | $98.9\%$ | $65.0\%$ | $4.1\%$ | $95.7\%$
HSE_1F_str | $98.4\%$ | $59.2\%$ | $7.1\%$ | $92.9\%$
HSE_4F_str | $96.9\%$ | $52.0\%$ | $16.7\%$ | $82.3\%$
HSE_9F_str | $93.8\%$ | $44.2\%$ | $24.1\%$ | $74.9\%$
HSE_20F_str | $77.3\%$ | $29.1\%$ | $35.1\%$ | $64.9\%$
HSE_1F_ld | $96.2\%$ | $52.9\%$ | $6.9\%$ | $93.1\%$
HSE_1F_cd | $96.7\%$ | $53.7\%$ | $6.4\%$ | $93.6\%$
HSE_1F_hd | $96.7\%$ | $52.7\%$ | $8.3\%$ | $91.7\%$
HSE_1F_hd_aniso | $96.8\%$ | $52.6\%$ | $8.3\%$ | $91.7\%$
HSE_1F_hd_ns | $99.9\%$ | $27.9\%$ | $100.00\%$ | -b
HSE_1F_ld_ns | $99.9\%$ | $18.6\%$ | $100.00\%$ | -b

  * a

The energy and momentum exchange efficiency is averaged over
$t=20.0t_{0}-25.0t_{0}$. The time interval is as when the simulation reaches
quasi-steady exchange rate. The simulations in this table correspond to the
$L_{y}$ reported in Table 1 to avoid significant CR momentum and energy
outflow.

  * b

Because streaming absent, there is no direct CR heating.

The first column of Table 2 lists the momentum exchange efficiencies. In most
of the simulations, $\gtrsim 95\%$ of momentum carried by the injected CR flux
is transferred to gas within the simulation box, whether the flux is
transported primarily by streaming or diffusion. The exception is HSE_20F_str,
in which the momentum efficiency is $\sim 77\%$. The relatively low efficiency
is due to rapid acceleration and subsequent outflow produced in this
simulation. Figure 10 shows that small boxes cannot adequately capture the
acceleration region so the expectation is that the CR flux escaping through
the upper boundary would continue to accelerate the gas if the simulation
domain were extended.

Other effects that would limit coupling, such as preferential escape along low
impedance channels appears to be absent. This contrasts with radiatively
accelerated outflows, where the impact of the radiative Rayleigh-Taylor
instability likely limits the coupling between radiation and dusty gas by
creating low density channels through with radiation can preferentially escape
(Krumholz & Thompson, 2012; Davis et al., 2014). It is not entirely clear
whether similar effects would be expected to be present in CR driven outflows,
but the non-linear development of instabilities such as the CR acoustic
instability (Begelman & Zweibel, 1994) or Parker instability (Parker, 1966;
Ryu et al., 2003; Heintz & Zweibel, 2018) could plausibly act in a fashion
similar to Rayleigh-Taylor instability in radiatively driven flows. Indeed, we
attribute the shocks seen in our simulations to the CR acoustic instability,
as discussed further in section 4.3 below. In these simulations, however, the
growth of the instability does not seem to be modifying the momentum coupling,
although it concentrates heating near the pressure jumps. Nor do we see any
evidence of the Parker instability, although this is likely impacted by our
adoption of a nearly vertical field in the acceleration region. Other
simulation results suggest this instability may play an important role in
overall structure of the gas disk (Heintz et al., 2020).

The transfer of energy between the CRs and gas is more complex. The third
column in Table 2 lists the time-averaged total energy transfer efficiency,
and the fourth and fifth column list the fractions of kinetic energy and CR
heating respectively. As with momentum transfer, we find that the CR energy
transfer is sensitive to our domain size, with outflow from the top boundary
carrying away a significant CR energy flux that is expect to generate further
heating in larger boxes. Hence, it is difficult to judge the overall level of
coupling, which will inevitably depend on the global structure of the flow and
galaxy potential. Nevertheless, we can draw conclusions about the relative
importance of CR heating to CR acceleration.

Comparison of the third and fourth columns shows that for streaming runs, the
majority of the energy transfer is due to CR heating. This is true even in
hybrid simulations where diffusion is significant. In contrast, there is no
explicit CR heating in the diffusion runs, due to the assumptions of our CR
transport model and all energy transfer increases the kinetic energy of the
gas.

Figure 13: The energy transfer efficiency as a function of time for
simulations with diffusion. The total energy transferred to the gas from CR
includes heating (the solid line) and kinetic energy (the dashed line).
HSE_1F_str (red) is the streaming simulation with $\kappa_{\rm diff}=10^{-8}$.
HSE_1F_ld (black), HSE_1F_cd (light blue) and HSE_1F_hd (green) are streaming-
diffusion simulations with $\kappa_{\rm diff}=1.0,~{}10.0,~{}100.0$
respectively. HSE_1F_ld_str (brown) and HSE_1F_str (orange) are pure diffusion
simulations with $\kappa_{\rm diff}=1.0,~{}100.0$ respectively.  Figure 14:
The energy transfer efficiency as a function of time for simulations with
diffusion. The total energy transferred to the gas from CR includes heating
(the solid line) and kinetic energy (the dashed line). These are all streaming
simulations with different CR flux: HSE_hF_str (green), HSE_1F_str (red),
HSE_4F_str (orange) and HSE_9F_str (blue) and HSE_20F_str (purple).

We show the time evolution of energy transfer components in Figure 13 for
varying CR transport mechanisms. After an initial transient stage on the order
of the Alfvén wave crossing time, the efficiency is relatively constant. The
two non-streaming simulations (brown and yellow dashed curves) have only a
kinetic energy term, and share similar transfer efficiency corresponding to
about $20\%-30\%$ of total injected CR energy. The CR heating is weakly
sensitive to the diffusivity, showing a trend towards increased heating
efficiency when CR diffusivity is lowered.

When streaming dominates the energy transfer, the overall efficiency and
partition between kinetic energy and heating depends on the injecting CR flux
as shown in Figure 14. A higher injected CR flux leads to smaller ratio
between heating and kinetic energy. The kinetic energy term is proportional to
fluid velocity $v_{\rm x}$, while the heating term is proportional to local
Alfvén velocity. The gas is accelerated to higher velocity in the larger CR
flux simulations, so the kinetic energy term rises. In our highest CR flux run
(HSE_20F_str), the kinetic energy and heating are roughly comparable to each
other.

The end result is that CR acceleration of outflows is most efficient when
heating is limited. This is trivially achieved in our non-streaming runs where
heating is absent. If streaming operates, increasing the injected CR flux
accelerates the gas more rapidly, lowering the ratio of gas velocity to the
streaming (Alfvén) velocity. This highlights the importance of understanding
the detailed dynamics of CR transport in star-forming galaxies and, if
streaming is important, introduces stronger dependence on CR flux than might
be anticipated by a simple comparison of the CR flux to our derived Eddington
flux.

###  4.2 Implications for CR Feedback in Galaxies

As discussed in section 2.3.1, the non-radiative simulations are scaled to
physical units by choosing $\rho_{0}$, $l_{0}$ and $T_{0}$. In this section,
we vary the scaling to infer properties of CR outflows in different galactic
environments.

Following Socrates et al. (2008), we estimate the CR flux $F_{\rm CR}$ from
galaxies based on the surface star formation rate $\dot{\Sigma}_{\star}$ and
scaling the CR production as a fraction of the expected supernova rate. Here
we adopt the scaling of Crocker et al. (2021):

$\displaystyle F_{\rm CR}$ $\displaystyle\approx$ $\displaystyle 5.6\times
10^{47}\frac{\rm erg}{M_{\odot}}\Big{(}\frac{E_{\rm SN}}{10^{51}\rm
erg}\Big{)}\Big{(}\frac{90M_{\odot}}{M_{\rm\star,SN}}\Big{)}\Big{(}\frac{\dot{\Sigma}_{\star}}{M_{\odot}~{}\rm
pc^{-2}\rm Myr^{-1}}\Big{)}$ (12) $\displaystyle=$ $\displaystyle 1.85\times
10^{-3}\Big{(}\frac{\dot{\Sigma}_{\star}}{M_{\odot}~{}\rm pc^{-2}\rm
Myr^{-1}}\Big{)}\rm erg~{}cm^{-2}~{}s^{-1},$

This assumes a Chabrier (2005) initial mass function so that, on average,
there is a supernova every $90M_{\odot}$ stellar mass. Each supernova event
roughly converts $\eta_{\rm SN}\sim 10\%$ of kinetic energy $E_{\rm SN}$ into
CR energy. Assuming these assumptions are correct, this flux is an upper limit
because it ignores any CR losses that may present, such as those resulting
from pion production after collision with ISM atoms or molecules.

In Section 2.2, we derived an estimate for the CR Eddington flux by balancing
the CR force with gravity. We estimate the gravity $g$ in Equation 5 according
to disk surface density as $g=2\pi G\Sigma_{\rm gas}/f_{\rm gas}$, where
$f_{\rm gas}$ is gas fraction. We also define a length scale $z_{\rm disk}$
that estimates the average density at atmosphere base by disk surface density
$\rho\sim\Sigma_{\rm gas}/z_{\rm disk}$. Then we can re-write the CR Eddington
flux as:

$\displaystyle F_{\rm Edd,str}$ $\displaystyle\approx$ $\displaystyle
2.67\times 10^{-7}\Big{(}\frac{z_{\rm disk}}{100\rm
pc}\Big{)}^{\frac{1}{2}}\Big{(}\frac{\Sigma_{\rm gas}}{M_{\odot}\rm
pc^{-2}}\Big{)}^{\frac{3}{2}}$ (13) $\displaystyle\times$
$\displaystyle\Big{(}\frac{B}{1\mu G}\Big{)}\Big{(}\frac{0.5}{f_{\rm
gas}}\Big{)}\rm erg~{}cm^{-2}~{}s^{-1}$ $\displaystyle F_{\rm Edd,diff}$
$\displaystyle\approx$ $\displaystyle 5.93\times
10^{-8}\Big{(}\frac{\Sigma_{\rm gas}}{M_{\odot}\rm
pc^{-2}}\Big{)}\Big{(}\frac{z_{\rm disk}}{100\rm pc}\Big{)}$ (14)
$\displaystyle\times$ $\displaystyle\Big{(}\frac{0.5}{f_{\rm
gas}}\Big{)}\Big{(}\frac{\kappa_{\rm diff}}{10^{26}\rm cm^{-2}s}\Big{)}\rm
erg~{}cm^{-2}~{}s^{-1}$

We can estimate the surface density and star formation rates where a galaxy
might be expected to exceed the CR Eddington limit and drive outflow by
setting Equation 12 equal to either Equation 13 or Equation 14 when streaming
or diffusion (respectively) dominates. The results are shown as magenta and
green solid lines in the upper panel of Figure 15. The magenta solid line is
for streaming with $f_{\rm gas}=0.5$, $z_{\rm disk}=100$pc and $B=1\mu G$. The
green solid line is for diffusion with $f_{\rm gas}=0.5$, $z_{\rm disk}=100$pc
and $\kappa_{\rm diff}=10^{26}\rm cm^{-1}s$. The upper left region of each
line is where the CR flux is above the Eddington limit, and vice versa. We
also show normal galaxies measured in Kennicutt Jr & De Los Reyes (2021) and
de los Reyes & Kennicutt Jr (2019) in Figure 15.

We see that many galaxies lie near or above these estimates for the Eddington
limit. This could be an indication that Equation 13 is overestimate the CR
flux in the wind launching region due to CR losses during transport through
the gas disk. Alternatively, it may mean that these galaxies are all driving
outflows due to CRs but that the mass loading factors are not large enough to
completely quench star formation. Therefore, it is useful to estimate how much
mass outflow is consistent with our simulations.

Figure 15: Upper panel: Applying CR Eddington flux to normal galaxies. The
magenta solid line gives the estimated limit where the CR flux from the galaxy
(Equation 12) equals to CR Eddington flux in streaming limit (Equation 13),
and the green solid line is in diffusion limit (Equation 14). The black dashed
line shows the classical Kennicutt-Schmidt relation (Kennicutt Jr, 1998). The
non-filled circles are observational data from de los Reyes & Kennicutt Jr
(2019) and Kennicutt Jr & De Los Reyes (2021). The filled data points are
streaming simulations with different scaling units. The shape separates data
groups with different surface star formation rate and surface gas density.
Each group includes four data points corresponding to four different injecting
CR flux, red, orange, blue and purple are $F_{\rm CR}=15,~{}60,~{}135,~{}300$
respectively. Lower panel: The mass loading factor $\eta$ as a function of gas
surface density. For each column of data points connected by a vertical line,
from top to bottom corresponds to different temperature scaling
$T_{0}=300K,~{}10^{3}K,~{}5\times 10^{3}K,~{}10^{4}\rm K$. The filled symbols
satisfies geometrical constrains in our model, and the non-filled symbols do
not satisfies the geometrical constrains .

In order to place our simulation results in this plot, we scale our non-
radiative simulations to different units. Once we fix the injecting CR flux,
each data point in the $\Sigma_{\rm gas}-\dot{\Sigma}_{\star}$ plane provides
two constraints on the three free parameters: density scale $\rho_{0}$, length
scale $l_{0}$, and temperature scale $T_{0}$. We adopt the convention of
choosing the temperature and using $\Sigma_{\rm gas}$ and
$\dot{\Sigma}_{\star}$ to fix $l_{0}$ and $\rho_{0}$ for given $T_{0}$.
Additionally, we consider two geometrical constrains in our simulation and
models that limit $l_{0}$ and $\rho_{0}$. First, our initial atmosphere with
scale height of $2l_{0}$ describes a fraction of disk vertical profile, so
$2l_{0}<z_{\rm disk}$. Second, the density at the base atmosphere should not
be exceed the average disk density $\rho_{0}\lesssim\Sigma_{\rm gas}/z_{\rm
disk}$.

We also define the mass loading factor
$\eta\equiv\dot{M}_{\rm{out}}/\Sigma_{\rm SFR}$, where $\dot{M}_{\rm{out}}$ is
the outflow mass per area. We show a series of simulations with the same
$\rho_{0}$ and $l_{0}$ but different $T_{0}$ in the lower panel of Figure 15.
We tested four $T_{0}=300K,~{}10^{3}K,~{}5\times 10^{3}K,~{}10^{4}\rm K$, and
show them from top to bottom in each column. As discussed below, lowering
$T_{0}$ yield higher $\eta$, but some of lowest $T_{0}$ models do not
satisfies the geometry constrains, so we show them as non-filled symbols.

Comparing the filled symbols, the Eddington limit seems to provide an
approximate estimation the of CR’s ability to drive a significant mass-loaded
outflow in galaxies. Galaxies below the CR Eddington limit mostly do not yield
$\eta\gtrsim 1$ outflow, while the super-Eddington systems tend to be more
efficient. If we focus on the internally consistent filled symbols we see that
it is challenging to drive outflows with large mass loading efficiency in the
higher surface density galaxies for the assumed $z_{\rm disk}=100$pc galaxy
model.

Figure 16: Similar as Figure 15, but modelled for compact star bursting
regions.

For the higher surface density star forming galaxies, with more compact star
bursting regions we show results from a model with $z_{\rm disk}$ in Figure
16, comparing to observed starburst regions from Kennicutt Jr & De Los Reyes
(2021). Comparing to the less compact galaxy model in Figure 15, the streaming
CR Eddington limit (the magenta line) is slightly lower in high surface
density region. The diamond data points in both groups shares similar
$\Sigma_{\rm gas}$, but the compact model implies $\eta\gtrsim 1$ outflows are
possible at lower star formations rates.

It is important to emphasise that the density at the base of the atmosphere
has a large effect on the implied mass loading in all models. The
dimensionless CR flux in our simulations are scaled by
$\rho_{0}v_{0}^{3}\propto\rho_{0}T_{0}^{3/2}$ so matching the model galactic
CR fluxes fixes this scaling. The highest densities correspond to the lowest
temperatures, with the highest mass loading $\eta$ corresponding to
$T_{0}\lesssim 10^{3}$K. This is because a lower $T_{0}$ corresponds to a
relatively low characteristic velocity. The increase in density more than
offsets the decrease in velocity ($v_{0}\propto T_{0}^{1/2}$), leading to
larger implied mass outflow. The assumption of lower temperatures is broadly
consistent with the larger fractions of dense molecular gas in LIRGS and
ULIRGs but poses a self-consistency problem for the CR transport models
adopted here, which assume the gas is sufficiently ionised so that CRs stream
at the Alfvén speed and have diffusivities similar to the Milky Way. If
ionisation fractions are low, the effective Alfvén speed may be significantly
higher. Furthermore, the streaming instability or turbulent fluctuations that
underly diffusive models may be strongly damped in molecular gas (Bustard &
Zweibel, 2020). This may increase the escape of CRs from highest density
regions and reduce the impact of destruction processes like pion production,
but could also lower the coupling between CRs and outflowing gas unless
sufficiently high ionisation fractions are present in the wind launching
regions of the galaxy, which are presumably closer to the disk surface.
Resolving these questions is beyond the scope of this work, but is central to
addressing the question of how much CR feedback can impact galaxies with high
molecular gas fractions.

###  4.3 Acoustic Instability

In all our simulations with streaming, a series of shocks grow and saturate
while gas is being accelerated by CRs. For example, the shocks in the gas
density and velocity in Figure 2, or the horizontal bands in the 2D gas
density snapshots in Figure 3 and Figure 6. These shocks do not form in our
diffusive only simulations. We find that the growth of these shocks are
consistent cosmic rays acoustic instability, and our present whenever the
resolution is sufficiently high.

Begelman & Zweibel (1994) studied the CR acoustic instability in the streaming
limit with a uniform background. They found that the growth rate varies with
$m\equiv v_{\rm A}/c_{\rm s}$ and $c\equiv c_{\rm CR}/c_{\rm s}$, where
$v_{\rm A}$, $c_{\rm s}$ and $c_{\rm CR}$ are Alfvén speed, sound speed and CR
acoustic speed (see definitions in Appendix A). Tsung et al. (2021 in prep)
extends this work and discusses the physical origin of the shock structures
and various relevant microphysical processes. In this section, we refer to the
gradient in gas density, gas pressure and CR pressure as the “background
gradient”, and derive the dispersion relation in Appendix A.

The emergence and growth of shocks in our simulation is spontaneous and
dynamical. They usually appear after the gas density and pressure are
redistributed to the larger scale height compared to the initial condition.
The shocks are approximately even-spaced and have similar amplitude. We
measure the average shock amplitude increment over time to estimate the
acoustic instability growth rate.

We first filter the high spatial frequency shocks and fit a smooth curve to
the $y$ direction averaged gas density using the LOWESS (Locally Weighted
Scatterplot Smoothing) method. Then we take the weighted volume average
residual as the shock amplitude $\delta\rho^{\prime}$. We show the resulting
amplitude evolution as the solid curves in Figure 17 for three simulations.
The shocks experience an initial linear growth, which we fit with the dashed
lines, before dropping at later time. The drop in the average amplitude at
later time reflects that the shocks begin to be advected out of the domain
with the fluid.

Figure 17: The growth rate of the shocks caused by acoustic instability for
HSE_1F_str (blue), HSE_4F_str (orange) and HSE_1F_ld (green). The dashed line
is the fitted growth rate, with slopes noted in the legend. The estimated
slopes from solving dispersion relation are 0.2, 0.5 and 0.8 respectively.

We solve dispersion relation of the CR acoustic instability (see Appendix A)
to obtain theoretical growth rates and compare them with the above
measurements. Our calculation is based on the sampling relevant $y$ direction
averaged quantities at the time and location where we first see the shocks.
The calculation for HSE_1F_str is based on the profiles at $t=6.7t_{0}$ and
$x=33.0l_{0}$, HSE_4F_str is based on the profiles at $t=3.6t_{0}$ and
$x=25.0l_{0}$, and HSE_1F_ld uses profile at $t=3.2t_{0}$ and $x=25.0l_{0}$.
The estimated growth rates are $0.2$, $0.5$ and $0.8$ for HSE_1F_str,
HSE_4F_str and HSE_1F_ld respectively (the corresponding slopes in Figure 17).
The good correspondence between estimated growth rates and measured growth
rates is consistent with shocks being seeded by the acoustic instability.

The expected growth rates in our simulation are relatively low due to the low
ratio between CR pressure and gas pressure in the diffuse region. Figure 19
shows the growth rate for the $k=1$ mode, the yellow circles label some
typical $c_{\rm CR}/c_{\rm{s}}-v_{\rm A}/c_{\rm{s}}$ pairs in our simulations
near the position and time the shocks starts to grow. The “bottleneck” slows
down CR propagation in dense gas and yields a low CR pressure in diffuse
region, while the gas is heated by the initial CR shock. As a result, the low
$c_{\rm CR}/c_{\rm cs}$ and high $v_{\rm A}/c_{\rm s}$ constrains the growth
rate to low level in the simulations.

The instability eventually saturates and the shock amplitudes are relatively
constant over time. The shocks spacing increase as the pressure scale heights
increase (Tsung et al. 2021 in prep)), and advects with the gas. The
instability creates multiple “stairs” in CR pressure. The steep pressure
gradients at shock fronts dominate the CR heating at later time when the over
CR pressure gradient is shallow.

##  5 Summary and Conclusions

Many studies have identified CR pressure and heating as promising mechanisms
for driving galactic outflows and ISM turbulence. In this work, we studied the
dynamical interaction due to an imposed CR flux at the base of an initially
hydrostatic atmosphere to investigate the outflow launching process. We
derived estimates for the CR Eddington flux in the streaming limit (equations
5 and 13) and diffusion limit (equations 6 and 14). We performed simulations
with different CR transport models including streaming, streaming-diffusion
and diffusion. In the set of streaming simulations, we vary the injected CR
flux from sub-Eddington to super-Eddington to quantify how this impacts the
atmosphere and outflow that develops. We measured the outflow rate from
simulations and scale simulated systems to realistic star forming
environments. We also consider simulations with radiative cooling, showing
that cooling modifies the wind density and pressure profile, resulting in
generally more inhomogeneous and cooler outflow.

We measured the momentum and energy transfer efficiency from CR to gas, and
summarised how it changes with the injected CR flux and the relative
importance between streaming and diffusion. We summarise our results as
follows:

  * •

Star forming systems near or above the classical Kennicutt-Schmidt relation
are likely to host super-Eddington CR flux (Figure 15 and Figure 16) if a
large fraction of the CRs survive transfer through the gas to reach the
surface. This suggests that CRs potentially contribute to driving outflows in
broad range of environments and star-formation rates. Higher CR flux usually
results in more mass outflow. The typical mass-loading of outflow from the
initially warm atmosphere ($\sim 10^{4}K$) is $\sim 0.1$. Increasing the
density at the base of the atmosphere while lowering the gas temperature leads
to higher mass-loading.

  * •

We identify the growth and saturation of CR acoustic instability in our
streaming simulations. However, the presence of instability does not limit the
momentum transfer from CR to gas. The momentum transfer efficiency is usually
$\gtrsim 95\%$ regardless the CR flux and does not strongly depend on whether
streaming or diffusion dominates transport, suggesting that CR and gas are
well-coupled in our simulations.

  * •

The energy transfer between CR and gas is sensitive to CR flux and the
dominant CR transport mechanism. The fraction of total energy transferred to
gas from CRs increases with CR flux. When streaming is included, a substantial
fraction of transferred energy is in the form of Alfvén wave heating,
resulting in larger pressure scale height and outflowing gas hotter than the
initial temperature by more than one order of magnitude. Typically $\lesssim
10\%$ of the total CR energy is imparted to gas kinetic energy. Pure diffusion
without the Alfvén wave heating leads to roughly twice as much energy being
transferred to the kinetic energy of the outflow.

  * •

For the transport scheme used here, streaming tends to limit the CR force
applied to the gas for a given CR flux. In absence of streaming, a lower
diffusivity requires a larger CR pressure gradient to transport the same CR
flux, leading to a large CR force on the gas. In contrast, when streaming
operates with a low diffusivity and a relatively high streaming speed, it will
dominate the transport of CRs relative to diffusion, providing a weaker CR
pressure gradient and a lower CR force. For a sufficiently high diffusivity,
CR diffusion can again dominate the transport but will then produce a weaker
CR pressure gradient and subsequently weaker force.

  * •

Radiative cooling is potentially important especially in sub-Eddington and
near-Eddington systems in terms of setting the gas temperature and pressure
profile of the extended atmosphere and outflow. When the CR flux is highly
super-Eddington, CR pressure is more important than gas pressure and cooling
has a weaker effect on mass loss. Cooling significantly lowers the wind
temperature, promoting gas condensation near the density irregularities seeded
by acoustic instability, resulting generally cooler and less homogeneous
outflow.

## Acknowledgements

We would like to thank Yan-fei Jiang for sharing the cosmic rays module
implemented in Athena++ as well as for support and throughout the project. We
would also thank Loreto Barcos-Muñoz, S. Peng Oh, Yiqing Song, Stephanie
Tonnesen, Tsun Hin Navin Tsung, and Dong Zhang for helpful discussions and
insights. This work used the computational resources provided by the Advanced
Research Computing Services (ARCS) at the University of Virginia. We also used
the Extreme Science and Engineering Discovery Environment (XSEDE), which is
supported by National Science Foundation (NSF) grant No. ACI-1053575. This
work was supported by the NSF under grant AST-1616171. X. Huang was supported
by the Simons Foundation through the Flatiron Institute’s Pre-Doctoral
Research Fellowship program while a portion of this work was completed.

## Data Availability

The Athena++ code is publicly available at
https://github.com/PrincetonUniversity/athena. The cosmic ray module used here
is not yet part of the public release but was obtained from Dr. Yan-Fei Jiang
who made it available upon request. The initialisation files, simulation
outputs, and data analysis scripts for this study will be shared upon request
to the corresponding author.

##  Appendix A Acoustic Instability Dispersion Relation

We briefly summarise the acoustic instability dispersion relation based on CR
transport equations for the conditions relevant to our simulations. We refer
the reader to Tsung et al. 2021 in prep) for a more thorough discussion of the
instability . To keep the problem tractable, we assume a one-dimensional
geometry with a 1-moment approximation similar to that used in previous work
(Begelman & Zweibel, 1994). The 1-moment and 2-moment scheme gives similar
results in the derived growth rate (Tsung et al. 2021 in prep). We first
define the scale height for the background CR pressure $L_{\rm c}$, gas
pressure $L_{\rm g}$ and gas density $L_{\rho}$ as:

$\displaystyle L_{\rm{c}}^{-1}\equiv-\frac{\partial\ln P_{\rm c}}{\partial
x},\quad L_{\rm{g}}^{-1}\equiv-\frac{\partial\ln P_{\rm g}}{\partial x},\quad
L_{\rho}^{-1}\equiv-\frac{\partial\ln\rho}{\partial x}$ (15)

We also define the scaled diffusivity $l_{\rm c}$

$l_{\rm c}\equiv\frac{\kappa_{\rm diff}}{c_{\rm s}}$ (16)

Using WKB approximation with perturbations $\propto\exp(i\omega t-ikx)$, when
$kL\gg 1$, we can rewrite the perturbed equations as:

$\displaystyle\omega\delta\rho$ $\displaystyle=$
$\displaystyle\bigg{(}1-\frac{i}{kL_{\rho}}\bigg{)}k\rho\delta v$
$\displaystyle\omega\delta v$ $\displaystyle=$
$\displaystyle\frac{k}{\rho}\delta P_{\rm g}+\frac{k}{\rho}\delta P_{\rm
c}-i\frac{\delta\rho}{\rho}g$ $\displaystyle\omega\delta P_{\rm g}$
$\displaystyle=$ $\displaystyle\bigg{(}1-\frac{i}{\gamma_{\rm g}kL_{\rm
g}}\bigg{)}\gamma_{\rm{g}}kP_{\rm g}\delta v$ $\displaystyle+$ $\displaystyle
i(\gamma_{\rm{g}}-1)\frac{P_{\rm c}v_{\rm{A}}}{L_{\rm
c}}\frac{\delta\rho}{2\rho}+(\gamma_{\rm{g}}-1)kv_{\rm{A}}\delta P_{\rm c}$
$\displaystyle\omega\delta P_{\rm c}$ $\displaystyle=$
$\displaystyle\bigg{(}1-\frac{i}{\gamma_{\rm{c}}kL_{\rm
c}}\bigg{)}\gamma_{\rm{c}}kP_{\rm c}\delta v$ (18) $\displaystyle+$
$\displaystyle(kv_{\rm{A}}+i\kappa_{\rm
diff}k^{2}+\frac{i\gamma_{\rm{c}}v_{\rm A}}{2L_{\rho}})\delta P_{\rm c}$ (19)
$\displaystyle-$
$\displaystyle(1+\frac{3i}{2kL_{\rho}}-\frac{i}{\gamma_{\rm{c}}kL_{\rm
c}})\gamma_{\rm{c}}kv_{\rm{A}}P_{\rm c}\frac{\delta\rho}{\rho}$ (20)

The quantities with $\delta$ are the perturbed quantities, $\omega$ is the
frequency and $k$ is the wave number. To compare with Begelman & Zweibel
(1994) and other previous work, these equations differ slightly from the
equations evolved by Athena++ in two ways. Equation (1f) is replaced by a
relation for the CR flux of the form

$\displaystyle F_{\rm c}=(v+v_{\rm A})(E_{\rm c}+P_{\rm c})-\kappa_{\rm
diff}\frac{\partial E_{\rm c}}{\partial x}.$ (21)

This is generally a good approximation since the $\partial F_{\rm c}/\partial
t$ term in equation (1f) is negligible. Furthermore, previous treatments
replace $\partial E_{\rm c}/\partial t$ with $\partial P_{\rm c}/\partial t$
in equation (1e).

We adopt Begelman & Zweibel (1994) notation and define:

$\nu\equiv\frac{\omega}{kc_{\rm{s}}},\quad
m\equiv\frac{v_{\rm{A}}}{c_{\rm{s}}},\quad
c\equiv\frac{c_{\rm{CR}}}{c_{\rm{s}}},$ (22)

where $c_{\rm{CR}}=\sqrt{\gamma_{\rm{c}}P_{\rm c}/\rho}$ is the CR acoustic
speed with $\gamma_{\rm{c}}=4/3$, and $c_{\rm{s}}$ is the sound speed. Solving
the equations gives the dispersion relation:

$\displaystyle(-ikl_{\rm{c}}-m+\nu-\frac{i\gamma_{\rm{c}}}{2kL_{\rho}}m)\bigg{[}\bigg{(}1-\frac{i}{kL_{\rho}}\bigg{)}\bigg{(}1+\frac{g}{kc_{\rm{s}}^{2}\nu^{2}}\bigg{)}\nu^{3}$
$\displaystyle-\nu\bigg{(}1-\frac{i}{kL_{\rm
g}\gamma_{\rm{g}}}\bigg{)}-i\frac{(1-\frac{i}{kL_{\rho}})c^{2}m(\gamma_{\rm{g}}-1)}{2L_{\rm
c}\gamma_{\rm{c}}k}\bigg{]}$ $\displaystyle-c^{2}\nu^{2}(1-\frac{i}{kL_{\rm
c}\gamma_{\rm{c}}})(1+\frac{m}{\nu}(\gamma_{\rm{g}}-1))$
$\displaystyle+\frac{1}{2}m^{2}c^{2}(1-\frac{i}{kL_{\rho}})\bigg{(}1+\frac{3i}{2kL_{\rho}}-\frac{i}{kL_{\rm
c}\gamma_{\rm{c}}}\bigg{(}\frac{\nu}{m}+(\gamma_{\rm{g}}-1)\bigg{)}$
$\displaystyle=0$ (23)

When $kl_{\rm{c}}\gg 1$ and the background is uniform (drop the $L_{\rho}$,
$L_{\rm c}$ and $L_{\rm g}$ terms), without gravity $g$, Equation A is equal
to the Eq 3.18 in Begelman & Zweibel (1994) in 1D form.

Figure 18: The re-scaled growth rate -Im($\nu$) as a function of $v_{\rm
A}/c_{\rm s}$ and $c_{\rm CR}/c_{\rm s}$ for different sets of scale height
$L_{\rm c}=L_{\rho}=L_{\rm g}=50$ (left) and $=10$ (right). In both plots, the
gravity $g=-0.826$, the magnetic field is assumed to be uniform and static.
The white solid lines represent solutions with streaming transport, the yellow
dashed curves are contours represent streaming with moderate CR diffusion.
Figure 19: The re-scaled growth rate -Im($\nu$) for two of the cases shown in
Figure 17: HSE_1F_str (left) and HSE_1F_ld (right). In both plots, the white
contours label the -Im($\nu$)=0.1, 1.0, 2.0 and 3.0. The yellow circles are
typical $m-c$ pair samples relevant to the instability growth, averaged over
$y$ direction.

Figure 18 shows the re-scaled growth rate $-\rm{Im}\nu$ for $k=1$ mode (white
contour) for two sets of scale heights $L_{\rm g}=L_{\rm c}=L_{\rho}=50$
(left) and $L_{\rm g}=L_{\rm c}=L_{\rho}=10$ (right) when streaming dominates,
corresponding to $l_{\rm{c}}=10^{-8}$. The yellow curves are contours for
growth rate calculated with the same set of parameter, but with moderate
diffusion $l_{\rm{c}}=1$. In this parameter space, adding moderate diffusion
has limited impact to the growth rate.

Similarly, Figure 19 shows the re-scaled growth rate solution with two sets of
parameters that approximate the simulation HSE_1F_str (left) and HSE_1F_ld
(right). The yellow circles label $c_{\rm CR}/c_{\rm s}-v_{\rm A}/c_{\rm s}$
pairs that sampled from the diffuse region in the atmosphere at early time,
where the instability emergent. The “bottleneck” limits CR pressure while
strong CR heating preserves gas pressure. As a results, the low $c_{\rm
CR}/c_{\rm cs}$ and high $v_{\rm A}/c_{\rm s}$ constrains the growth rate to
relatively low level in the simulations.

In the left panel, we sample data from $t=6.8t_{0}$ to $t=7.2t_{0}$ with step
$\Delta t=0.1t_{0}$. In each time stamp, the data points are collected at
location of $x=25l_{0},30l_{0},35.0l_{0}$ averaging over the $y$ direction. In
the right panel, the sampling period is $t=3.0t_{0}$ to $t=3.5t_{0}$ with step
$\Delta t=0.1t_{0}$. In each time stamp, the data points are sampled at
$x=20l_{0},25l_{0},30.0l_{0}$, averaging over the $y$ direction.

## References

  * Bai et al. (2019) Bai X.-N., Ostriker E. C., Plotnikov I., Stone J. M., 2019, The Astrophysical Journal, 876, 60
  * Bambic et al. (2021) Bambic C. J., Bai X.-N., Ostriker E. C., 2021, arXiv preprint arXiv:2102.11877
  * Begelman & Zweibel (1994) Begelman M. C., Zweibel E. G., 1994, ApJ, 431, 689
  * Booth et al. (2013) Booth C., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, The Astrophysical Journal Letters, 777, L16
  * Boulares & Cox (1990) Boulares A., Cox D. P., 1990, The Astrophysical Journal, 365, 544
  * Breitschwerdt et al. (1991a) Breitschwerdt D., McKenzie J., Völk H., 1991a, Astronomy and Astrophysics, 245, 79
  * Breitschwerdt et al. (1991b) Breitschwerdt D., McKenzie J. F., Voelk H. J., 1991b, A&A, 245, 79
  * Breitschwerdt et al. (1993) Breitschwerdt D., McKenzie J., Völk H., 1993, Astronomy and Astrophysics, 269, 54
  * Brüggen & Scannapieco (2020) Brüggen M., Scannapieco E., 2020, The Astrophysical Journal, 905, 19
  * Buck et al. (2020) Buck T., Pfrommer C., Pakmor R., Grand R. J., Springel V., 2020, Monthly Notices of the Royal Astronomical Society, 497, 1712
  * Bustard & Zweibel (2020) Bustard C., Zweibel E. G., 2020, arXiv preprint arXiv:2012.06585
  * Bustard et al. (2020) Bustard C., Zweibel E. G., D’Onghia E., Gallagher III J., Farber R., 2020, The Astrophysical Journal, 893, 29
  * Butsky & Quinn (2018) Butsky I. S., Quinn T. R., 2018, The Astrophysical Journal, 868, 108
  * Butsky et al. (2020) Butsky I. S., Fielding D. B., Hayward C. C., Hummels C. B., Quinn T. R., Werk J. K., 2020, The Astrophysical Journal, 903, 77
  * Chabrier (2005) Chabrier G., 2005, The Initial Mass Function 50 Years Later, pp 41–50
  * Chevalier & Clegg (1985) Chevalier R. A., Clegg A. W., 1985, Nature, 317, 44
  * Crocker et al. (2021) Crocker R. M., Krumholz M. R., Thompson T. A., 2021, Monthly Notices of the Royal Astronomical Society, 502, 1312
  * Davis et al. (2014) Davis S. W., Jiang Y.-F., Stone J. M., Murray N., 2014, ApJ, 796, 107
  * Everett et al. (2008) Everett J. E., Zweibel E. G., Benjamin R. A., McCammon D., Rocks L., Gallagher III J. S., 2008, The Astrophysical Journal, 674, 258
  * Ferriere (2001) Ferriere K. M., 2001, Reviews of Modern Physics, 73, 1031
  * Fielding et al. (2018) Fielding D., Quataert E., Martizzi D., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3325
  * Girichidis et al. (2018) Girichidis P., Naab T., Hanasz M., Walch S., 2018, Monthly Notices of the Royal Astronomical Society, 479, 3042
  * Gronke & Oh (2020) Gronke M., Oh S. P., 2020, Monthly Notices of the Royal Astronomical Society, 492, 1970
  * Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, The Astrophysical Journal Letters, 777, L38
  * Heintz & Zweibel (2018) Heintz E., Zweibel E. G., 2018, ApJ, 860, 97
  * Heintz et al. (2020) Heintz E., Bustard C., Zweibel E. G., 2020, ApJ, 891, 157
  * Holguin et al. (2019) Holguin F., Ruszkowski M., Lazarian A., Farber R., Yang H. K., 2019, Monthly Notices of the Royal Astronomical Society, 490, 1271
  * Hopkins et al. (2020) Hopkins P. F., et al., 2020, Monthly Notices of the Royal Astronomical Society, 492, 3465
  * Huang et al. (2020) Huang X., Davis S. W., Zhang D., 2020, The Astrophysical Journal, 893, 50
  * Ipavich (1975) Ipavich F. M., 1975, The Astrophysical Journal, 196, 107
  * Jiang & Oh (2018) Jiang Y.-F., Oh S. P., 2018, ApJ, 854, 5
  * Jubelgas et al. (2008) Jubelgas M., Springel V., Enßlin T., Pfrommer C., 2008, Astronomy & Astrophysics, 481, 33
  * Kempski & Quataert (2020) Kempski P., Quataert E., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1801
  * Kennicutt Jr (1998) Kennicutt Jr R. C., 1998, The Astrophysical Journal, 498, 541
  * Kennicutt Jr & De Los Reyes (2021) Kennicutt Jr R. C., De Los Reyes M. A., 2021, The Astrophysical Journal, 908, 61
  * Kennicutt Jr & Evans (2012) Kennicutt Jr R. C., Evans N. J., 2012, Annual Review of Astronomy and Astrophysics, 50, 531
  * Kim et al. (2021) Kim J.-G., Ostriker E. C., Filippova N., 2021, The Astrophysical Journal, 911, 128
  * Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, The Astrophysical Journal, 654, 304
  * Krumholz & Thompson (2012) Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155
  * Mao & Ostriker (2018) Mao S. A., Ostriker E. C., 2018, The Astrophysical Journal, 854, 89
  * Martin (2005) Martin C. L., 2005, The Astrophysical Journal, 621, 227
  * Parker (1966) Parker E. N., 1966, ApJ, 145, 811
  * Quataert et al. (2021) Quataert E., Thompson T. A., Jiang Y.-F., 2021, arXiv preprint arXiv:2102.05696
  * Ramzan et al. (2020) Ramzan B., Ko C., Chernyshov D., 2020, The Astrophysical Journal, 905, 117
  * Ruszkowski et al. (2017) Ruszkowski M., Yang H.-Y. K., Zweibel E., 2017, The Astrophysical Journal, 834, 208
  * Ryu et al. (2003) Ryu D., Kim J., Hong S. S., Jones T. W., 2003, ApJ, 589, 338
  * Salem & Bryan (2014) Salem M., Bryan G. L., 2014, Monthly Notices of the Royal Astronomical Society, 437, 3312
  * Scannapieco & Brüggen (2015) Scannapieco E., Brüggen M., 2015, The Astrophysical Journal, 805, 158
  * Semenov et al. (2021) Semenov V. A., Kravtsov A. V., Caprioli D., 2021, The Astrophysical Journal, 910, 126
  * Skilling (1971) Skilling J., 1971, ApJ, 170, 265
  * Skinner & Ostriker (2013) Skinner M. A., Ostriker E. C., 2013, ApJS, 206, 21
  * Socrates et al. (2008) Socrates A., Davis S. W., Ramirez-Ruiz E., 2008, The Astrophysical Journal, 687, 202
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, Annual Review of Astronomy and Astrophysics, 53, 51
  * Spilker et al. (2020) Spilker J. S., et al., 2020, The Astrophysical Journal, 905, 85
  * Strickland & Heckman (2009) Strickland D. K., Heckman T. M., 2009, The Astrophysical Journal, 697, 2030
  * Uhlig et al. (2012) Uhlig M., Pfrommer C., Sharma M., Nath B. B., Enßlin T., Springel V., 2012, Monthly Notices of the Royal Astronomical Society, 423, 2374
  * Veilleux et al. (2005) Veilleux S., Cecil G., Bland-Hawthorn J., 2005, Annu. Rev. Astron. Astrophys., 43, 769
  * Veilleux et al. (2020) Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, The Astronomy and Astrophysics Review, 28, 1
  * Wiener et al. (2013) Wiener J., Oh S. P., Guo F., 2013, Monthly Notices of the Royal Astronomical Society, 434, 2209
  * Wiener et al. (2017) Wiener J., Pfrommer C., Oh S. P., 2017, Monthly Notices of the Royal Astronomical Society, 467, 906
  * Wiener et al. (2019) Wiener J., Zweibel E. G., Ruszkowski M., 2019, MNRAS, 489, 205
  * Wiersma et al. (2009) Wiersma R. P., Schaye J., Smith B. D., 2009, Monthly Notices of the Royal Astronomical Society, 393, 99
  * Zhang (2018) Zhang D., 2018, Galaxies, 6, 114
  * Zhang et al. (2018) Zhang D., Davis S. W., Jiang Y.-F., Stone J. M., 2018, The Astrophysical Journal, 854, 110
  * Zweibel (2020) Zweibel E. G., 2020, The Astrophysical Journal, 890, 67
  * de los Reyes & Kennicutt Jr (2019) de los Reyes M. A., Kennicutt Jr R. C., 2019, The Astrophysical Journal, 872, 16