# The Physics of Galactic Winds Driven by Cosmic Rays II: Isothermal Streaming
Solutions

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

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

###### Abstract

We use analytic calculations and time-dependent spherically-symmetric
simulations to study the properties of isothermal galactic winds driven by
cosmic-rays (CRs) streaming at the Alfvén velocity. The simulations produce
time-dependent flows permeated by strong shocks; we identify a new linear
instability of sound waves that sources these shocks. The shocks substantially
modify the wind dynamics, invalidating previous steady state models: the CR
pressure $p_{c}$ has a staircase-like structure with $dp_{c}/dr\simeq 0$ in
most of the volume, and the time-averaged CR energetics are in many cases
better approximated by $p_{c}\propto\rho^{1/2}$, rather than the canonical
$p_{c}\propto\rho^{2/3}$. Accounting for this change in CR energetics, we
analytically derive new expressions for the mass-loss rate, momentum flux,
wind speed, and wind kinetic power in galactic winds driven by CR streaming.
We show that streaming CRs are ineffective at directly driving cold gas out of
galaxies, though CR-driven winds in hotter ISM phases may entrain cool gas.
For the same physical conditions, diffusive CR transport (Paper I) yields
mass-loss rates that are a few-100 times larger than streaming transport, and
asymptotic wind powers that are a factor of $\simeq 4$ larger. We discuss the
implications of our results for galactic wind theory and observations; strong
shocks driven by CR-streaming-induced instabilities produce gas with a wide
range of densities and temperatures, consistent with the multiphase nature of
observed winds. We also quantify the applicability of the isothermal gas
approximation for modeling streaming CRs and highlight the need for
calculations with more realistic thermodynamics.

###### keywords:

Galaxies: Winds – Cosmic Rays

††pubyear: 2021††pagerange: The Physics of Galactic Winds Driven by Cosmic
Rays II: Isothermal Streaming Solutions–B

##  1 Introduction

A significant fraction of the mechanical energy supplied by supernovae to the
interstellar medium (ISM) goes into cosmic-ray (CR) protons (e.g., Blandford &
Eichler 1987). Those cosmic-rays may in turn play a number of important roles
in galaxy formation. Cosmic-rays regulate the ionization state of the dense
interstellar medium (e.g., Dalgarno 2006), can contribute to driving galactic
winds from star-forming galaxies (e.g., Ipavich 1975), and may be an important
source of heating in low-density phases of the interstellar, circumgalactic,
and intracluster medium (e.g., Guo & Oh 2008; Wiener et al. 2013b). Despite
their potential importance in galaxy formation, our understanding of the
impacts of CRs is still relatively rudimentary. The primary theoretical and
observational challenge is that the physical processes regulating CR transport
are still not fully understood (see, e.g., Amato & Blasi 2018 for a review in
the context of the Milky Way). Empirically, we know that CRs have a short
mean-free path and thus do not leave galaxies on a light-crossing time,
despite their relativistic energies. Theoretically, CRs scatter off of small-
scale magnetic fluctuations that can either be the small-scale tail of a
turbulent cascade (e.g., Yan & Lazarian 2002) or fluctuations generated by the
CRs themselves (e.g., Kulsrud & Pearce 1969). If CRs are not efficiently
scattered by ambient turbulence, as is plausibly the case in many physical
conditions, any net drift of the cosmic-rays exceeding the Alfvén speed
$v_{A}$ will excite the gyro-resonant streaming instability (Lerche, 1967);
short-wavelength Alfvén waves can then grow to the point that they can scatter
the cosmic-rays and limit the resulting CR streaming speed (e.g., Bai et al.
2019).

The mechanism of CR transport can have a large effect on their broader
astrophysical impacts. For example, if CRs stream relative to the thermal gas,
they inevitably heat the gas (mediated via the streaming-excited waves) at a
rate $|{\bf v_{A}\cdot\nabla}p_{c}|$, where $p_{c}$ is the CR pressure
(Wentzel 1971). This heating is absent in the case of pure diffusive
transport. Previous numerical work has also shown that the properties of
galactic winds driven by CRs change significantly depending on the mechanism
of CR transport (e.g., Wiener et al. 2017), as does the impact of CRs on the
circumgalactic medium (e.g., Butsky & Quinn 2018; Hopkins et al. 2021).

In this paper and a companion (Quataert et al. 2021; hereafter Paper I) we
study the physical properties of galactic winds driven by CRs and their
dependence on galaxy properties and the mechanism of cosmic-ray transport.
There is a large body of previous analytic and numerical work on galactic
winds driven by CRs with either diffusive or streaming transport (e.g.,
Ipavich 1975; Breitschwerdt et al. 1991; Everett et al. 2008; Booth et al.
2013; Recchia et al. 2016; Chan et al. 2019; see the introduction to Paper I
for a more comprehensive discussion of previous work). Our work provides an
analytic framework for understanding these previous results and provides
estimates of wind properties suitable for use in cosmological simulations or
semi-analytic models of galaxy formation. Paper I considered the case of
cosmic-ray transport by diffusion. In this paper, we consider the case of CR
transport by streaming at the Alfvén speed. We also directly compare these
solutions to their counterparts with diffusion. In both papers we assume that
the gas is isothermal, as is plausible for photoionized gas or when cooling is
rapid. In addition to our analytic estimates, we carry out spherically
symmetric time-dependent simulations of CR-driven winds using the two-moment
CR transport scheme of Jiang & Oh (2018). In the case of diffusive CR
transport in Paper I, the simulations largely validated the analytic
estimates. As we show in this paper, the case of CR transport by streaming is
more interesting: the time-dependent simulations show that steady-state wind
models with streaming are linearly unstable. This invalidates previous steady-
state calculations because the winds are intrinsically time-dependent and have
a time-averaged structure that is different from their steady-state
counterparts. The simulations motivate a revised theory of CR-driven winds
that accounts for the effects of this instability and the non-linear
structures it produces. We also show that the assumption of isothermal gas
breaks down as a result of CR heating in important regimes of parameter space
where earlier isothermal CR streaming wind models have been applied.

The remainder of this paper is organized as follows. In §2, we present
analytic estimates of the mass-loss rate and terminal velocity in galactic
winds driven by CRs streaming at the Alfvén speed using standard theoretical
assumptions. In §3, we present time-dependent numerical simulations, which
show that in most cases the solutions are highly time-dependent, with strong
shocks permeating the flow. These shocks significantly modify the time-
averaged energetics of the CRs relative to standard CR wind models in the
literature. In §4, we develop a modified CR-driven wind model that accounts
for the energetics seen in the simulations. §5 synthesizes our numerical
results compared to the analytics (§5.1), compares the properties of galactic
winds driven by streaming CRs to those driven by diffusing CRs (§5.2; see
Paper I), quantitatively assesses the validity of the isothermal gas
approximation used throughout this work (§5.3), and discusses the
observational implications of our results (§5.4). In §6, we summarize our main
results. Appendix A presents linear stability analyses relevant for
interpreting our numerical simulations. We identify two new (to the best of
our knowledge) instabilities driven by CR streaming, which are the origin of
the time dependence seen in the simulations. Appendix B discusses a few
aspects of the time-averaged CR energetics not addressed in the main text.

##  2 Analytic Approximations for Galactic Winds Driven by Streaming Cosmic
Rays

The equation for the CR energy density $E_{c}$ in the absence of CR sources
and pionic losses, and including diffusion along magnetic field lines and
streaming at the Alfvén velocity down the CR pressure gradient, is given by

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

where $p_{c}=E_{c}/3$ is the CR pressure, ${\bf v_{s}}=-{\bf v_{A}}|\nabla
p_{c}|/\nabla p_{c}$ with ${\bf v_{A}}={\bf B}/(4\pi\rho)^{1/2}$, and 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 2 only when time variations are sufficiently slow (see eq. 42). 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).$ (2)

Here $\kappa$ is the diffusion coefficient and ${\bf n}={\bf v_{\rm A}}/|\bf
v_{\rm A}|$. Equation 1, and indeed a scalar CR pressure in the momentum
equation for the gas, 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.

In this paper we focus primarily on the properties of galactic winds in which
CRs stream at the Alfvén velocity. We also compare those results to the case
of CR diffusion discussed in Paper I. We consider the simplified model problem
of a spherical CR-driven wind in an isothermal gravitational potential. We
further assume an isothermal gas of sound speed $c_{i}$, as is a priori
plausible for rapidly cooling or photoionized gas. In §5.3 we assess the
isothermal gas assumption and show that it breaks down in important regimes of
parameter space because of the inevitable presence of CR heating when CRs
stream at the Alfvén speed. The assumption of spherical symmetry also implies
that the magnetic field is assumed to be a split-monopole configuration, and
only shows up dynamically in setting the Alfvén speed (and hence the streaming
speed) as a function of radius.

In this section we derive steady-state analytic approximations to the wind
solutions by considering hydrostatic equilibrium below the sonic point and
assuming rapid streaming (high $v_{A}$); the latter simplification is
analogous to the assumption of rapid diffusion made in the analytic estimates
of Paper I. In Section 3, we then treat the same problem numerically by
solving the time-dependent CR-driven wind problem using the numerical scheme
of Jiang & Oh (2018). These numerical calculations invalidate the analytic
solutions for many physical parameters, because the numerical solutions are
time-dependent (due to instabilities discussed in §3.4.1 & Appendix A).
Nevertheless, the standard streaming case presented here provides important
pedagogical insights, makes contact with the previous literature (e.g.,
Ipavich 1975; Mao & Ostriker 2018), and sets the stage for §4, in which we
present a modification to the streaming wind solutions that better captures
the dynamics found in our numerical simulations.

###  2.1 The Density Profile and Mass-loss Rate

Under the approximation of time-steady spherical flow, subject to gas pressure
$p$, CR pressure $p_{c}$, and an isothermal gravitational potential
characterized by velocity dispersion $V_{g}$, the equations describing mass
and momentum conservation are

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

Assuming $dp_{c}/dr<0$, the steady-state CR streaming flux is given by

${\bf F_{c}}=4p_{c}({\bf v}+{\bf v_{A}}).$ (5)

Neglecting CR source terms and hadronic losses, the CR energy equation is
given by

$\nabla\cdot\left[4\,\left({\bf v}+{\bf v}_{\rm
A}\right)p_{c}\right]=\left({\bf v}+{\bf v_{\rm A}}\right)\cdot\nabla p_{c}.$
(6)

Noting that in spherical symmetry

$\nabla\cdot{\bf v_{\rm A}}=-\frac{1}{2}v_{\rm A}\frac{d\ln\rho}{dr},$ (7)

the CR energy equation becomes

$\frac{dp_{c}}{dr}=\frac{4}{3}\frac{p_{c}}{\rho}\frac{d\rho}{dr}\left(\frac{v_{\rm
A}/2+v}{v_{\rm A}+v}\right).$ (8)

Combining this expression with equations 3 and 4, we obtain the wind equation
for the fluid in the streaming limit:

$\frac{1}{v}\frac{dv}{dr}=\frac{2}{r}\left(\frac{c_{i}^{2}+c_{\rm
eff}^{2}-V_{g}^{2}}{v^{2}-c_{i}^{2}-c_{\rm
eff}^{2}}\right)\equiv\frac{N_{st}}{D_{st}}$ (9)

where we have defined the numerator $N_{st}$ and denominator $D_{st}$ of the
wind equation, and where

$c_{\rm eff}^{2}=\frac{4}{3}\frac{p_{c}}{\rho}\left(\frac{v_{\rm
A}/2+v}{v_{\rm A}+v}\right),$ (10)

is the effective CR sound speed. We will use several different definitions of
the CR sound speed in this paper, depending on the exact context; see Table
LABEL:tab:cs.

Equation (9) closely follows the relations derived in Ipavich (1975). Setting
the numerator and denominator in the wind equation equal to zero at the sonic
point $r_{s}$ implies that

$v(r_{s})=V_{g}\ \ \ {\rm and\ \ \ }c_{\rm
eff}^{2}(r_{s})=V_{g}^{2}-c_{i}^{2}\equiv{\rm V^{2}_{g,eff}}$ (11)

where we have defined an effective gravitational speed ${\rm V_{g,eff}}$ for
future use.

Table 1: Summary of different definitions of cosmic-ray sound speed used in
this paper. In Paper I we exclusively used $c_{c}$. The value of the CR sound
speed at the base of the wind appears in many of our analytic estimates, and
is denoted $c_{\rm eff,0}$ or $c_{c,0}$.

Quantity Definition Context Use in Paper General $c_{\rm eff}$ eq. 10 full
streaming wind equation §2 $c_{\rm eff}$ for ${\rm\gamma_{eff}}=2/3$
$\sqrt{\frac{2}{3}\frac{p_{c}}{\rho}}$ (eq. 14) high $v_{A}$ limit of standard
streaming analytics §2 Analytics $c_{\rm eff}$ for general ${\rm\gamma_{eff}}$
$\sqrt{{\rm\gamma_{eff}}\frac{p_{c}}{\rho}}$ (eq. 47) high $v_{A}$ limit of
modified streaming analytics §4.1 Analytics $c_{\rm eff}$ with
${\rm\gamma_{eff}}=1/2$ $\sqrt{\frac{1}{2}\frac{p_{c}}{\rho}}$ (eq. 47 with
${\rm\gamma_{eff}}=1/2$) high $v_{A}$ limit of modified streaming analytics
§4.2 Analytics $c_{\rm eff}$ with ${\rm\gamma_{eff}}=4/3$
$\sqrt{\frac{4}{3}\frac{p_{c}}{\rho}}$ (eq. 78) adiabatic CR sound speed §A.1
Stability Calculation $c_{c}$ $\sqrt{\frac{p_{c}}{\rho}}$ ‘Isothermal’ CR
Sound Speed §3 Simulations and §A.2 Stability Calculation

We now seek an analytic expression for the mass loss rate in CR-driven winds
in the streaming limit for comparison with the more complete numerical
calculations presented in Section 3. To do so, we employ the commonly-used
strategy of assuming that the system can be approximated as in hydrostatic
equilibrium with $v\rightarrow 0$ in the region below the sonic point
$r<r_{s}$. This approximation is justified deep in the atmosphere of the wind,
but becomes increasingly suspect as the sonic point is approached.
Nevertheless, as in discussions of the isothermal and polytropic Parker-type
thermal winds, this approximation yields the scalings of how the mass-loss
rate depends on physical parameters and an estimate for the normalization of
the mass loss rate $\dot{M}_{w}$.

To make analytic progress we assume that in the hydrostatic portion of the
solution, $v\ll v_{A}$. We check below a posteriori when this assumption is
valid (eq. 28) and note that our numerical simulations in §3 do not employ the
high $v_{A}$ assumption. Taking the limit $v\ll v_{\rm A}$, equation (8)
becomes

$\frac{d\ln
p_{c}}{dr}=\frac{2}{3}\frac{d\ln\rho}{dr}\Longrightarrow\left(\frac{p_{c}}{p_{c,0}}\right)=\left(\frac{\rho}{\rho_{0}}\right)^{2/3},$
(12)

where $p_{c,0}$ and $\rho_{0}$ represent parameters of the problem at the base
of the outflow, or at any other reference point. Note that equation (12) can
also be written as

$c_{\rm eff}^{2}=c_{\rm
eff,0}^{2}\left(\frac{\rho_{0}}{\rho}\right)^{1/3}=\frac{2}{3}\frac{p_{c,0}}{\rho_{0}}\left(\frac{\rho_{0}}{\rho}\right)^{1/3},$
(13)

where $c_{\rm eff}^{2}=(2/3)p_{c}/\rho$ in the limit that $v\ll v_{\rm A}$.
These expressions show that the CR sound speed increases as the gas density
decreases and the wind accelerates, allowing the flow to eventually reach the
sonic conditions given in equation 11. If, however, $v\gg v_{A}$, $c_{\rm
eff}^{2}\simeq(4/3)p_{c}/\rho$ and equation 8 becomes
$p_{c}\propto\rho^{4/3}$, i.e., the CR are adiabatic. In this case $c_{\rm
eff}^{2}\propto\rho^{1/3}$ and the CR sound speed decreases with decreasing
density. This shows that there cannot be an isothermal CR-driven outflow in
the limit of small $v_{A}$, because the CR sound speed can never reach $V_{g}$
if $c_{\rm eff,0}\lesssim V_{g}$, i.e., if the CR are initially bound. The
condition $v\lesssim v_{A}$ must thus apply over a significant range of radii
in order for the outflow to accelerate and reach the escape speed. Equations
12 & 13 are an approximation to this requirement that enables considerable
analytic progress. We will use these for the remainder of this section.

We now calculate the density profile in the hydrostatic portion of the flow
assuming $p_{c}\propto\rho^{2/3}$ (eq. 12), via

$\frac{c_{i}^{2}}{\rho}\frac{d\rho}{dr}+\frac{c_{\rm
eff}^{2}}{\rho}\frac{d\rho}{dr}=\frac{c_{i}^{2}}{\rho}\frac{d\rho}{dr}+\frac{c_{\rm
eff,0}^{2}}{\rho}\left(\frac{\rho_{0}}{\rho}\right)^{1/3}\frac{d\rho}{dr}=-\frac{2V_{g}^{2}}{r},$
(14)

which has the solution

$c_{i}^{2}\ln\left(\frac{\rho}{\rho_{0}}\right)-3c_{\rm
eff,0}^{2}\left(\frac{\rho_{0}^{1/3}}{\rho^{1/3}}-1\right)=-2V_{g}^{2}\ln\left(\frac{r}{r_{0}}\right).$
(15)

This expression yields a simple implicit relation for the density profile that
is easily solved numerically for $\rho(r)$. Using the numerical solutions to
equation 15 we can also estimate the mass-loss rate in galactic winds driven
by streaming CRs in the high $v_{A}$ limit as follows: using equation 13,
equation 11 directly implies the density at the sonic point, namely

$\rho(r_{s})\simeq\rho_{0}\left(\frac{c_{\rm eff,0}}{{\rm
V_{g,eff}}}\right)^{6}.$ (16)

The sonic radius $r_{s}$ is then determined numerically using the solution to
equation 15 as the radius where $\rho(r_{s})$ is reached; we derive analytic
approximations for $r_{s}$ below. The mass-loss rate is then $\dot{M}_{\rm
w}\simeq 4\pi r_{s}^{2}\rho({r_{s}})V_{g}$ (eq. 11).

Figure 1 shows the mass-loss rate as a function of the two key dimensionless
parameters in the problem, namely the strength of gravity $(V_{g}/c_{i})$ and
the base CR pressure ($p_{c,0}/\rho_{0}c_{i}^{2}$). As in Paper I, we express
the mass-loss rate 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)$ (17)

where $n_{0}\equiv\rho_{0}/m_{p}$. Figure 1 shows that the mass-loss rate due
to streaming CRs is a particularly strong function of the strength of gravity
relative to the sound speed in the galactic disk; as we show explicitly below,
the dominant dependence is actually on $V_{g}/c_{\rm eff,0}$, i.e., on the CR
sound speed rather than the gas sound speed. It is useful to also think of
different values of $c_{i}$ and $c_{\rm eff,0}$ in Figure 1 as corresponding
to different phases of the ISM. Figure 1 implies that most of the mass loss
will be from the warmer phases of the ISM, which have larger $c_{i}$ and
$c_{\rm eff,0}$, unless $p_{c,0}$ is much lower in those phases.

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

We can also find analytic approximations to the numerical solutions in Figure
1. We first consider the limit where $c_{i}^{2}\ll 3c_{\rm eff,0}^{2}$. The
analytic solution to equation (15) is then

$\frac{\rho}{\rho_{0}}\simeq\left[1+\frac{2}{3}\frac{V_{g}^{2}}{c_{\rm
eff,0}^{2}}\,\ln\left(\frac{r}{r_{0}}\right)\right]^{\,-3}.$ (18)

Equation 18 is in fact a good approximation to the numerical solution of 15
even for the case of $c_{\rm eff,0}=c_{i}$. Given equation 18, the mass loss
rate can be estimated using equation (3), but with all quantities evaluated at
the sonic point. Comparing equation 16 and 18 we derive an approximation for
the sonic point radius

$\frac{r_{s}}{r_{0}}\simeq\exp\left[\frac{3}{2}\frac{c_{\rm
eff,0}^{2}}{V_{g}^{2}}\left(\frac{{\rm V^{2}_{g,eff}}}{c_{\rm
eff,0}^{2}}-1\right)\right].$ (19)

The general expression for the wind mass loss rate in the streaming limit,
assuming a hydrostatic density profile with $c_{i}^{2}\ll 3c_{\rm eff,0}^{2}$,
is then

$\dot{M}_{w}\simeq 4\pi r_{0}^{2}\,\rho_{0}\,V_{g}\,\left(\frac{\,c_{\rm
eff,0}}{{\rm V_{g,eff}}}\right)^{6}\,\exp\left[\frac{3c_{\rm
eff,0}^{2}}{V_{g}^{2}}\left(\frac{{\rm V^{2}_{g,eff}}}{c_{\rm
eff,0}^{2}}-1\right)\right].$ (20)

This expression for $\dot{M}_{w}$ is independent of the gas sound speed
$c_{i}$. The apparent dependence on $c_{i}$ in Figure 1 is primarily because
$c_{i}$ is used to non-dimensionalize $\dot{M}_{w}$, $c_{\rm eff,0}$ and
$V_{g}$; there is also a weak dependence on $c_{i}$ at low $c_{\rm
eff,0}/c_{i}$ derived below (eq. 27). We retain this choice of dimensionless
variables for consistency with Paper I.

For case of massive galaxies like the Milky Way, where one expects
$V_{g}^{2}\gg c_{\rm eff,0}^{2}$ and $V_{g}^{2}\gg c_{i}^{2}$, the expression
for the sonic point radius reduces to the remarkably simple form

$\frac{r_{s}}{r_{0}}=e^{3/2},$ (21)

and equation (20) becomes

$\begin{split}\dot{M}_{w}=&\,4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{3}\,\left(\frac{\,c_{\rm
eff,0}}{V_{g}}\right)^{6}\simeq 6\times 10^{-4}\,{\rm
M_{\odot}\,yr^{-1}}\left(\frac{r_{0}}{1\,{\rm kpc}}\right)^{2}\\\
&\times\,\left(\frac{n_{0}}{1\,{\rm cm^{-3}}}\right)\left(\frac{c_{\rm
eff,0}}{10\,\,{\rm km\,\,s}^{-1}}\right)^{6}\left(\frac{V_{g}}{100\,\,{\rm
km\,\,s}^{-1}}\right)^{-5}\end{split}$ (22)

For reference, if we scale for parameters appropriate to a galaxy like the
Milky Way with mass-averaged ISM conditions $V_{g}\simeq 150$ km s-1, $c_{\rm
eff,0}\simeq 10$ km s-1, $n_{0}\simeq 1$ cm-3, and $r_{0}\sim 5$ kpc, this
expression yields $\dot{M}_{w}\simeq 2\times 10^{-3}$ M⊙ yr-1, much less than
the star formation rate, with the most important factor being the strong
suppression $(c_{\rm eff,0}/V_{g})^{6}\simeq 10^{-7}$. However, the mass-loss
rate in equation 22 can also be expressed as $\dot{M}_{w}\propto p_{c,0}c_{\rm
eff,0}^{4}/V_{g}^{5}$. Assuming that the hot ISM has a similar CR pressure,
its mass-loss rate will be much larger than that of the cold ISM, with, e.g.,
$\dot{M}_{w}\simeq 1\,{\rm M_{\odot}\,yr^{-1}}$ for $c_{\rm eff,0}\simeq
50\,{\rm km\,\,s}^{-1}$, comparable to the star formation rate. Most of the
mass loss is thus likely to originate from the hotter phases of the ISM. Cold
gas can still be unbound by such an outflow if it is entrained in the hot CR
wind. Equation 22 and Figure 1 show, however, that direct CR driving of winds
from $\sim 10^{4}$ K gas is very inefficient.

An instructive expression for the wind mass-loss rate can be obtained by
comparing the mass-loss rate estimated here to the star formation rate. To do
so, we note that feedback-regulated models of star formation in galaxies
predict that (e.g., Thompson et al. 2005; Ostriker & Shetty 2011)
$\dot{M}_{*}\approx\pi r_{0}^{2}\dot{\Sigma}_{*}$ with $\dot{\Sigma}_{*}\simeq
2\sqrt{2}\pi G\Sigma_{g}^{2}\phi/(p_{*}/m_{*})$ where $p_{*}/m_{*}\approx
3000\,{\rm km\,\,s}^{-1}$ is the momentum per unit mass associated with
stellar feedback, which supports the disk against its own self-gravity, and
$\phi$ quantifies the stellar and dark matter contribution to the gravity of
the disk. Equation 22 for the mass-loss rate in the streaming limit takes the
form

$\begin{split}\frac{\dot{M}_{\rm
w}}{\dot{M}_{*}}\simeq&\,0.06\left(\frac{c_{\rm eff,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)^{4}\left(\frac{p_{*}/m_{*}}{3000\,\,{\rm
km\,\,s}^{-1}}\right)\\\ &\times\left(\frac{V_{g}}{100\,\,{\rm
km\,\,s}^{-1}}\right)^{-5}\frac{p_{c,0}}{\pi G\Sigma_{g}^{2}\phi}\end{split}$
(23)

For Milky-way like conditions in which the CR pressure is comparable to that
needed for hydrostatic equilibrium in the galactic disk, the wind mass-loss
rate due to streaming CRs from the ‘mass-average’ ISM (with $c_{\rm
eff,0}\simeq 10\,{\rm km\,\,s}^{-1}$) is much less than the star formation
rate. As noted previously, however, volume filling lower density gas with
larger $c_{\rm eff,0}$ will have a significantly higher mass-loss rate,
plausibly comparable to the star formation rate. For example, taking
$V_{g}=150$ km s-1 with $c_{\rm eff,\,0}=30$ km s-1 and $50$ km s-1 gives
$\dot{M}_{w}/\dot{M}_{*}\simeq 0.2$ and $\simeq 2$, respectively, if
$p_{c,0}\simeq 0.3\pi G\Sigma_{g}^{2}\phi$ in the warmer phases of the ISM.

We now consider analytic approximations for the mass-loss rate in the limit of
weak CR pressure compared to gas pressure at the base in the galactic disk,
i.e., $c_{\rm eff,0}\ll c_{i}$. In this limit the gas density profile is
initially set by gas pressure, with

$\frac{\rho}{\rho_{0}}=\left(\frac{r_{0}}{r}\right)^{2V_{g}^{2}/c_{i}^{2}}\ \
\ \ (r<r_{tr})$ (24)

However, as the density drops, the CR pressure increases in importance
relative to the gas pressure. So long as $V_{g}\gtrsim c_{i}$,222This is the
only interesting regime for a wind solution since otherwise the gas is
effectively unbound in the galactic disk even without CRs. the sonic point
condition (eq. 11) requires that the pressure be CR dominated at the sonic
point. The transition between gas pressure and CR pressure support happens at
a radius

$r_{tr}\simeq
r_{0}\left(\frac{\rho_{0}c^{2}_{i}}{p_{c,0}}\right)^{\frac{3c_{i}^{2}}{2V_{g}^{2}}}.$
(25)

The density profile exterior to this transition radius is then like equation
18 but with a different boundary condition set by continuity at $r_{tr}$. This
yields

$\frac{\rho}{\rho_{0}}\simeq\left(\frac{p_{c,0}}{\rho_{0}c^{2}_{i}}\right)^{3}\left[1+\frac{V_{g}^{2}}{c_{i}^{2}}\,\ln\left(\frac{r}{r_{tr}}\right)\right]^{\,-3}\
\ \ \ (r>r_{tr})$ (26)

The sonic point is then located at $r_{s}\simeq e^{3/2}r_{tr}$ and the mass-
loss rate is

$\dot{M}_{w}\simeq 4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{3}\,\left(\frac{\rho_{0}c^{2}_{i}}{p_{c,0}}\right)^{\frac{3\,c_{i}^{2}}{V_{g}^{2}}}\left(\frac{c_{\rm
eff,0}}{V_{g}}\right)^{6}.$ (27)

Note that equation 27 (valid for $c_{\rm eff,0}<c_{i}$) differs from equation
22 (valid for $c_{\rm eff,0}>c_{i}$) only by a factor of
$(\rho_{0}c^{2}_{i}/p_{c,0})^{3c_{i}^{2}/V_{g}^{2}}$, which is not that
different from one for massive galaxies with $V_{g}\gg c_{i}$; the latter is
the condition that $r_{tr}\simeq r_{0}$ (eq. 25). Physically, the sonic point
conditions uniquely determine the density at the sonic point (eq. 16). So long
as $r_{tr}\simeq r_{0}$, then the radius of the sonic point can be estimated
entirely neglecting the gas-dominated region near the base of the wind, and is
still given by equation 21.

###  2.2 Validity of the High Alfvén Velocity Approximation

The analytic estimate of the mass-loss rate in Figure 1 assumes that $v_{A}\gg
v$. Since the mass-loss rate is set by conditions at the sonic point, our
estimate is applicable when $v_{A}(r_{s})\gtrsim v(r_{s})=V_{g}$ so that
equation 12 is valid out to the sonic point. Using equations 11 and 13 this
constraint can be rewritten as

$v_{A,0}\gtrsim v_{A,crit}\equiv{c_{\rm eff,0}}\,\frac{V_{g}\,c^{2}_{\rm
eff,0}}{V^{3}_{\rm
g,eff}}\left(\frac{r_{s}}{r_{0}}\right)^{2}\simeq\frac{e^{3}c_{\rm
eff,0}^{3}}{V_{g}^{2}}\ \ {\rm High\ v_{A}\ Limit}$ (28)

where the last expression is appropriate for massive galaxies with
$V_{g}^{2}\gg c_{\rm eff,0}^{2},c_{i}^{2}$. For sufficiently small $v_{A,0}\ll
v_{A,crit}$ there cannot be an unbound wind because the CRs are effectively
adiabatic (see the discussion after eq. 13). A plausible conjecture is in fact
that if $v_{A,0}\lesssim v_{A,crit}$ the solution will not be able to pass
through the sonic point and will thus formally be a ‘breeze’ rather than a
wind. We shall see that this is borne out by our numerical solutions in §3.6.

Comparing equation 23 and equation 28 we see that for a given base Alfvén
speed $v_{A,0}$, a larger base CR sound speed $c_{\rm eff,0}$ will both
increase the mass-loss rate and make the high $v_{A}$ approximation used to
derive equation 23 suspect. One might thus anticipate that the largest mass-
loss rate for a given set of base conditions would be obtained for conditions
that just satisfy $v_{A}\sim v_{A,crit}$. We shall see below in §2.4 that this
is indeed correct (eq. 38).

###  2.3 Energetics, Momentum Flux, & Terminal Velocity

There is no conserved total energy flux for our model problem of isothermal
gas with streaming CRs. The reason is that the gas is heated by the CRs (the
standard streaming term $\sim v_{A}dp_{c}/dr$) but this energy is assumed to
be instantaneously radiated away to keep the gas isothermal. Instead, the
steady state CR energy equation takes the form

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}F_{c}\right)=\left(v+v_{A}\right)\frac{dp_{c}}{dr}$
(29)

An effective energy equation for the gas can be derived using the gas momentum
equation, which can be rewritten as

$\frac{d}{dr}\left(\frac{1}{2}v^{2}+c_{i}^{2}\ln\rho+2V_{g}^{2}\ln
r\right)=-\frac{1}{\rho}\frac{dp_{c}}{dr}$ (30)

Multiplying by $r^{2}v$ and combining with equation 29 yields

$\frac{1}{r^{2}}\,\frac{d}{dr}\left(\dot{M}\left[\frac{1}{2}v^{2}+c_{i}^{2}\ln\rho+2V_{g}^{2}\ln
r\right]+\dot{E}_{c}\right)=4\pi v_{A}\frac{dp_{c}}{dr}$ (31)

where $\dot{E}_{c}=4\pi r^{2}F_{c}$ is the CR ‘luminosity/power.’ Because
$dp_{c}/dr<0$, the right-hand-side of equation 31 is negative and so the total
power carried by the wind decreases with increasing radius.

To the extent that $p_{c}\propto\rho^{2/3}$ (eq. 8), the CR pressure gradient
term in equation 30 can be rewritten as a CR enthalpy, leading to a conserved
Bernoulli-like constant (e.g., Mao & Ostriker 2018)

$\frac{d}{dr}\left(\frac{1}{2}v^{2}+c_{i}^{2}\ln\rho+2V_{g}^{2}\ln r-3c_{\rm
eff,0}^{2}\left[\frac{\rho}{\rho_{0}}\right]^{-1/3}\right)=0.$ (32)

We can roughly estimate the terminal speed of CR-driven winds in the streaming
limit using equation 32 as follows. The velocity of the gas will increase with
increasing radius so long as $v_{A}>v$. This is because the CR enthalpy $\sim
c_{\rm eff}^{2}\propto\rho^{-1/3}$ increases to arbitrarily large values with
decreasing density so long as $p_{c}\propto\rho^{2/3}$. The acceleration of
the flow ceases, however, when $v_{A}<v$ because then the gas becomes
adiabatic and $p_{c}\propto\rho^{4/3}$, i.e., $c_{\rm
eff}^{2}\propto\rho^{1/3}$ (eq. 10), which decreases with decreasing density.
There is thus a key radius in the flow where $v\sim v_{A}$ which determines
where the acceleration ceases. We call this the Alfvén radius333The Alfvén
radius in magnetocentrifugal winds is conceptually different from that defined
here. $r_{A}$ and estimate the terminal speed using $v_{\infty}\simeq
v_{A}(r_{A})$. We can estimate $r_{A}$ and $v_{\infty}$ by equating the
Bernoulli-like constant in equation 32 at the sonic point $r_{s}$ and $r_{A}$
and using equations 11. In doing so we neglect the gravity term which makes
only a small logarithmic correction to the result (and for any potential with
a finite escape speed the correction is even less important). This yields

$v_{\infty}^{2}\simeq V_{g}^{2}+6{\rm
V^{2}_{g,eff}}\left(\left[\frac{\rho(r_{A})}{\rho(r_{s})}\right]^{-1/3}-1\right)$
(33)

Given $\rho(r_{s})$ from the sonic conditions (eq. 16), equation 33 has 3
unknowns: $v_{\infty}$, $r_{A}$, $\rho(r_{A})$. The two additional equations
are the definition $v_{A}(r_{A})=v_{\infty}$, which yields
$\rho(r_{A})=B_{0}^{2}/(4\pi v_{\infty}^{2})\,(r_{0}/r_{A})^{4}$ and
conservation of mass ($\dot{M}$ = const) between $r_{s}$ and $r_{A}$, which
yields $\rho(r_{A})r_{A}^{2}v_{\infty}\simeq r_{s}^{2}\rho(r_{s})V_{g}$. These
two equations combine to yield

$\frac{\rho(r_{A})}{\rho(r_{s})}=\left(\frac{r_{s}}{r_{0}}\right)^{4}\left(\frac{V_{g}\,c_{\rm
eff,0}}{{\rm V_{g,eff}}\,v_{A,0}}\right)^{2}\left(\frac{c_{\rm eff,0}}{{\rm
V_{g,eff}}}\right)^{4}=\left(\frac{v_{A,crit}}{v_{A,0}}\right)^{2}$ (34)

where the second equality uses equation 28. Recall that $v_{A,0}>v_{A,crit}$
is required for the validity of our high $v_{A}$ approximation, in which case
$v_{\infty}\gtrsim V_{g}$ from equation 33. Equation 34 can be substituted
into equation 33 to estimate $v_{\infty}$. The expression is particularly
simple in the high $v_{A}$ limit of massive galaxies with ${\rm
V_{g,eff}}\simeq V_{g}\gg c_{\rm eff,0}\gtrsim c_{i}$ and $r_{s}\simeq
e^{3/2}r_{0}$ (eq. 21):

$\begin{split}v_{\infty}\simeq&\
\sqrt{6}\frac{V_{g}}{e}\left(\frac{V_{A,0}}{c_{\rm
eff,0}}\right)^{1/3}\left(\frac{V_{g}}{c_{\rm eff,0}}\right)^{2/3}\simeq
400\,{\rm km\,\,s}^{-1}\,\times\\\ &\,\left(\frac{V_{g}}{100\,{\rm
km\,\,s}^{-1}}\right)^{5/3}\left(\frac{v_{A,0}}{10\,{\rm
km\,\,s}^{-1}}\right)^{1/3}\left(\frac{c_{\rm eff,0}}{10\,{\rm
km\,\,s}^{-1}}\right)^{-1}\end{split}$ (35)

where we again note that equation 35 requires that equation 28 is satisfied.
When $v_{A,0}\lesssim v_{A,crit}$, either there is no wind (for very small
$v_{A,0}$) or $v_{\infty}\lesssim V_{g}$. As discussed in the context of
equations 22 and 23, the largest mass-loss rates generally arise in phases of
the ISM with larger $c_{\rm eff,0}$. Equation 35 shows that a consequence of
this larger mass-loading is that the velocity of the outflow is significantly
smaller at larger $c_{\rm eff,0}$, probably at most $\sim V_{g}$.

For the case of massive galaxies with $v_{A,0}>v_{A,crit}$, equations 22 and
35 can be combined to yield the total energy carried by the wind. We express
this wind ‘luminosity’ in terms of the CR ‘luminosity’ at the base of the
wind, $\dot{E}_{c}(r_{0})=16\pi r_{0}^{2}p_{c,0}v_{A,0}$, yielding

$\frac{0.5\dot{M}_{w}v_{\infty}^{2}}{\dot{E}_{c}(r_{0})}\simeq\frac{e}{2}\frac{c_{\rm
eff,0}^{2}}{v_{A,0}^{1/3}V_{g}^{5/3}}\simeq\frac{1}{2}\left(\frac{v_{A,crit}}{v_{A,0}}\right)^{1/3}\left(\frac{c_{\rm
eff,0}}{V_{g}}\right)\ll 1$ (36)

where the second equality uses the definition of $v_{A,crit}$ valid for
massive galaxies. The final inequality in equation 36 follows from requiring
$c_{\rm eff,0}<V_{g}$ and $v_{A,0}>v_{A,crit}$ (the latter for the validity of
our analytics) and implies that in the streaming limit, the total energy
carried by the gas at large radii is less than that supplied to CRs at the
base of the wind. This is because for an isothermal gas equation of state,
streaming losses $\sim v_{A}dp_{c}/dr$ remove energy from the CRs and are
assumed to be rapidly radiated away by the gas.

Finally, we can derive an expression for the asymptotic momentum flux in the
wind $\dot{p}_{w}=\dot{M}_{w}v_{\infty}$. This quantity can be compared with
the total momentum rate carried by photons from star formation
$\dot{p}_{ph}=L_{\rm bol}/c=\epsilon_{ph}\dot{M}_{*}c$, where
$\epsilon_{ph,\,-3.3}=\epsilon_{ph}/5\times 10^{-4}$ for steady-state star
formation and a standard IMF. Using equations 23 and 35, we then have that

$\begin{split}\frac{\dot{p}_{w}}{\dot{p}_{ph}}&\simeq\frac{\dot{M}_{w}}{\dot{M}_{*}}\frac{v_{\infty}}{\epsilon_{ph}c}\\\
&\simeq 0.16\,\epsilon_{ph,\,-3.3}^{-1}\,\left(\frac{c_{\rm eff,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)^{3}\left(\frac{p_{*}/m_{*}}{3000\,\,{\rm
km\,\,s}^{-1}}\right)\\\ &\times\left(\frac{V_{g}}{100\,\,{\rm
km\,\,s}^{-1}}\right)^{-10/3}\left(\frac{v_{A,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)^{1/3}\frac{p_{c,0}}{\pi G\Sigma_{g}^{2}}\end{split}$
(37)

Equation 37 shows that the momentum flux in isothermal CR-driven winds is in
general less than or at most comparable to that in the stellar radiation
field. For comparison, IMF-averaged SNe-ejecta have momentum fluxes $\sim L/c$
and the momentum can increase by a factor of $\sim 10$ or more in the Sedov-
Taylor phase (Ostriker & Shetty 2011).

###  2.4 Maximum Mass-Loss Rate

There is an upper limit to the mass-loss rate associated with CR-driven winds
that is set by energy conservation. This upper limit can be rigorously derived
for diffusive CR transport because in that regime there is a conserved wind
power that is independent of radius (eq. 34 of Paper I). This is not the case
for isothermal winds with CR streaming because streaming losses remove energy
from the cosmic-rays (and that energy is assumed to be radiated away). It is
nonetheless instructive to provide a rough estimate of the maximum mass-loss
rate for a given set of base conditions.

The cosmic-ray power supplied at the base of the wind is
$\dot{E}_{c}(r_{0})=16\pi r_{0}^{2}p_{c,0}v_{A,0}$. Assume that a fraction
$\zeta$ of this energy goes into lifting gas out to larger radii (the
remaining energy is lost radiatively). The maximum mass-loss rate is realized
when the asymptotic kinetic energy vanishes, so that $\dot{M}_{\rm max}\simeq
2\zeta\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.444$v_{\rm esc}(r_{0})$ is not
defined for a $\ln(r)$ potential, 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).

Equation 36 provides an estimate of the energy in the wind at large radii,
i.e., $\zeta$. However, this result only holds when the asymptotic kinetic
energy is finite, which is not the case as
$\dot{M}_{w}\rightarrow\dot{M}_{max}$. We do not have an analogous
approximation for $\zeta$ valid as $\dot{M}_{w}\rightarrow\dot{M}_{max}$. To
account for this we use equation 36 but multiply the result by a dimensionless
factor $\zeta_{0}$. Doing so, it is straightforward to combine equation 22,
28, & 36 to find

$\frac{\dot{M}_{w}}{\dot{M}_{\rm
max}}\simeq\zeta_{0}^{-1}\,\left(\frac{v_{A,crit}}{v_{A,0}}\right)^{2/3}$ (38)

Equation 38 shows that in the high base Alfvén speed limit, the predicted
mass-loss rate is below $\dot{M}_{max}$, consistent with the finite terminal
speed predicted by equation 35. For lower base Alfvén speeds our analytic
estimates are less applicable, but extrapolation of our results suggests that
$\dot{M}_{w}\simeq\dot{M}_{\rm max}$ will be realized for $v_{A,0}\lesssim
v_{A,crit}$. We shall see in §3 that this extrapolation is borne out: our
numerical solutions for low $v_{A,0}$ produce outflows but these outflows have
$\dot{M}_{w}\simeq\dot{M}_{\rm max}$ and never reach the critical point at
which $v=V_{g}$ (Fig. 9); they are formally ‘breezes’ rather than transonic
winds.

We can also express the maximum mass-loss rate $\dot{M}_{\rm max}$ in a form
that is easier to compare to observational constraints. To do so, we write the
CR energy injection rate at the base of the outflow as

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

where $\dot{M}_{*}$ is the star formation rate and $\epsilon_{c}\equiv
10^{-6.3}\epsilon_{c,-6.3}$ is set by 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. Our expression for the maximum mass-loss rate can thus be written as

$\frac{\dot{M}_{\rm
max}}{\dot{M}_{*}}\simeq\,0.2\,\zeta_{0}\,\epsilon_{c,-6.3}\,\left(\frac{100\,\,{\rm
km\,\,s}^{-1}}{V_{g}}\right)^{7/3}\,\left(\frac{v_{A,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)^{1/3}$ (40)

where we have assumed $v^{2}_{\rm esc}(r_{0})\simeq 4V_{g}^{2}$ and have
maximized the CR energy flux by taking $v_{A,0}\simeq v_{A,crit}$ and hence
$c_{\rm eff,0}/V_{g}\simeq e^{-1}(v_{A,0}/V_{g})^{1/3}$ (eq. 28) and
$\zeta\simeq 0.18\zeta_{0}\,(v_{A,0}/V_{g})^{1/3}$ (eq. 36).

Table 2: Parameters and properties of our numerical simulations. Columns are
diffusion coefficient $\kappa$, base Alfvén speed $v_{A,0}$, 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}$, mass-loss rate in
units of $\dot{M}_{0}=4\pi r_{0}^{2}\rho_{0}c_{i}$ (eq. 17), velocity of the
gas at the outer radius $v_{(}r_{out})$, and ratio of the wind power at the
outer radius to the CR power at the base
$\dot{E}_{w}(r_{out})/\dot{E}_{c}(r_{0})$. Solutions labeled as ’Breeze’ are
not transonic and never reach velocities comparable to $V_{g}$. Numerical
resolution is $\delta r/r=5.25\times 10^{-4}$ unless noted otherwise.

Transport | $\kappa$ | $v_{A,0}$ | $V_{g}$ | $p_{c,0}$ | $r_{out}$ | $V_{m}$ | $\dot{M}_{sim}$ | $v(r_{out})$ | $\frac{\dot{E}_{w}(r_{out})}{\dot{E}_{c}(r_{0})}$ | Notes
---|---|---|---|---|---|---|---|---|---|---
| ($r_{0}c_{i}$) | ($c_{i}$) | ($c_{i}$) | ($\rho_{0}c_{i}^{2}$) | ($r_{0}$) | ($c_{i}$) | ($\dot{M}_{0}$) | ($V_{g}$) | – |
Streaming |  |  |  |  |  |  |  |  |  |
| – | 1 | 10 | 1 | 10 | 3000 | $5\times 10^{-4}$ | 1.9 | 0.16 |
| – | 1 | 10 | 1 | 10 | $3\times 10^{4}$ | $4.6\times 10^{-4}$ | 2.2 | 0.15 | High $V_{m}$
| – | 3 | 10 | 1 | 5 | 3000 | $9\times 10^{-4}$ | 2.2 | 0.17 |
| – | 10 | 10 | 1 | 5 | 6000 | $1.2\times 10^{-3}$ | 2.8 | 0.16 |
| – | 10 | 10 | 1 | 5 | 3000 | $2.1\times 10^{-3}$ | 2.9 | 0.19 | $\delta r/r=2.46\times 10^{-5}$
| – | 10 | 10 | 1 | 10 | $2\times 10^{4}$ | $1.1\times 10^{-3}$ | 4.1 | 0.14 | High $V_{m}$a
| – | 1 | 10 | 0.3 | 5 | 3000 | $4.1\times 10^{-5}$ | 2.4 | 0.1 |
| – | 0.1 | 6 | 1 | 55 | 3000 | $1.3\times 10^{-4}$ | 0.56 | 0.12 | Breeze
| – | 0.3 | 6 | 1 | 55 | 3000 | $4.3\times 10^{-4}$ | 0.65 | 0.14 | Breeze
| – | 1 | 6 | 1 | 15 | 3000 | $2\times 10^{-3}$ | 0.77 | 0.19 |
| – | 3 | 6 | 1 | 15 | 3000 | $4.2\times 10^{-3}$ | 2.6 | 0.21 |
| – | 10 | 6 | 1 | 15 | 3000 | $6\times 10^{-3}$ | 4.1 | 0.19 |
| – | 1 | 3 | 1 | 15 | 3000 | $0.014$ | 0.75 | 0.29 |
| – | 3 | 3 | 1 | 15 | 3000 | $0.04$ | 1.1 | 0.33 |
| – | 10 | 3 | 1 | 5 | 3000 | $0.065$ | 3.0 | 0.3 |
Streaming & Diffusion |  |  |  |  |  |  |  |  |  |
| 0.03 | 1 | 10 | 1 | 5 | 3000 | $2.6\times 10^{-3}$ | 0.93 | 0.33 |
| 0.33 | 1 | 10 | 1 | 5 | 3000 | 0.012 | 1.0 | 0.5 |
| 3.3 | 1 | 10 | 1 | 5 | 3000 | 0.033 | 4.3 | 0.72 |
Diffusion |  |  |  |  |  |  |  |  |  |
| $3.3$ | 0 | 10 | 1 | 15 | 3000 | 0.03 | 4.6 | 1 |

a Differences in velocity and $\dot{E}_{w}/\dot{E}_{c}$ relative to
$V_{m}=6000$ are primarily due to the larger box: $v(r=5)\simeq 2.6V_{g}$ &
$\dot{E}_{w}(r=5)/\dot{E}_{c}(r_{0})\simeq 0.19$.

##  3 Numerical Simulations

###  3.1 Equations

We solve the cosmic ray hydrodynamic equations in one dimensional spherical
polar coordinates using the two moment approach developed by Jiang & Oh (2018)
. 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}\log r$.

The magnetic field is assumed to be a split-monopole, with

$\displaystyle B(r)=B_{0}(r_{0}/r)^{2},$ (41)

where $B_{0}$ is the magnetic field at the bottom boundary of our simulation
box with radius $r_{0}$. We do not evolve the magnetic field in our 1D
calculations. It is only used to calculate the Alfvén velocity, which is
needed for cosmic ray streaming.

The full set of equations for gas density $\rho$, flow velocity $v$, CR energy
density $E_{c}$, CR pressure $p_{c}=E_{c}/3$ 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+v_{s})\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].$ (42)

The streaming velocity is $v_{s}=-v_{A}\nabla p_{c}/|\nabla p_{c}|$, where the
Alfvén velocity is calculated based on the assumed magnetic field and gas
density $\rho$ as $v_{A}=B/(4\pi\rho)^{1/2}$. The reduced speed of light is
$V_{m}$, which is chosen to be much larger than $v$ and $v_{A}$ in the whole
simulation box. We will primarily carry out simulations with cosmic ray
transport by streaming; in a few cases, diffusion will be included as well.
For the case of streaming,

$\displaystyle\sigma_{c}=\frac{|\nabla p_{c}|}{v_{A}(E_{c}+p_{c})}.$ (43)

while for cases with both diffusion555As in the radiation transfer literature,
Jiang & Oh (2018) write the diffusive flux as $F_{c}=-\sigma_{c}^{-1}\nabla
p_{c}$. In the rest of this paper, and in Paper I, we follow the CR literature
and define the diffusive flux as $F_{c}=-\kappa\nabla E_{c}$ (eq. 1). This
accounts for the factor of 3 relating $\kappa$ and $\sigma_{c}^{-1}$ in
equation 44. and streaming

$\displaystyle\sigma_{c}^{-1}=3\kappa+\frac{v_{A}(E_{c}+p_{c})}{|\nabla
p_{c}|}$ (44)

Substituting the right-hand-side of the fourth of equations 42 into the 3rd
term on the right-hand-side of the 2nd of equations 42 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
(this latter term is small in our simulations, as we have verified post
facto). Note also that multiplying the fourth of equations 42 by equation 44
yields equation 2, along with another term related to the time variation of
the CR flux. The time dependent term in the CR flux equation is dynamically
important in our simulations in some regions as we will describe more below.
Physically this means that the assumption that the CR energy flux is given by
equation 2 is not well-satisfied in many locations in our simulations.

A few additional comments about the physical content of the two-moment model
are in order. In the presence of CR diffusion, so that $\kappa\neq 0$,
$\sigma_{c}$ is always finite. This is not the case for transport by streaming
alone with $\kappa\rightarrow 0$. In that case, $\sigma_{c}\rightarrow 0$ when
$\nabla p_{c}\rightarrow 0$ and the two-moment CR equations reduce to a
hyperbolic system in which the intrinsic transport speed is the speed of light
$V_{m}$. Physically, this captures the fact that if $\nabla p_{c}=0$ there is
nothing driving the streaming instability and thus nothing to limit the CR
transport to the Alfvén speed. This is also precisely the regime in the which
the time-dependent term in the fourth of equations 42 is important. A
transport speed of $V_{m}$ does not, however, imply $F_{c}\sim V_{m}p_{c}$.
Rather, the relation between CR flux and energy density is set by boundary
conditions. This is analogous to the fact that for photons, $F\sim pc$ in the
optically thin limit for a single source but not if one is measuring the flux
in between two equal sources of photons, in which case $F\sim 0$.

In what follows we refer to the CRs as “coupled" to the gas when there is a
finite CR pressure gradient such that the CR energy flux due to streaming is
close to the one-moment value of $4p_{c}(v_{A}+v)$. Formally, this holds when
the time dependent term in the CR flux equation (eq. 42) can be neglected.

###  3.2 Initial and Boundary Conditions

For each simulation, we pick gas density $\rho_{0}$, magnetic field $B_{0}$
(parameterized by $v_{A,0}=v_{A}(r_{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}(r_{0})\left(r/r_{0}\right)^{-V_{g}^{2}/c_{i}^{2}}.$ (45)

We apply floor values of $10^{-4}$ for $\rho/\rho_{0}$ and
$p_{c}(r)/p_{c}(r_{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, we fix the cosmic ray energy density
$E_{c}$ and gas density $\rho$ to the desired values at $r_{0}$. Then we
determine the gradients of $E_{c}$ and $\rho$ in the ghost zones by assuming
steady state. 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. The
cosmic ray flux in the ghost zones is set to be $(v+v_{A})(E_{c}+p_{c})$ for
the streaming simulations.

For the boundary condition at the top of the simulation domain, we keep the
gradient of $\rho$ and $E_{c}$ continuous across the boundary but find that it
is necessary to set the cosmic ray flux to be $V_{m}E_{c}/3^{1/2}$ (so that
$F_{c}/E_{c}>v_{A}$) at the top boundary for the streaming case to get an
outflow solution (for the diffusion-only case it is sufficient to simply keep
the gradient of $F_{c}$ continuous across the outer boundary). The flow
velocity at the top boundary is set by requiring the mass flux $\rho vr^{2}$
to be continuous. The boundary conditions for the simulations with both
streaming and diffusion are the same as those with streaming alone.

###  3.3 Simulation Suite

Table 3: Summary of Units for Numerical Simulations. $c_{i}$ is the gas
isothermal sound speed, $\rho_{0}$ is the gas density at the base of the wind
(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}$ Alfvén Speed
$v_{A}$ $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}$

Table LABEL:tab:compare summarizes our 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 base
Alfvén speed $v_{A,0}$ (in units of $c_{i}$) and, for a few of the
simulations, the diffusion coefficient $\kappa$ (in units of $r_{0}c_{i}$).
The key numerical parameters are the resolution, reduced speed of light
$V_{m}$, and box size. Table LABEL:tab:compare also includes one of our
diffusion only simulations from Paper I, for comparison to the analogous
simulations with both streaming and diffusion and streaming alone. In §5.2 we
present a detailed comparison of the streaming results of this paper and the
diffusion results from Paper I.

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}$.

###  3.4 Streaming-Only Simulations

Figure 2: Left: Instantaneous CR pressure and density profiles in several of
our $V_{g}=10$ simulations. The isothermal gas sound speed $c_{i}=1$ in our
units so gas density is also gas pressure. The flow is permeated by strong
shocks, which are sourced by a linear instability near the base of the wind
(§3.4.1). The density increase in the vicinity of the shock triggers the
‘bottleneck’ effect which locally flattens the CR pressure profile leading to
the staircase like structure in $p_{c}$. Figure 4 shows a zoom-in on the
structure of one of the shocks. Right: Density profile near the base of the
wind where the instability and shocks originate; they grow more quickly at
higher $v_{A,0}$ and higher resolution (HR); see §3.4.1 and Appendix A.

In this section we summarize the properties of our numerical models of
galactic winds with cosmic-ray transport via streaming at the Alfvén speed.
The models explored are summarized in Table LABEL:tab:compare. We also compare
the numerical solutions to the analytic solutions summarized in §2 (and, in a
few cases, to the modified analytic models in §4 that are motivated by the
numerical solutions). We begin by analyzing simulations with $V_{g}/c_{i}=10$
and $v_{A,0}/c_{i}=1,10$ because these most dramatically highlight some of the
new features revealed by our time-dependent simulations. We then consider
different values of $V_{g}/c_{i}$ in §3.5 and a wider range of $v_{A,0}$ in
§3.6.

A key feature of nearly all of the streaming simulations in this paper is that
they are time dependent due to an instability that develops near the base of
the wind that steepens into shocks. This is shown in Figure 2, which plots the
CR pressure and gas density as a function of radius once the simulations have
evolved to a statistical steady state.666The full domain is larger than that
shown in Figures 2-3 (see Table LABEL:tab:compare), but we focus on
$r/r_{0}=1-2$ because otherwise the number of shocks in the Figures is large
and obscures legibility. Recall that in our units the isothermal sound speed
is 1 so that the gas density in Figure 2 and subsequent Figures is also the
gas pressure. The right panel of Figure 2 zooms in further to the inner region
near the base of the wind where the instability and shocks originate. We
derive the physical origin of the instability in Appendix A and summarize the
results below in §3.4.1. As a point of contrast, we note that the diffusion
simulations in Paper I nearly all reached a laminar steady state with no
analogous instability or shocks.

Additional features of the time-dependent streaming solution for $v_{A,0}=1$
are shown in Figures 3 & 4: the left panel of Figure 3 shows the gas velocity,
Alfvén velocity, and the ‘isothermal’ CR sound speed $(p_{c}/\rho)^{1/2}$;
while the right panel shows the true CR flux from our time dependent solution
compared to the steady state flux assumed in one-moment CR transport models.
Figure 4 zooms in and shows the detailed structure of one of the shocks.

The most striking features of Figure 2 are (1) the outflows are permeated by
strong shocks in which the gas density increases by at least a factor of $\sim
10$, and (2) the CR pressure profile is not continuous, but is instead largely
flat with significant changes only near the shocks. These two features are
related because of the bottleneck effect: if CRs are well-coupled to the gas
by scattering due to the streaming instability, they stream at the Alfvén
speed relative to the gas and the solution to the steady state CR energy
equation for $v_{A}\gg v$ is $p_{c}\propto\rho^{2/3}$ (see §2). In regions
where $\rho$ increases outwards, which occupy a large fraction of the volume
in our wind solutions (Fig. 2), $p_{c}$ thus needs to increase outwards if the
CRs are well-coupled. This conclusion is, however, inconsistent with the fact
that CRs stream down the CR pressure gradient. As a result, the only physical
solution when the gas density increases outwards is that the CR pressure is
flat. The absence of a CR pressure gradient implies that there is no driving
of the streaming instability and hence no scattering to pin the CR streaming
speed at $v_{A}$. A corollary of the fact that $p_{c}\sim$ constant over much
of the domain is that the CR flux deviates from the canonical streaming value
of $F_{c}=4p_{c}(v+v_{A})$. The reason is that when $p_{c}\sim$ constant, the
steady state CR energy equation only requires $r^{2}F_{c}\sim$ constant, with
no direct constraint on its value. Figure 3 (right panel) shows that in
general the CR flux is a factor of few less than the standard streaming
solution, with $F_{c}=4p_{c}(v+v_{A})$ only near the base of the wind prior to
the shocks developing and in the vicinity of the shocks where the CR pressure
changes significantly.

Figure 3: Left: Alfvén velocity, flow velocity, and ‘isothermal’ CR sound
speed for our $V_{g}=10,v_{A,0}=1$ simulation; gas sound speed $c_{i}=1$. Note
the hierarchy of velocities at large radii with
$v_{A}>\sqrt{p_{c}/\rho}>|v|>c_{i}$. Right: CR flux $F_{c}$ evolved in our two
moment formalism and ‘equilibrium’ CR flux $4p_{c}(v+v_{A})$ assumed in one-
moment CR models. The equilibrium CR flux is a good approximation near the the
base of the wind but not exterior to that once the strong shocks and CR
bottleneck develop.

Figure 4 shows the properties of the shocks in detail. From Figure 3 (left
panel) there is a rough hierarchy of velocities at the shocks with
$v_{A}\gg(p_{c}/\rho)^{1/2}\sim v\gtrsim c_{i}$. The gas density jump at the
shock is roughly consistent with an isothermal shock in which the density
increases by a factor of $(V_{s}/c_{i})^{2}$ where $V_{s}$ is the speed of the
shock. This is apriori somewhat surprising given that $p_{c}\sim\rho
V_{s}^{2}$ so one might expect CR pressure to limit how compressive the shocks
are. Indeed, if the CRs were fully coupled, $p_{c}\propto\rho^{2/3}$ would
imply that the density could only jump by order unity for a shock in which
$p_{c}\sim\rho V_{s}^{2}$. The reason this is not realized in our numerical
solutions is that the CRs are in fact not well-coupled throughout most of the
shock. Figure 4 shows this explicitly: the CR pressure is flat in much of the
region where the density changes appreciably and $p_{c}/\rho^{2/3}$ varies
strongly as a result. Another way to see this is that only in a very small
sub-region of the shock is the CR energy flux equal to the coupled one-moment
value of $4p_{c}(v_{A}+v)$.

Figure 4: Zoom-in on the structure of one of the shocks in the
$V_{g}=10,v_{A,0}=1$ simulation. The CRs are well coupled, i.e.,
$dp_{c}/dr\neq 0$ and $F_{c}\approx 4p_{c}(v_{A}+v)$ only in a small subregion
of the shock where the density is larger than average and $v_{A}$ is smaller
than average. This minimizes the overall magnitude of CR streaming losses
($v_{A}dp_{c}/dr$) in the wind. The large density compression is roughly
consistent with that expected for an isothermal shock (recall that $c_{i}=1$
in our units so this is a high Mach number shock).

We interpret the fact that the CRs are not well coupled throughout the shock
as a consequence of a global constraint on the CR pressure profile. This
constraint regulates the change in the CR pressure at each shock. Note that
per Figure 3 most of the flow has a velocity less than the CR sound speed and
so is all in causal contact. As we show in §3.4.2 below, the time-averaged CR
pressure profile roughly satisfies $p_{c}\propto\rho^{1/2}$ so that on average
the CR pressure scale-height is only a bit larger than the gas scale-height.
The CR pressure only changes significantly at the shocks so this constraint on
the CR pressure scale-height implies that on average $N_{shock}\Delta
p_{c}/p_{c}\sim 1/2$ where $N_{shock}$ is the number of shocks per density
scale-height and $\Delta p_{c}/p_{c}$ is the fractional change in CR pressure
at each shock. The instabilities that seed the shock initially grow on a wide
range of length-scales but the longest wavelength modes saturate at the
highest amplitudes, and some shocks merger with each other, leading to
$N_{shock}\sim 1$, i.e., modes with wavelengths comparable to the local
density scale-height dominate. This prediction of $\Delta p_{c}/p_{c}\sim 1/2$
is consistent with the modest change in CR pressure per shock in the numerical
simulations (e.g., Fig. 2).

####  3.4.1 Physical Origin of the Instability

Begelman & Zweibel (1994) showed that sound waves are driven unstable by CR
streaming when the magnetic energy density is larger than the gas energy
density. This result only applies, however, in the absence of strong radiative
cooling. For an isothermal gas equation of state as we use here, sound waves
are stable given the physics included in Begelman & Zweibel (1994)’s analysis.
To understand the origin of the instability in our simulations, in Appendix A
we calculate the local linear stability of an isothermal hydrostatic
atmosphere with CRs. We identify two new instabilities in the presence of CR
streaming. The first is present in one-moment CR models and is driven by
background gradients in density and CR pressure. This instability is a
streaming analogue of the acoustic instability with CR diffusion and a
background CR pressure gradient studied by Drury & Falle (1986). It leads to
sound waves being amplified after propagating a few scale-heights. We believe
that this is the dominant instability in most of our simulations since the
numerical and analytic growth rates are reasonably consistent with each other
(Fig. 17) and the numerical growth rates are largely independent of the
reduced speed of light, consistent with an instability in the well-coupled
one-moment regime. We also identify a second instability driven by CR
streaming in the two-moment formulation that is not present in one-moment CR
transport models, i.e., the instability relies on the finite speed of light.
The unstable mode is a sound wave at high $c_{i}/v_{A}$ but the entropy mode
at low $c_{i}/v_{A}$; the latter is the most relevant to CR-driven wind models
which tend to be magnetically dominated. Physically, the finite speed of light
introduces a phase shift between the CR flux and the CR pressure (eq. 42).
This phase shift effectively produces a negative diffusion coefficient that
amplifies linear waves (eq. 86). We believe that this instability causes the
very rapid short wavelength modes to grow in our highest resolution simulation
in Figure 2.

Galactic winds are likely subject to other CR-driven instabilities that might
affect their dynamics in a way qualitatively akin to that found here. We
return to this point in §6.

####  3.4.2 Quasi-Steady Structure of Streaming Solutions

Figure 5: Time-averaged $V_{g}=10$ solutions in the $p_{c}-\rho$ plane, along
with one instantaneous snapshot. At high densities near the base of the wind,
$p_{c}\propto\rho^{2/3}$. However, once the shocks set in, the effective
equation of state is closer to $p_{c}\propto\rho^{1/2}$. This occurs at
smaller radii for higher $v_{A,0}$ and/or higher resolution (Fig. 2).

In this section we show that despite the time-dependent nature of our galactic
wind simulations with CR streaming, the time-averaged properties of the winds
can be understood using a modified version of CR-driven wind theory. The key
modification is that the inhomogeneous nature of the flow leads to an average
relation between CR pressure and gas density that differs from the standard
$p_{c}\propto\rho^{2/3}$ result derived for CR streaming at high $v_{A}$ (eq.
12). To see this, Figure 5 shows the time-averaged CR pressure vs. time-
averaged CR density for three of our simulations. For comparison, we also show
the instantaneous $p_{c}(\rho)$ profile for the $v_{A,0}=1$ simulation. Near
the base of the wind $p_{c}\propto\rho^{2/3}$. Once the shocks develop,
however, the time-averaged relation between $p_{c}-\rho$ flattens to be closer
to an effective adiabatic index $\gamma_{\rm eff}\sim 1/2$, i.e.,
$p_{c}\propto\rho^{1/2}$. This transition happens closer to the base of the
wind at higher $v_{A,0}$ and higher resolution, consistent with where the
shocks first develop in these simulations in Figure 2.

Figure 6: Time-averaged properties of galactic winds with CR streaming for
$V_{g}=10$ and $v_{A,0}=1,10$. Upper Left:. Time-averaged density profiles
compared to our analytics; standard high $v_{A}$ (${\rm\gamma_{eff}}=2/3$)
models underestimate the density at larger radii and thus the mass-loss rate
while the lower ${\rm\gamma_{eff}}$ models from §4 do much better. Upper
Right:. Time-averaged mass and energy loss-rates for $V_{g}=10$ & $v_{A,0}=1$.
The total energy flux is dominated by streaming even at large radii; the gas
kinetic energy flux is subdominant. The green solid line quantifies how the
inhomogeneous flow, with strong shocks surrounded by bottleneck regions in
which $p_{c}=$ constant (Fig. 2), minimizes streaming losses by concentrating
the regions with a large CR pressure gradient where the density is high and
the Alfvén speed is low. This suppression of streaming losses makes
${\rm\gamma_{eff}}\sim 0.5$ a better approximation than the canonical
${\rm\gamma_{eff}}=2/3$ solution (see Fig. 5). Lower Left: Time-averaged
velocity profiles. The green vertical dotted line is the analytic prediction
of the sonic point from eq. 59. Lower Right: Numerator and Denominator of the
steady state wind equation (eq. 11). Because the flow is time dependent it
does not strictly satisfy the steady state wind equation. Nonetheless the
analytic prediction of the sonic point for ${\rm\gamma_{eff}}=0.5$ (vertical
dotted line) is a reasonable estimate of the radius at which $v\sim c_{c}\sim
V_{g}$ and $N_{st}\sim D_{st}\sim 0$.

To understand why $p_{c}(\rho)$ becomes shallower than
$p_{c}\propto\rho^{2/3}$, it is helpful to remember the origin of the latter:
for a steady solution with $v_{A}\gg v$, the CR energy equation becomes
$r^{-2}d(r^{2}F_{c})/dr=v_{A}dp_{c}/dr$ with $F_{c}=4p_{c}v_{A}$, which
implies $p_{c}\propto\rho^{2/3}$ for a split-monopole. Physically, the change
in CR power $\sim r^{2}F_{c}$ is due to energy transfer from the CRs to the
gas at a rate $v_{A}dp_{c}/dr$. In the time dependent simulations, however,
over most of the volume $dp_{c}/dr=0$ and the CRs and gas are not well-coupled
and do not exchange energy with each other. With $dp_{c}/dr=0$,
$r^{2}F_{c}\sim$ const so that $F_{c}\sim p_{c}v_{A}$ implies
$p_{c}\propto\rho^{1/2}$. In more detail, the exchange of energy between CRs
and the gas, i.e., $dp_{c}/dr\neq 0$, occurs only in the vicinity of the
shocks which are regions of higher than average density and thus lower than
average $v_{A}$ (see Fig. 3 and 4). This inhomogeneity causes $\langle
v_{A}dp_{c}/dr\rangle<\langle v_{A}\rangle\langle dp_{c}/dr\rangle$, where
$\langle\rangle$ denotes a time average, as we show explicitly in the upper
right panel of Figure 6 (green line). As a result, the time-averaged energy
equation for the CRs is better modeled as
$r^{-2}d(r^{2}F_{c})/dr=fv_{A}dp_{c}/dr$ with $f\sim 1/3$ or so. In the limit
$f\ll 1$ and $F_{c}\sim p_{c}v_{A}$, the time-averaged CR energy equation
again reduces simply to $r^{2}F_{c}\sim$ constant, i.e.,
$p_{c}\propto\rho^{1/2}$. This is roughly consistent with the time-averaged
$p_{c}(\rho)$ in Figure 5 and the time-averaged CR luminosity shown in Figure
6 (solid line in the upper right panel), which does not change much with
radius away from the base of the wind.

Figure 6 shows additional properties of our time-averaged galactic wind
solutions with CR streaming. In the upper left panel, we show that the density
profiles near the base of the wind, where the solution is roughly hydrostatic,
are well approximated by the analytic theory developed in §4, in which the CRs
are modeled as a fluid with an effective adiabatic index $\gamma_{\rm
eff}<2/3$, consistent with the modified CR energetics summarized in Figure 5.
The change from ${\rm\gamma_{eff}}=2/3$ to ${\rm\gamma_{eff}}=1/2$ leads to
significantly higher gas densities (factor of $\sim 10$) in the nearly
hydrostatic portion of the wind. As a result, the mass-loss rate in the wind
is larger than standard wind theory would predict: in §4 we show analytically
that the seemingly modest of change of $\gamma_{\rm eff}=2/3\rightarrow 0.5$
in fact can change the mass-loss rate in CR-driven winds by a factor of a
few-100. This is driven primarily by the change in the gas density profile
seen in Figure 6.

The lower left panel of Figure 6 shows the time-averaged velocity profiles of
the wind for $v_{A,0}=1,10$ (in this plot, and those that follow, we define
the time-averaged velocity as $\langle\rho v\rangle/\langle\rho\rangle$). The
gas accelerates to a speed somewhat larger than both the CR sound speed and
the circular velocity of the potential. The lower right panel of Figure 6
shows the numerator and denominator of the steady state wind equation (eq. 9).
Both pass through 0 as expected for a steady wind. However, the numerator and
denominator do not vanish at the same radius, as they would for a true steady
state wind (e.g., Fig. 3 of Paper I). The reason is that the wind equation
does not formally hold for the time average of a time dependent
solution.777For example, interpreting eq. 9 in terms of time-averaged
variables requires assuming $\langle\rho v\rangle=\langle\rho\rangle\langle
v\rangle$, which is not in general true for inhomogeneous flows like those
found here. The vertical dotted line in the lower panels of Figure 6 is the
analytic prediction of the sonic point for our modified solutions from §4,
$r_{s}\simeq er_{0}\simeq 2.7r_{0}$ (see eq. 49); the standard CR-wind theory
in §2 predicts a larger sonic point radius of $r_{s}=e^{3/2}r_{0}\simeq
4.5r_{0}$ (eq. 11). Despite the significant time dependence of the solutions,
the analytic estimate does a reasonable job of capturing the location where
the flow reaches $v\sim c_{c}\sim V_{g}$, even though the solution does not
rigorously satisfy the time-steady critical point conditions.

Finally, the upper right panel of Figure 6 quantifies the energy fluxes in the
steady state solution. For the steady wind problem with CR streaming and
isothermal gas there is no conserved total energy flux because streaming
losses from the CRs to the gas are assumed to be instantaneously radiated
away. Nor is there formally a separately conserved gas or CR energy flux. As
noted in our discussion of the origin of $p_{c}\propto\rho^{1/2}$, however, CR
streaming losses are less in our inhomogeneous wind solutions than in standard
wind theory, so that $\propto r^{2}F_{c}$ is only a weak function of radius,
particularly away from the base of the wind. Figure 6 also shows that the
steady state CR energy flux is a factor of few lower than the canonical value
of $4(v+v_{A})p_{c}$, consistent with the individual snapshot in Figure 3.
Finally, the gas kinetic energy flux is a factor of $\sim 10$ less than the CR
energy flux, even at large radii where the gas is supersonic.

###  3.5 Solutions For Different $V_{g}$

Figure 7: Left: Instantaneous density profiles for different $V_{g}$, with two
different $v_{A,0}$ for $V_{g}=3$. The instability and shocks that
significantly modify the wind dynamics are present in all cases except for
$V_{g}=3,v_{A,0}=1$ (see §3.5 for a physical explanation). Right:. Time-
averaged radial velocity and Alfvén velocity for the same four simulations.
The kinematics of the four different simulations are similar scaled by $V_{g}$
as we have done here.

In this section we present wind solutions for different values of
$V_{g}/c_{i}=3,6,10$. Decreasing $V_{g}/c_{i}$ corresponds to a shallower
potential for fixed gas thermodynamics or hotter gas for a given escape speed.
The latter, i.e., varying $c_{i}$, can physically be thought of as describing
cosmic-rays coupling primarily to different phases of the ISM (e.g., warm
ionized medium for lower $c_{i}$ vs. the hot ISM for higher $c_{i}$). For
these three values of $V_{g}/c_{i}$, the left panel of Figure 7 shows
instantaneous density profiles (at randomly chosen times in the statistical
steady state) while the right panel shows time-averaged velocity and Alfvén
speed profiles. We normalize the velocity profile using $v/V_{g}$ rather than
$v/c_{i}$ because this highlights the similarity of the kinematics across a
range of $V_{g}$.

Figure 8: Time-averaged solutions for different $V_{g}$ in the $p_{c}-\rho$
plane. Solutions with $V_{g}=6$ show the transition to
$p_{c}\propto\rho^{1/2}$ just as the $V_{g}=10$ solutions do in Figure 5. By
contrast, solutions with $V_{g}=3$ are more consistent with
$p_{c}\propto\rho^{2/3}$ (until at lower densities $v>v_{A}$ and
$p_{c}\propto\rho^{4/3}$). This is due to the smaller density fluctuations
introduced by shocks when the gas temperature is effectively larger (smaller
$V_{g}/c_{i}$); see Fig. 7.

The left panel of Figure 7 shows that for our fiducial $v_{A,0}=c_{i}$
magnetic field, the strong shocks are present for $V_{g}=6$ but absent for
$V_{g}=3$. The reason for the latter is that the sound wave instability
identified in Appendix A operates most effectively when $v_{A}\gg v$; in the
opposite regime of $v<v_{A}$ there is no instability because the gas is nearly
adiabatic. Smaller values of $V_{g}$ correspond to significantly larger
density scale-heights (left panel of Fig. 7) and thus smaller values of
$v_{A}$. This inhibits the acoustic instability that is driven by rapid
streaming. Additional evidence for this interpretation is that a simulation
with $V_{g}=3$ and a larger base Alfvén speed, $v_{A,0}=3$ (also shown in
Figure 7), does show the instability and shocks with similar properties to the
higher $V_{g}$ simulations.

A second striking result of Figure 7 is that simulations with different values
of $V_{g}$ have similar kinematics, with terminal velocities that are within a
factor of $\lesssim 2$ of $V_{g}$. We compare this to our analytic predictions
in §5.1.

Figure 8 shows the time-averaged solutions from Figure 7 in the $p_{c}-\rho$
plane (analogous to Figure 5). The $V_{g}=6$ simulation is similar to the
$V_{g}=10$ simulation, with $p_{c}\propto\rho^{1/2}$ a better approximation to
the effective equation of state at intermediate densities
$10^{-2}-10^{-3.5}\rho_{0}$ where the shocks are present and $v_{A}/v$ is
large. This is not really the case for the $V_{g}=3$ simulations, including
the $v_{A,0}=3$ simulation in which strong shocks are present just as in the
$V_{g}=6$ and 10 simulations. Our interpretation is that this is because for
higher $c_{i}/V_{g}$, the density jump at the shocks is smaller, as is evident
in the left panel of Figure 7. The smaller density contrast between the shocks
and the rest of the solution minimizes the inhomogeneous nature of the flow,
and thus the differences relative to canonical streaming solutions. A smaller
value of $v_{A}/v$ also increases the effective adiabatic index of the gas,
with ${\rm\gamma_{eff}}\rightarrow 4/3$ for $v\ll v_{A}$ (this indeed applies
at the lowest densities in Fig. 8). However, the $V_{g}=6,v_{A,0}=1$ and
$V_{g}=3,v_{A,0}=3$ solutions in Figure 7 have very similar values of
$v_{A}/v$ as a function of radius, and yet different ${\rm\gamma_{eff}}$ in
Figure 8. This points to the magnitude of the density jump at the shocks as
the primary reason that the standard ${\rm\gamma_{eff}}=2/3$ solution is a
better approximation for $V_{g}=3$ (and both $v_{A,0}$) than it is for
$V_{g}=6$ and higher.

###  3.6 Solutions For Different $v_{A,0}$

Figure 9 shows time-averaged radial velocity and Alfvén velocity profiles for
different values of $v_{A,0}=0.1,0.3,1$, and 3, all for $V_{g}=6$ and
$p_{c,0}=\rho_{0}c_{i}^{2}$. The most striking result in Figure 9 is that the
two low base Alfvén velocity solutions ($v_{A,0}=0.1,0.3$) never reach $v\sim
V_{g}$, i.e., they never reach the sonic point. These solutions are formally
‘breezes’ rather than transonic winds. Note also that since $v_{A}\lesssim v$
at large radii for these solutions, the effective adiabatic index of the gas
${\rm\gamma_{eff}}$ is increasing towards the adiabatic value of $4/3$ (eq.
14). Once ${\rm\gamma_{eff}}\gtrsim 1$, such that the CR sound speed
$\propto\rho^{{\rm\gamma_{eff}}-1}$ decreases with decreasing density, the
flow can no long accelerate by tapping into the increasing CR sound speed.
Thus the low $v_{A,0}$ solutions in Figure 9 will not be able to accelerate to
$v\sim V_{g}$ at yet larger radii. By contrast, the $v_{A,0}=3$ solution in
Figure 9 is a robust transonic wind and $v_{A,0}=1$ solution is on the border
with $v\sim V_{g}$.

The transition between transonic winds and breezes in Figure 9 at $v_{A,0}\sim
1$ is reasonably well captured by our analytic estimate of $v_{A,crit}$, the
base Alfvén velocity required to maintain $v_{A}>v$ out to the sonic point
(§2.2). For the parameters of Figure 9, equation 28 yields $v_{A,crit}\simeq
0.5$ while equation 64 (derived in the next section for the modified CR
energetics found in our simulations) yields $v_{A,crit}\simeq 0.6$. This
supports our conjecture in §2.2 that $v_{A,0}\gtrsim v_{A,crit}$ is an
approximate condition for a supersonic wind.

It is also instructive to compare the mass-loss rate for the solutions with
varying $v_{A,0}$ to the maximum mass-loss rate allowed by energy conservation
(§2.4). For our simulations with a $\ln(r)$ potential, we can define the
escape speed as the speed at $r_{0}$ needed to just reach the outer radius of
the box, i.e., $v_{esc}=2\sqrt{\ln(r_{out})}V_{g}$. Given the energy flux in
the wind at $r_{out}$ in Table LABEL:tab:compare, we find that the simulations
in Figure 9 with $v_{A,0}=0.1,0.3,1,\&\,3$ have
$\dot{M}_{w}/\dot{M}_{max}=0.78,0.74,0.51,\&\,0.33$, respectively. The
$V_{g}=6$ simulation with $v_{A,0}=10$ has $\dot{M}_{w}/\dot{M}_{max}=0.15$.
These results confirm the conjecture in §2.2 & 2.4 that solutions with smaller
values of $v_{A,0}$ (relative to $v_{A,crit}$; eq. 28) have mass-loss rates
approaching the maximum possible given the CR energy available in the wind.
$\dot{M}_{w}\simeq\dot{M}_{\rm max}$ is also consistent with the low speeds
$\lesssim V_{g}$ of the lower $v_{A,0}$ models at all radii in Figure 9.

Figure 9: Time-averaged velocity (solid) and Alfvén velocity (dotted) profiles
for $V_{g}=6$, $p_{c,0}=\rho_{0}c_{i}^{2}$ and four different values of the
base Alfvén velocity. The solutions with lower base Alfvén speed,
$v_{A,0}=0.1,0.3$, never reach $v\sim V_{g}$; these are outflows but ‘breezes’
rather than transonic winds.

###  3.7 Solutions with Streaming and Diffusion

Figure 10: Time-averaged properties of galactic winds compared for CR
streaming alone ($v_{A,0}=1$), diffusion alone ($\kappa=3.3$) and both
streaming and diffusion for several diffusion coefficients; all simulations
are for $V_{g}=10$. With increasing $\kappa$, the solutions with streaming and
diffusion approach the pure diffusion solutions; even small diffusion
coefficients $\kappa\sim 0.033-0.33$ (which on their own do not drive fast
winds), significantly modify the streaming solution. The kinks in the time-
averaged density and velocity profile for $\kappa=3.3$ at $r\sim 1.7$ are
associated with the onset of shocks at this radius (see Fig. 11).  Figure 11:
Instantaneous density profiles for solutions with streaming ($v_{A,0}=1$;
$\kappa=0$) and diffusion; all simulations take $V_{g}=10$. For increasing
diffusion coefficient $\kappa$, the instability giving rise to shocks is
suppressed and shows up at progressively larger radii. For $\kappa=3.3$, the
solution is stable throughout.

CR transport is unlikely to be well-modeled as either pure streaming or pure
diffusion: even in models of self-confinement of CRs by the streaming
instability, there is in general a diffusive-like correction to streaming
whose magnitude is larger when damping processes counteract the streaming
instability and suppress the amplitude of the resulting Alfvénic fluctuations
(e.g., Wiener et al. 2013a). To understand the implications of this property
of CR transport for galactic winds, we briefly consider the properties of
winds with both streaming and diffusion. We find that even a modest diffusion
coefficient significantly modifies the properties of streaming-only solutions.

Figure 10 compares the steady state density, pressure, and velocity profiles
of winds with pure streaming (a fiducial $v_{A,0}=1$ simulation), pure
diffusion ($\kappa=3.3r_{0}c_{i}$) and streaming plus diffusion for a range of
diffusion coefficients (all in units of $r_{0}c_{i}$). The mass loss rates for
these simulations are given in Table LABEL:tab:compare. Even a small diffusion
coefficient of $\kappa\sim 0.033-0.33$, which absent streaming does not
produce a high speed wind (see Paper I), changes the mass-loss rate of the
wind by a factor of $\sim 10-20$ relative to streaming alone. Indeed, by
$\kappa=1/3$, the mass-loss rate in the simulation with streaming and
diffusion is within a factor of a few of the pure diffusion solutions.
Likewise Figure 10 shows that the density and CR pressure profiles in the
simulations with both streaming and diffusion are within a factor of a few of
the pure diffusion result by $\kappa\sim 1/3$. The velocity profiles of the
simulations with both transport processes are somewhat less sensitive to
$\kappa$ and only accelerate as rapidly as the pure diffusion solution once
$\kappa\sim 3$ (Fig. 10).

The strong sensitivity of the galactic wind solutions driven by CR streaming
to a modest diffusion coefficient can be understood as follows. We can
estimate the magnitude of the diffusive CR flux in a calculation with
streaming alone. The CR pressure gradient is then set by the density scale
height which is $H\simeq r_{0}(c_{i}/V_{g})^{2}$ for our spherical model. Thus

$\frac{F_{\rm diff}}{F_{\rm
str}}\simeq\frac{\kappa}{r_{0}c_{i}}\frac{V_{g}^{2}}{v_{A,0}c_{i}}\simeq
100\frac{\kappa}{r_{0}c_{i}}\left(\frac{V_{g}}{100\,{\rm
km\,\,s}^{-1}}\right)^{2}\left(\frac{(10\,{\rm
km\,\,s}^{-1})^{2}}{v_{A,0}c_{i}}\right)$ (46)

As a result if $\kappa/(r_{0}c_{i})\gtrsim v_{A,0}c_{i}/V_{g}^{2}$ the
diffusive CR flux will be larger than the streaming flux near the base of the
wind. For the simulation in Figure 10 this condition becomes
$\kappa/(r_{0}c_{i})\gtrsim 0.01$, which is satisfied by all of the
simulations in Figure 10. Physically, solutions with streaming alone have
smaller density and CR pressure scale-heights than diffusion solutions, so
that only a small diffusion coefficient is needed for diffusive CR transport
to be energetically important. The latter tends to establish the much flatter
CR pressure profiles seen in Figure 10, which in turn accelerates the gas more
efficiently and increases the mass-loss rate (see Paper I).

Figure 11 highlights an additional feature of the simulations with both CR
streaming and diffusion: even modest diffusion coefficients suppress the
linear instability and shocks described in §3.4. In particular, with
increasing diffusion coefficient, the strong shocks in the outflow only show
up at larger radii and by $\kappa=3.3$ the solution is laminar and steady
state with no shocks.

##  4 Analytic Approximations with Modified Cosmic Ray Energetics

The results in §3.4 show that the standard assumptions of models of galactic
winds driven by streaming CRs, in particular equations 10 & 12, do not in fact
apply to the time-averaged properties of many of our numerical wind solutions.
A self-consistent generalization of CR wind theory with streaming that
accounts for this would require developing in detail a new closure model for
the time-averaged CR energy equation. Appendix B briefly discusses some of the
subtleties of doing so. Here we restrict ourselves to the simpler task of
using an effective equation of state in which
$p_{c}\propto\rho^{\rm\gamma_{eff}}$ with ${\rm\gamma_{eff}}<2/3$. This
generalizes the $v_{A}\gg v$ analytics of §2 to account for the different CR
energetics found in the simulations (Fig. 5). We first consider general
${\rm\gamma_{eff}}$ and then specialize to ${\rm\gamma_{eff}}=1/2$.

###  4.1 General ${\rm\gamma_{eff}}$

If we model the CRs as a fluid with $p_{c}\propto\rho^{\rm\gamma_{eff}}$, the
effective CR sound speed is

$c_{\rm eff}^{2}={\rm\gamma_{eff}}\frac{p_{c}}{\rho}$ (47)

and the steady state sonic point equations become

$\frac{1}{v}\frac{dv}{dr}=\frac{N_{\rm eff}}{D_{\rm eff}}$ (48)

with

$D_{\rm eff}=v^{2}-c_{i}^{2}-c_{\rm eff}^{2}\ \ \ \ \ \ N_{\rm
eff}=\frac{2}{r}\left(c_{i}^{2}+c_{\rm eff}^{2}-V_{g}^{2}\right)$ (49)

We proceed as in §2 by considering the nearly hydrostatic portion of the flow
interior to the sonic point, for which

$c_{i}^{2}\frac{d\ln\rho}{dr}+c_{\rm
eff,0}^{2}\left(\frac{\rho_{0}}{\rho}\right)^{1-{\rm\gamma_{eff}}}\frac{d\ln\rho}{dr}=-\frac{2V_{g}^{2}}{r}.$
(50)

Equation 50 can be solved numerically for $\rho(r)$ in the hydrostatic portion
of the wind. Figure 6 (top left panel) shows that these solutions reproduce
the density profiles in our time dependent simulations much better than
standard ${\rm\gamma_{eff}}=2/3$ CR wind theory.

The critical point conditions (eqs 49) determine the density and flow velocity
at the sonic point, namely

$v(r_{s})=V_{g}\ \ \ \ {\rm and}\ \ \ \ \rho(r_{s})=\rho_{0}\left(\frac{c_{\rm
eff,0}^{2}}{{\rm V^{2}_{g,eff}}}\right)^{\frac{1}{1-{\rm\gamma_{eff}}}}$ (51)

The mass-loss rate in the wind can then be estimated as

$\dot{M}_{w}=4\pi r_{s}^{2}\rho(r_{s})V_{g}$ (52)

Note that $\rho(r_{s})$ in equation 51 is analytic in terms of the base
properties of the wind, but the location of the sonic point $r_{s}$ can in
general only be determined numerically given the solution to equation 50.
Figure 12 shows the resulting mass-loss rate (in units of $\dot{M}_{0}$; eq.
17) for ${\rm\gamma_{eff}}=1/2$ as a function of $p_{c,0}/\rho_{0}c_{i}^{2}$
and $V_{g}/c_{i}$, as well as the ratio of the mass-loss rate predicted by
${\rm\gamma_{eff}}=1/2$ to that predicted by ${\rm\gamma_{eff}}=2/3$. The
larger CR pressure predicted by ${\rm\gamma_{eff}}=1/2$ significantly
increases the mass-loss rate, particularly for larger $V_{g}/c_{\rm eff,0}$.

Figure 12: Left:. Analytic mass-loss rate for CR driven galactic winds in the
limit of rapid streaming (large $v_{A}$) for ${\rm\gamma_{eff}}=1/2$ (see eq.
52 and associated text), as a function of the strength of gravity relative to
the gas sound speed in the disk ($V_{g}/c_{i}$) and the base CR pressure
(${\rm p_{c,0}/\rho_{0}c_{i}^{2}}$). The mass-loss rate is normalized by eq.
17. Right:. Ratio of the predicted mass-loss rate for ${\rm\gamma_{eff}}=1/2$
to that for ${\rm\gamma_{eff}}=2/3$ (again using eq. 52 with the appropriate
${\rm\gamma_{eff}}$ in eqs. 50 and 51). The modified CR thermodynamics due to
strong shocks leading to ${\rm\gamma_{eff}}=1/2$ can significantly increase
the mass-loss rate, particularly for larger $V_{g}/c_{i}$ or smaller
$p_{c,0}/\rho_{0}c_{i}^{2}$. For both panels, labeled values of $V_{g}/c_{i}$
on the color bar are logarithmically distributed and correspond to the curves
on the plot.

As in §2 we can analytically approximate the density profile in Figure 6 and
the mass-loss rate in Figure 12 if we focus on either of the limits $c_{\rm
eff,0}\gg c_{i}$ or $c_{\rm eff,0}\ll c_{i}$. Consider first the case of
massive galaxies with $V_{g}\gg c_{\rm eff,0}\gtrsim c_{i}$. In this limit,
the solution of equation 50 is

$\rho(r)\simeq\rho_{0}\left(1+2(1-{\rm\gamma_{eff}})\frac{V_{g}^{2}}{c_{\rm
eff,0}^{2}}\ln[r/r_{0}]\right)^{-\frac{1}{1-{\rm\gamma_{eff}}}},$ (53)

the sonic point is located at

$r_{s}\simeq r_{0}e^{\frac{1}{2(1-{\rm\gamma_{eff}})}}$ (54)

and the mass loss rate is

$\dot{M}_{w}\simeq
4\pi\rho_{0}r_{0}^{2}V_{g}e^{\frac{1}{1-{\rm\gamma_{eff}}}}\left(\frac{c_{\rm
eff,0}^{2}}{V_{g}^{2}}\right)^{\frac{1}{1-{\rm\gamma_{eff}}}}.$ (55)

In the opposite limit of weak CR pressure compared to gas pressure at the base
in the galactic disk, i.e., $c_{\rm eff,0}\ll c_{i}$, the gas density profile
is initially set by gas pressure, and is given by equation 24. As the density
drops, the CR pressure increases in importance relative to the gas pressure.
As in §2, the sonic point condition (eq. 49) requires that the pressure be CR
dominated at the sonic point (assuming $c_{i}<V_{g}$). The transition between
gas pressure and CR pressure support happens at a radius

$r_{tr}\simeq
r_{0}\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{\frac{c_{i}^{2}}{2V_{g}^{2}(1-{\rm\gamma_{eff}})}}.$
(56)

The density profile exterior to this transition radius is then like equation
53 but with a different boundary condition set by continuity at $r_{tr}$. This
yields

$\frac{\rho}{\rho_{0}}\simeq\left(\frac{p_{c,0}}{\rho_{0}c_{i}^{2}}\right)^{\frac{1}{1-{\rm\gamma_{eff}}}}\left[1+\frac{2(1-{\rm\gamma_{eff}})}{{\rm\gamma_{eff}}}\frac{V_{g}^{2}}{c_{i}^{2}}\,\ln\left(\frac{r}{r_{tr}}\right)\right]^{\frac{-1}{1-{\rm\gamma_{eff}}}}\
\ (r>r_{tr})$ (57)

The sonic point is then located at $r_{s}\simeq
e^{0.5/(1-{\rm\gamma_{eff}})}r_{tr}$ and the mass-loss rate is

$\dot{M}_{w}\simeq 4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{\frac{1}{1-{\rm\gamma_{eff}}}}\,\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{\frac{c_{i}^{2}}{V_{g}^{2}(1-{\rm\gamma_{eff}})}}\left(\frac{c_{\rm
eff,0}^{2}}{V_{g}^{2}}\right)^{\frac{1}{1-{\rm\gamma_{eff}}}}.$ (58)

###  4.2 Specialization to ${\rm\gamma_{eff}}=1/2$, i.e.,
$p_{c}\propto\rho^{1/2}$

To obtain somewhat simpler expressions, we now focus on the case of
${\rm\gamma_{eff}}=1/2$ motivated by the simulations in §3.4. In this case,
the expressions for the sonic point and mass-loss rate for massive galaxies
with $V_{g}\gg c_{\rm eff,0}\gtrsim c_{i}$ become

$r_{s}\simeq e\,r_{0}$ (59)

and

$\begin{split}\dot{M}_{w}\simeq&\,4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{2}\,\left(\frac{\,c_{\rm
eff,0}}{V_{g}}\right)^{4}\simeq 0.023\,{\rm
M_{\odot}\,yr^{-1}}\left(\frac{r_{0}}{1\,{\rm kpc}}\right)^{2}\\\
&\times\,\left(\frac{n_{0}}{1\,{\rm cm^{-3}}}\right)\left(\frac{c_{\rm
eff,0}}{10\,\,{\rm km\,\,s}^{-1}}\right)^{4}\left(\frac{V_{g}}{100\,\,{\rm
km\,\,s}^{-1}}\right)^{-3}.\end{split}$ (60)

while for $c_{\rm eff,0}\lesssim c_{i}\ll V_{g}$ they become

$r_{s}\simeq
e\,r_{0}\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{\frac{c_{i}^{2}}{V_{g}^{2}}}.$
(61)

and

$\dot{M}_{w}\simeq\,4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{2}\,\left(\frac{\,c_{\rm
eff,0}}{V_{g}}\right)^{4}\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{\frac{2c_{i}^{2}}{V_{g}^{2}}}$
(62)

Note that equation 62 for high base CR sound speed differs from equation 60
for low base CR sound speed only by a factor of
$\left(\rho_{0}c_{i}^{2}/p_{c,0}\right)^{2c_{i}^{2}/V_{g}^{2}}>1$ (this is
similar to the standard theory in §2; see the discussion around eq. 27). As a
result, the two can be approximately combined to yield

$\dot{M}_{w}\simeq\,4\pi
r_{0}^{2}\,\rho_{0}\,V_{g}\,e^{2}\,\left(\frac{\,c_{\rm
eff,0}}{V_{g}}\right)^{4}{\rm
Max}\left[1,\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{\frac{2c_{i}^{2}}{V_{g}^{2}}}\right]$
(63)

Perhaps most importantly, the mass-loss rate in equation 60 is a factor of
$\sim(V_{g}/c_{\rm eff,0})^{2}\gg 1$ larger than the standard CR streaming
mass-loss rate (eq. 22). This highlights how the seemingly modest modification
to the CR energetics introduced by the strong shocks in our numerical
simulations (${\rm\gamma_{eff}}<2/3$) in fact significantly modifies the
properties of the resulting winds. Indeed, if we scale for parameters
appropriate to the ‘average’ ISM of the Milky Way with $V_{g}\simeq 150$ km
s-1 and $c_{\rm eff,0}\simeq 10$ km s-1, the CR driven mass loss rate
estimated in equation 55 is a factor of $\sim 100$ larger than the standard
theory.

The expression for the mass-loss rate in equation 60 can be rewritten as
$\dot{M}_{w}\propto p_{c,0}\,c_{\rm eff,0}^{2}/V_{g}^{4}$. This shows that, to
the (uncertain) extent that the CR pressure is roughly similar in different
phases of the ISM, the mass-loss rate will be dominated by the hotter ISM
phases, which have larger values of $c_{\rm eff,0}$. As an example, if we take
$p_{c,0}\simeq 10^{-12}\,{\rm dyne\,cm^{-2}}$, $V_{g}\simeq 150\,{\rm
km\,\,s}^{-1}$, and $r_{0}\sim 3$ kpc as representative of the MW, equation 60
implies $\dot{M}_{w}\simeq 0.06$, 0.5, and 1.5 $\,{\rm M_{\odot}\,yr^{-1}}$
for $c_{\rm eff,0}\simeq 10$, 30, and 50 $\,{\rm km\,\,s}^{-1}$, respectively.
The latter two estimates of the mass-loss rate, representing CR-driven winds
from the warm-hot ISM, are dynamically important in the sense of being of
order the star formation rate, while the direct mass loss from the colder ISM
is negligible. The same qualitative conclusion was true for standard CR
streaming wind theory in §2.

It is also instructive to generalize other properties of CR-driven wind theory
to ${\rm\gamma_{eff}}=1/2$. In particular, the condition for the validity of
the high Alfvén speed limit becomes

$v_{A,0}\gtrsim v_{A,crit}\equiv{c_{\rm eff,0}}\,\frac{V_{g}\,c_{\rm
eff,0}}{V^{2}_{\rm g,eff}}\left(\frac{r_{s}}{r_{0}}\right)^{2}\simeq c_{\rm
eff,0}\,e^{2}\frac{c_{\rm eff,0}}{V_{g}}$ (64)

where the second equality is for massive galaxies with $V_{g}\gg c_{\rm
eff,0}$. As before, the high $v_{A}$ approximation is most applicable for
massive galaxies with large $V_{g}/c_{\rm eff,0}$ (physically, larger
$V_{g}/c_{\rm eff,0}$ corresponds to lower mass-loss rates and thus lower
densities and larger $v_{A}$ for a given base magnetic field strength).

The energetics of the wind can be understood by noting that so long as
$p_{c}\propto\rho^{1/2}$, the CR pressure gradient in equation 30 can again be
rewritten as a CR enthalpy, leading to a conserved Bernoulli-like constant

$\frac{d}{dr}\left(\frac{1}{2}v^{2}+c_{i}^{2}\ln\rho+2V_{g}^{2}\ln r-2c_{\rm
eff,0}^{2}\left[\frac{\rho}{\rho_{0}}\right]^{-1/2}\right)=0.$ (65)

The wind terminal velocity, power, and momentum loss-rate can then be
estimated using arguments that parallel those given in §2. The generalization
of equations 33, 35, 36, and 37 to ${\rm\gamma_{eff}}=1/2$ are

$v_{\infty}^{2}\simeq V_{g}^{2}+4{\rm
V^{2}_{g,eff}}\left(\left[\frac{\rho(r_{A})}{\rho(r_{s})}\right]^{-1/2}-1\right),$
(66)

$\begin{split}v_{\infty}\simeq&\
\frac{2V_{g}}{e}\frac{\sqrt{v_{A,0}V_{g}}}{c_{\rm eff,0}}\simeq 230\,{\rm
km\,\,s}^{-1}\,\times\\\ &\,\left(\frac{V_{g}}{100\,{\rm
km\,\,s}^{-1}}\right)^{3/2}\left(\frac{v_{A,0}}{10\,{\rm
km\,\,s}^{-1}}\right)^{1/2}\left(\frac{c_{\rm eff,0}}{10\,{\rm
km\,\,s}^{-1}}\right)^{-1}\end{split}$ (67)
$\frac{0.5\dot{M}_{w}v_{\infty}^{2}}{\dot{E}_{c}(r_{0})}\sim 0.25,$ (68)

and

$\begin{split}\frac{\dot{p}_{w}}{\dot{p}_{ph}}&\simeq\frac{\dot{M}_{w}}{\dot{M}_{*}}\frac{v_{\infty}}{\epsilon_{ph}c}\\\
&\simeq 2.4\,\epsilon_{ph,\,-3.3}^{-1}\,\left(\frac{c_{\rm eff,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)\left(\frac{p_{*}/m_{*}}{3000\,\,{\rm
km\,\,s}^{-1}}\right)\\\ &\times\left(\frac{V_{g}}{100\,\,{\rm
km\,\,s}^{-1}}\right)^{-3/2}\left(\frac{v_{A,0}}{10\,\,{\rm
km\,\,s}^{-1}}\right)^{1/2}\frac{p_{c,0}}{\pi G\Sigma_{g}^{2}}\end{split}$
(69)

Finally, the maximum mass-loss rate in CR-driven winds set by energy
conservation (§2.4) is

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

where we have assumed that a fraction $0.25\,\zeta_{0.25}$ of the base CR
power is available to drive gas to larger radii. This normalization of $\zeta$
is reasonably consistent with eq. 68 and Figure 13 (bottom panel) discussed in
the next section.

Equation 69 shows that the momentum loss-rate in CR-driven winds with
${\rm\gamma_{eff}}\sim 1/2$ is substantially larger than the standard theory
(eq. 37), with values that are comparable to or exceed the momentum carried by
the stellar radiation field. Equation 68 is particularly striking compared to
the standard theoretical result in equation 36 and highlights how the modified
CR energetics found in our simulations implies a much larger fraction of the
CR energy being transferred to kinetic energy of the gas at large radii.
Equation 68 can be understood physically by noting that
${\rm\gamma_{eff}}=1/2$ is motivated in the first place by the suppression of
streaming losses due to the effects of the strong shocks in the simulations
(e.g., Fig 6). Absent streaming losses $r^{2}p_{c}v_{A}\sim{\rm const}$ so
that a significant fraction of the CR energy at the base of the wind is
available to convert into wind energy at large radii.

##  5 Discussion

In this section we first synthesize the results of our streaming analytics and
simulations (§5.1) and then compare the properties of CR-driven winds with
streaming and diffusive CR transport (§5.2). §5.3 discusses the validity of
the isothermal gas approximation used in this work. We conclude this section
by summarizing some of the observational implications of our results (§5.4).

###  5.1 Synthesis of Streaming Results

Figure 13: Top: Ratio of the analytically predicted-mass loss rate in the
limit of rapid streaming, for both the canonical ${\rm\gamma_{eff}}=2/3$ and
our modified ${\rm\gamma_{eff}}=1/2$, to the simulated mass-loss rate for all
of our streaming simulations; ${\rm\gamma_{eff}}=1/2$ (eq. 60) is the better
predictor of the simulated mass-loss rates. Middle: Ratio of the analytically
predicted terminal velocity in the limit of rapid streaming, for both the
canonical ${\rm\gamma_{eff}}=2/3$ and our modified ${\rm\gamma_{eff}}=1/2$, to
the simulated velocity at the edge of the box for all of our streaming
simulations. The analytics somewhat overestimate the terminal speed, largely
because the acceleration in streaming solutions is slow (Fig. 6) and the flow
has often not reached the terminal velocity on the domain. Bottom: Energy flux
at the top of the domain relative to that input in CR at the base, for
streaming simulations and the diffusion simulations from Paper I. CR transport
via diffusion is energy-conserving while with streaming, the terminal energy
flux is only $\sim 20\%$ of the input energy flux because of transfer of
energy to the gas and our assumed isothermal equation of state. Dotted lines
connect simulations with streaming and diffusion to their pure streaming or
pure diffusion counterparts.

Figure 13 brings together our various numerical and analytic solutions for the
properties of CR-driven galactic winds with streaming at the Alfvén speed. In
the top panel we compare our analytic estimates of the mass-loss rate for
streaming CRs with ${\rm\gamma_{eff}}=0.5$ and $2/3$ (in particular, eqs. 60 &
23, respectively) to the numerical simulations. Perhaps most strikingly, the
high $v_{A}$ simulations (red circles and $\times$) agree very well with the
${\rm\gamma_{eff}}=1/2$ analytics over a factor of $\sim 100$ in mass-loss
rate while the standard CR wind theory is significantly in error at low mass-
loss rates (which are the higher $V_{g}/c_{i}$ simulations in the lower left
of the plot). For lower $v_{A}$, the modified ${\rm\gamma_{eff}}=1/2$
analytics overestimates $\dot{M}_{w}$ by a factor of few; this is not so
surprising since our high $v_{A}$ analytics is then somewhat less applicable.
Overall, we find that equation 60 with a numerical pre-factor of $f\sim
0.3f_{0.3}$ captures the numerical mass-loss rates reasonably well.

The middle panel of Figure 13 compares the analytic estimate of the wind
terminal velocity to the speed in the simulations at the top of the box. There
is reasonably good agreement, though the analytic speeds tend to be somewhat
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. This is more true of the streaming simulations in this
paper than the diffusion simulations in Paper I because the acceleration of
the flow is significantly slower with CR streaming than diffusion (see §5.2
below). Finally, the bottom panel of Figure 13 shows the energy flux at large
radii in our simulated winds relative to the energy flux in CRs at the base.
For the case of CR diffusion from Paper I, energy is globally conserved so
that the energy flux at large radii is essentially identical to that supplied
by the CRs at the base of the wind. For streaming on the other hand, energy is
lost because of the $v_{A}dp_{c}/dr$ work done by the CRs on the (isothermal)
gas: the net energy flux at large radii is $\sim 0.1-0.3$ of that supplied by
the CRs at the base of the wind, roughly consistent with equation 68.

###  5.2 Comparison of Winds Driven by CRs in the Streaming and Diffusive
Regimes

In this section we briefly compare some of the properties of galactic winds
driven by streaming CRs (this paper) with those produced by diffusing CRs
(Paper I). The top panel of Figure 14 shows the density, velocity, and CR
pressure profiles for our $V_{g}=10c_{i}$, $p_{c,0}=\rho_{0}c_{i}^{2}$
simulations with streaming transport ($v_{A,0}=c_{i}$) and diffusive transport
($\kappa=3.3r_{0}c_{i}$). The latter was discussed in detail in Paper I and is
included here in Table LABEL:tab:compare as well. Note that the mass-loss rate
is a factor of $\simeq 60$ times larger for the diffusion solution than for
the streaming solution.

The key difference between the streaming and diffusion solutions is the CR
pressure profile. With diffusion, the CR scale-height is large and the
solution approaches $p_{c}\propto 1/r$ in the limit of rapid diffusion. By
contrast, with CR streaming and $v_{A}\gg v$,
$p_{c}\propto\rho^{{\rm\gamma_{eff}}}$ with ${\rm\gamma_{eff}}=2/3$ for
canonical CR-wind theory (eq. 13) and ${\rm\gamma_{eff}}\sim 0.5$ for our
modified theory accounting for strong shocks due to instabilities in the flow
(Fig. 5). In either of the latter two regimes, the CR pressure scale-height is
thus of order the gas scale-height. In general, this implies a smaller CR
scale-height for CR streaming than for CR diffusion. The larger CR pressure
(gradient) in the presence of CR diffusion leads to the other key features of
the solution shown in Figure 14: winds with CR diffusion in general have much
larger mass-loss rates, larger terminal speeds, and retain a larger fraction
of the CR energy supplied at the base of the wind (a factor of $\simeq 60$,
2.5, and 6 larger, respectively, for the comparison in Figure 14).

To quantify the properties of winds with CR diffusion vs. streaming over a
wide parameter regime, we turn to our analytic estimates validated by
simulations. Figure 6 of Paper I shows that for $\kappa\gtrsim r_{0}c_{i}$,
winds driven by CR diffusion have mass-loss rates

$\dot{M}_{\rm diff}\simeq\frac{2\pi r_{0}^{2}\,p_{c,0}c_{i}}{h_{c}{\rm
V^{2}_{g,eff}}}$ (71)

where

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

is the base CR scale-height in units of the base radius of the flow.

The bottom panel of Figure 14 compares the diffusion and streaming estimates
of the mass-loss rate for a wide range of base CR pressures $p_{c,0}$ and
strength of gravity relative to the isothermal gas sound speed $V_{g}/c_{i}$.
For streaming we estimate the mass-loss rate by numerically solving for the
location of the sonic point as described below equation 49. We multiply the
resulting $\dot{M}_{\rm stream}$ by a factor of $0.3$ motivated by Fig. 13
(top panel). The key result of Figure 14 (bottom panel) is that mass-loss
rates due to diffusing CRs are a factor of a few-1000 larger than those
produced by streaming CRs for the same base conditions. It is worth stressing
that this is true even though we have taken into account that our simulations
(and the modified analytics in §4) have larger mass-loss rates driven by
streaming CRs than canonical CR wind theory (i.e., the diffusion vs. streaming
mass-loss rates would have been yet more different using the standard theory
of §2).

A useful analytic expression for the relative mass-loss rates in the diffusive
and streaming limits can be derived comparing equations 63 & 71, which yields

$\frac{\dot{M}_{\rm diff}}{\dot{M}_{\rm
str}}\simeq\frac{0.5}{f_{0.3}}\frac{V_{g}}{c_{\rm eff,0}}\frac{{\rm
Max}\left(\frac{c_{i}}{c_{\rm eff,0}},\frac{V_{g}}{c_{\rm
eff,0}}\sqrt{\frac{r_{0}c_{i}}{\kappa}}\right)}{{\rm
Max}\left[1,\left(\frac{\rho_{0}c_{i}^{2}}{p_{c,0}}\right)^{2c_{i}^{2}/V_{g}^{2}}\right]}$
(73)

where $0.3f_{0.3}$ is the dimensionless number we multiply equation 63 by,
motivated by Figure 13. Equation 73 reproduces the results of Figure 14 to
better than a factor of 2, with differences of this magnitude only occurring
for the smallest $V_{g}/c_{i}$.

Figure 14: Top: Comparison of the time-averaged velocity, CR pressure, and gas
density profiles for streaming ($v_{A,0}=c_{i}$) and diffusion
($\kappa=3.3r_{0}c_{i}$) solutions for CR-driven galactic winds with
$V_{g}=10c_{i}$ and $p_{c,0}=\rho_{0}c_{i}^{2}$. With CR diffusion, the CR
pressure declines significantly more slowly near the base of the wind,
increasing the mass-loss rate and causing the flow to accelerate much more
rapidly than with streaming transport. Bottom: Mass-loss rate predicted by CR
diffusion relative to CR streaming as a function of the base CR pressure
$p_{c,0}$ and the depth of the gravitational potential $V_{g}$; mass-loss
rates are analytic estimates validated by our simulations (see §5.2 for
details). For most parameters, diffusive transport predicts a mass-loss rate
significantly larger than streaming transport. This is a consequence of the
larger CR pressure (and pressure gradient) in the subsonic portion of the wind
at small radii (top panel).

###  5.3 Validity of the Isothermal Gas Approximation

Figure 15: Ratio of the time-averaged gas heating rate due to CR streaming
$q_{A}$ to the optically thin cooling rate ${q_{\rm cool}}$ for our solutions
with $v_{A,0}=c_{i}$ and $p_{c,0}=\rho_{0}c_{i}^{2}$. Cooling is calculated
post facto and is not included in the simulations. Self-consistency of our
isothermal solutions requires rapid cooling, i.e., ${q_{\rm cool}}>q_{A}$. We
assign physical units to our simulations by assuming $v_{A,0}=10\,{\rm
km\,\,s}^{-1}$, $p_{c,0}=10^{-12}$ erg cm-3, $r_{0}=1\,{\rm kpc}$, $n_{0}=1$
cm-3, and $\Lambda=10^{-21}$ erg s-1 cm3; for other parameter values the
curves in this Figure can be scaled up or down using equation 74.

In his original treatment of galactic winds driven by streaming CRs, Ipavich
(1975) separately studied solutions that included Alfvénic heating of the gas
($\sim v_{A}dp_{c}/dr$), but neglected radiative losses, and "zero-
temperature" solutions in which the heating of the gas is assumed to be
rapidly radiated away. The subsequent literature has considered variants of
one or both of these limits. For example, Breitschwerdt et al. (1991)
neglected Alfvénic heating of the gas but also evolved the gas energy equation
without radiative cooling.888 Breitschwerdt et al. (1991) effectively assumed
that energy is transferred from CRs to Alfvén waves by the streaming
instability, but that the waves are not damped, and hence that there is no
heating of the gas. However, all physical damping mechanisms for Alfvén waves
excited by the streaming instability (e.g., ion-neutral, turbulent damping,
etc.; see Squire et al. 2021 for a recent discussion) lead to the Alfvén waves
rapidly heating the gas, so we view this model variant as relatively
unphysical. Mao & Ostriker (2018) assumed an isothermal equation of state, as
we have in this paper (our analysis in §2 is most analogous to their work).

In this section we provide estimates of when the isothermal gas approximation
is plausibly self-consistent in the sense that the cooling time of the gas in
the wind is shorter than the heating timescale associated with energy transfer
from the CRs to the gas, mediated by the Alfvén waves. This is a necessary but
not sufficient condition for the plausibility of the isothermal gas
approximation – necessary because the Alfvénic heating term is necessarily
present in models of CR transport mediated by the streaming instability, but
not sufficient because it focuses only on Alfvénic heating neglecting other
heating mechanisms (e.g., shocks).

The optically thin cooling rate per unit volume given a cooling function
$\Lambda(T)$ is ${q_{\rm cool}}=n^{2}\Lambda(T)$ while the Alfvénic heating
rate per unit volume is given by $q_{A}=v_{A}|dp_{c}/dr|$. The ratio of the
heating and cooling rates can thus be written as

$\frac{q_{A}}{{q_{\rm
cool}}}=\left(\frac{v_{A,0}p_{c,0}}{n_{0}^{2}r_{0}\Lambda}\right)\left(\frac{\rho_{0}}{\rho}\right)^{2}\frac{v_{A}}{v_{A,0}}\frac{r_{0}}{p_{c,0}}\bigg{|}\frac{dp_{c}}{dr}\bigg{|}$
(74)

where we have written quantities in terms of dimensionless variables such that
the pre-factor in front is roughly the ratio of the heating time to the
cooling time at the base of the wind.

To evaluate equation 74 in our simulations, we need to assume physical values
for the quantities in the first parentheses. We take $v_{A,0}=10\,{\rm
km\,\,s}^{-1}$, $p_{c,0}=10^{-12}$ erg cm-3, $r_{0}=1\,{\rm kpc}$, $n_{0}=1$
cm-3, and $\Lambda=10^{-21}$ erg s-1 cm3, in which case the first term in
parentheses in equation 74 is $\simeq 3\times 10^{-7}$. For these parameters,
Figure 15 shows the ratio of the time-averaged heating and cooling rates
$q_{A}/{q_{\rm cool}}$ for our $v_{A,0}=c_{i}$, $p_{c,0}=\rho_{0}c_{i}^{2}$,
and $V_{g}/c_{i}=3$, 6, $10$ simulations. In evaluating $q_{A}$ and ${q_{\rm
cool}}$, we explicitly calculate $\langle v_{A}dp_{c}/dr\rangle$ and
$\langle\rho^{2}\rangle$ to best capture the volume averaged heating and
cooling rates in the inhomogeneous flow, though this only makes a factor of
few difference for the results in Figure 15. The curves in Figure 15 can be
trivially scaled to other choices of $c_{i}$, $r_{0}$, $\Lambda$, & $n_{0}$
using equation 74. For example, a higher base density of $n_{0}\simeq 10-1000$
cm-3 and a somewhat smaller physical scale of the system, $r_{0}\sim 0.1-0.3$
kpc, would be appropriate for starbursts, which would decrease $q_{A}/{q_{\rm
cool}}$ in Figure 15 by a factor of few-100s.

The sonic point in our analytic models is $r_{s}\simeq 2.7r_{0}$ (eq. 49) and
indeed this is roughly the radius at which the flow reaches $v\sim V_{g}\sim
v_{\infty}$ in the simulations (Fig. 6). As a result, the isothermal
approximation is well-motivated if $q_{A}\lesssim{q_{\rm cool}}$ out to
$r/r_{0}\sim$ few, so that cooling can regulate the temperature and maintain a
roughly constant value of $c_{i}$ at the radii within which the mass-loss rate
and wind velocity are determined.999In our isothermal solutions, gas pressure
is only important at much smaller radii than $r\sim r_{s}\sim$ few $r_{0}$,
namely $r\lesssim r_{tr}$ (eq. 25). One might thus be tempted to evaluate the
importance of cooling only at those radii. This is not correct, however,
because if $q_{A}/{q_{\rm cool}}\gtrsim 1$ at $r_{tr}\lesssim r\lesssim
r_{s}$, the gas will heat up appreciably, significantly increasing the
importance of gas pressure at radii interior to the sonic point in the
isothermal solutions, thus fundamentally modifying the structure of the flow.
A second comparison that is important is the ratio of the heating/cooling
timescale and the flow time in the wind $\sim H/v$. For our CR-driven wind
solutions, this ratio is $t_{\rm adv}/t_{\rm
heat}\simeq(v_{A}/v)(p_{c}/p_{gas})$. Since $v\lesssim v_{A}$ and
$p_{c}\gtrsim p_{gas}$ near the sonic point, $t_{\rm heat}\lesssim t_{\rm
adv}$: the short heating timescale implies that it is indeed the comparison of
the heating and cooling rates (captured by equation 74 and Fig. 15) that sets
the thermodynamics of the gas.

Figure 15 shows that near the base of the wind ${q_{\rm cool}}\gg q_{A}$ for
all of our models so that cooling is important at small radii.101010The
exception to this is if the gas is very hot to begin with, e.g., if
$c_{i}\simeq 300\,{\rm km\,\,s}^{-1}$, $\Lambda\simeq 10^{-22.5}$ erg s-1 cm3,
in which case the curves in Figure 15 would all be scaled up by $\sim 10^{6}$.
This is essentially a thermally-driven wind in which CRs do not play a
dominant role. By $r\sim 2r_{0}$, however the importance of cooling depends
strongly on the depth of the gravitational potential $V_{g}$, with cooling
dominant for $V_{g}=3c_{i}$, marginally important for $V_{g}=6c_{i}$, and
negligible $V_{g}=10c_{i}$. This implies that the $V_{g}=10c_{i}$ isothermal
solutions are only plausibly self-consistent if the pre-factor in equation 74
is significantly smaller than assumed in Figure 15, e.g., if the gas density
at the base is $n_{0}\gg 1$ cm-3, as might be the case in a starburst or at
higher redshift. For $V_{g}=3c_{i}$, the mass-loss rate is sufficiently larger
that cooling of the gas is dynamically important at all radii and so the
isothermal approximation is reasonable. $V_{g}=6c_{i}$ lies in between the
$V_{g}=3$ & $10c_{i}$ cases, with comparable heating and cooling rates for
$r\sim 2-5r_{0}$.

The strong dependence of $q_{A}/{q_{\rm cool}}$ on $V_{g}$ in Figure 15
reflects the strong dependence of the mass-loss rate on $V_{g}$ (e.g., Fig.
12). These results can be understood analytically as follows. We evaluate the
importance of cooling at the analytic sonic point using our modified streaming
models with general ${\rm\gamma_{eff}}$. The density at the sonic point is
then given by equation 51 and the location of the sonic point by equation 59.
Assuming a split monopole magnetic field and taking $|d\ln p/d\ln
r|_{r_{s}}\sim 1$, we find

$\begin{split}&\frac{q_{A}}{{q_{\rm
cool}}}\bigg{|}_{r_{s}}\simeq\frac{v_{A,0}p_{c,0}}{n_{0}^{2}r_{0}\Lambda}\,e^{-3/(2-2{\rm\gamma_{eff}})}\,\left(\frac{c_{\rm
eff,0}}{{\rm
V_{g,eff}}}\right)^{-\frac{2(2.5-{\rm\gamma_{eff}})}{1-{\rm\gamma_{eff}}}}\\\
&\simeq\frac{5}{n_{0}\Lambda_{-21}}\,\frac{v_{A,0}\,c_{\rm
eff,0}^{2}}{(10\,{\rm km\,\,s}^{-1})^{3}}\frac{\rm
kpc}{r_{0}}\,\left(\frac{c_{\rm eff,0}}{0.1\,{\rm V_{g,eff}}}\right)^{-8}\
({\rm\gamma_{eff}}=1/2)\end{split}$ (75)

where $\Lambda=10^{-21}\Lambda_{-21}$ erg cm3 s-1 and we have assumed
${\rm\gamma_{eff}}=1/2$ in the last equality. Recall that optically thin
atomic cooling of solar metallicity gas has $\Lambda\sim 10^{-22}-10^{-21}$
erg s-1 cm3 for $T\sim 10^{4-6}$K, with slower cooling (lower $\Lambda$) at
yet higher T. For comparison, for ${\rm\gamma_{eff}}=2/3$ the result is
similar to equation 75 but $\propto(c_{\rm eff,0}/{\rm V_{g,eff}})^{-11}$.

Equation 75 describes well the numerical results in Figure 15, both the
normalization of $q_{A}/{q_{\rm cool}}$ at $r\sim{\rm few}\,r_{0}$ and the
strong dependence on $V_{g}$. Overall, we conclude that the validity of the
isothermal gas approximation is determined primarily by the ratio of the base
CR sound speed in the wind to the galaxy escape speed, via the strong
dependence of $q_{A}/{q_{\rm cool}}$ on $c_{\rm eff,0}/V_{g}$ in equation 75.
Indeed, the constraint $q_{c}\gtrsim q_{A}$ at $r\sim r_{s}$ requires

$\frac{c_{\rm eff,0}}{{\rm V_{g,eff}}}\gtrsim
0.13\left(\frac{v_{A,0}}{10\,{\rm km\,\,s}^{-1}}\frac{\rm
kpc}{r_{0}}n_{0}^{-1}\Lambda_{-21}^{-1}\right)^{1/6}\left(\frac{{\rm
V_{g,eff}}}{100\,{\rm km\,\,s}^{-1}}\right)^{1/3}.$ (76)

As noted in §2 & 4, the mass-loss rate increases significantly with increasing
$c_{\rm eff,0}$ such that warmer phases of the ISM, which likely have larger
$c_{\rm eff,0}$, will probably dominate the mass-loss rate. Equation 76 shows
that precisely because of the larger mass-loss rate, these phases will also
have stronger cooling and are more likely to be in the (roughly) isothermal
limit studied in this paper. That being said, many of the isothermal CR-
streaming-driven solutions in the literature with lower values of $c_{\rm
eff,0}$, including several in this paper, are not self-consistent because CR
heating of the gas is likely to overwhelm radiative cooling, invalidating the
isothermal gas approximation.

One important caveat to our analysis in this section is that we have compared
the time-averaged CR heating and radiative cooling rates in the wind. Figure 2
shows, however, that CR heating is negligible over most of the volume because
$dp_{c}/dr\simeq 0$. Cooling will still be important in those regions, driving
the gas towards photoionization equilibrium. Likewise, CR heating and
radiative cooling are both enhanced in the vicinity of the strong shocks that
permeate the flow (Fig. 4). This is not accounted for in our analysis of time-
averaged heating and cooling rates in Figure 15. Time-dependent simulations
with both streaming heating and cooling are needed to fully understand when
(and, indeed, if) the isothermal approximation used here and in previous work
is appropriate.

###  5.4 Observational Implications

In Paper I, we highlighted the importance of empirically constraining CR
transport models in other galaxies given the theoretical uncertainties in the
mechanism(s) of CR transport. We argued that GeV gamma-ray observations by
FERMI (Abdo et al., 2010b, a; Ackermann et al., 2012), though modest in
number, are particularly valuable because pion decay directly constrains the
energetics and transport of the GeV protons that dominate the CR energy
density (e.g., Lacki et al. 2010; Lacki et al. 2011; Crocker et al. 2021). We
defer to a future paper detailed models of streaming-driven galactic winds
calibrated to gamma-ray observations. It is worth stressing, however, that the
gamma-ray calibrated diffusion coefficients in Paper I of $\kappa\sim 10^{29}$
cm2 s-1 correspond to CR escape times of $H_{c}^{2}/\kappa\sim 30$ Myrs for
$H_{c}\sim 3$ kpc. For streaming transport to produce a similar escape
timescale with $H_{c}\sim 3$ kpc, we require an Alfvén speed of $v_{A,0}\sim
100\,{\rm km\,\,s}^{-1}$. This requires strong magnetic fields in the warm
ionized or hot phases of the ISM.

Our general expression for the mass-loss rate validated by the simulations
(eq. 60) can be written as $\dot{M}_{w}\propto p_{c,0}\,c_{\rm
eff,0}^{2}/V_{g}^{4}$, where $p_{c,0}$ is the base CR pressure in the wind.
This shows that the mass-loss rate will in general be dominated by the hotter
ISM phases, which have larger values of $c_{\rm eff,0}$; see also Fig. 12.
Importantly, streaming CRs are thus ineffective at directly driving cold gas
out of galaxies, though CR-driven winds in hotter phases of the ISM may
entrain cooler gas along the way (e.g., Wiener et al. 2019; Brüggen &
Scannapieco 2020). The mass-loss rate in the warm-hot ISM phases can in
principle approach the maximum mass-loss rate allowed by energy conservation,
$\dot{M}_{\rm max}/\dot{M}_{*}\simeq\,2\,\epsilon_{c,-6.3}\,(100\,\,{\rm
km\,\,s}^{-1}/v_{\rm esc}(r_{0}))^{2}$ (eq. 70). This maximum mass-loss rate
is, however, itself already modest compared to the mass-loss rates typically
required in cosmological models to reconcile the stellar mass and dark matter
halo mass-functions (e.g., Somerville & Davé 2015). The reason is that
isothermal winds with CR streaming can only supply at most $\sim 2.5$% of the
energy of SNe to galactic winds: $\sim 10$% of the SNe energy goes into CRs
but the CRs only retain $\sim 1/4$ of this energy in the wind when cooling is
efficient (Fig. 13). This suggests that streaming CRs are an energetically
inefficient source of feedback for driving galactic winds. Pionic losses,
which are likely important in luminous starbursts given their gamma-ray
luminosities (e.g., Lacki et al. 2010), would strengthen this conclusion by
further decreasing the CR energy available to drive galactic winds.

Our simulations show that the terminal speed in CR streaming-driven galactic
winds is typically within a factor of few of the galaxy escape speed set by
$V_{g}$ (see Table LABEL:tab:compare). This property of our wind models is
intriguingly similar to some observational inferences (e.g., Martin 2005;
Weiner et al. 2009) and to the class of ‘momentum-conserving’ wind models
(Murray et al., 2005) that have been successful at reproducing the galaxy
stellar mass function and the galaxy mass-metallicity relationship (e.g.,
Oppenheimer & Davé 2006; Oppenheimer et al. 2010). However, as noted in the
previous paragraph, the somewhat low (maximum) mass-loss rates in our models
are in tension with what is required in cosmological models of galaxy
formation. The wind momentum fluxes we find are similar to those expected from
supernova-heated hot gas wind models and those driven by radiation pressure
(eq. 69; Lochhaas et al. 2020).

One of the most important new results presented in this paper vis-a-vis
observations of galactic winds is the presence of strong shocks that permeate
the flow (Fig. 2 & Fig. 7). These shocks significantly increase the mass-loss
rate relative to previous steady state calculations (Fig. 12). It is likely
that such shocks would also produce a wide range of observational signatures.
In particular, shocks are likely to give rise to multi-phase gas with hot
post-shock gas cooling rapidly and producing X-ray emission, optical line
emission, and even (re-)forming molecules and dust in the post shock region if
there is ambient dense gas in the flow (Hollenbach & McKee, 1989). Strong in-
situ shocks could also accelerate electrons producing synchrotron emission off
the plane of the galactic disk. This could be important for interpreting
gradients in synchrotron emission and synchrotron spectral indices above
galaxy disks and/or in galaxy outflows (e.g., Dumke et al. 1995; Buckman et
al. 2020), since models typically focus on the interplay between cooling and
escape (e.g., Seaquist & Odegard 1991), rather than re-acceleration. Future
work extending our results to models with heating, cooling, and magnetic
fields would be able to directly predict the observational signatures of the
strong shocks found here, for comparison to observations.

##  6 Summary

In this paper we have used analytic estimates and idealized time-dependent
spherically symmetric simulations to study the properties of isothermal
galactic winds driven by CRs streaming at the Alfvén velocity. This
complements Paper I (Quataert et al., 2021) in which we carried out a similar
analysis for diffusive CR transport. These two mechanisms of CR transport
differ dramatically in their predictions for how the CR pressure decreases
away from a galaxy (see the top panel of Fig. 14): in the limit of rapid CR
streaming, $p_{c}\propto\rho^{2/3}$ according to standard CR-driven wind
theory (eq. 12) and so the CR pressure scale-height is tied to that of the
gas, while for rapid CR diffusion, $p_{c}\propto r^{-1}$, i.e., the CR
pressure scale-height is of order the size of the system. One aim of this
paper and its companion is to understand how these differences in the CR
pressure profiles translate into differences in the predicted properties of
galactic winds.

To start, we analytically estimated galactic wind properties using standard
theoretical assumptions for CR streaming (§2). In particular, when the Alfvén
velocity is large compared to the flow velocity, the standard treatment of CR
streaming implies $p_{c}\propto\rho^{2/3}$ (Ipavich 1975; see eq. 12). This
effective equation of state applies at small radii in the wind, before the
flow has accelerated to large velocities $\gtrsim v_{A}$. The assumption
$p_{c}\propto\rho^{2/3}$ is sufficient to allow analytic estimates of the wind
properties: the mass-loss rate $\dot{M}_{w}$ (eq. 23), terminal speed
$v_{\infty}$ (eq. 35), asymptotic wind power $0.5\dot{M}_{w}v_{\infty}^{2}$
(eq. 36), and wind momentum loss-rate $\dot{p}_{w}$ (eq. 37).

To test these models, we carried out time-dependent simulations for the same
idealized spherically-symmetric galactic wind problem that we studied
analytically (isothermal gas and an isothermal gravitational potential). The
simulations use the two-moment numerical method for CR transport developed by
Jiang & Oh (2018). The two-moment method is particularly important for this
study because we require an accurate treatment of CR streaming. For our
spherically-symmetric problem, the magnetic field is not evolved dynamically
but is only present via the Alfvén velocity, which sets the CR streaming
speed.

For most of the parameter space we surveyed (Table LABEL:tab:compare), the
simulations yield solutions that are intrinsically time-dependent,
invalidating previous steady-state wind calculations. Small amplitude waves
grow near the base of the wind and then steepen into strong shocks that
permeate the flow (Fig. 2-4). To interpret the instabilities found in the
simulations, we identified two new linear instabilities driven by CR streaming
(§3.4.1 and Appendix A). The first is an acoustic instability at high Alfvén
speed due to the finite background CR pressure gradient present in hydrostatic
equilibrium; this instability, which dominates in most of our simulations, is
a streaming analogue of the acoustic instability discussed by Drury & Falle
(1986). The second instability is a short-wavelength instability of the CR
entropy mode due to the finite speed of light retained in the two-moment CR
equations (eqs. 42). We note that the sound wave instability driven by CR
streaming discovered by Begelman & Zweibel (1994) is not present in our
simulations because of our assumption of isothermal gas, but it may be present
in winds where cooling is inefficient (e.g., Huang & Davis 2021).

The shocks seeded by the streaming-generated linear instabilities have at
least two particularly important dynamical consequences. (1) As the waves grow
in amplitude they eventually reverse the sign of the local pressure gradient
$dp_{c}/dr$. Since CRs cannot stream up the pressure gradient, this leads to
the CR pressure profile having a ’staircase’ like structure (Fig. 2): most of
the volume is filled with $dp_{c}/dr=0$ while the CR pressure changes
appreciably only in small regions near the shocks. (2) The latter implies that
the CR pressure changes appreciably only in denser than average regions and
thus regions of smaller than average Alfvén velocity (Fig. 4). This leads to
the streaming-induced transfer of energy from the CRs to the gas $|{\bf
v_{A}\cdot\nabla}p_{c}|$ being less efficient than would be the case in a
laminar solution. This is particularly true for lower gas sound speeds $c_{i}$
(i.e., larger $V_{g}/c_{i}$ where the circular velocity of the potential
$\sqrt{2}V_{g}$ sets the velocity scale for shocks), because shock-induced
density compressions are then larger. The net result is that the time-averaged
relation between CR pressure and gas density can be shallower than the
canonical $p_{c}\propto\rho^{2/3}$, with $\langle
p_{c}\rangle\propto\langle\rho\rangle^{{\rm\gamma_{eff}}}$ and
${\rm\gamma_{eff}}\simeq 1/2$ rather than ${\rm\gamma_{eff}}=2/3$ (Figs. 5 &
8). The scaling $p_{c}\propto\rho^{1/2}$ is equivalent to
${\bf\nabla}\cdot(p_{c}{\bf v_{A})}\sim$ constant, i.e., the CR power is
roughly independent of radius. This reflects the fact that the inhomogeneous
time-dependent solutions minimize the energetic importance of streaming losses
which otherwise sap energy out of the CRs. It is worth noting that our fluid
model for CRs neglects the fact that CRs can be efficiently accelerated at
shocks, which would produce an additional source of CRs in the vicinity of
shocks. This could modify how CRs behave in shocks from what we find in our
calculations (Fig. 4).

Motivated by the simulation results we carried out analytic estimates of the
wind mass-loss rate, kinetic power, terminal velocity, and momentum loss-rate
for general ${\rm\gamma_{eff}}$ (§4). The case of ${\rm\gamma_{eff}}\simeq
1/2$ does a good job of reproducing our simulation results (see §5.1 and Fig.
13), including the density profiles in the inner hydrostatic portion of the