# The Physics of Galactic Winds Driven by Cosmic Rays I: Diffusion

Eliot Quataert1, Todd A. Thompson2,3, and Yan-Fei Jiang(姜燕飞)4
1Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
2 Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus, OH 43210, USA
3 Center for Cosmology and Astro-Particle Physics (CCAPP), The Ohio State
University, 191 West Woodruff Ave., Columbus, OH 43210, USA
4 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY, 10010, USA

###### Abstract

The physics of Cosmic ray (CR) transport remains a key uncertainty in
assessing whether CRs can produce galaxy-scale outflows consistent with
observations. In this paper, we elucidate the physics of CR-driven galactic
winds for CR transport dominated by diffusion. A companion paper considers CR
streaming. We use analytic estimates validated by time-dependent spherically-
symmetric simulations to derive expressions for the mass-loss rate, momentum
flux, and speed of CR-driven galactic winds, suitable for cosmological-scale
or semi-analytic models of galaxy formation. For CR diffusion coefficients
$\kappa\gtrsim r_{0}c_{i}$ where $r_{0}$ is the base radius of the wind and
$c_{i}$ is the isothermal gas sound speed, the asymptotic wind energy flux is
comparable to that supplied to CRs, and the outflow rapidly accelerates to
supersonic speeds. By contrast, for $\kappa\lesssim r_{0}c_{i}$, CR-driven
winds accelerate more slowly and lose most of their energy to gravity, a CR
analogue of photon-tired stellar winds. Given CR diffusion coefficients
estimated using Fermi gamma-ray observations of pion decay, we predict mass-
loss rates in CR-driven galactic winds of order the star formation rate for
dwarf and disc galaxies. The dwarf galaxy mass-loss rates are small compared
to the mass-loadings needed to reconcile the stellar and dark matter halo mass
functions. For nuclear starbursts (e.g., M82, Arp 220), CR diffusion and pion
losses suppress the CR pressure in the galaxy and the strength of CR-driven
winds. We discuss the implications of our results for interpreting
observations of galactic winds and for the role of CRs in galaxy formation.

###### keywords:

Galaxies; Winds; Cosmic Rays

##  1 Introduction

Galactic winds play a key role in galaxy evolution, shaping the galaxy stellar
mass function, the mass-metallicity relation, and affecting the morphology of
galaxies over cosmic time (e.g., Somerville & Davé 2015). Despite their
importance, many puzzles persist. These include the acceleration mechanism for
the cool atomic and molecular gas seen in emission and absorption from rapidly
star-forming galaxies across the universe (e.g., Veilleux et al. 2020), the
very large mass fluxes and low star formation efficiencies inferred for dwarf
galaxies, and the apparent need for massive outflows from normal star-forming
galaxies in the local universe in order to match some models of Galactic
chemical evolution (e.g., Andrews et al. 2017).

Among the physical mechanisms suggested for driving large amounts of cool gas
from star-forming galaxies, cosmic rays (CRs) are of particular interest
because their pressure is dynamically important with respect to gravity in the
Milky Way disc (Boulares & Cox, 1990). In this picture, CRs are injected into
the disc of the galaxy by supernovae and stellar processes, and scatter off of
magnetic fluctuations in the ISM, slowly diffusing (or streaming) away from
the disc. This process sets up a pressure gradient that can in principle
accelerate gas away from the host galaxy.

Many treatments of CR driven winds exist in the literature, starting with the
foundational work of Ipavich (1975) who computed time-steady solutions for CR
driven winds assuming streaming transport at the Alfvén velocity. More
detailed models including CR streaming and hydromagnetic wave pressure were
computed by Breitschwerdt et al. (1991), deriving a mass-loss rate of order M⊙
yr-1 from the Galaxy. Everett et al. (2008) further developed a hybrid thermal
and CR driven wind model for the Galaxy. More generally, on the basis of
momentum conservation, Socrates et al. (2008) argued that CRs could drive
significant mass fluxes from star-forming galaxies.

In parallel to these more analytic and steady-state treatments, there is a
large and growing body of work on multi-dimensional simulations of galaxy
formation with CRs, which suggest that they can play a number of important
roles: e.g., driving cold-gas outflows from galaxies (e.g., Booth et al.
2013), modifying the thermal and ionization state of the CGM (e.g., Ji et al.
2020), and heating gas in galaxy groups and clusters, suppressing cooling
flows (e.g., Ruszkowski et al. 2017). These numerical models vary in their
treatment of CR transport, including isotropic diffusion (e.g, Uhlig et al.
2012; Simpson et al. 2016), anisotropic diffusion along magnetic fields (e.g.,
Pakmor et al. 2016; Chan et al. 2019; Hopkins et al. 2020), and/or CR
streaming along magnetic fields (e.g., Ruszkowski et al. 2017; Chan et al.
2019; Hopkins et al. 2020).

As suggested by the diversity of approaches in the literature, one of the
significant uncertainties in modeling the properties of CRs in galaxy
formation is that the physical mechanism(s) coupling CRs to the thermal plasma
are not fully understood (e.g., Amato & Blasi 2018). Small-scale fluctuations
in the magnetic field scatter CRs, setting their mean free path. On scales
larger than this mean free path, the CRs can be approximated as a fluid
(Skilling, 1971). However, it is not clear whether the fluctuations that
scatter CRs are self-excited by the CRs themselves (e.g., the streaming
instability; Kulsrud & Pearce 1969) or produced by background turbulent
fluctuations in the interstellar/circumgalactic-medium (ISM/CGM) cascading to
small-scales where they couple to the CRs (e.g., Yan & Lazarian 2002). Nor is
it clear whether the dominant scattering mechanism depends on the
thermodynamic phase of the ISM/CGM, the local magnetic field strength, or
other physical properties. The distinction between self-excited versus
background turbulence as the source of scattering that sets the CR mean free
path is at the core of the streaming/diffusion dichotomy that dominates CR
transport modeling.

Wiener et al. (2017) explored the difference between galactic winds driven
with CR diffusion and CR streaming in simulations, finding that diffusive CR
transport results in much larger overall mass-loss rates than in CR streaming
models. Chan et al. (2019) and Hopkins et al. (2020) found similar results and
further argued that models with CR transport dominated by diffusion were
required for CRs to efficiently escape Milky-way, M31, and Magellanic Cloud-
like galaxies, as is required by gamma-ray observations of pion decay produced
by CRs interacting with the ISM (Abdo et al., 2010b, a; Ackermann et al.,
2012). These results highlight that a better understanding of CR microphysics
is critical for understanding the role of CRs in galaxy formation.

Together with the numerical and analytic treatments of CR-driven winds in the
literature, there is also a body of phenomenological work interpreting the
non-thermal radio and gamma-ray emission from star-forming galaxies. This
directly informs our understanding of the underlying CR population and their
role in driving galactic winds. The far-infrared-radio correlation and the
gamma-ray emission from star-forming galaxies, including those with strong
winds like M82, NGC 253, and Arp 220 can be used to constrain the average CR
injection rate per unit star formation, the gas density seen by CRs as they
propagate, the magnetic field strength, and the CR escape timescale (e.g.,
Pavlidou & Fields 2001; Torres 2004; Lacki et al. 2010; Lacki & Thompson 2013;
Yoast-Hull et al. 2013, 2014; Buckman et al. 2020; Crocker et al. 2021, 2020).
These works thus provide benchmarks for understanding how the inner or “base”
boundary conditions for a putative CR-driven wind vary as a function of galaxy
properties.

This paper is the first in a series that aims to elucidate the physics of
galactic winds driven by CRs using a combination of analytic calculations and
idealized numerical simulations. In this paper we focus on the case of CR
diffusion. A companion paper discusses the case of CR streaming, which turns
out to be much physically richer than the diffusion limit considered here. In
§2 we present analytic estimates of the mass-loss rate and terminal velocity
of galactic winds driven by CRs in the diffusion approximation. In §3 we
validate these estimates with time-dependent numerical simulations, which
further show that the winds reach a laminar steady state. §4 calibrates CR
diffusion coefficients and CR pressures in galaxies using Fermi gamma-ray
observations of pion decay, and calculates the resulting implications for CR-
driven galactic winds. We summarize and discuss the implications of our
results in §5. Appendix A presents a linear stability analysis and shows that
although sound waves are formally unstable in CR-driven winds with diffusion
(Drury & Falle, 1986), the growth rates are too slow for the instability to be
dynamically important in almost all cases (the exception is extremely low gas
sound speeds). We also derive the linear WKB dispersion relation for the two-
moment cosmic-ray transport model simulated in §3 and show that entropy and
sound waves are linearly stable in the WKB limit for CR transport by diffusion
in the two-moment CR transport model.

##  2 Analytic Approximations for Cosmic Ray Driven Winds with Diffusion

In this section we focus on analytic approximations to the spherical steady
state galactic wind problem in the presence of cosmic rays. Following these
analytical approximations, in Section 3 we treat the more complete numerical
problem and solve the time-dependent spherically symmetric CR-driven wind
problem using the numerical scheme of Jiang & Oh (2018). We show that the
analytic approximations in this section do an excellent job of reproducing the
more complete numerical results.

The general equation of motion, including gas pressure $p$, CR pressure
$p_{c}$, a general gravitational potential $\Phi$, and magnetic forces in the
limit of magnetohydrodynamics can be written as

$\frac{\partial{\bf v}}{\partial t}+{\bf v}\cdot\nabla{\bf
v}=-\frac{1}{\rho}\nabla p-\frac{1}{\rho}\nabla
p_{c}-\nabla\Phi+\frac{1}{4\pi\rho}(\nabla\times{\bf B})\times{\bf B}.$ (1)

The energy equation for the CRs in the absence of CR sources and pionic
losses, and including diffusion and streaming at the Alfvén velocity down the
CR pressure gradient, i.e., ${\bf v_{s}}=-{\bf v_{A}}|\nabla p_{c}|/\nabla
p_{c}$ with ${\bf v_{A}}={\bf B}/(4\pi\rho)^{1/2}$, can be written as

$\frac{\partial E_{c}}{\partial t}+\nabla\cdot{\bf F}_{c}=\left({\bf v}+{\bf
v_{\rm s}}\right)\cdot\nabla p_{c},$ (2)

where the “equilibrium” CR flux111The two-moment CR model we solve numerically
in §3 evolves ${\bf F_{c}}$ as an independent variable and the flux reduces to
equation 3 only when time variations are sufficiently slow (see eq. 43). is

${\bf F}_{c}=4p_{c}\left({\bf v}+{\bf v_{\rm s}}\right)-\kappa\,{\bf
n}\left({\bf n}\cdot\nabla E_{c}\right),$ (3)

$E_{c}$ is the CR energy density, $p_{c}=E_{c}/3$, $\kappa$ is the diffusion
coefficient, and ${\bf n}={\bf v_{\rm A}}/|\bf v_{\rm A}|$. Equation 2, and
indeed the scalar CR pressure in equation 1, is formally valid only on scales
larger than the mean-free-path of the $\sim$ GeV energy CRs that dominate the
total energy of the CR population (e.g., Skilling 1971).222Even given this
restrictive assumption, the diffusion term in equation 3 in general depends on
local plasma conditions. It can also depend on the cosmic ray energy density
itself, in which case the ‘diffusion’ term in eq. 3 is not even truly
diffusive (e.g., Skilling 1971; Wiener et al. 2013). We do not consider this
complication in the present work.

Throughout this paper we focus exclusively on the case of CR diffusion, which
corresponds to ${\bf v_{s}}=0$ in equation 3. In a companion paper, we
consider the case of CR streaming. We further assume a hydrodynamic model
which drops the field-aligned diffusive flux in equation 3 in favor of an
effective isotropic diffusion equation with a spatially constant diffusion
coefficient $\kappa$. We also consider a simple model for the gravity of a
galaxy and its host dark matter potential: an isothermal sphere for which
$\Phi=2V_{g}^{2}\ln r$ where $\sqrt{2}V_{g}$ is the circular velocity of the
potential. Finally, we simplify the thermodynamics of the gas by using an
isothermal equation of state with sound speed $c_{i}$, as would be appropriate
for a warm gas in ionization equilibrium, or which might approximate the
effects of turbulence in the atmosphere of the host galaxy. In addition, we
consider a single phase flow and do not consider possible variation of CR
transport with, e.g., the ionization state of the gas.

With these approximations, the steady state equations of motion are

$\dot{M}_{w}=4\pi r^{2}\rho v={\rm const},$ (4)
$v\frac{dv}{dr}=-\frac{1}{\rho}\frac{dp}{dr}-\frac{1}{\rho}\frac{dp_{c}}{dr}-\frac{2V_{g}^{2}}{r},$
(5)

and

$v\frac{dp_{c}}{dr}=-\frac{4p_{c}}{3r^{2}}\frac{dr^{2}v}{dr}+\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\kappa\frac{dp_{c}}{dr}\right).$
(6)

A key simplification can be made for the purposes of analytic estimates by
noting that the order of magnitude of the diffusion term in equation 6
relative to both of the other terms (which describe CR advection and adiabatic
energy changes) is

$\frac{\rm Diffusion}{\rm
Advection}\sim\frac{\kappa}{rv}\sim\left(\frac{\kappa}{10^{28}\,{\rm
cm^{2}\,s^{-1}}}\right)\left(\frac{30\,{\rm km\,s^{-1}\,kpc}}{v\,r}\right),$
(7)

which suggests that in the limit of rapid diffusion, and near the base of the
outflow where the gas velocity is small, the advective and adiabatic terms in
the CR energy equation can simply be neglected. For reference, the diffusion
coefficient in the Milky Way is estimated to be $\sim 10^{29}$ cm2 s-1 (e.g.,
Trotta et al. 2011), although there is significant uncertainty in this
estimate because of a degeneracy between the diffusion coefficient and the
size of the CR propagation region (the ‘halo;’ Linden et al. 2010; Trotta et
al. 2011). In §4, we return to constraints on the diffusion coefficient in
galaxies using observations of gamma-ray emission from pion decay. For now, we
proceed under the assumption that for sufficiently large diffusion
coefficients $\kappa$, one can neglect the advective and adiabatic terms in
equation 6 and focus solely on the diffusion term, for which the steady state
solution is

$p_{c}=p_{c,0}+\frac{\dot{E}_{c}}{12\pi\kappa}\left(\frac{1}{r}-\frac{1}{r_{0}}\right)$
(8)

where $p_{c,0}$ is the base CR pressure at radius $r_{0}$ and

$\dot{E}_{c}=-12\pi r^{2}\kappa\frac{dp_{c}}{dr}=-12\pi
r_{0}^{2}\kappa\frac{dp_{c}}{dr}\bigg{|}_{r_{0}}$ (9)

is the energy per unit time supplied to the CRs by supernovae and other
processes in the galaxy. The assumption that diffusion is much faster than
advection near the base of the wind implies that $\dot{E}_{c}$ is independent
of radius. In our numerical wind solutions in §3 we shall see that this
assumption is valid at small radii near the sonic point (as assumed here to
estimate $\dot{M}$), but that advection then takes over as the dominant energy
transport mechanism at larger radii (see Fig. 4 discussed below).

We will make frequent use of the CR scale height near the base of the wind:

$H_{c}=\left(\frac{-d\ln p_{c}}{dr}\bigg{|}_{r_{0}}\right)^{-1}\equiv
r_{0}\,h_{c}$ (10)

where we have defined $h_{c}$ as the dimensionless CR pressure scale-height at
$r_{0}$. With this definition, $\dot{E}_{c}=12\pi r_{0}h_{c}^{-1}\kappa
p_{c,0}$ and equation 8 becomes $p_{c}=p_{c,0}[1+h_{c}^{-1}(r_{0}/r-1)]$. Note
that solutions with $h_{c}<1$ have $p_{c}\rightarrow 0$ at a finite radius
$r=r_{0}/(1-h_{c})$. If CR diffusion were the only energy transport mechanism,
$h_{c}\simeq 1$ would be the physical solution extending to large radii.
However, we will see below that $h_{c}\lesssim 1$ is typically the appropriate
approximate solution at small radii (where our analysis here applies), because
advection of CR energy eventually takes over from diffusion as the dominant
transport mechanism.

Equations 4 and 5 can be combined to yield a wind equation

$\frac{1}{v}\frac{dv}{dr}\left(v^{2}-c_{i}^{2}\right)=\frac{2c_{i}^{2}}{r}-\frac{1}{\rho}\frac{dp_{c}}{dr}-\frac{2V_{g}^{2}}{r}.$
(11)

With equation 8, i.e., assuming $\dot{E}_{c}$ is a constant because diffusion
is rapid, the wind equation becomes

$\frac{1}{v}\frac{dv}{dr}\left(v^{2}-c_{i}^{2}\right)=\frac{2c_{i}^{2}}{r}+\frac{\dot{E}_{c}}{12\pi
r^{2}\rho\kappa}-\frac{2V_{g}^{2}}{r}.$ (12)

In the portion of the flow where $v<c_{i}$, an approximation to the density
profile can be derived by assuming hydrostatic equilibrium so that

$\frac{c_{i}^{2}}{\rho}\frac{d\rho}{dr}-\frac{\dot{E}_{c}}{12\pi
r^{2}\rho\kappa}=-\frac{2V_{g}^{2}}{r}.$ (13)

Equation 13 has the solution

$\rho(r)=\rho_{0}\left[x^{-\xi}+A(x^{-1}-x^{-\xi})\right]$ (14)

where $x=r/r_{0}$ and $\xi$ and $A$ are the two key dimensionless parameters
in this problem, which measure the strength of gravity relative to the gas and
CR pressures at the outflow’s base:

$\xi\equiv\frac{2V_{g}^{2}}{c_{i}^{2}}$ (15)

and

$\begin{split}A&\equiv\frac{\dot{E}_{c}}{12\pi\kappa
c_{i}^{2}\rho_{0}r_{0}(\xi-1)}\simeq\frac{\dot{E}_{c}}{24\pi\kappa
V_{g}^{2}\rho_{0}r_{0}}\\\
&\simeq\frac{p_{c,0}}{2h_{c}\rho_{0}V_{g}^{2}}=\frac{1}{2h_{c}}\frac{c_{c,0}^{2}}{V_{g}^{2}}\lesssim
1,\end{split}$ (16)

where we assume $\xi\gg 1$ in the approximations after the first equality, and
where the last equality defines the base CR sound speed
$c_{c,0}=(p_{c,0}/\rho_{0})^{1/2}$. The conclusion that $A\lesssim 1$ follows
from the fact that $A$ is roughly the ratio of the CR pressure force to the
gravitational force at the base of the outflow. If that ratio is $\gtrsim 1$,
then the initial assumption of HE is invalid and the gas distribution would
expand out until $A\lesssim 1$.

The two power-laws in equation 14, namely $\rho\sim r^{-\xi}$ and $\rho\sim
r^{-1}$, correspond to the gas pressure and CR pressure dominated phases of
the solution, respectively. The transition between the two happens at a radius
$r_{tr}\simeq r_{0}\,A^{-1/(\xi-1)}$. The CR-dominated hydrostatic
$\rho\propto p_{c,0}/r$ solution in eq. 14 was derived by Ji et al. 2020 for
the CGM in their cosmological zoom-in simulations with CRs. Hopkins et al.
(2021a) further used this solution to estimate the properties outflows on CGM
scales driven by CRs. The key difference between the solutions here and their
estimates is that we self-consistently match our CR-dominated solution onto
the gas-pressure dominated solution near the galaxy in order to calculate the
properties of galaxy-scale winds. We do not include a CGM at larger radii in
our calculations.

The mass-loss rate in the wind is set by the conditions at the critical
(sonic) point of equation 12, which we denote as $r_{s}$. Setting the
numerator and denominator of the wind equation to zero simultaneously yields
two conditions:

$v(r_{s})=c_{i}\ \ \ {\rm and}\ \ \
\frac{r_{s}\,\rho(r_{s})}{r_{0}\rho_{0}}=\frac{c_{c,0}^{2}}{2h_{c}{\rm
V^{2}_{g,eff}}}.$ (17)

where

${\rm V^{2}_{g,eff}}=V_{g}^{2}-c_{i}^{2},$ (18)

These conditions allow for a compact expression for the wind mass-loss rate:

$\dot{M}_{w}=4\pi r_{s}^{2}\rho(r_{s})c_{i}=\frac{2\pi
r_{s}r_{0}\rho_{0}c_{i}}{h_{c}}\frac{c_{c,0}^{2}}{{\rm V^{2}_{g,eff}}}.$ (19)

To estimate the radius of the sonic point $r_{s}$, one can numerically solving
the 2nd critical point equation in eq. 17 using the analytic density profile
in equation 13, i.e.,

$\frac{c_{c,0}^{2}}{2h_{c}(V_{g}^{2}-c_{i}^{2})}\simeq
x_{s}^{-\xi+1}+\frac{c_{c,0}^{2}}{h_{c}(2V_{g}^{2}-c_{i}^{2})}\left(1-x_{s}^{-\xi+1}\right)$
(20)

A good approximate solution to equation 20 can be found by expanding for
$c_{i}^{2}\ll 2V_{g}^{2}$ which yields

$r_{s}\simeq
r_{0}\,\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})}$
(21)

As an example, if $c_{i}\simeq c_{c,0}\simeq 10\,{\rm km\,\,s}^{-1}$ (as in
the Milky Way), equation 21 predicts $r_{s}/r_{0}\simeq 1.38,$ 1.13, and 1.05
for $h_{c}=1$ and $V_{g}=30$, 60, and 100 $\,{\rm km\,\,s}^{-1}$, while
numerical solution of equation 20 yields nearly identical results of
$r_{s}/r_{0}\simeq 1.39$, 1.13, and 1.05, respectively. The sonic point is
quite close to the base of the wind, i.e., $r_{s}\simeq r_{0}$, unless
$c_{i}\sim V_{g}$.

Using equation (21) for the sonic point in equation (19), we obtain an
expression for the mass-loss rate:

$\dot{M_{w}}\simeq\frac{2\pi
r_{0}^{2}\rho_{0}c_{i}}{h_{c}}\,\frac{c_{c,0}^{2}}{{\rm
V^{2}_{g,eff}}}\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})}.$
(22)

In the limit that $2V_{g}\gg c_{i}$,

$\dot{M}_{w}\sim\frac{2\pi r_{0}^{2}c_{i}p_{c,0}}{h_{c}V_{g}^{2}}=\frac{2\pi
r_{0}^{2}\rho_{0}c_{i}}{h_{c}}\left(\frac{c_{c,0}}{V_{g}}\right)^{2}.$ (23)

We stress that for equation 23 to be applicable, CR diffusion must dominate
over advection out to at least the sonic point since the mass-loss rate is set
by the flow properties at and interior to the sonic point. This, together with
equation 7, implies that the relevant criterion for the validity of our
analytics is roughly $\kappa\gtrsim r_{0}c_{i}$. We shall see that this is
borne out by the numerical simulations in §3. Diffusion coefficients
satisfying this constraint are also strongly suggested by gamma-ray data on
pion decay in nearby star-forming galaxies, as we show in §4.

Note that equations 22 and 23 imply that the mass-loss rate does not
explicitly depend on the diffusion coefficient $\kappa$ at fixed base CR
pressure $p_{c,0}$ (though there is an implicit dependence via $h_{c}$ as we
will see below). However, for a fixed CR injection power $\dot{E}_{c}$,
equations 22 and 23 imply $\dot{M}_{w}\propto 1/\kappa$ because the base CR
pressure itself scales as
$p_{c,0}=\rho_{0}c_{c,\,0}^{2}\propto\dot{E}_{c}/\kappa$ (see eq. 8).
Substituting into equation (22), we find that

$\dot{M}_{w}\simeq\,\left(\frac{r_{0}\,c_{i}}{6\kappa}\right)\left(\frac{\dot{E}_{c}}{{\rm
V^{2}_{g,eff}}}\right)\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})},$
(24)

which makes the $\kappa$ dependence, and the competition between diffusion and
advection at the base of the outflow, explicit.

We now present an order of magnitude estimate of the scale-height $h_{c}$ by
determining the radius $r_{\rm adv}\simeq r_{0}(1+h_{c})$ at which the
advective flux $F_{\rm adv}=4p_{c}v=\dot{M}_{w}p_{c}/(\pi r^{2}\rho)$ is
comparable to the diffusive flux $F_{\rm diff}=-3\kappa dp_{c}/dr$. This
occurs when the velocity of the outflow is given by

$v(r_{\rm adv})\approx-\kappa\frac{d\ln
p_{c}}{dr}\approx\frac{c_{i}}{h_{c}}\frac{\kappa}{c_{i}r_{0}}$ (25)

where we have dropped a factor of $3/4$ consistent with the rough nature of
the estimates that follow. Since $h_{c}\lesssim 1$ and our analytics assumes
$\kappa\gtrsim c_{i}r_{0}$, equation 25 implies that the transition from
diffusion to advection happens exterior to the sonic point. Gas pressure is
then negligible and the momentum equation becomes $\rho vdv/dr\simeq-
dp_{c}/dr-2\rho V_{g}^{2}/r$. At the sonic point, the two terms on the right-
hand-side of $dv/dr$ are comparable (see eq. 17). Since the gas density scale-
height is smaller than the CR scale-height, somewhat exterior to the sonic
point $\rho vdv/dr\simeq-dp_{c}/dr$. Multiplying by $4\pi r^{2}$ gives

$\dot{M}\frac{dv}{dr}\bigg{|}_{\sim r_{\rm adv}}\sim 4\pi r_{\rm
adv}^{2}\frac{p_{c,0}}{h_{c}}\rightarrow v(r_{\rm
adv})\sim\frac{h_{c}V_{g}^{2}}{c_{i}}$ (26)

where in the second expression we have used the approximate version of
$\dot{M}$ from equation 23 and have taken $r_{\rm adv}\approx r_{0}$,
consistent with $h_{c}\lesssim 1$. Equating equations 25 and 26 then yields
$h_{c}\sim(c_{i}/V_{g})(\kappa/r_{0}c_{i})^{1/2}$. For
$\kappa\rightarrow\infty$ this gives $h_{c}>>1$, inconsistent with the local
approximation used here; the solution should be $p_{c}\propto 1/r$, i.e.,
$h_{c}\simeq 1$, so that

$h_{c}\sim\min\bigg{(}1,\frac{c_{i}}{V_{g}}\sqrt{\frac{\kappa}{r_{0}c_{i}}}\bigg{)}.$
(27)

As an example, if $\kappa\sim 10^{29}$ cm2 s${}^{-1}\simeq
10\,(r_{0}c_{i}/30\,{\rm km\,\,s}^{-1}\,{\rm kpc})$, $c_{i}\sim 10\,\,{\rm
km\,\,s}^{-1}$, and $V_{g}\sim 100\,\,{\rm km\,\,s}^{-1}$, $h_{c}\sim 1/3$. In
our analytic scalings that follow in this section we primarily normalize
$h_{c}$ to a value of $1/4$ motivated by these fiducial parameters, but in
several plots we will use equation 27. In §3, we also compare equation 27 to
our numerical solutions and find good agreement.

Our final expression for the mass-loss rate in equation 22 can be written as

$\begin{split}\dot{M}_{w}\simeq&\frac{4\pi
r_{0}^{2}\,\rho_{0}\,c_{i}}{h_{c}}\,\left(\frac{c_{c,0}^{2}}{2{\rm
V^{2}_{g,eff}}}\right)\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})}\\\
\simeq&\,\,0.06\,\,{\rm M_{\odot}\,yr^{-1}}\left(\frac{r_{0}}{1\,\,{\rm
kpc}}\right)^{2}\left(\frac{1/4}{h_{c}}\right)\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})}\\\
&\times\,\left(\frac{n_{0}}{1\,{\rm cm^{-3}}}\right)\left(\frac{c_{i}\
c^{2}_{c,0}}{[10\,\,{\rm km\,\,s}^{-1}]^{3}}\right)\left(\frac{{\rm
V_{g,eff}}}{100\,\,{\rm km\,\,s}^{-1}}\right)^{-2}\end{split}$ (28)

where $n_{0}=\rho_{0}/m_{p}$. For reference, if we scale for parameters
appropriate to a galaxy like the Milky Way (with $V_{g}\simeq 150$ km s-1,
$c_{i}=c_{c,0}\simeq 10$ km s-1, $n_{0}\simeq 1$ cm-3, and $r_{0}\sim 5$ kpc),
equation 28 yields $\dot{M}_{w}\simeq 1\,{\rm M_{\odot}\,yr^{-1}}$, comparable
to the star formation rate. Equation 28 also predicts $\dot{M}_{w}\propto
p_{c,0}c_{i}/h_{c}$, which scales $\propto p_{c,0}c_{i}$ or $\propto
p_{c,0}c_{i}^{1/2}$, depending on which regime of equation 27 is appropriate
(these scalings assume $r_{s}\simeq r_{0}$ for simplicity). To the (uncertain)
extent that the CR pressure is comparable in different phases of the ISM, the
outflow is thus likely to be dominated by the warmer ISM phases, though only
by a factor of a few.

Figure 1 shows the mass-loss rate (eq. 28 with $h_{c}$ from eq. 27 taking
$\kappa=10r_{0}c_{i}$) as a function of the two key dimensionless parameters
in the problem, namely the strength of gravity ($V_{g}/c_{i}=(\xi/2)^{1/2}$;
see eq. 15) and the base CR sound speed
($p_{c,0}/\rho_{0}c_{i}^{2}=c_{c,0}^{2}/c_{i}^{2}\propto A/\xi$); see eq. 16).
The mass-loss rate is given in units of

$\dot{M}_{0}=4\pi r_{0}^{2}\rho_{0}c_{i}\simeq 3.2\,\,{\rm
M_{\odot}\,yr^{-1}}\,\left(\frac{r_{0}}{\,{\rm
kpc}}\right)^{2}\left(\frac{n_{0}}{1\,{\rm
cm^{-3}}}\right)\left(\frac{c_{i}}{10\,\,{\rm km\,\,s}^{-1}}\right).$ (29)

To express the mass-loss rate in CR driven winds in a more intuitive form we
take $V_{g}\gg c_{i}$ and use the simplified expression for the mass-loss rate
in equation 23. We then write the CR energy injection rate at the base of the
outflow as

$\dot{E}_{c}=\epsilon_{c}\dot{M}_{*}c^{2},$ (30)

where $\dot{M}_{*}$ is the star formation rate and $\epsilon_{c}\equiv
10^{-6.3}\epsilon_{c,-6.3}$ is related to the fraction of SNe energy that goes
into CRs: for $10^{51}$ ergs per SNe and 1 SNe per 100 $M_{\odot}$ of stars
formed, $\epsilon_{c}=10^{-6.3}$ if $10\%$ of the SNe energy goes into primary
CRs. Equation 23 can then be approximately rewritten as (taking $r_{s}\simeq
r_{0}$ to simplify eq. 21)

$\displaystyle\frac{\dot{M}_{w}}{\dot{M}_{*}}$ $\displaystyle\simeq$
$\displaystyle 0.8\
\epsilon_{c,-6.3}\left(\frac{c_{i}r_{0}}{\kappa}\right)\left(\frac{100\,\,{\rm
km\,\,s}^{-1}}{V_{\rm g,\,eff}}\right)^{2}$ (31) $\displaystyle\simeq$
$\displaystyle 0.08\ \epsilon_{c,-6.3}\,\left(\frac{c_{i}r_{0}}{{30\,\rm
km\,\,s^{-1}\,kpc}}\right)$ $\displaystyle\times\left(\frac{10^{29}\,{\rm
cm^{2}\,s^{-1}}}{\kappa}\right)\left(\frac{100\,\,{\rm km\,\,s}^{-1}}{V_{\rm
g,\,eff}}\right)^{2}.$

Per equation 7 and the associated discussion, $\kappa\sim\,10c_{i}r_{0}$ is
plausible; we have scaled to representative values in the second line. This
expression shows that significant wind mass loading relative to the global
star formation rate is in principle possible, and that it should grow strongly
with decreasing $V_{\rm g,\,eff}$. However, if $\kappa\sim 10^{29}$ cm2 s-1 is
appropriate, the mass-loading factors in dwarf galaxies are relatively modest
compared to the values $\dot{M}_{w}\gg\dot{M}_{*}$ needed to reconcile the
stellar and dark matter mass functions. We return in Section 4 to an
observational calibration of the diffusion coefficients in other star-forming
galaxies.

A second instructive expression for the wind mass-loss rate can be obtained by
comparing the mass-loss rate estimated here to the star formation rate
predicted by feedback-regulated models of star formation in galaxies, namely
$\dot{M}_{*}\approx\pi r_{0}^{2}\dot{\Sigma}_{*}$ where the surface density of
star formation is (e.g., Thompson et al. 2005)

$\dot{\Sigma}_{*}\simeq\frac{\sqrt{8}\pi G\Sigma_{g}^{2}\phi}{v_{*}}$ (32)

where $\phi=1+\Sigma_{*}/\Sigma_{g}$ describes the contribution of the stellar
disc with surface density $\Sigma_{*}$ to the local gravitational potential
and $v_{*}=p_{*}/m_{*}\approx 3000\,{\rm km\,\,s}^{-1}$ is the momentum per
unit mass of star formed associated with stellar feedback, which supports the
disc against its own self-gravity (Ostriker & Shetty, 2011). Equation 23 for
the mass-loss rate can then be recast as

$\frac{\dot{M}_{\rm w}}{\dot{M}_{*}}\simeq
8\left(\frac{1/4}{h_{c}}\right)\left(\frac{c_{i}v_{*}}{3\times 10^{4}\,{\rm
km^{2}\,s^{-2}}}\right)\left(\frac{100\,\,{\rm
km\,\,s}^{-1}}{V_{g}}\right)^{2}\frac{p_{c,0}}{\pi G\Sigma_{g}^{2}\phi}$ (33)

Equation 33 again shows that for Milky-way like conditions in which the CR
pressure is comparable to that needed for hydrostatic equilibrium in the
galactic disc ($p_{c,0}\simeq 1/3\,\pi G\Sigma_{g}^{2}\phi$; Boulares & Cox
1990), the wind mass-loss rate driven by diffusing CRs can be of order or
larger than the star formation rate. As we discuss in Section 4, equation (33)
also shows that for dense starburst galaxies, in which $p_{c,0}\ll\pi
G\Sigma_{g}^{2}\phi$ (Lacki et al., 2010; Lacki et al., 2011), the mass-loss
rate in CR driven winds is significantly reduced.

Figure 1: Analytic mass-loss rates for CR driven galactic winds in the limit
of rapid diffusion (eq. 28 with $h_{c}$ from eq. 27), as a function of the
strength of gravity relative to the gas sound speed in the disc
($V_{g}/c_{i}$) and the base CR pressure (${\rm p_{c,0}/\rho_{0}c_{i}^{2}}$).
The mass-loss rates are normalized by equation (29). Labeled values of
$V_{g}/c_{i}$ on the color bar are logarithmically distributed and correspond
to the curves on the plot.

###  2.1 Wind Energetics, Momentum Flux, and Velocity

For a steady state solution, the momentum equation (eq. 5) and CR energy
equation with diffusion (eq. 6) can be combined to yield a total conserved
energy outflow rate, namely

$\dot{E}_{w}=\dot{M}_{w}\left(\frac{1}{2}v^{2}+c_{i}^{2}\ln\rho+4c_{c}^{2}+\Phi\right)+\dot{E}_{\rm
c}={\rm constant}$ (34)

where $\dot{E}_{\rm c}=-12\pi r^{2}\kappa\,dp_{c}/dr$ (as before) and where
the four terms in parentheses in equation 34 correspond to the gas kinetic
energy flux ($\equiv\dot{E}_{k}$), the gas enthalpy/advective flux (assuming
our isothermal equation of state for the gas), the CR enthalpy/advective flux
($\equiv\dot{E}_{c,\,\rm adv}$), and the gravitational energy flux
($\equiv\dot{E}_{\rm grav}$), respectively. Note that equation 34 no longer
assumes that diffusion dominates over advection at all radii.

Table 1: Parameters and properties of our numerical simulations. Numerical
resolution is $\delta r/r=5.25\times 10^{-4}$ unless otherwise noted. Columns
are diffusion coefficient $\kappa$, velocity $V_{g}$ of the isothermal
gravitational potential, base CR pressure $p_{c,0}$, outer radius of domain
$r_{out}$, reduced speed of light $V_{m}$, simulation mass0loss rate
$\dot{M}_{sim}$ in units of $\dot{M}_{0}$ (eq. 29), simulation mass-loss rate
in units of $\dot{M}_{max}$ (the maximum mass-loss rate allowed by energy
conservation; eq. 47), the dimensionless base CR scale-height $h_{c}$ (eq.
10), and the kinetic energy and CR enthalpy fluxes at the top of the domain
relative to the diffusive CR flux at the base
$(\dot{E}_{k}(r_{out})+\dot{E}_{c,adv}(r_{out}))/\dot{E}_{c}(r_{0})$. When the
latter is $\simeq 1$, most of the CR energy flux at the base of the wind goes
into kinetic energy at large radii, while when it is $\lesssim 1$, most of the
CR energy goes into work against gravity, and $\dot{M}_{w}\simeq\dot{M}_{max}$
(eq. 40 & 47).

$\kappa$ | $V_{g}$ | $p_{c,0}$ | $r_{out}$ | $V_{m}$ | $\dot{M}_{sim}$ | $\dot{M}_{sim}$ | $v(r_{out})$ | $h_{c}$ | $\dot{E}_{k}(r_{out})+\dot{E}_{c,adv}(r_{out})$
---|---|---|---|---|---|---|---|---|---
($r_{0}c_{i}$) | ($c_{i}$) | ($\rho_{0}c_{i}^{2}$) | ($r_{0}$) | ($c_{i}$) | ($\dot{M}_{0}$) | ($\dot{M}_{max}$) | ($V_{g}$) | – | $(\dot{E}_{c}[r_{0}])$
$20^{a}$ | 200 | 400 | 7.6 | 60000 | 0.21 | 0.032 | 13.8 | 0.023 | 0.95
$10^{b}$ | 100 | 100 | 7.6 | 30000 | 0.16 | 0.067 | 9.5 | 0.03 | 0.94
$3.3^{c}$ | 33 | 11 | 7.6 | 10000 | 0.09 | 0.19 | 5 | 0.053 | 0.8
$30$ | 10 | 1 | 5 | 3000 | 0.016 | 0.018 | 13.6 | 0.32 | 0.93
$10$ | 10 | 1 | 5 | 3000 | 0.022 | 0.054 | 8.2 | 0.23 | 0.93
$3.3$ | 10 | 1 | 5 | 3000 | 0.03 | 0.15 | 4.6 | 0.16 | 0.83
$2$ | 10 | 1 | 5 | 3000 | 0.03 | 0.42 | 3.3 | 0.13 | 0.57
$1$ | 10 | 1 | 15 | 3000 | 0.03 | 0.75 | 0.5 | 0.14 | 0.23
$0.33$ | 10 | 1 | 15 | 3000 | 0.007 | 0.78 | 0.4 | 0.21 | 0.16
$0.11$ | 10 | 1 | 25 | 3000 | 0.0019 | 0.85 | 0.43 | 0.23 | 0.08
$10$ | 10 | 0.3 | 5 | 3000 | 0.007 | 0.058 | 8.2 | 0.23 | 0.99
$10$ | 10 | 0.1 | 5 | 3000 | 0.0027 | 0.056 | 8.1 | 0.24 | 0.95
$10$ | 10 | 3 | 15 | 3000 | 0.065 | 0.089 | 9.4 | 0.23 | 0.9
$10$ | 6 | 1 | 5 | 3000 | 0.043 | 0.055 | 7.5 | 0.33 | 0.89
$3.3$ | 6 | 1 | 15 | 3000 | 0.055 | 0.26 | 4.9 | 0.24 | 0.67
$10$ | 3 | 1 | 15 | 3000 | 0.136 | 0.11 | 7.6 | 0.5 | 0.9
$3.3$ | 3 | 1 | 15 | 3000 | 0.156 | 0.3 | 4.2 | 0.4 | 0.74

a $\delta r/r=8.25\times 10^{-6}$ for $1\leq r\leq 1.14$; $c_{i}$ is 1/20 of
$\kappa=1$, $V_{g}=10$ sim
b $\delta r/r=2.5\times 10^{-5}$ for $1\leq r\leq 1.11$; $c_{i}$ is 1/10 of
$\kappa=1$, $V_{g}=10$ sim
c $\delta r/r=6.6\times 10^{-5}$ for $1\leq r\leq 1.11$; $c_{i}$ is 1/3 of
$\kappa=1$, $V_{g}=10$ sim

The total energy flux and terminal velocity of the wind can be estimated as
follows. Under our assumption of $\kappa\gtrsim r_{0}c_{i}$, at small radii
near the base of the wind, the energy flux at small radii is almost entirely
due to CR diffusion, so that $\dot{E}_{w}\simeq\dot{E}_{\rm c}$. This neglects
gravity near the base of the wind, which we return to below. The CR diffusive
flux near the base can be related to the mass-loss flux using equation 24
where in this expression we now identify $\dot{E}_{c}$ near the base as the
total wind luminosity $\dot{E}_{w}$ in equation 34. At large radii, the energy
flux will be dominated by the gas with $\dot{E}_{w}\simeq
0.5\dot{M}_{w}v_{\infty}^{2}$. Equating our expressions for the total wind
luminosity at small and large radii yields the terminal velocity of the wind

$\begin{split}&v_{\infty}\simeq 2{\rm
V_{g,eff}}\left(\frac{3\kappa}{r_{0}c_{i}}\right)^{1/2}\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{-c_{i}^{2}/(4V_{g}^{2})}\\\
&\simeq 10^{3}\,\,{\rm km\,\,s}^{-1}\left(\frac{V_{g}}{100\,{\rm
km\,s^{-1}}}\right)\left(\frac{\kappa}{10^{29}{\rm cm^{2}\,s^{-1}}}\frac{\rm
30\,kpc\,{\rm km\,s^{-1}}}{r_{0}c_{i}}\right)^{1/2}\end{split}$ (35)

where the second expression assumes $V_{g}\gg c_{i}$. Recall that we require
$\kappa\gtrsim r_{0}c_{i}$ given our assumptions used in deriving the
approximate wind solutions here. In this case, equation 35 shows that
$v_{\infty}\gtrsim{\rm V_{g,eff}}$, with $v_{\infty}\sim 1-10V_{g}$ plausible.
We show below (Fig. 7) that equation 35 agrees very well with our numerical
simulations.

Using equations (28) and (35) we can write down an expression for the
asymptotic total momentum loss rate carried by the wind,

$\dot{p}_{w}=\dot{M}_{w}v_{\infty}\simeq\frac{4\pi
r_{0}^{2}c_{i}p_{c,0}}{h_{c}V_{\rm
g}}\left(\frac{3\kappa}{r_{0}c_{i}}\right)^{1/2}.$ (36)

where we assume $V_{g}\gg c_{i}$ to simplify the expressions and highlight the
key scalings. Assuming no pionic losses in the host galaxy, we can write this
expression in terms of the star formation rate using equations (24) & (30) as

$\dot{p}_{w}=\dot{M}_{w}v_{\infty}\simeq\frac{\epsilon_{c}\dot{M}_{*}c^{2}}{V_{\rm
g}}\left(\frac{r_{0}c_{i}}{3\kappa}\right)^{1/2}.$ (37)

This quantity can be compared with the total momentum rate carried by photons
from star formation $\dot{p}_{*}=L_{*}/c=\epsilon_{*}\dot{M}_{*}c$, where
$\epsilon_{*,\,-3.3}=\epsilon_{*}/5\times 10^{-4}$ for steady-state star
formation and a standard IMF. We then have that

$\displaystyle\frac{\dot{p}_{w}}{\dot{p}_{*}}$ $\displaystyle\simeq$
$\displaystyle\frac{\epsilon_{c}}{\epsilon_{*}}\frac{c}{V_{\rm
g,\,eff}}\left(\frac{r_{0}c_{i}}{3\kappa}\right)^{1/2}$ (38)
$\displaystyle\simeq$ $\displaystyle
0.5\,\frac{\epsilon_{c,-6.3}}{\epsilon_{*,\,-3.3}}\left(\frac{100\,\rm
km\,\,s^{-1}}{V_{\rm g}}\right)$ $\displaystyle\times$
$\displaystyle\left(\frac{r_{0}\,c_{i}}{30\,\rm km\,\,s^{-1}\,\,kpc}\frac{\rm
10^{29}\,cm^{2}\,\,s^{-1}}{\kappa}\right)^{1/2}.$

Note that the asymptotic kinetic energy loss rate
$\dot{E}_{w}=0.5\dot{M}_{w}v_{\infty}^{2}=\dot{E}_{c}=\epsilon_{c}\dot{M}_{*}c^{2}$.
CR-driven winds in the rapid diffusion approximation are thus energy-
conserving in the sense that the wind kinetic energy power at large radii is
of order that supplied to the CRs at small radii. As in the discussions
following equations (31) and (33), we reiterate that this conclusion assumes
that there are no pionic losses in the host galaxy (see §4).

Figure 2: Left: CR pressure and gas density as a function of radius for
$V_{g}=10$ and several diffusion coefficients (in units of $r_{0}c_{i}$; see
Table LABEL:tab:units). The dotted profile is the analytic hydrostatic
solution in equation 14 for $h_{c}=1/4$. The analytic density profile is
reasonable interior to the sonic point (and thus for estimating $\dot{M}$) but
not at larger radii for the higher $\kappa$ solutions. Right: Flow velocity
and CR sound speed for the same solutions, in units of gas sound speed
$c_{i}$. The vertical green dotted line in both panels is the analytic
estimate of the location of the critical point (eq. 21). For the $\kappa=3-30$
solutions, the flow accelerates very rapidly and the velocity exceeds the
escape velocity just exterior to the base of the wind. For $\kappa=1$ the flow
is significantly slower, and the density profile is closer to hydrostatic. For
$\kappa=0.33$ the velocities are subsonic and well below the escape speed at
all radii shown here. We return to the $\kappa=0.33$ solution in Figure 6 and
show that it does drive a wind but one whose properties are very different
from the high $\kappa$ simulations.

###  2.2 Maximum Mass-Loss Rate

There is a strict upper limit to the mass-loss rate associated with CR-driven
winds that is set by energy conservation. This is the regime in which gravity
significantly modifies the energetics of the outflow and cannot be neglected
as we did in deriving equation 35. If the asymptotic speed of the wind
vanishes then $\dot{E}_{c}(r_{0})+\dot{M}_{\rm max}\Phi(r_{0})\simeq 0$, i.e.,
$\dot{M}_{\rm max}\simeq 2\dot{E}_{c}(r_{0})/v^{2}_{\rm esc}(r_{0})$ where
$v_{\rm esc}(r_{0})$ is the escape speed from the base of the wind.333$v_{\rm
esc}(r_{0})$ is formally not defined for the $\ln(r)$ potential focused on
here, but it is $\sim 2V_{g}$ for a more realistic potential which deviates
from isothermal at larger radii (and/or for a computational domain of
reasonable size even with a $\ln(r)$ potential). This regime is analogous to
photon-tired winds in stellar wind theory (Owocki & Gayley, 1997). Using
$\dot{E}_{c}=\epsilon_{c}\dot{M}_{*}c^{2}$ (eq. 30), we can write

$\frac{\dot{M}_{\rm
max}}{\dot{M}_{*}}\simeq\,9\,\epsilon_{c,-6.3}\,\left(\frac{100\,\,{\rm
km\,\,s}^{-1}}{v_{\rm esc}(r_{0})}\right)^{2},$ (39)

which has an “energy-driven” velocity scaling (Murray et al., 2005).
Alternatively, using $\dot{E}_{c}(r_{0})\simeq 12\pi r_{0}\kappa
p_{c,0}/h_{c}$ we find

$\displaystyle\dot{M}_{\rm max}$ $\displaystyle\simeq$
$\displaystyle\frac{24\pi r_{0}\kappa p_{c,\,0}}{v_{\rm
esc}^{2}(r_{0})h_{c}}\simeq\frac{24\pi
r_{0}\kappa\rho_{0}}{h_{c}}\left(\frac{c_{c,\,0}}{v_{\rm
esc}(r_{0})}\right)^{2}$ (40) $\displaystyle\simeq$ $\displaystyle 23\,{\rm
M_{\odot}\,\,yr^{-1}}\left(\frac{r_{0}}{\rm
kpc}\right)\left(\frac{\kappa}{10^{29}\rm
cm^{2}\,\,s^{-1}}\right)\left(\frac{n_{0}}{\rm cm^{-3}}\right)$
$\displaystyle\times$
$\displaystyle\left(\frac{1/4}{h_{c}}\right)\left(\frac{c_{c,\,0}}{10\,{\rm
km\,\,s^{-1}}}\right)^{2}\left(\frac{100\,{\rm km\,\,s^{-1}}}{v_{\rm
esc}(r_{0})}\right)^{2}$

The maximum cosmic ray-driven wind momentum will be realized for mass-loss
rates a bit below $\dot{M}_{\rm max}$ when $v_{\infty}\sim v_{\rm esc}(r_{0})$
(at $\dot{M}_{\rm max}$, $v_{\infty}\equiv 0$ and so $\dot{p}_{w}=0$). A rough
estimate is

$\frac{\dot{p}_{w,\,{\rm max}}}{\dot{p}_{*}}\sim
5\frac{\epsilon_{c,-6.3}}{\epsilon_{*,\,-3.3}}\left(\frac{100\,\rm
km\,\,s^{-1}}{v_{\rm esc}(r_{0})}\right).$ (41)

but more detailed calculations are required to precisely determine the
numerical pre-factor in this expression.

It is straightforward to show that the mass-loss rate estimated in equation 28
is

$\dot{M}_{w}\simeq\dot{M}_{\rm max}\,\left(\frac{r_{0}c_{i}}{3\kappa}\right),$
(42)

where we have assumed $v^{2}_{\rm esc}(r_{0})\simeq 4V_{g}^{2}$. Equation 42
shows that the wind mass-loss rate is in general $<\dot{M}_{\rm max}$ given
our restriction to $\kappa\gtrsim r_{0}c_{i}$. Our analysis so far does not
preclude that winds exist for $\kappa$ smaller than $r_{0}c_{i}$, but our
analytic approximations break down in this regime. Equation 42 shows, however,
that if one extrapolates our estimate of $\dot{M}_{w}$ to $\kappa<r_{0}c_{i}$,
the solutions reach the maximum mass-loss rate allowed by energy conservation
for $\kappa\sim 1/3r_{0}c_{i}$. Thus a plausible conjecture is that any wind
with $\kappa\lesssim r_{0}c_{i}$ will be a CR analogue of photon-tired winds
(Owocki & Gayley, 1997). In the next section we show that this conjecture is
correct: numerical simulations with $\kappa\lesssim r_{0}c_{i}$ do produce a
wind, but one whose character is very different from the high $\kappa$
simulations. See Section 3.5. In particular, winds with $\kappa\lesssim
r_{0}c_{i}$ have $\dot{M}_{w}\simeq\dot{M}_{\rm max}$ (and thus lose most of
their energy to work against gravity) and accelerate much more slowly, with
the sonic point at radii many times larger than $r_{0}$. Finally, we note that
for $\kappa\equiv 0$, there cannot be a wind so long as the base CR and gas
sound speeds are below the escape speed (as assumed here) since both fluids
are then adiabatic.

Equations 31, 39, & 42 highlight the importance of the value of the diffusion
coefficient for assessing the implications of our results for CR-driven
galactic winds. The maximum mass-loss rate a given CR energy flux can sustain
is significantly larger than the galaxies’ star formation rate (eq. 39).
However, this is only realizable for relatively small diffusion coefficients,
$\kappa\lesssim r_{0}c_{i}$. If, on the other hand, $\kappa\gg r_{0}c_{i}$,
$\dot{M}_{w}/\dot{M}_{*}\lesssim 1$ (eq. 31) and thus less dynamically
important. In §4 we estimate the diffusion coefficient in star-forming
galaxies using gamma-ray constraints on pion losses.

##  3 Numerical Simulations

###  3.1 Equations

We solve the time-dependent cosmic ray hydrodynamic equations based on the two
moment approach as developed by Jiang & Oh (2018) in one dimensional spherical
polar coordinates.444The two-moment approach enables a much more accurate
treatment of CR streaming, which is not essential for the present paper but
will be in later publications in this series. This algorithm has been
implemented in the magneto-hydrodynamic code Athena++ (Stone et al., 2020). As
in §2, we use an isothermal equation of state for the gas with isothermal
sound speed $c_{i}$ and take the gravitational potential to be
$\Phi=2V_{g}^{2}\ln r$.

Neglecting CR streaming, the full set of equations for gas density $\rho$,
flow velocity $v$, comic ray energy density $E_{c}$ and flux $F_{c}$ in 1D
spherical polar coordinates are

$\displaystyle\frac{\partial\rho}{\partial
t}+\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\rho v\right)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\frac{\partial\left(\rho
v\right)}{\partial t}+\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\rho v^{2}\right)$ $\displaystyle=$
$\displaystyle-\rho\frac{\partial\Phi}{\partial
r}-c_{i}^{2}\frac{\partial\rho}{\partial
r}+\sigma_{c}\left[F_{c}-v(E_{c}+p_{c})\right],$ $\displaystyle\frac{\partial
E_{c}}{\partial
t}+\frac{1}{r^{2}}\frac{\partial\left(r^{2}F_{c}\right)}{\partial r}$
$\displaystyle=$ $\displaystyle-v\,\sigma_{c}[F_{c}-v(E_{c}+p_{c})],$
$\displaystyle\frac{1}{V_{m}^{2}}\frac{\partial F_{c}}{\partial
t}+\frac{\partial p_{c}}{\partial r}$ $\displaystyle=$
$\displaystyle-\sigma_{c}\left[F_{c}-v(E_{c}+p_{c})\right].$ (43)

Here the cosmic ray pressure is $p_{c}=E_{c}/3$. The reduced speed of light is
$V_{m}$, which is chosen to be much larger than $v$ in the whole simulation
box. For CR transport by spatially independent diffusion, $\sigma_{c}$ is a
constant and is related to the diffusion coefficient used elsewhere in this
paper by $\kappa\equiv 1/(3\sigma_{c})$.555Jiang & Oh (2018) follow the
conventions of the radiation transfer literature in which the diffusive flux
is $F_{c}=-\sigma_{c}^{-1}\nabla p_{c}$. In the rest of this paper, we follow
the conventions of the CR literature and define the diffusive flux as
$F_{c}=-\kappa\nabla E_{c}$ (eq. 2). This accounts for the factor of 3
relating $\kappa$ and $\sigma_{c}^{-1}$.

Substituting the right-hand-side of the fourth of equations 43 into the 3rd
term on the right-hand-side of the 2nd of equations 43 produces the usual
$\partial p_{c}/\partial r$ cosmic-ray pressure gradient in the momentum
equation, along with another term related to the time variation of the CR
flux. Likewise, the term on the right-hand-side of the 3rd of equations 43
becomes the usual $v\partial p_{c}/\partial r$ term related to CR pdV work,
again with another term related to the time variation of the CR flux. The
steady state versions of equations 43 are thus equivalent to the steady state
equations 4-6 solved in §2.

Table 2: Summary of Units for Numerical Simulations. $c_{i}$ is the gas
isothermal sound speed (assumed to be a constant), $\rho_{0}$ is the gas
density at the base of the wind (at radius $r=r_{0}$), and $p_{c,0}$ is the
base CR pressure.

Quantity Symbol Units Radial Velocity $v$ $c_{i}$ ‘Isothermal’ CR Sound Speed
$c_{c}\equiv\left(\frac{p_{c}}{\rho}\right)^{1/2}$ $c_{i}$ Gravitational
Velocity $V_{g}$ $c_{i}$ Density $\rho$ $\rho_{0}$ CR Pressure $p_{c}$
$p_{c,0}$ CR Flux $F_{c}$ $c_{i}p_{c,0}$ CR Diffusion Coefficient $\kappa$
$c_{i}r_{0}$

###  3.2 Initial and Boundary Conditions

For each simulation, we pick gas density $\rho_{0}$ and cosmic ray pressure
$p_{c,0}$ at the bottom boundary $r_{0}$ and then initialize the gas density
and cosmic ray energy density at each radius $r$ as

$\displaystyle\rho(r)$ $\displaystyle=$
$\displaystyle\rho_{0}\left(r/r_{0}\right)^{-V_{g}^{2}/c_{i}^{2}},$
$\displaystyle E_{c}(r)$ $\displaystyle=$ $\displaystyle
3p_{c,0}\left(r/r_{0}\right)^{-V_{g}^{2}/c_{i}^{2}}.$ (44)

We apply floor values of $10^{-4}$ for $\rho/\rho_{0}$ and $p_{c}(r)/p_{c,0}$
in the initial condition. The flow velocity and cosmic ray flux are set to be
zero in the whole simulation domain initially.

For the bottom boundary condition at $r_{0}$, the cosmic ray energy density
$E_{c}$ and gas density $\rho$ are fixed to be the initial values. The flow
velocity in the ghost zones is set such that mass flux $\rho vr^{2}$ is
continuous from the last active zone to the ghost zones. We keep the gradient
of the cosmic ray flux continuous across the bottom boundary. For the boundary
condition at the top of the simulation domain, we keep the gradient of
$\rho,E_{c}$ and $F_{c}$ continuous across the boundary. The flow velocity at
the top boundary is set by requiring the mass flux $\rho vr^{2}$ to be
continuous.

###  3.3 Overview of Simulations

Table LABEL:tab:compare summarizes our suite of simulations. The key physical
parameters of each simulation are $V_{g}/c_{i}$ and
$p_{c,0}/\rho_{0}c_{i}^{2}$, as in the analytics of §2, as well as the
diffusion coefficient $\kappa$ (in units of $r_{0}c_{i}$). The key numerical
parameters are the resolution, reduced speed of light, and box size.

The units for the results of our numerical simulations are summarized in Table
LABEL:tab:units: gas density is in units of the base density $\rho_{0}$,
speeds are in units of the gas isothermal sound speed $c_{i}$, CR pressure is
in units of the base CR pressure $p_{c,0}$, and CR fluxes are in units of
$p_{c,0}c_{i}$.

Figure 3: Numerator and denominator (eqs 45) of the steady state wind equation
evaluated for our $\kappa=10$, $V_{g}=10$ simulation. $N_{2}$ is the
approximation to the numerator of the wind equation used in our analytics (eq.
& 46). The vertical dotted line is the analytic estimate of the location of
the critical point (eq. 21). The numerical simulation reaches a steady state
that passes through a critical point close to that predicted by our
(approximate) analytics. Figure 4: Energetics of the wind evaluated for our
$\kappa=10$, $V_{g}=10$ simulation. The total energy flux ($\dot{E}_{w}$; see
eq. 34) and mass flux ($\dot{M}_{w}$) are independent of radius. At the base
of the wind, nearly all of the energy is carried by the diffusing CRs
($\dot{E}_{c}$). Advection of CR enthalpy by the motion of the gas
($\dot{E}_{c,adv}$) is important at intermediate radii while at large radii
most of the energy is in the kinetic energy of the outflowing gas
($\dot{E}_{k}$). For the logarithmic potential used in our simulations, we
define $\dot{E}_{grav}$ using $\Phi=2V_{g}^{2}\ln(r/r_{0})$ so that
$\dot{E}_{grav}>0$.

Nearly all of our numerical solutions reach a steady state state that does not
evolve significantly in time once we run for a few box crossing times. Because
the solutions reach a steady state, they are independent of the reduced speed
of light $V_{m}$, which only enters into the time-dependent terms in equations
43.

Figure 5: Dependence of density and pressure profiles (Left) and velocity and
CR sound speed profiles (Right) on the depth of the gravitational potential
$V_{g}$ (in units of $c_{i}$, the isothermal gas sound speed). All solutions
are for $\kappa=10$ and equal gas and CR pressures at the base of the wind
($p_{c,0}=\rho_{0}c_{i}^{2}$). The CR pressure profiles are relatively
independent of the gravitational potential, being largely set by rapid CR
diffusion. By contrast, the density falls off more gradually, and the flow
accelerates more gradually, for weaker gravity.

Drury & Falle (1986) showed that rapid CR diffusion drives sound waves
unstable in the presence of a background CR pressure gradient, so it is not a
priori obvious that our wind solutions should reach a steady state. We show in
Appendix A that the growth rate of this instability is too slow to grow
significantly in galactic winds, unless the gas sound speed is extremely small
(the $V_{g}=200$ simulation, which is our lowest $c_{i}$ simulation, is the
only one to show significant time dependence; see Fig. 10). We also derive the
linear WKB dispersion relation for sound waves and entropy modes in the two-
moment CR system of equations, and show that the linear modes are stable,
consistent with the simulations.

###  3.4 Numerical Solutions For $\kappa\gtrsim r_{0}c_{i}$

Figure 2 shows the resulting density, CR pressure, fluid velocity, and CR
sound speed profiles for $V_{g}/c_{i}=10$, varying $\kappa$ from $0.33-30$ (in
units of $r_{0}c_{i}$; Table LABEL:tab:units). The vertical dotted line in Fig
2 (as well as Figs. 3 & 4 below) is the critical point predicted in equation
21, namely $r_{s}\simeq 1.047r_{0}$.

The left panel of Figure 2 shows that for all of the solutions, the CR
pressure falls off much more slowly than the density near the base. This is
because the gas density scale-height is initially quite small, $\sim
r_{0}(c_{i}/V_{g})^{2}\sim 10^{-2}r_{0}$ while the CR pressure scale-height is
significantly larger due to CR diffusion $\sim h_{c}r_{0}\sim 0.25r_{0}$ (eq.
27). For all of the solutions, the hydrostatic equilibrium approximation for
the density profile (eq. 14) is reasonable at very small radii. There is,
however, a bifurcation in the density profiles as a function of $\kappa$. The
solutions with $\kappa$ = 0.33 and 1 are closer to hydrostatic over the entire
radial range while for $\kappa\gtrsim 1$ the density falls off much more
rapidly at large radii. The right panel of Figure 2 shows that this
bifurcation in density profiles corresponds to a bifurcation in velocity
profiles. The solutions with $\kappa\gtrsim 3$ drive strong winds that
accelerate to a velocity larger than the escape speed by $r\sim 1.2r_{0}$. By
contrast, for $\kappa=1$ the velocities are much lower and for $\kappa=0.33$
the solution is subsonic for all of the radii shown in Figure 2. We return to
the lowest $\kappa$ simulation in §3.5 and show that it does produce a wind,
but one whose character is very different from the higher $\kappa$ simulations
focused on here.

Figure 3 quantifies the extent to which the numerical solution satisfies the
time steady critical point conditions, for the $\kappa=10$ numerical
simulation. Specifically, from equation 11 we define the denominator and
numerator of the critical point equation to be

$D\equiv v^{2}-c_{i}^{2}\ \ \ {\rm and}\ \ \
N\equiv\frac{2c_{i}^{2}}{r}-\frac{1}{\rho}\frac{dp_{c}}{dr}-\frac{2V_{g}^{2}}{r}$
(45)

and we further define

$N_{2}\equiv\frac{2c_{i}^{2}}{r}+\frac{p_{c,0}r_{0}}{\rho
r^{2}h_{c}}-\frac{2V_{g}^{2}}{r}$ (46)

to be the numerator of the critical point equation under the approximation
that $\dot{E}_{c}$ is constant as a function of radius, which was used in our
analytic derivations in §2. Figure 3 shows that $N$ and $D$ pass through zero
at the same point, as expected for a time steady solution that passes through
the critical point. The location of the critical point is close to the
prediction in equation 21, which is shown by the vertical dotted line in Fig.
3; the deviation is because the analytics makes the approximation that the
hydrostatic density profile (eq. 14) is valid all the way to the sonic point.
Equation 46 (used in our analytics) is also a reasonable, though not perfect,
approximation to the true numerator in the critical point equation.

Figure 4 quantifies the radial profiles of the mass-loss rate $\dot{M}$ and
various contributions to the wind energy flux in the numerical solution for
$\kappa=10$. The mass-loss rate $\dot{M}$ is independent of radius as expected
for a steady state solution. The contributions to the energy flux shown in
Figure 4 are defined in and below equation 34 (we do not show the gas
enthalpy/advective flux, which is negligible). We note that for the
logarithmic potential used here, we define $\Phi=2V_{g}^{2}\ln(r/r_{0})$. The
zero point of the potential is thus at the base of the wind ($r=r_{0}$),
rather than at infinity, as is more common; the latter is not possible given
the divergence of the logarithmic potential as $r\rightarrow\infty$. As a
result of this choice of zero point, $\dot{E}_{grav}>0$.

Figure 4 shows that interior to the critical point, the energy flux is
dominated by CR diffusion, as assumed in the analytic calculation in §2. At
intermediate radii the CR advective flux dominates the energy transport while
eventually at large radii, when the CR pressure has accelerated the gas to
well above the CR sound speed and escape speeds, the gas kinetic energy flux
dominates. Note that the gravitational contribution to the energy flux is
small at all radii, consistent with $\dot{M}_{w}\ll\dot{M}_{\rm max}$ and
$v>V_{g}$. The qualitative features of Figure 4 are the same for
$\kappa=3-30$. For $\kappa=1$, however, most of the energy at the outer edge
of the domain is in CR enthalpy rather than gas kinetic energy.666This implies
that $v(r_{out})=0.5V_{g}=5c_{i}$ in Table LABEL:tab:compare is a lower limit.
If all of the CR enthalpy is converted to kinetic energy, the final velocity
would be $\simeq 1.8V_{g}$, close to the analytic estimate in equation 35. For
$\kappa=0.33$, the energetics of the flow is yet more different from that in
Figure 4, as we discuss in §3.5.

Figure 5 shows how the $\kappa=10$ solutions depend on the depth of the
gravitational potential $V_{g}$ (in units of the isothermal sound speed
$c_{i}$). The solutions are qualitatively very similar. The primary difference
is that weaker gravity (lower $V_{g}$) implies a larger density scale-height
at the base of the wind. Correspondingly, the wind accelerates significantly
more slowly.

###  3.5 Low $\kappa\lesssim r_{0}c_{i}$ Solutions

The analytic derivations in §2 primarily assumed $\kappa\gtrsim r_{0}c_{i}$.
Our numerical simulations with $\kappa\gtrsim r_{0}c_{i}$ are consistent with
the properties of these analytic solutions, as we discuss further in the next
section. Here we present numerical wind models with $\kappa\lesssim
r_{0}c_{i}$, which differ dramatically from their high $\kappa$ counterparts.
We stress that the higher $\kappa$ solutions are likely the most
astrophysically relevant, given estimates of diffusion coefficients from MW CR
data (§2) and gamma-ray data from pion decay (§4).

Figure 6: Kinematics (top) and energetics (bottom) of the wind for our
$\kappa=0.33$, $V_{g}=10$ simulation, which differs dramatically from the
higher $\kappa$ simulations. The flow accelerates very gradually, reaching the
sonic point ($v/c_{i}=1$) only at $r\sim 8r_{0}$, compared to $r\simeq
1.02r_{0}$ for high $\kappa$ solutions. The flow velocity is also
significantly lower ($\lesssim V_{g}$) than for the high $\kappa$ solutions
(Fig. 2). At the base of the wind, nearly all of the energy is carried by the
diffusing CRs ($\dot{E}_{c}$). However, almost all of this energy is used to
overcome the gravitational potential, leading to
$\dot{E}_{grav}\simeq\dot{E}_{w}$ at large radii (note that $\dot{E}_{grav}>0$
for our logarithmic potential with $\Phi=2V_{g}^{2}\ln(r/r_{0})$). This
corresponds to $\dot{M}_{w}\simeq\dot{M}_{\rm max}$ (eq. 40) and is a CR
analogue of photon-tired winds in stellar wind theory (Owocki & Gayley, 1997).
For larger $\kappa$ the energy lost to gravity is negligible and nearly all of
the CR energy at the base goes into CR enthalpy and kinetic energy at larger
radii (see Fig. 4).

Figure 6 shows the kinematics (top) and energetics (bottom) of the wind for
our $\kappa=0.33$, $V_{g}=10$ simulation, which differs significantly from the
higher $\kappa$ simulations. The velocity profile in Figure 6 is shown over a
larger range of radii than was plotted in Figure 2. The flow accelerates very
gradually reaching the sonic point only at $r_{s}\sim 8r_{0}$777Note that
there are small oscillations in the solution near the sonic point; these
remain at late times even once the solution settles into a steady state. We do
not have a definitive explanation for these oscillations but they are not
present in any of our $\kappa\gtrsim r_{0}c_{i}$ solutions (the majority in
Table LABEL:tab:compare)., compared to $r_{s}\simeq r_{0}$ for the high
$\kappa$ solutions. The most striking property of the energetics is that
$\dot{E}_{grav}\simeq\dot{E}_{w}$ (see eq. 34) at large radii while at the
base of the wind, nearly all of the energy is carried by CR diffusion
($\dot{E}_{c}$); recall that for our logarithmic potential
$\Phi=2V_{g}^{2}\ln(r/r_{0})$ and so $\dot{E}_{grav}>0$. The fact that
$\dot{E}_{grav}(r_{out})\simeq\dot{E}_{c}(r_{0})$ implies that most of the
energy supplied to CRs at the base goes into lifting (almost hydrostatically)
the gas out to large radii. As a result, the gas kinetic energy is a minor
contribution to the energy flux at large radii, in marked contrast to the high
$\kappa$ simulation shown in Figure 4. The condition
$\dot{E}_{c}\simeq\dot{E}_{grav}$ is precisely the condition for the CR-
analogue of photon-tired winds; this is also the regime in which
$\dot{M}_{w}\sim\dot{M}_{\rm max}$ (eq. 40). For our simulations, we can
define the escape velocity from the base of the wind to be the velocity needed
to reach the outer edge of the box, i.e, $v_{\rm
esc}(r_{0})=2V_{g}\sqrt{\ln(r_{out}/r_{0})}$. Equation 40 then becomes

$\frac{\dot{M}_{max}}{\dot{M}_{0}}\simeq\frac{3}{2h_{c}\ln(r_{\rm
out}/r_{0})}\frac{\kappa}{r_{0}c_{i}}\frac{c_{c,0}^{2}}{V_{g}^{2}}$ (47)

Table 1 gives the simulation mass-loss rates in units of $\dot{M}_{max}$. For
$\kappa=0.33r_{0}c_{i}$, $\dot{M}_{w}\simeq 0.78\dot{M}_{max}$, close to the
maximum value allowed by energy conservation. This is consistent with the
energetics in Figure 6. For comparison, we note that $\dot{M}_{w}/\dot{M}_{\rm
max}=0.85,0.78,0.75,0.42,0.15,0.054,0.018$ for our simulations with
$\kappa/r_{0}c_{i}=0.11,0.33,1,2,3.3,10,\&30$, respectively. This is
consistent with the expectation from equation 42 that
$\dot{M}_{w}/\dot{M}_{max}\propto 1/\kappa$ for $\kappa\gtrsim r_{0}c_{i}$.
The bifurcation in density and velocity profiles in Figure 2 at $\kappa\simeq
r_{0}c_{i}$ thus corresponds to solutions with $\dot{M}_{w}\sim\dot{M}_{max}$
($\kappa\lesssim r_{0}c_{i}$, nearly hydrostatic, slower acceleration) vs.
those with $\dot{M}_{w}\ll\dot{M}_{max}$ ($\kappa\gtrsim r_{0}c_{i}$, rapid
acceleration). More generally, the last column in Table LABEL:tab:compare
shows that for solutions with $\kappa>r_{0}c_{i}$ nearly all of the CR energy
supplied at the base goes into kinetic energy or enthalpy of the wind at large
radii. By contrast, for solutions with $\kappa<r_{0}c_{i}$, the asymptotic
kinetic and enthalpy flux is small compared to the CR energy flux at the base
because most of the energy is lost escaping the gravitational potential.

The importance of gravity for the low $\kappa$ solution in Figure 6 implies
that some of the details of the solution – though not the fact that
$\dot{M}_{w}\simeq\dot{M}_{max}$ – will likely be sensitive to the form of the
gravitational potential. In particular, we suspect that the exact acceleration
profile for the gas, including the final terminal speed (which depends on the
small residual energy not lost to work against gravity), will depend on the
details of the potential.

###  3.6 Synthesis of Analytics vs. Numerics

Figure 7: Ratio of the analytically predicted mass-loss rate (eq. 28 with
$h_{c}$ from eq. 27), terminal velocity (eq. 35), and base CR scale-height
(eq. 27) compared to the simulation results. The agreement is good for the
full range of simulated parameters. Multiplying $h_{c}$ in eq. 27 by $\simeq
0.75$ would remove the small systematic offset in $h_{c}$ and $\dot{M}_{w}$.

Figure 7 compares our numerical and analytic solutions for the mass-loss rate,
terminal velocity, and base CR scale-height. We focus on the
$\kappa>r_{0}c_{i}$ solutions for which we have the most detailed analytic
estimates. For $\kappa<r_{0}c_{i}$, $\dot{M}_{w}\simeq\dot{M}_{max}$ (eq. 47),
but we do not have a prediction for the terminal velocity in this regime.

Overall, Figure 7 shows that the analytics and numerics agree well using our
analytically estimated base CR scale-height (eq. 27). The latter in turn
agrees well with the full numerical simulations (bottom panel of Fig. 7).
There is a slight systematic offset in the mean analytical estimates of
$h_{c}$ and $\dot{M}$ that could be removed by multiplying equation 27 by
$\simeq 0.75$.

The agreement in Figure 7 holds over a factor of $\sim 30$ in CR diffusion
coefficient, $\sim 100$ in base CR pressure, and $\sim 10^{4}$ in the ratio of
gravity to gas pressure ($\propto V_{g}^{2}$) which is also a proxy for galaxy
mass; see Table LABEL:tab:compare for the full range of simulations. Figure 7
also compares the analytic estimate of the wind terminal velocity (eq. 35) to
the speed in the simulations at the top of the box. There is again reasonably
good agreement, though the analytic speeds tend to be somewhat ($\sim 30\%$)
larger than the simulations. This is primarily because some of the energy
remains in the CRs in the simulations given the finite outer radius of the
computational domain, while the analytic estimate assumes that all of the CR
energy has been transferred to the gas.

##  4 Models Calibrated to Gamma-ray Observations

The analytic and numerical calculations in §2 and 3 show how the properties of
galactic winds driven by CR diffusion depend on the diffusion coefficient. For
example, winds lose most of their energy to gravity if $\kappa\lesssim
r_{0}c_{i}$ (Table LABEL:tab:compare), while if $\kappa\gtrsim r_{0}c_{i}$,
the terminal speed of the wind $\propto\kappa^{1/2}$ (eq. 35) and the base CR
pressure and the wind mass-loss rate depend on $\kappa$ for a given star
formation rate (eq. 24). To assess the implications of our results for CR
driven galactic winds we must thus have a handle on $\kappa$ in other
galaxies. Unfortunately, however, there is still sufficient uncertainty in CR
transport that it is non-trivial to determine from first principles if the
transport is indeed diffusive (vs. streaming), let alone the value of the
diffusion coefficient in different phases of the ISM and for the full range of
conditions realized in galaxy formation. As a result, we appeal instead to
observations to calibrate physically reasonable diffusion coefficients (see
also Chan et al. 2019; Hopkins et al. 2020) and then use those, together with
the results of §2 and 3, to quantify the implications for CR-driven galactic
wind models.

The non-thermal radio and gamma-ray emission from galaxies provide direct
observational constraints on the properties of CRs in external galaxies. These
in turn inform the CR pressure and diffusion coefficient that are critical for
setting the strength of CR-driven galactic winds. In this section, we provide
simple estimates to elucidate these constraints. We focus on the gamma-ray
emission from pion decay in star forming galaxies observed by Fermi (Abdo et
al., 2010b, a; Ackermann et al., 2012), rather than non-thermal radio
emission, despite the fact that there are far more observations available for
the latter (see Lacki et al. 2010). The reason is that the Fermi data directly
constrains the properties of the GeV protons that dominate the CR energy
density.

We consider a simple one-zone model in which a galaxy is characterized by its
size $r_{0}$, gas surface density $\Sigma_{g}$, and isothermal/turbulent
velocity $c_{i}$. The cosmic-ray scale-height is $H_{c}$ while the gas scale-
height is $H_{g}$. For CR diffusion we expect $H_{c}>H_{g}$, as is indeed the
case in our simulations in §3. The gamma-ray luminosity is set by the ratio of
the timescale for pion losses $t_{\pi}\simeq t_{0}\,(n_{\pi}/n_{eff})$ (where
$t_{0}=5\times 10^{7}$ yrs and $n_{\pi}=1$ cm-3 are set by the cross-section
for pion production and the effective density $n_{eff}$ is defined below) to
the timescale for CRs to escape by diffusion, $t_{\rm diff}\simeq
H_{c}^{2}/\kappa$. The effective density $n_{eff}$ is the density averaged
over the region the CRs are diffusing through and so is given by $\simeq
nH_{g}/H_{c}$ where $n$ is the mid-plane density of the galaxy. For reasons
that will become clear, we choose to express $t_{\pi}$ in terms of gas surface
density by using $\Sigma_{g}\simeq 2\rho H_{g}$. The ratio of the diffusion
timescale to the pion loss timescale is then

$\frac{t_{\rm
diff}}{t_{\pi}}\simeq\frac{\Sigma_{g}H_{c}}{2\kappa\rho_{\pi}t_{0}}$ (48)

where $\rho_{\pi}=n_{\pi}m_{p}\simeq 1.67\times 10^{-24}$ cm-3. Equation 48
shows that for a fixed ratio of hadronic losses ($t_{\pi}$) to CR escape
($t_{\rm diff}$) there is a degeneracy between the CR scale-height and
diffusion coefficient, with $\kappa\propto H_{c}$. This is consistent with the
known degeneracy between CR diffusion coefficient and ‘halo size’ in the
literature (e.g., Fig. 3 of Trotta et al. 2011).

The gamma-ray emission from pion decay in a star-forming galaxy is given by
$L_{\gamma}\simeq 1/3\dot{E}_{\pi}$ where
$\dot{E}_{\pi}\simeq\dot{E}_{c}\min(t_{\rm diff}/t_{\pi},1)$ is the rate of
energy loss to pion decay and the factor of $1/3$ quantifies the fraction of
pion decay in neutral vs. charged pions. Defining
$\dot{E}_{c}=E_{cr}\dot{M}_{*}/m_{*}$ with888Note that $\epsilon_{c}$ defined
in equation 30 is given by $\epsilon_{c}=E_{cr}/m_{*}c^{2}$. In this section
it is convenient to separate out $\epsilon_{c}$ into a part that depends on
the initial mass function ($m_{*}$) and a part that depends on the CR energy
supplied per SN ($E_{cr}$). $E_{cr}=10^{50}E_{cr,50}$ the energy per SN going
into CRs and $m_{*}\simeq 100M_{\odot}$ the total stellar mass formed per
core-collapse SNe, we find (Lacki et al., 2011)

$L_{\gamma}\simeq A_{\gamma}\,L_{*}\min\left(\frac{t_{\rm
diff}}{t_{\pi}},1\right)$ (49)

where $L_{*}\simeq\epsilon_{*}\dot{M}_{*}c^{2}$ is the total luminosity
produced by star formation. The factor $A_{\gamma}\simeq 3.3\times
10^{-4}E_{cr,50}/(\epsilon_{*}m_{*}17M_{\odot}^{-1})$ quantifies the maximum
possible gamma-ray luminosity from pion decay, which is realized in the limit
$t_{\pi}<t_{\rm diff}$ when all of the CR proton energy is lost to pion decay
before the CRs can escape (Thompson et al., 2007; Lacki et al., 2011). Note
that the factor $\epsilon_{*}m_{*}17M_{\odot}^{-1}\simeq 1$ is relatively
independent of the stellar initial mass function because massive stars set the
SN rate and the luminosity of a star-forming population. To good accuracy, we
can thus take $A_{\gamma}\simeq 3.3\times 10^{-4}E_{cr,50}$.

Figure 8: Empirically constrained cosmic-ray properties in star-forming
galaxies inferred from gamma-ray observations of pion decay by Fermi, as a
function of galaxy star formation rate. All calculations assume gas isothermal
sound speed $c_{i}=10\,{\rm km\,\,s}^{-1}$; scalings to other values of
$c_{i}$ are shown in Figure 9 and discussed in the text. Star formation rates
of example galaxies are indicated near the x-axis. Three examples are
considered corresponding to nuclear starbursts (black lines; $r_{0}=$ 0.3 kpc,
$V_{g}=150\,{\rm km\,\,s}^{-1}$), dwarf galaxies (red lines; $r_{0}=1$ kpc,
$V_{g}=50\,{\rm km\,\,s}^{-1}$), and star-forming disc galaxies (blue lines;
$r_{0}=3$ kpc, $V_{g}=150\,{\rm km\,\,s}^{-1}$). Left: $t_{\rm
diff}/t_{\pi}\propto L_{\gamma}/L_{*}$ is the ratio of the CR diffusion time
to the pion loss time and sets the gamma-ray luminosity of the galaxy (eq.
49). We infer $t_{\rm diff}/t_{\pi}$ for different galaxy star formation rates
using the Fermi correlation between $L_{\gamma}$ and star formation rate. We
then calculate the CR diffusion coefficient $\kappa$ and the dimensionless CR
scale-height $h_{c}=H_{c}/r_{0}$ using eq. 48, eq. 32 and eq. 27. Middle: The
fractional contribution of CR pressure to the pressure in the galactic disc
($p_{HE}=\pi G\Sigma_{g}^{2}\phi)$ depends on the assumed size of the star-
forming disc $r_{0}$, with nuclear starbursts having suppressed $p_{c}/p_{HE}$
because of rapid pion losses and CR diffusion in the dense nuclear regions.
The predicted mass-loss rate relative to the star formation rate
$\dot{M}_{w}/\dot{M}_{*}$ is largest for the disc galaxy model. Right:
Asymptotic wind speed $v_{\infty}$ and momentum flux in the wind $\dot{p}_{w}$
relative to the stellar radiation field ($\dot{p}_{*}=L_{*}/c$)

Observations of star-forming galaxies by Fermi show that there is a
correlation between gamma-ray luminosity and star formation rate (or infrared
luminosity; Ackermann et al. 2012; Griffin et al. 2016; Linden 2017) and a
correlation between gamma-ray luminosity and gas surface density (Lacki et
al., 2011). Both correlations imply that high star formation rate and high gas
surface density systems approach the ‘proton calorimeter’ limit (Pohl, 1994;
Thompson et al., 2007) in which most of the CR proton energy is lost to pion
decay before the CRs escape the galaxy. Because the correlation between gamma-
ray luminosity and infrared luminosity is better constrained than the
correlation between gamma-ray luminosity and gas surface density, we use the
former to calibrate the diffusion coefficients in our models:
$L_{\gamma}\simeq 2.3\times 10^{4}L_{\odot}(L_{TIR}/10^{9}L_{\odot})^{1.25}$
(Griffin et al., 2016) where $L_{TIR}$ is the total infrared luminosity from
$0.1-1000\,\mu m$.

At high infrared luminosities, the infrared luminosity is linearly
proportional to the star formation rate but this is not true for lower
infrared luminosities where the ultraviolet radiation makes an increasingly
important contribution to the total radiated starlight from galaxies. We
correct for this using Bell (2003) who finds

$\dot{M}_{*}\simeq 0.12\,{\rm
M_{\odot}\,yr^{-1}}\left(\frac{L_{TIR}}{10^{9}L_{\odot}}\right)\left(1+\sqrt{\frac{10^{9}L_{\odot}}{L_{TIR}}}\right)$
(50)

Given $L_{\gamma}(L_{TIR})$ and $L_{TIR}(\dot{M}_{*})$, equation 49 allows us
to infer $t_{\rm diff}/t_{\pi}$ as a function of star formation rate. This is
shown with the purple line in the left panel Figure 8: $t_{\rm diff}\simeq
t_{\pi}$ at $\dot{M}_{*}\sim 10^{3}\,{\rm M_{\odot}\,yr^{-1}}$ while $t_{\rm
diff}\ll t_{\pi}$ for much lower star formation rates. The latter is a direct
consequence of the fact that $L_{\gamma}\ll A_{\gamma}L_{*}$ in systems like
the Milky Way, M31, and the Magellanic clouds (e.g., Ackermann et al. 2012) so
that CR protons escape before losing most of their energy to pion decay.

Figure 9: Empirically constrained mass-loss rates, terminal velocities, and
momentum fluxes in galactic winds, as a function of galaxy star formation
rate, for different values of the gas isothermal sound speed $c_{i}$. The
nuclear starburst model (top) assumes $r_{0}=$ 0.3 kpc, $V_{g}=150\,{\rm
km\,\,s}^{-1}$, and $\phi=1$ while the dwarf galaxy model assumes $r_{0}=1$
kpc, $V_{g}=50\,{\rm km\,\,s}^{-1}$, and $\phi=5$. Note the different x-axis
range and normalization of $v_{\infty}$ for the two panels. The results for
the star-forming disc galaxy model shown in Figure 8 are nearly independent of
$c_{i}$ (see text) and so are not plotted here. Starburst mass-loss rates are
$\ll\dot{M}_{*}$ independent of $c_{i}$ while dwarf galaxy mass-loss rates can
reach $\sim\dot{M}_{*}$ for larger values of $c_{i}$. For starbursts with high
star formation rates, the terminal speed and momentum flux are best
interpreted as upper limits (see §4 for details).

Given $t_{\rm diff}/t_{\pi}$ as a function of galaxy star formation rate
(purple line Fig. 8), we now use equation 48 to constrain the CR diffusion
coefficient. To do so, however, we need an estimate of the CR scale-height and
gas surface density as a function of star formation rate. We model
$h_{c}=H_{c}/r_{0}$ using equation 27 which assumes that the CR scale-height
is set by diffusion in a CR-driven galactic wind. We model the gas surface
density $\Sigma_{g}$ using equation 32 which yields $\dot{M}_{*}\simeq\pi
r_{0}^{2}\sqrt{8}\pi G\Sigma_{g}^{2}\phi/v_{*}$ where $r_{0}$ is the size of
the galactic disc. The free parameters of our model are thus the size of the
star-forming galactic disc $r_{0}$, the galaxy circular velocity set by
$V_{g}$, the gas isothermal sound speed $c_{i}$, and the dimensionless stellar
contribution to the gravitational potential $\phi$ (as well as several
‘microphysics’ parameters such as the CR energy per supernovae). Given choices
for these parameters, as well as the observed $L_{\gamma}-L_{*}$ correlation,
we solve equation 27, 32, and 48 for $\kappa$, $h_{c}$, and $\Sigma_{g}$. In
what follows, we initially assume $c_{i}=10\,{\rm km\,\,s}^{-1}$ and consider
parameters approximating 3 different classes of galaxies. Results for these
classes are given in Figure 8:

  1. 1.

Nuclear starbursts: black lines in Fig. 8; $r_{0}=$ 0.3 kpc, $V_{g}=150\,{\rm
km\,\,s}^{-1}$, $\phi=1$. This model is meant to approximate the nuclear
starbursts M82 and NGC 253 with $\dot{M}_{*}\simeq 10$ M⊙ yr-1, as well as
ultra-luminous starbursts like Arp 220 with $\dot{M_{*}}\sim 10^{2}$ M⊙ yr-1.

  2. 2.

Dwarf galaxies: red lines in Fig. 8; $r_{0}=1$ kpc, $V_{g}=50\,{\rm
km\,\,s}^{-1}$, $\phi=5$. This model is meant to approximate normal star-
forming dwarf galaxies like the LMC or SMC and their potentially much more
rapidly star-forming counterparts.

  3. 3.

Star-forming disc galaxies: blue lines in Fig. 8; $r_{0}=3$ kpc,
$V_{g}=150\,{\rm km\,\,s}^{-1}$, $\phi=5$. This model is meant to approximate
the Milky Way, M31, and other star-forming spirals in the local
($\dot{M}_{*}\sim\,{\rm M_{\odot}\,yr^{-1}}$) and high-redshift
($\dot{M}_{*}\sim 10-100\,{\rm M_{\odot}\,yr^{-1}}$) universe.

We discuss below (Fig. 9 and eqs. 55-65) the scaling of our results to other
gas sound speeds $c_{i}$.

The dotted lines in the left panel of Figure 8 show the CR diffusion
coefficient in units of $10^{29}$ cm2 s-1 required to explain the gamma-ray
luminosities of star-forming galaxies, per the method described in the
previous paragraph. Figure 8 also shows the dimensionless CR scale-height
$h_{c}$ we derive (dashed lines). For the dwarf and starburst models in Figure
8, $h_{c}\simeq 1$ and the diffusion coefficient is $\sim 2\times 10^{29}$ cm2
s-1 $\gg r_{0}c_{i}$ with only a modest factor of few variation with star
formation rate. For the star forming disc model, however, $h_{c}\sim 0.1$ and
$\kappa\sim r_{0}c_{i}\sim 10^{28}$ cm2 s-1. The dimensionless scale-height
for our ‘disc’ model is comparable to that inferred from synchrotron emission
in nearby star-forming disc galaxies (e.g., Krause et al. 2018), though the
latter depends on both the magnetic field and CR scale-heights. Our inferred
values of $h_{c}\sim 0.1$ and $\kappa\sim 10^{28}$ cm2 s-1 for disc galaxies
in Figure 8 are, however, factors of $\sim 5$ smaller than preferred in
phenomenological models constrained by Milky Way CR data (e.g., Trotta et al.
2011). However, the estimated gamma-ray luminosity of the Milky-Way $\simeq
8\times 10^{38}$ erg s-1 (Strong et al., 2010) is a factor of $\simeq 3$ lower
than what would be predicted by its infrared luminosity given the correlations
from Ackermann et al. (2012) and Griffin et al. (2016) used here. $h_{c}$ and
$\kappa$ would then be 3 and 10 times larger, respectively (see eqs 55 & 61
below), in better agreement with detailed Milky-Way modeling. In addition, we
show below that for $h_{c}\lesssim 1$, $h_{c}\propto c_{i}$ and $\kappa\propto
c_{i}$ (eqs. 55 & 61). Our disc galaxy model in Figure 8 would thus have a
larger scale-height and diffusion coefficient if we assumed that the CRs
primarily coupled to the warm-hot phase of the ISM, as is quite plausible.

The key dimensionless number from §3 that determines the properties of CR
driven galactic winds is $\kappa/(r_{0}c_{i})$; winds accelerate rapidly and
reach speeds significantly larger than $V_{g}$ if $\kappa/(r_{0}c_{i})\gtrsim
1$ (Fig. 2). For our 3 galaxy models in Figure 8, we note that
$r_{0}c_{i}\simeq 10^{27}$ cm2 s-1 (starburst; black lines), $r_{0}c_{i}\simeq
3\times 10^{27}$ cm2 s-1 (dwarf; red lines), and $r_{0}c_{i}\simeq 10^{28}$
cm2 s-1 (disc; blue lines). A key conclusion from Figure 8 is that the gamma-
ray data on star-forming galaxies are consistent with diffusion coefficients
that are in the regime of $\kappa\gtrsim r_{0}c_{i}$, though only marginally
so for our star-forming disc model. This does not prove that such diffusion
coefficients are correct, at a minimum because it is possible that CR escape
is not set by diffusion but rather by streaming or advection in winds driven
by other mechanisms (e.g., supernovae), in which case our constraint on the CR
diffusion coefficient using equation 49 would not apply. But the diffusion
coefficients inferred in Figure 8 are nonetheless a useful and instructive
observational check on diffusive CR-driven galactic wind models.

Given the diffusion coefficients inferred in the left panel of Figure 8, we
can estimate the CR energy density in the galactic disc - which sets the base
conditions for the wind - as follows: if the star formation rate per unit area
is $\dot{\Sigma}_{*}$, the CR pressure in the galaxy is given by (e.g.,
Thompson & Lacki 2013)999Note that this expression for $p_{c,0}$ is larger
than its spherical counterpart used in §2 (eq. 9) by a factor of 2. The
difference is the surface area of a sphere ($4\pi r_{0}^{2}$) vs. that of the
bottom and top of the disk ($2\times\pi r_{0}^{2}$). We use the disk version
in this section but our results are not sensitive to factor of 2 changes in
$p_{c,0}$.

$p_{c,0}\simeq\frac{\dot{\Sigma}_{*}E_{cr}}{6H_{c}m_{*}}{\rm min}(t_{\rm
diff},t_{\pi})$ (51)

where ${\rm min}(t_{\rm diff},t_{\pi})$ determines the effective loss/escape
time for CRs in the galaxy. Using equation 32 and that fact that $t_{\rm
diff}\lesssim t_{\pi}$ even at the highest star formation rates in Figure 8,
we find:

$\frac{p_{c,0}}{\pi
G\Sigma^{2}_{g}\phi}\simeq\frac{\sqrt{2}E_{cr}h_{c}r_{0}}{3\kappa p_{*}}$ (52)

The dashed lines in Figure 8 (middle panel) shows the resulting base CR
pressures for our 3 galaxy classes/models. For Milky Way-like “star-forming
disc” conditions with $r_{0}\sim 3$ kpc, we find that $p_{c,0}\sim 0.1\pi
G\Sigma_{g}^{2}$, a bit smaller than local measurements in the Milky Way
(Boulares & Cox, 1990), though we stress that our model is not intended to
reproduce solar-circle measurements, but rather approximate the disc-averaged
physical conditions. Figure 8 predicts that Milky Way-like galaxies have
roughly the largest fraction of disc pressure support from CRs, though for a
given value of $r_{0}$, the base CR pressure only varies by a factor of a few
over a factor of $\sim 10^{4}$ in star formation rate. The base CR pressure
is, however, sensitive to $r_{0}$, and is significantly smaller, $\sim$ few
$10^{-3}$, for nuclear starbursts like M82 and Arp 220. Physically this is
because for a given star formation rate, a smaller size for the star-forming
disc implies higher gas densities and thus more rapid pion losses. In order to
be compatible with the gamma-ray observations, the diffusion time must be
correspondingly shorter as well and thus the CR pressure cannot build to as
large a value. These conclusions are qualitatively similar to those of Lacki
et al. (2010); Lacki et al. (2011) (see also Crocker et al. 2021, 2020) who
developed one-zone galaxy models with CRs in order to reproduce the far-
infrared radio correlation. They concluded that $p_{c,0}\sim\pi
G\Sigma_{g}^{2}$ at low gas-surface densities, but that $p_{c,0}\ll\pi
G\Sigma_{g}^{2}$ at the high gas densities of nuclear starbursts (see Fig. 15
of Lacki et al. 2010 and section 6.3 of Lacki et al. 2011).

Finally, if we combine equation 52 and equation 28 we arrive at a simple
expression for the mass-loss rate in CR-driven galactic winds:

$\frac{\dot{M}_{w}}{\dot{M}_{*}}\simeq\frac{E_{cr}/m_{*}}{{\rm
V^{2}_{g,eff}}}\frac{c_{i}r_{0}}{3\kappa}\left(\frac{4h_{c}V_{g}^{4}}{c_{c,0}^{2}c_{i}^{2}}\right)^{c_{i}^{2}/(2V_{g}^{2})}$
(53)

where $\sqrt{E_{cr}/m_{*}}$ is a velocity scale associated with CR feedback,
which is $\simeq 220\,{\rm km\,\,s}^{-1}$ for $E_{cr}=10^{50}$ erg and
$m_{*}=100M_{\odot}$.

The solid lines in the middle panel of Figure 8 shows the mass-loading of CR-
driven galactic winds $\dot{M}_{w}/\dot{M}_{*}$ from equation 53 using
$\kappa$ in the left panel of Figure 8 calibrated to gamma-ray observations.
We again show results for our 3 fiducial galaxy models. The mass-loss rates in
CR-driven winds are significant for a wide range of normal disc galaxy
conditions with $\dot{M}_{w}\sim\dot{M}_{*}$, but are strongly suppressed
($\dot{M}_{w}\sim 10^{-3}\dot{M}_{*}$) in nuclear starbursts like M82, NGC
253, and Arp 220 because rapid CR diffusion and pion losses suppress the base
CR pressure (Fig. 8) and thus the mass-loss rate in the wind. Dwarf galaxies
lie somewhere in between with $M_{w}\sim 0.2\dot{M}_{*}$.

The right panel of Figure 8 shows the terminal velocity (eq. 35) and momentum
flux $\dot{p}_{w}=\dot{M}_{w}v_{\infty}$ (in units of the photon momentum flux
$\dot{p}_{*}=L/c$) for the same three galaxy models (compare with Lochhaas et
al. 2020). For the disc and dwarf models $v_{\infty}\sim 1000\,{\rm
km\,\,s}^{-1}$ and $\dot{p}_{w}\sim\dot{p}_{*}$ while for the starburst model
$v_{\infty}\sim 10^{4}\,{\rm km\,\,s}^{-1}$ and $\dot{p}_{w}\sim
0.1\dot{p}_{*}$. The large velocities and correspondingly lower momentum
fluxes for the starburst model are a consequence of $\kappa\gg r_{0}c_{i}$
needed to avoid overproducing the gamma-ray luminosities.

One of the uncertainties in assessing the implications of our results for
observations is the appropriate value of the gas isothermal sound speed. Most
of the mass in the ISM is in cooler phases but most of the volume is in warmer
phases. As a result, it is plausible, though not guaranteed, that the warmer
phases set the scattering rate and diffusion coefficient for the cosmic-rays.
Figure 9 shows how the mass-loss rates, terminal velocities, and momentum
fluxes we infer from gamma-ray data in our starburst and dwarf models depend
on the assumed value of $c_{i}$. We do not show similar results for the normal
star-forming disc model because that model is in the regime $h_{c}\lesssim 1$
where the mass-loss rate, terminal velocity, and momentum flux given gamma-ray
inferred diffusion coefficients are a very weak function of $c_{i}$ (e.g., the
results in Figure 8 for the ‘disc’ model apply to better than a factor of 2
accuracy for $c_{i}\lesssim 100\,{\rm km\,\,s}^{-1}$); this is derived
analytically below. Figure 9 shows that the mass-loss rate increases notably
with increasing $c_{i}$ for both our starburst and dwarf galaxy models. The
momentum flux increases more slowly with $c_{i}$ and the terminal velocity of
the wind decreases due to the larger mass-loadings. However, the qualitative
conclusions drawn from Figure 8 remain robust. Namely, for the starburst
models $\dot{M}_{w}\ll M_{*}$ and for the dwarf models $\dot{M}_{w}$ is at
most $\sim{\rm few\times}\,\dot{M}_{*}$. The latter is, however, still smaller
than the large mass-loadings in dwarf galaxies typically needed to reconcile
the galaxy stellar mass and dark matter halo mass functions (e.g., Somerville
& Davé 2015).

Luminous starbursts (e.g., Arp 220, and to a lesser extent M82) have gamma-ray
luminosities approaching the calorimeter limit $L_{\gamma}\simeq
A_{\gamma}L_{*}$ due to $t_{\pi}\lesssim t_{diff}$ (e.g., Lacki et al. 2011;
Ackermann et al. 2012; Griffin et al. 2016). In the calorimeter limit, the
constraints on $\kappa$ in Figure 8 are best interpreted as upper limits,
since the gamma-ray emission is roughly independent of $\kappa$ for
$t_{\pi}\lesssim t_{diff}$. In this regime, the estimated mass-loss rate is
independent of $\kappa$ because the base CR pressure is set by $t_{\pi}$ in
equation 52 rather than $t_{diff}$. However, because the asymptotic wind speed
is $\propto\kappa^{1/2}$ (eq. 59), in the calorimeter limit, we can only
empirically derive an upper limit on $v_{\infty}$ and the asymptotic wind
kinetic energy and momentum flux. Thus, particularly for the more luminous
starbursts in Figures 8 & 9, our results may be better interpreted as upper
limits on the wind terminal velocity and momentum flux. This further
strengthens our conclusion that winds due to cosmic-rays alone are weak in
starburst galaxies and cannot drive their exceptional outflows (e.g., see
Barcos-Muñoz et al. 2018)

We now derive analytic approximations to the results in Figures 8 and 9. These
are valuable because they show how the results depend on all of the physical
parameters of the problem. In the analytics we assume that $r_{s}\sim r_{0}$
(see eq. 21), i.e., that the factor
$\left(4h_{c}V_{g}^{4}/[c_{c,0}^{2}c_{i}^{2}]\right)^{c_{i}^{2}/(2V_{g}^{2})}\sim
1$. This is an excellent approximation for Figure 8 in which $c_{i}=10\,{\rm
km\,\,s}^{-1}$, but is less applicable for the largest values of $c_{i}$ in
Figure 9. In our analytic estimates, we also use a fit to our observational
calibration of the CR diffusion timescale. An approximate fit to the results
in Figure 8 is given by

$\frac{t_{\rm diff}}{t_{\pi}}\simeq\alpha\ E_{cr,50}^{-1}\
\left(\frac{\dot{M}_{*}}{1\,\,{\rm M_{\odot}\,yr^{-1}}}\right)^{0.5}$ (54)

with $\alpha\equiv 0.07\alpha_{0.07}\simeq 0.07$. Equation 54 is accurate to
better than a factor of 2 over the entire range of $\dot{M}_{*}$ shown in
Figure 8. It is not asymptotically correct, however, at either high or low
star formation rates. For high star formation rates, $L_{\gamma}\propto
L_{TIR}^{1.25}\propto\dot{M}_{*}^{1.25}$. Thus
$L_{\gamma}/L_{*}\propto\dot{M}_{*}^{0.25}$ instead of
$\propto\dot{M}_{*}^{0.5}$. For low star formation rates, however,
$L_{TIR}\propto\dot{M}_{*}^{2}$ so that
$L_{\gamma}/L_{*}\propto\dot{M}_{*}^{1.5}$. In practice, we find that equation
54 is a good compromise, particularly given its simplicity. In particular,
following through the derivations of $h_{c}$, $\kappa$, $p_{c,0}$, and
$\dot{M}_{w}/\dot{M}_{*}$ per equations 27, 48, 52, and 53 we can show how the
results in Figure 8 depend on the various micro (CR and star formation) and
macro (global galaxy) parameters in the problem. Combining equations 27, 48, &
54 we find

$\begin{split}h_{c}\simeq&\min\bigg{[}1,\frac{0.9\,E_{cr,50}v_{*,8.5}^{1/2}}{\alpha_{0.07}\,\phi^{1/2}}\left(\frac{c_{i}}{10\,{\rm
km\,\,s}^{-1}}\right)\\\ &\times\,\left(\frac{\rm
kpc}{r_{0}}\right)\left(\frac{100\,{\rm
km\,\,s}^{-1}}{V_{g}}\right)^{2}\bigg{]}\end{split}$ (55)

where $v_{*,8.5}=v_{*}/3000\,{\rm km\,\,s}^{-1}$.

There are two regimes depending on whether $h_{c}<1$ or $h_{c}=1$ in equation
55. For $h_{c}\sim 1$, which is the regime appropriate for dwarf galaxies and
nuclear starbursts in our models in Figs 8 & 9,

$\kappa\simeq 3\times 10^{29}\,{\rm
cm^{2}\,s^{-1}}\,\left(\frac{h_{c}E_{cr,50}v_{*,8.5}^{1/2}}{\alpha_{0.07}\,\phi^{1/2}}\right),$
(56) $\frac{p_{c,0}}{\pi G\Sigma_{g}^{2}\phi}\simeq
0.03\,\alpha_{0.07}\,\phi^{1/2}v_{*,8.5}^{-3/2}\,m_{*,2}^{-1}\,\left(\frac{r_{0}}{3\,{\rm
kpc}}\right)$ (57) $\begin{split}\frac{\dot{M}_{w}}{\dot{M}_{*}}&\simeq
0.1\,\alpha_{0.07}\,\phi^{1/2}\,v_{*,8.5}^{-1/2}\,m_{*,2}^{-1}\\\
&\times\left(\frac{100\,{\rm km\,\,s}^{-1}}{{\rm
V_{g,eff}}}\right)^{2}\,\left(\frac{c_{i}r_{0}}{h_{c}30\,{\rm
km\,\,s}^{-1}\,{\rm kpc}}\right),\end{split}$ (58)

where $m_{*,2}=m_{*}/100M_{\odot}$. Equations 56 and 35 can also be combined
to estimate the terminal velocity of the wind:

$\begin{split}v_{\infty}&\simeq 1600\,\,{\rm
km\,\,s}^{-1}\,\alpha_{0.07}^{-1/2}\,E_{cr,50}^{1/2}\,\phi^{-1/4}\,v_{*,8.5}^{1/4}\\\
&\times\,\left(\frac{{\rm V_{g,eff}}}{100\,{\rm
kms^{-1}}}\right)\left(\frac{\rm 30\,h_{c}\,kpc\,{\rm
km\,s^{-1}}}{r_{0}c_{i}}\right)^{1/2}\end{split}$ (59)

Equations 58 and 59 then imply

$\begin{split}\frac{\dot{p}_{w}}{\dot{p}_{*}}&\simeq
1\,\frac{\alpha_{0.07}^{1/2}\,E_{cr,50}^{1/2}\,\phi^{1/4}}{m_{*,2}\,\epsilon_{*,-3.3}\,v_{*,8.5}^{1/4}}\\\
&\times\,\left(\frac{100\,{\rm km\,s^{-1}}}{{\rm
V_{g,eff}}}\right)\left(\frac{r_{0}c_{i}}{\rm 30\,h_{c}\,kpc\,{\rm
km\,s^{-1}}}\right)^{1/2}\end{split}$ (60)

In the opposite regime of $h_{c}\lesssim 1$, which is the regime of star-
forming disc galaxies in Figure 8, we find

$\begin{split}\kappa&\simeq 10^{29}\,{\rm
cm^{2}\,s^{-1}}\left(\frac{E_{cr,50}^{2}v_{*,8.5}}{\alpha_{0.07}^{2}\,\phi}\right)\,\left(\frac{c_{i}}{10\,{\rm
km\,\,s}^{-1}}\right)\\\ &\times\,\left(\frac{100\,{\rm
km\,\,s}^{-1}}{V_{g}}\right)^{2}\,\left(\frac{3\,\,{\rm
kpc}}{r_{0}}\right),\end{split}$ (61) $\frac{p_{c,0}}{\pi
G\Sigma_{g}^{2}\phi}\simeq
0.03\,\alpha_{0.07}\,\phi^{1/2}v_{*,8.5}^{-3/2}\,m_{*,2}^{-1}\,\left(\frac{r_{0}}{3\,{\rm
kpc}}\right)$ (62) $\frac{\dot{M}_{w}}{\dot{M}_{*}}\simeq
0.2\,\alpha^{2}_{0.07}\,\phi\,v_{*,8.5}\,m_{*,2}^{-1}\,\left(\frac{r_{0}}{3\,{\rm
kpc}}\right)^{2},$ (63) $v_{\infty}\simeq 1000\,\,{\rm
km\,\,s}^{-1}\,\alpha_{0.07}^{-1}\,E_{cr,50}\,\phi^{-1/2}\,v_{*,8.5}^{1/2}\,\left(\frac{3\,{\rm
kpc}}{r_{0}}\right)$ (64)

and

$\frac{\dot{p}_{w}}{\dot{p}_{*}}\simeq
1.3\frac{\alpha_{0.07}\,E_{cr,50}\,\phi^{1/2}\,v_{*,8.5}^{3/2}}{m_{*,2}\,\epsilon_{*,-3.3}}\,\left(\frac{r_{0}}{3\,{\rm
kpc}}\right)$ (65)

We reiterate that $\phi$ in equations 55-65 quantifies the contribution of
stars and dark matter to the gravitational acceleration with $\phi\sim 1$ for
gas-dominated systems and $\phi\sim 3-10$ for typical star-forming galaxies in
the local Universe. Also note that equations 57 and 62 are identical, i.e. the
same for $h_{c}<1$ and $h_{c}=1$. This is because $p_{c}$ depends only on the
ratio $H_{c}/\kappa$ (eq. 52) which is uniquely determined by our calibration
of $t_{\pi}/t_{\rm diff}$ (eqs. 48 & 54).

Equations 55-65 do a good job of reproducing many key results shown in Figure
8: (1) the characteristic diffusion coefficient $\sim 10^{28-29}$ cm2 s-1
relatively independent of galaxy properties required to reproduce the observed
gamma-ray-star formation rate correlation, and (2) a base CR pressure
$p_{c,0}$ and wind mass-loading $\dot{M}_{w}/\dot{M}_{*}$ that are relatively
independent of star formation rate, but a strong function of the size of the
star-forming disc $r_{0}$ at fixed star formation rate (because diffusion is
much more rapid for smaller $r_{0}$ and decreases $p_{c,0}$ and
$\dot{M}_{w}$). .

Equations 55-65 also elucidate how our gamma-ray inferences depend on the gas
isothermal sound speed. Figure 8 takes $c_{i}=10\,{\rm km\,\,s}^{-1}$. Our
star-forming disc model in Figure 8 is in the regime with $h_{c}\lesssim 1$,
so that equations 61-65 apply. In this regime $\kappa\propto c_{i}$ so that
$\kappa/r_{0}c_{i}\sim 1$ from Figure 8 is in fact true relatively independent
of $c_{i}$. Moreover, for $h_{c}\lesssim 1$, $\dot{M}_{w}$ is independent of
$c_{i}$ (eq. 63). The results in the middle and right panels of Figure 8 are
thus applicable over a wide range of $c_{i}$ (this is the reason that we do
not plot the ‘disc’ model in Fig. 9). Our conclusion that CR-driven winds
constrained by gamma-ray observations of disc galaxies are dynamically
important with $\dot{M}_{w}\sim\dot{M}_{*}$ is thus robust to uncertainties in
$c_{i}$.

In contrast to the disc model, the dwarf and starburst models in Figure 8 are
in the regime where $h_{c}\sim 1$ so that equations 56-60 apply. In this case
the inferred $\kappa$ is independent of $c_{i}$ (eq. 56); however, even
$c_{i}\sim 100\,{\rm km\,\,s}^{-1}$ still corresponds to
$\kappa/r_{0}c_{i}\gtrsim 1$ so that the qualitative physics of the wind does
not change. Equations 58, 59, & 60 predict that the inferred mass-loss rate,
terminal velocity, and momentum flux given gamma-ray constraints scale
$\propto c_{i},c_{i}^{-1/2}$ and $c_{i}^{1/2}$, respectively. This is
consistent with, though weaker than, the trends in Figure 9. The stronger
dependence on $c_{i}$ in Figure 9 is because for larger values of $c_{i}$ the
approximation
$\left(4h_{c}V_{g}^{4}/[c_{c,0}^{2}c_{i}^{2}]\right)^{c_{i}^{2}/(2V_{g}^{2})}\sim
1$ used to derive equations 56-65 is no longer as accurate. Physically, this
is because for larger values of $c_{i}$ and smaller $c_{c,0}$ the sonic point
$r_{s}$ is at somewhat larger radii $\sim{\rm few\times}r_{0}$, which
increases $\dot{M}_{w}$ (eq. 28) and decreases $v_{\infty}$ (eq. 35).

We conclude this section by reiterating that our most general expressions for
the base CR pressure $p_{c,0}$, mass-loss rate $\dot{M}_{w}$, terminal
velocity $v_{\infty}$, and momentum flux $\dot{p}_{w}$ in CR-driven galactic
winds are given in §2. Those results depend, however, on a theoretically and
observationally uncertain CR diffusion coefficient. In this section (and in
Figs 8 & 9) we have estimated the CR diffusion coefficient $\kappa$ using the
existing (but limited) gamma-ray data from Fermi on pion-decay in star-forming
galaxies, thus enabling more concrete predictions of the properties of CR-
driven galactic winds.

##  5 Summary and Discussion

The physics of cosmic ray (CR) transport in galaxies and in the circumgalactic
medium remains a significant uncertainty in assessing the impact of CRs on
galaxy formation. A central question is what determines the scattering mean
free path of CRs, and how this depends on local plasma conditions (e.g., Amato
& Blasi 2018; Hopkins et al. 2021b). In this paper, we have assumed that CR
transport can be modeled by a spatially independent diffusion coefficient. The
diffusion approximation for CR transport is particularly appropriate if
ambient turbulence scatters the CRs (vs. scattering by fluctuations excited by
the CRs themselves). A companion paper will consider the case of CR transport
mediated by the streaming instability. These two mechanisms of CR transport
differ dramatically in their predictions for how the CR pressure decreases
away from a galaxy: in the limit of rapid CR diffusion, $p_{c}\propto r^{-1}$
(eq. 8), i.e., the CR pressure scale-height is of order the size of the
system, while in the limit of rapid CR streaming, $p_{c}\propto\rho^{2/3}$
(for a split-monopole field geometry; e.g., Mao & Ostriker 2018) and so the CR
pressure scale-height is tied to that of the gas. This difference in the
dynamics of the CRs in general leads to significantly different wind
properties for the two CR transport models, as has been highlighted previously
in numerical simulations (e.g., Wiener et al. 2017; Chan et al. 2019). One aim
of this paper and its companion is to understand these differences
analytically and using idealized time-dependent numerical simulations, thus
elucidating how the properties of CR-driven galactic winds depend on global
galaxy properties and the physics of CR transport.

In this paper, we analytically estimated the properties of galactic-winds
driven by diffusion by assuming that the CR diffusion timescale is short
compared to the flow time (or dynamical time) near the base of the wind; this
requires CR diffusion coefficients $\kappa\gtrsim r_{0}c_{i}$ where $r_{0}$ is
the size of the galaxy (i.e., the star-forming disc) and $c_{i}$ is the gas
sound speed. In this limit, the asymptotic kinetic energy flux carried by the
wind is comparable to that supplied to the CRs at the base of the wind, i.e.,
the wind is energy conserving. The mass-loss rate of CR driven winds has the
form $\dot{M}_{w}\sim 2\pi r_{0}^{2}\,\rho_{0}\,c_{i}\,(c_{c,0}/V_{g})^{2}\sim
2\pi r_{0}^{2}p_{c,0}c_{i}/V_{g}^{2}$ (eq. 28; see Fig. 1), and the asymptotic
wind speed is $V_{\infty}\simeq 2V_{g}(3\kappa/r_{0}c_{i})^{1/2}$ (eq. 35)
where $\rho_{0}$, $p_{c,0}$ and $c_{c,0}$ are the gas density, CR pressure,
and CR sound speed at the base of the outflow and $\sqrt{2}V_{g}$ is the
rotation velocity of the galaxy. Equation 31 compares this estimate of the
mass-loss rate in CR-driven winds to the galaxy star formation rate, with
$\dot{M}_{w}/\dot{M}_{*}\propto 1/\kappa$. Physically, for a given rate of CR
production, set by the star formation rate, the CR pressure in the galaxy, and
thus the strength of the wind, decreases with increasing diffusion coefficient
since the CRs escape the galaxy more rapidly.

In addition to our analytic estimates, we also carried out time-dependent
spherically symmetric simulations of CR-driven winds using the two-moment CR
transport scheme for Athena++ developed by Jiang & Oh (2018). The simulations
show that, for $\kappa\gtrsim r_{0}c_{i}$, the analytic estimates for the
mass-loss rate, terminal speed, and CR scale-height near the base of the wind
are accurate to $\sim 50\%$ over a factor of $\sim 30$ in CR diffusion
coefficient, $\sim 30$ in base CR pressure, and $\sim 100$ in the ratio of the
escape speed to the gas sound speed (Fig. 7; see Table LABEL:tab:compare for
the full range of simulations). In addition, the simulations show that there
is a critical value of the CR diffusion coefficient $\kappa\simeq r_{0}c_{i}$
below which the character of the solution changes considerably. For
$\kappa\lesssim r_{0}c_{i}$, CR-driven winds accelerate much more slowly and
are nearly hydrostatic over a very extended radial range. In this regime most
of the energy supplied to CRs at the base of the wind goes into work against
gravity expanding to large radii (Fig. 6): the asymptotic kinetic energy flux
in the wind is only a small fraction of that initially supplied to the CRs
(see the last column of Table LABEL:tab:compare). These low $\kappa$ solutions
are CR analogues of photon-tired stellar winds (Owocki & Gayley, 1997). The
mass-loss rate in this regime can be accurately estimated from global energy
conservation as $\dot{M}_{w}\simeq\dot{M}_{max}\simeq
2\dot{E}_{c}/v_{esc}^{2}$ (eqs. 40 & 47), where $\dot{E}_{c}$ is the energy
per unit time supplied to CRs at the base of the wind. This maximum possible
mass-loss rate in CR-driven winds is quite large,
$\simeq\dot{M}_{*}(300\,\,{\rm km\,\,s}^{-1}/v_{esc})^{2}$ (eq. 39). For
$\kappa>r_{0}c_{i}$, however, the actual outflow rate is much less than this
maximal value (eq. 42).

A key difference between our treatment of CR-driven winds in this paper and
analogous treatments of stellar winds driven by radiation in the diffusion
approximation (e.g., Owocki et al. 2017) is that stellar wind theory is
typically formulated in terms of a given photon-matter cross section $\sigma$,
which sets the Eddington luminosity. By contrast, here we are considering a
fixed CR diffusion coefficient, equivalent to a fixed value of the mean-free
path $1/(\sigma\rho)$. This difference means that many solutions in stellar
wind theory do not directly carry over to the CR problem, although many of the
important concepts do.

In our models with $\kappa\gtrsim r_{0}c_{i}$, the properties of CR-driven
winds are largely set close to the ‘base’ of the wind, i.e., near the galaxy.
In particular, the sonic point - which sets the mass-loss rate - is close to
the base of the wind (eq. 21) unless $c_{i}\sim V_{g}$ and the energy flux in
the wind - which sets the terminal velocity - is set by the CR diffusive flux
at the base (eq. 35 and associated discussion). As a result, we suspect that
the properties of these solutions are unlikely to be that sensitive to spatial
variation in the CR diffusion coefficient unless there are large variations at
small radii near the sonic point. By contrast, our solutions with
$\kappa\lesssim r_{0}c_{i}$ accelerate much more slowly (Fig. 6) and are
likely much more sensitive to spatial variation in the microphysics of CR
transport. In addition, because the low $\kappa$ solutions have a kinetic
power $\dot{E}_{k}$ at large radii that is small compared to the cosmic ray
power at the base of the wind, they are likely more sensitive to the ambient
pressure in the CGM, which could confine lower $\dot{E}_{k}$ outflows.

Our time-dependent simulations allow us to study the stability of the analytic
steady state wind solutions. Nearly all of our simulations reach a laminar
steady state with no evidence of instability. This is at first glance
surprising since Drury & Falle (1986) showed that CR diffusion in the presence
of a background CR pressure gradient renders sound waves linearly unstable. We
show in Appendix A, however, that the growth rate of the sound wave
instability is not fast enough compared to the flow time in the wind for the
instability to grow significantly; the one exception to this is our lowest gas
sound speed simulation (the $V_{g}=200c_{i}$ simulation in Table
LABEL:tab:compare; see Figure 10). Appendix A also carries out a WKB linear
stability calculation (neglecting the background cosmic-ray pressure gradient)
for the two-moment CR transport scheme used in our simulations, and shows that
sound waves and entropy modes are linearly stable in the presence of CR
diffusion, consistent with the steady state solutions found in the
simulations.

A key parameter that sets the strength of the galactic wind in our models is
the CR pressure in the bulk of the ISM (with $\dot{M}_{w}\propto p_{c,0}$). If
$p_{c,0}\sim\pi G\Sigma_{g}^{2}\phi$ (the pressure required for hydrostatic
equilibrium), CR-driven winds will have dynamically important mass-loss rates
with $\dot{M}_{w}\gtrsim\dot{M}_{*}$ (eq. 29). If, however, $p_{c,0}\ll\pi
G\Sigma_{g}^{2}\phi$, then since $\dot{M}_{w}\propto p_{c,0}$ (eq. 28), the
mass-loss rates will be significantly smaller. The equilibrium CR pressure
$p_{c,0}$ is in turn set by CR escape (i.e., the diffusion coefficient
$\kappa$) and/or hadronic losses (eq. 51). To assess the implications of our
results for the role of CRs in driving galactic winds, it is thus necessary to
estimate the CR diffusion coefficient in other galaxies. This remains a
daunting task from first principles, so we instead turned to observations (see
§4). In particular, observations of the non-thermal emission from CRs in other
galaxies provide direct constraints on CR diffusion coefficients and the CR
pressure in galaxies (e.g., Lacki et al. 2010; Lacki et al. 2011; Crocker et
al. 2021). The non-thermal gamma-ray emission from neutral pion decay is
particularly important in this regard because (1) it constrains the properties
of CR protons (vs. synchrotron emission), and (2) observations at GeV energies
by Fermi, though modest in number, directly constrain the CRs that dominate
the total CR pressure. In §4 we developed a simple analytic model interpreting
gamma-ray observations in the context of diffusive CR transport. This model
essentially derives the theoretically and observationally uncertain diffusion
coefficient as a function of the observed gamma-ray luminosity of galaxies. We
find that a model with a diffusion coefficient $\sim 10^{28-29}$ cm2 s-1 (Fig.
8 and eqs. 56 & 61) is consistent with the Fermi data on gamma-ray emission
from star-forming galaxies. This is consistent with similar estimates by Chan
et al. (2019) and Hopkins et al. (2020) and their more detailed numerical
calculations.

Our constraint on the diffusion coefficient in other galaxies also translates
into an estimate of the CR pressure in galactic discs. For typical star
forming galaxies with disc sizes $r_{0}\sim 3$ kpc, we find that the CR
pressure is of order 10% of the pressure required for vertical hydrostatic
equilibrium in the disc (Fig. 8 and eq. 57). This is reasonably consistent
with Milky Way measurements. However, the fractional contribution of CRs to
pressure support in the disc is $\propto r_{0}$ and is only $\sim 10^{-2.5}$
for typical nuclear starburst conditions (Fig. 8 and eq. 57). Physically, this
is because in more compact star-forming regions, the gas densities are higher
and thus pion losses are stronger. In addition, the CR diffusion time is
shorter. There is thus less time for the CR pressure to build up and so the
equilibrium CR pressure in the disc is lower. These conclusions are consistent
with the earlier work of Lacki et al. (2010) based on modeling the far
infrared-radio correlation.

Our results on the CR diffusion coefficient and CR pressure implied by gamma-
ray observations can be used to estimate the properties of CR-driven galactic
winds across a wide range of galaxies. The middle panel of Figure 8 plots the
resulting ratio of the CR-driven mass-loss rate to the star formation rate for
three fiducial galaxy models, while the right panel shows the terminal
velocity and momentum flux of the wind. We find that for massive star-forming
disc galaxies, the mass-loss rates are of order the star formation rate,
momentum fluxes are of order $\dot{p}_{*}=L/c$, and terminal velocities are
$\sim 500\,{\rm km\,\,s}^{-1}$ (a few times the circular velocity). For lower-
mass dwarf galaxies, however, we find that CRs are somewhat less efficient at
driving winds ($\dot{M}_{w}\sim 0.2\dot{M}_{*}$ and $\dot{p}_{w}\sim L/c$),
primarily because the CR diffusion time is so short (to explain the gamma-ray
data) that the CR pressure in the disc is comparatively low. This is even more
true in nuclear starbursts: CRs become much less efficient at driving winds
with $\dot{M}_{w}/\dot{M}_{*}\propto r_{0}$ (eq. 58) and
$\dot{M}_{w}/\dot{M}_{*}\sim 10^{-2}-10^{-3}$ for well-studied local
starbursts like M82 and Arp 220 (Fig. 8). This conclusion fundamentally rests
on our inference that $p_{c,0}\ll\pi G\Sigma_{g}^{2}\phi$ given CR diffusion
coefficients and pion loss timescales motivated by gamma-ray observations. An
independent observational probe of the CR proton pressure in other galaxies
would be a valuable test of our models.

One of the uncertain parameters in applying our results to observations is the
appropriate isothermal gas sound speed. This depends on the phase of the ISM
that the cosmic-rays most effectively couple to. Figure 8 assumes
$c_{i}=10\,{\rm km\,\,s}^{-1}$, which is an appropriate mass-averaged value in
the Milky Way. For typical star-forming disc galaxy parameters, we find that
the properties of the winds using gamma-ray constrained diffusion coefficients
are weakly dependent on $c_{i}$ (also derived analytically in equations
63-65). Our conclusion that cosmic-rays are a significant source of winds in
normal disc galaxies is thus reasonably robust to the uncertainty of the phase
of the ISM that primarily determines CR transport.

Figure 9 shows our gamma-ray inferred wind properties for dwarf galaxy and
nuclear starburst models for larger values of $c_{i}$, appropriate if cosmic-
rays primarily couple to volume filling warm-hot gas. For these galaxy models,
the mass-loss rate can increase significantly for larger values of $c_{i}$, as
does the momentum flux in the wind; the terminal speed of the wind is
correspondingly smaller for larger $c_{i}$. However, our general conclusions
are reasonably robust to uncertainties in $c_{i}$: the mass-loss rates due to
CRs alone in starburst galaxies are $\ll\dot{M}_{*}$ and in dwarf galaxies are
at most $\sim{\rm few}\times\dot{M}_{*}$. The latter is still below what is
typically needed to reconcile the stellar and dark matter halo mass functions
(see, e.g., Muratov et al. 2015 Table 3).

As noted earlier in the discussion, the maximum mass-loss rate in CR-driven
winds allowed by energy conservation is appreciable, $\dot{M}_{\rm
max}\simeq\dot{M}_{*}(300\,\,{\rm km\,\,s}^{-1}/v_{esc})^{2}$ (eq. 39). Mass-
loss rates $\sim\dot{M}_{\rm max}$ would be particularly important in dwarf
galaxies. However, these large mass-loss rates are only realized when the
outflow is very slow and most of the energy supplied to CRs by star formation
goes into work leaving the gravitational potential of the galaxy (Fig. 6).
This is turn requires low CR diffusion coefficients. Such slow outflows would
produce gamma-ray luminosities in dwarf galaxies and compact nuclear
starbursts larger than are observed. This is the fundamental observational
constraint that leads us to favor larger diffusion coefficients and modest
mass-loss rates in dwarf and starburst galaxies. A corollary of this result is
that in all of our models calibrated to explain gamma-ray luminosities well
below the proton-calorimeter value, most of the CR proton energy is vented
into the CGM. Even if the CR-driven mass-loadings on galactic scales are
modest, CRs may play an important ‘preventive’ feedback role on CGM scales
and/or may significantly modify the dynamics and thermodynamics of the CGM (as
was indeed found in the simulations of Ji et al. 2020).

It is instructive to compare our results to related results in the literature.
For example, Booth et al. (2013) assumed $\kappa=3\times 10^{27}$ cm2 s-1 in
their simulations of the impact of cosmic-rays on star-forming galaxies. By
contrast, Salem & Bryan (2014) considered a range of diffusion coefficients
$\kappa=3\times 10^{27}-10^{29}$ cm2 s-1 in a similar study. Neither work
compared to gamma-ray observations. Our results strongly disfavor the low
diffusion coefficient used by Booth et al. (2013) and favor the upper end of
the values modeled in Salem & Bryan (2014). Chan et al. (2019) studied three-
dimensional simulations of idealized galaxies with CRs and other forms of
stellar feedback, and directly compared to gamma-ray observations. They also
concluded that CR diffusion coefficients of $\sim 10^{29}$ cm2 s-1 were
required for consistency with gamma-ray observations. Hopkins et al. (2020)
reached similar conclusions using cosmological zoom-in simulations. Our
analytics help firm up the conclusions drawn from these simulations and show
how they depend on other stellar feedback parameters and the galaxy model
(see, in particular, our analytic scalings in equations 56-65). Both Chan et
al. (2019) and Hopkins et al. (2020) also found, as we do, that while CRs can
drive winds in Milky-way mass galaxies, CRs are not very important wind-
drivers in dwarf galaxies relative to other mechanisms.

A significant difference between our solutions and the cosmological zoom-in
simulations with CRs of Hopkins et al. (2020), Ji et al. (2020), and Hopkins
et al. (2021a) is that we find that advection of CR energy by the gas motion
becomes the dominant CR transport mechanism relatively close to the base of
the wind, with the gas kinetic energy flux taking over at yet larger radii
(see Fig. 4 and eq. 25). By contrast, Hopkins et al. (2020), Ji et al. (2020),
and Hopkins et al. (2021a) argue that diffusion sets up a $p_{c}\propto
r^{-1}$ profile throughout the CGM. A possible resolution of this difference
is that diffusion would likely again be the dominant CR transport mechanism
exterior to a termination shock between a galactic wind and the CGM, which is
not included in our calculations. It is also worth noting that our simulations
require high resolution to resolve the acceleration of the gas at small radii,
particularly for colder phases of the ISM, i.e., larger $V_{g}/c_{i}$ (see
Table LABEL:tab:compare). This is not achievable in cosmological simulations.
If we take our fiducial $\kappa=10$, $V_{g}=10$ simulation (Table
LABEL:tab:compare) and reduce the resolution to $dr/r=0.05$, the mass-loss
rate increases by a factor of $\sim 4$. This is, however, almost certainly
boundary condition dependent, and it is not clear how this result would change
for cosmological simulations which do not have any boundary in the galaxy.

Finally, we stress that our observational calibration of CR diffusion
coefficients using Fermi gamma-ray data is based on a limited sample of
galaxies, primarily those in the local group, M82, NGC 253, and Arp 220
(Ackermann et al., 2012; Griffin et al., 2016). It is thus entirely possible
that there are physical correlations of CR transport with galaxy properties
(gas density, metallicity, galaxy size, …) that are not revealed by the
current data. Despite this caveat, given the particularly large theoretical
uncertainties in the microphysics of CR transport, we believe that
observational calibration of the models is an important constraint, and one
that will hopefully improve in the coming years.

## Data Availability

The numerical simulation results used in this paper will be shared on
reasonable request to the corresponding author.

## Acknowledgments

We thank Andrea Antoni, Phil Hopkins, Philipp Kempski, S. Peng Oh, Eve
Ostriker, and Jono Squire for useful conversations. EQ thanks the Princeton
Astrophysical Sciences department and the theoretical astrophysics group and
Moore Distinguished Scholar program at Caltech for their hospitality and
support. EQ was supported in part by a Simons Investigator Award from the
Simons Foundation and by NSF grant AST-1715070. The Center for Computational
Astrophysics at the Flatiron Institute is supported by the Simons Foundation.
TAT is supported in part by NSF grant #1516967 and NASA grant #80NSSC18K0526.
TAT acknowledges support from a Simons Foundation Fellowship and an IBM
Einstein Fellowship from the Institute for Advanced Study, Princeton, while a
portion of this work was completed.This research made extensive use of
Matplotlib (Hunter, 2007) and Astropy,101010http://www.astropy.org a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013, 2018).

## References

  * Abdo et al. (2010a) Abdo A. A., et al., 2010a, http://dx.doi.org/10.1051/0004-6361/201015759 , https://ui.adsabs.harvard.edu/abs/2010A&A…523L…2A 523, L2
  * Abdo et al. (2010b) Abdo A. A., et al., 2010b, http://dx.doi.org/10.1088/2041-8205/709/2/L152 , https://ui.adsabs.harvard.edu/abs/2010ApJ…709L.152A 709, L152
  * Ackermann et al. (2012) Ackermann M., Ajello M., Allafort A., Baldini L., Ballet J., Bastieri D., et al. 2012, http://dx.doi.org/10.1088/0004-637X/755/2/164 , https://ui.adsabs.harvard.edu/abs/2012ApJ…755..164A 755, 164
  * Amato & Blasi (2018) Amato E., Blasi P., 2018, http://dx.doi.org/10.1016/j.asr.2017.04.019 Advances in Space Research, https://ui.adsabs.harvard.edu/abs/2018AdSpR..62.2731A 62, 2731
  * Andrews et al. (2017) Andrews B. H., Weinberg D. H., Schönrich R., Johnson J. A., 2017, http://dx.doi.org/10.3847/1538-4357/835/2/224 , https://ui.adsabs.harvard.edu/abs/2017ApJ…835..224A 835, 224
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, http://dx.doi.org/10.1051/0004-6361/201322068 , https://ui.adsabs.harvard.edu/abs/2013A&A…558A..33A 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, http://dx.doi.org/10.3847/1538-3881/aabc4f , https://ui.adsabs.harvard.edu/abs/2018AJ….156..123A 156, 123
  * Barcos-Muñoz et al. (2018) Barcos-Muñoz L., et al., 2018, http://dx.doi.org/10.3847/2041-8213/aaa28d , http://adsabs.harvard.edu/abs/2018ApJ…853L..28B 853, L28
  * Bell (2003) Bell E. F., 2003, http://dx.doi.org/10.1086/367829 , https://ui.adsabs.harvard.edu/abs/2003ApJ…586..794B 586, 794
  * Booth et al. (2013) Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, http://dx.doi.org/10.1088/2041-8205/777/1/L16 , https://ui.adsabs.harvard.edu/abs/2013ApJ…777L..16B 777, L16
  * Boulares & Cox (1990) Boulares A., Cox D. P., 1990, http://dx.doi.org/10.1086/169509 , https://ui.adsabs.harvard.edu/abs/1990ApJ…365..544B 365, 544
  * Breitschwerdt et al. (1991) Breitschwerdt D., McKenzie J. F., Voelk H. J., 1991, , http://adsabs.harvard.edu/abs/1991A
  * Buckman et al. (2020) Buckman B. J., Linden T., Thompson T. A., 2020, http://dx.doi.org/10.1093/mnras/staa875 , https://ui.adsabs.harvard.edu/abs/2020MNRAS.494.2679B 494, 2679
  * Chan et al. (2019) Chan T. K., Kereš D., Hopkins P. F., Quataert E., Su K. Y., Hayward C. C., Faucher-Giguère C. A., 2019, http://dx.doi.org/10.1093/mnras/stz1895 , https://ui.adsabs.harvard.edu/abs/2019MNRAS.488.3716C 488, 3716
  * Crocker et al. (2020) Crocker R. M., Krumholz M. R., Thompson T. A., 2020, arXiv e-prints, https://ui.adsabs.harvard.edu/abs/2020arXiv200615821C p. arXiv:2006.15821
  * Crocker et al. (2021) Crocker R. M., Krumholz M. R., Thompson T. A., 2021, http://dx.doi.org/10.1093/mnras/stab148 , https://ui.adsabs.harvard.edu/abs/2021MNRAS.tmp..195C
  * Drury & Falle (1986) Drury L. O., Falle S. A. E. G., 1986, http://dx.doi.org/10.1093/mnras/223.2.353 , https://ui.adsabs.harvard.edu/abs/1986MNRAS.223..353D 223, 353
  * Everett et al. (2008) Everett J. E., Zweibel E. G., Benjamin R. A., McCammon D., Rocks L., Gallagher III J. S., 2008, http://dx.doi.org/10.1086/524766 , http://adsabs.harvard.edu/abs/2008ApJ…674..258E 674, 258
  * Griffin et al. (2016) Griffin R. D., Dai X., Thompson T. A., 2016, http://dx.doi.org/10.3847/2041-8205/823/1/L17 , https://ui.adsabs.harvard.edu/abs/2016ApJ…823L..17G 823, L17
  * Hopkins et al. (2020) Hopkins P. F., et al., 2020, http://dx.doi.org/10.1093/mnras/stz3321 , https://ui.adsabs.harvard.edu/abs/2020MNRAS.492.3465H 492, 3465
  * Hopkins et al. (2021a) Hopkins P. F., Chan T. K., Ji S., Hummels C. B., Kereš D., Quataert E., Faucher-Giguère C.-A., 2021a, http://dx.doi.org/10.1093/mnras/staa3690 , https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.3640H 501, 3640
  * Hopkins et al. (2021b) Hopkins P. F., Squire J., Chan T. K., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021b, http://dx.doi.org/10.1093/mnras/staa3691 , https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.4184H 501, 4184
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Ipavich (1975) Ipavich F. M., 1975, http://dx.doi.org/10.1086/153397 , http://adsabs.harvard.edu/abs/1975ApJ…196..107I 196, 107
  * Ji et al. (2020) Ji S., et al., 2020, http://dx.doi.org/10.1093/mnras/staa1849 , https://ui.adsabs.harvard.edu/abs/2020MNRAS.496.4221J 496, 4221
  * Jiang & Oh (2018) Jiang Y.-F., Oh S. P., 2018, http://dx.doi.org/10.3847/1538-4357/aaa6ce , https://ui.adsabs.harvard.edu/abs/2018ApJ…854….5J 854, 5
  * Krause et al. (2018) Krause M., et al., 2018, http://dx.doi.org/10.1051/0004-6361/201731991 , https://ui.adsabs.harvard.edu/abs/2018A&A…611A..72K 611, A72
  * Kulsrud & Pearce (1969) Kulsrud R., Pearce W. P., 1969, http://dx.doi.org/10.1086/149981 , https://ui.adsabs.harvard.edu/abs/1969ApJ…156..445K 156, 445
  * Lacki & Thompson (2013) Lacki B. C., Thompson T. A., 2013, http://dx.doi.org/10.1088/0004-637X/762/1/29 , https://ui.adsabs.harvard.edu/abs/2013ApJ…762…29L 762, 29
  * Lacki et al. (2010) Lacki B. C., Thompson T. A., Quataert E., 2010, http://dx.doi.org/10.1088/0004-637X/717/1/1 , http://adsabs.harvard.edu/abs/2010ApJ…717….1L 717, 1
  * Lacki et al. (2011) Lacki B. C., Thompson T. A., Quataert E., Loeb A., Waxman E., 2011, http://dx.doi.org/10.1088/0004-637X/734/2/107 , http://adsabs.harvard.edu/abs/2011ApJ…734..107L 734, 107
  * Linden (2017) Linden T., 2017, http://dx.doi.org/10.1103/PhysRevD.96.083001 , https://ui.adsabs.harvard.edu/abs/2017PhRvD..96h3001L 96, 083001
  * Linden et al. (2010) Linden T., Profumo S., Anderson B., 2010, http://dx.doi.org/10.1103/PhysRevD.82.063529 , https://ui.adsabs.harvard.edu/abs/2010PhRvD..82f3529L 82, 063529
  * Lochhaas et al. (2020) Lochhaas C., Thompson T. A., Schneider E. E., 2020, arXiv e-prints, https://ui.adsabs.harvard.edu/abs/2020arXiv201106004L p. arXiv:2011.06004
  * Mao & Ostriker (2018) Mao S. A., Ostriker E. C., 2018, http://dx.doi.org/10.3847/1538-4357/aaa88e , https://ui.adsabs.harvard.edu/abs/2018ApJ…854…89M 854, 89
  * Muratov et al. (2015) Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, http://dx.doi.org/10.1093/mnras/stv2126 , https://ui.adsabs.harvard.edu/abs/2015MNRAS.454.2691M 454, 2691
  * Murray et al. (2005) Murray N., Quataert E., Thompson T. A., 2005, http://dx.doi.org/10.1086/426067 , http://adsabs.harvard.edu/abs/2005ApJ…618..569M 618, 569
  * Ostriker & Shetty (2011) Ostriker E. C., Shetty R., 2011, http://dx.doi.org/10.1088/0004-637X/731/1/41 , https://ui.adsabs.harvard.edu/abs/2011ApJ…731…41O 731, 41
  * Owocki & Gayley (1997) Owocki S. P., Gayley K. G., 1997, ThePhysics of Stellar Winds Near the Eddington Limit. p. 121
  * Owocki et al. (2017) Owocki S. P., Townsend R. H. D., Quataert E., 2017, http://dx.doi.org/10.1093/mnras/stx2251 , https://ui.adsabs.harvard.edu/abs/2017MNRAS.472.3749O 472, 3749
  * Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, http://dx.doi.org/10.3847/2041-8205/824/2/L30 , https://ui.adsabs.harvard.edu/#abs/2016ApJ…824L..30P 824, L30
  * Pavlidou & Fields (2001) Pavlidou V., Fields B. D., 2001, http://dx.doi.org/10.1086/322455 , https://ui.adsabs.harvard.edu/abs/2001ApJ…558…63P 558, 63
  * Pohl (1994) Pohl M., 1994, , https://ui.adsabs.harvard.edu/abs/1994A&A…287..453P 287, 453
  * Ruszkowski et al. (2017) Ruszkowski M., Yang H. Y. K., Zweibel E., 2017, http://dx.doi.org/10.3847/1538-4357/834/2/208 , https://ui.adsabs.harvard.edu/#abs/2017ApJ…834..208R 834, 208
  * Salem & Bryan (2014) Salem M., Bryan G. L., 2014, http://dx.doi.org/10.1093/mnras/stt2121 , https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.3312S 437, 3312
  * Simpson et al. (2016) Simpson C. M., Pakmor R., Marinacci F., Pfrommer C., Springel V., Glover S. C. O., Clark P. C., Smith R. J., 2016, http://dx.doi.org/10.3847/2041-8205/827/2/L29 , https://ui.adsabs.harvard.edu/#abs/2016ApJ…827L..29S 827, L29
  * Skilling (1971) Skilling J., 1971, http://dx.doi.org/10.1086/151210 , https://ui.adsabs.harvard.edu/abs/1971ApJ…170..265S 170, 265
  * Socrates et al. (2008) Socrates A., Davis S. W., Ramirez-Ruiz E., 2008, http://dx.doi.org/10.1086/590046 , http://adsabs.harvard.edu/abs/2008ApJ…687..202S 687, 202
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, http://dx.doi.org/10.1146/annurev-astro-082812-140951 , https://ui.adsabs.harvard.edu/abs/2015ARA&A..53…51S 53, 51
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, http://dx.doi.org/10.3847/1538-4365/ab929b , https://ui.adsabs.harvard.edu/abs/2020ApJS..249….4S 249, 4
  * Strong et al. (2010) Strong A. W., Porter T. A., Digel S. W., Jóhannesson G., Martin P., Moskalenko I. V., Murphy E. J., Orlando E., 2010, http://dx.doi.org/10.1088/2041-8205/722/1/L58 , https://ui.adsabs.harvard.edu/abs/2010ApJ…722L..58S 722, L58
  * Thompson & Lacki (2013) Thompson T. A., Lacki B. C., 2013, The FIR-Radio Correlation in Rapidly Star-Forming Galaxies: The Spectral Index Problem and Proton Calorimetry. p. 283, http://dx.doi.org/10.1007/978-3-642-35410-6˙20 doi:10.1007/978-3-642-35410-6˙20
  * Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, http://dx.doi.org/10.1086/431923 , http://adsabs.harvard.edu/abs/2005ApJ…630..167T 630, 167
  * Thompson et al. (2007) Thompson T. A., Quataert E., Waxman E., 2007, http://dx.doi.org/10.1086/509068 , http://adsabs.harvard.edu/abs/2007ApJ…654..219T 654, 219
  * Torres (2004) Torres D. F., 2004, http://dx.doi.org/10.1086/425415 , https://ui.adsabs.harvard.edu/abs/2004ApJ…617..966T 617, 966
  * Trotta et al. (2011) Trotta R., Jóhannesson G., Moskalenko I. V., Porter T. A., Ruiz de Austri R., Strong A. W., 2011, http://dx.doi.org/10.1088/0004-637X/729/2/106 , https://ui.adsabs.harvard.edu/abs/2011ApJ…729..106T 729, 106
  * Uhlig et al. (2012) Uhlig M., Pfrommer C., Sharma M., Nath B. B., Enßlin T. A., Springel V., 2012, http://dx.doi.org/10.1111/j.1365-2966.2012.21045.x , http://adsabs.harvard.edu/abs/2012MNRAS.423.2374U 423, 2374
  * Veilleux et al. (2020) Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, http://dx.doi.org/10.1007/s00159-019-0121-9 Astronomy and Astrophysics Reviews, https://ui.adsabs.harvard.edu/abs/2020A&ARv..28….2V 28, 2
  * Wiener et al. (2013) Wiener J., Oh S. P., Guo F., 2013, http://dx.doi.org/10.1093/mnras/stt1163 , https://ui.adsabs.harvard.edu/abs/2013MNRAS.434.2209W 434, 2209
  * Wiener et al. (2017) Wiener J., Pfrommer C., Oh S. P., 2017, http://dx.doi.org/10.1093/mnras/stx127 , https://ui.adsabs.harvard.edu/abs/2017MNRAS.467..906W 467, 906