The physics of galactic winds driven by cosmic rays I: Diffusion

Quataert, Eliot; Thompson, Todd A; Jiang (姜燕飞), Yan-Fei

doi:10.1093/mnras/stab3273

ABSTRACT

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

cosmic rays, galaxies: dwarf, galaxies: starburst

1 INTRODUCTION

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

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

Many treatments of CR-driven winds exist in the literature, starting with the foundational work of Ipavich (1975) who computed time-steady solutions for CR-driven winds assuming streaming transport at the Alfvén velocity. More detailed models including CR streaming and hydromagnetic wave pressure were computed by Breitschwerdt, McKenzie & Voelk (1991), deriving a mass-loss rate of the order of M_⊙ yr⁻¹ from the Galaxy. Everett et al. (2008) further developed a hybrid thermal and CR-driven wind model for the Galaxy while Recchia, Blasi & Morlino (2016) solved for the Milky Way’s CR spectrum and wind dynamics self-consistently using the CR kinetic equations. More generally, on the basis of momentum conservation, Socrates, Davis & Ramirez-Ruiz (2008) argued that CRs could drive significant mass fluxes from star-forming galaxies.

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

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

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

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

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

2 ANALYTIC APPROXIMATIONS FOR COSMIC RAY-DRIVEN WINDS WITH DIFFUSION

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

The general equation of motion, including gas pressure p, CR pressure p_c, a general gravitational potential Φ, and magnetic forces in the limit of magnetohydrodynamics can be written as

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

(1)

The energy equation for the CRs in the absence of CR sources and pionic losses, and including diffusion and streaming at the Alfvén velocity down the CR pressure gradient, i.e. v_s = −v_A|∇p_c|/∇p_c with v_A = B/(4πρ)^1/2, can be written as

$$\begin{eqnarray*} \frac{\partial E_c}{\partial t}+\nabla \cdot {\bf F}_c=\left({\bf v}+{\bf v_{\rm s}}\right)\cdot \nabla p_{c}, \end{eqnarray*}$$

(2)

where the ‘equilibrium’ CR flux¹ is

$$\begin{eqnarray*} {\bf F}_c=4 p_c \left({\bf v}+{\bf v_{\rm s}}\right)-\kappa \, {\bf n}\left({\bf n}\cdot \nabla E_c\right), \end{eqnarray*}$$

(3)

E_c is the CR energy density, p_c = E_c/3, κ is the diffusion coefficient, and |${\bf n}={\bf v_{\rm A}}/|\bf v_{\rm A}|$|⁠. Equation (2), and indeed the scalar CR pressure in equation (1), is formally valid only on scales larger than the mean-free-path of the ∼GeV energy CRs that dominate the total energy of the CR population (e.g. Skilling 1971).²

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

With these approximations, the steady state equations of motion are

$$\begin{eqnarray*} \dot{M}_w = 4 \pi r^2 \rho v = {\rm const}, \end{eqnarray*}$$

(4)

$$\begin{eqnarray*} v \frac{\mathrm{ d}v}{\mathrm{ d}r} = -\frac{1}{\rho }\frac{\mathrm{ d}p}{\mathrm{ d}r} - \frac{1}{\rho }\frac{\mathrm{ d}p_c}{\mathrm{ d}r} - \frac{2 V_g^2}{r}, \end{eqnarray*}$$

(5)

and

$$\begin{eqnarray*} v\frac{\mathrm{ d}p_c}{\mathrm{ d}r} = -\frac{4 p_c}{3 r^2} \frac{\mathrm{ d} r^2 v}{\mathrm{ d}r} + \frac{1}{r^2}\frac{\rm d}{\mathrm{ d}r}\left(r^2 \kappa \frac{\mathrm{ d} p_c}{\mathrm{ d}r}\right). \end{eqnarray*}$$

(6)

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

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

(7)

which suggests that in the limit of rapid diffusion, and near the base of the outflow where the gas velocity is small, the advective and adiabatic terms in the CR energy equation can simply be neglected. For reference, the diffusion coefficient in the Milky Way is estimated to be ∼10²⁹ cm² s⁻¹ (e.g. Trotta et al. 2011), although there is significant uncertainty in this estimate because of a degeneracy between the diffusion coefficient and the size of the CR propagation region (the ‘halo;’ Linden, Profumo & Anderson 2010; Trotta et al. 2011); this degeneracy can be partially broken using observations of radioactive secondaries such as ¹⁰Be (Evoli et al. 2020). In Section 4, we return to constraints on the diffusion coefficient in galaxies using observations of gamma-ray emission from pion decay. For now, we proceed under the assumption that for sufficiently large diffusion coefficients κ, one can neglect the advective and adiabatic terms in equation (6) and focus solely on the diffusion term, for which the steady state solution is

$$\begin{eqnarray*} p_c = p_{c,0} + \frac{\dot{E}_c}{12 \pi \kappa } \left(\frac{1}{r} - \frac{1}{r_0}\right) , \end{eqnarray*}$$

(8)

where p_{c, 0} is the base CR pressure at radius r₀ and

$$\begin{eqnarray*} \dot{E}_c = -12 \pi r^2 \kappa \frac{\mathrm{ d}p_c}{\mathrm{ d}r} = -12 \pi r_0^2 \kappa \frac{\mathrm{ d}p_c}{\mathrm{ d}r}\bigg |_{r_0} \end{eqnarray*}$$

(9)

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

Throughout this paper, we interpret the base radius r₀ that appears in equations (8) and (9) to be roughly the size of the star-forming region within which most of the cosmic-rays are produced. In a nuclear starburst, this radius can be significantly smaller than the total size of the galactic disc. We will also make frequent use of the CR scale height near the base of the wind:

$$\begin{eqnarray*} H_c = \left(\frac{- \mathrm{ d} \ln p_c}{\mathrm{ d} r}\bigg |_{r_0}\right)^{-1} \equiv r_0 \, h_c , \end{eqnarray*}$$

(10)

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

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

$$\begin{eqnarray*} \frac{1}{v}\frac{\mathrm{ d}v}{\mathrm{ d}r}\left(v^2 - c_i^2 \right) = \frac{2 c_i^2}{r} -\frac{1}{\rho }\frac{\mathrm{ d}p_c}{\mathrm{ d}r} - \frac{2 V_g^2}{r}. \end{eqnarray*}$$

(11)

With equation (8), i.e. assuming |$\dot{E}_c$| is a constant because diffusion is rapid, the wind equation becomes

$$\begin{eqnarray*} \frac{1}{v}\frac{\mathrm{ d}v}{\mathrm{ d}r}\left(v^2 - c_i^2 \right) = \frac{2 c_i^2}{r} + \frac{\dot{E}_c}{12 \pi r^2 \rho \kappa } - \frac{2 V_g^2}{r}. \end{eqnarray*}$$

(12)

In the portion of the flow where v < c_i, an approximation to the density profile can be derived by assuming hydrostatic equilibrium so that

$$\begin{eqnarray*} \frac{c_i^2}{\rho }\frac{\mathrm{ d}\rho }{\mathrm{ d}r} - \frac{\dot{E}_c}{12 \pi r^2 \rho \kappa } = - \frac{2 V_g^2}{r}. \end{eqnarray*}$$

(13)

Equation (13) has the solution

$$\begin{eqnarray*} \rho (r) = \rho _0 \left[x^{-\xi } + A(x^{-1} - x^{-\xi })\right] , \end{eqnarray*}$$

(14)

where x = r/r₀ and ξ and A are the two key dimensionless parameters in this problem, which measure the strength of gravity relative to the gas and CR pressures at the outflow’s base:

$$\begin{eqnarray*} \xi \equiv \frac{2 V_g^2}{c_i^2} \end{eqnarray*}$$

(15)

and

$$\begin{eqnarray*} A & \equiv& \frac{\dot{E}_c}{12 \pi \kappa c_i^2 \rho _0 r_0 (\xi - 1)} \simeq \frac{\dot{E}_c}{24 \pi \kappa V_g^2 \rho _0 r_0} \nonumber\\ &\simeq& \frac{p_{c,0}}{2 h_c \rho _0 V_g^2} = \frac{1}{2 h_c}\frac{c_{c,0}^2}{V_g^2}\lesssim 1, \end{eqnarray*}$$

(16)

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

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

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

$$\begin{eqnarray*} v(r_s) = c_i \ \ \ {\rm and} \ \ \ \frac{r_s \, \rho (r_s)}{r_0 \rho _0} = \frac{c_{c,0}^2}{2 h_c {\rm V^2_{g,eff}}}, \end{eqnarray*}$$

(17)

where

$$\begin{eqnarray*} {\rm V^2_{g,eff}}=V_g^2-c_i^2. \end{eqnarray*}$$

(18)

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

$$\begin{eqnarray*} \dot{M}_w = 4 \pi r_s^2 \rho (r_s) c_i = \frac{2 \pi r_s r_0 \rho _0 c_i}{h_c} \frac{c_{c,0}^2}{{\rm V^2_{g,eff}}}. \end{eqnarray*}$$

(19)

To estimate the radius of the sonic point r_s, one can numerically solving the second critical point equation in equation (17) using the analytic density profile in equation (13), i.e.

$$\begin{eqnarray*} \frac{c_{c,0}^2}{2 h_c \left(V_g^2-c_i^2\right)} \simeq x_s^{-\xi +1} + \frac{c_{c,0}^2}{h_c \left(2V_g^2-c_i^2\right)}\left(1 - x_s^{-\xi +1}\right). \end{eqnarray*}$$

(20)

A good approximate solution to equation (20) can be found by expanding for |$c_i^2 \ll 2 V_g^2$| which yields

$$\begin{eqnarray*} r_{s} \simeq r_0 \, \left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/(2V_g^2)} . \end{eqnarray*}$$

(21)

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

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

$$\begin{eqnarray*} \dot{M_w}\simeq \frac{2 \pi r_0^2\rho _0 c_i}{h_c} \, \frac{c_{c,0}^2}{{\rm V^2_{g,eff}}} \left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/(2V_g^2)}. \end{eqnarray*}$$

(22)

In the limit that 2V_g ≫ c_i,

$$\begin{eqnarray*} \dot{M}_w \sim \frac{2 \pi r_0^2 c_i p_{c,0}}{h_c V_g^2}=\frac{2 \pi r_0^2 \rho _0 c_i}{h_c} \left(\frac{c_{c,0}}{V_g}\right)^2. \end{eqnarray*}$$

(23)

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

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

$$\begin{eqnarray*} \dot{M}_w\simeq \, \left(\frac{r_0\, c_i}{6\kappa }\right)\left(\frac{\dot{E}_c}{{\rm V^2_{g,eff}}}\right)\left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/(2V_g^2)}, \end{eqnarray*}$$

(24)

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

We now present an order of magnitude estimate of the scale-height h_c by determining the radius r_adv ≃ r₀(1 + h_c) at which the advective flux |$F_{\rm adv} = 4 p_c v = \dot{M}_w p_c/(\pi r^2 \rho)$| is comparable to the diffusive flux F_diff = −3κdp_c/dr. This occurs when the velocity of the outflow is given by

$$\begin{eqnarray*} v(r_{\rm adv}) \approx - \kappa \frac{\mathrm{ d} \ln p_c}{\mathrm{ d}r} \approx \frac{c_i}{h_c} \frac{\kappa }{c_i r_0} , \end{eqnarray*}$$

(25)

where we have dropped a factor of 3/4 consistent with the rough nature of the estimates that follow. Since h_c ≲ 1 and our analytics assumes κ ≳ c_ir₀, equation (25) implies that the transition from diffusion to advection happens exterior to the sonic point. Gas pressure is then negligible and the momentum equation becomes |$\rho v \mathrm{ d}v/\mathrm{ d}r \simeq -\mathrm{ d}p_c/\mathrm{ d}r - 2\rho V_g^2/r$|⁠. At the sonic point, the two terms on the right-hand-side of dv/dr are comparable (see equation 17). Since the gas density scale-height is smaller than the CR scale-height, somewhat exterior to the sonic point ρvdv/dr ≃ −dp_c/dr. Multiplying by 4πr² gives

$$\begin{eqnarray*} \dot{M} \frac{\mathrm{ d}v}{\mathrm{ d}r}\bigg |_{\sim r_{\rm adv}} \sim 4 \pi r_{\rm adv}^2 \frac{p_{c,0}}{h_c} \rightarrow v(r_{\rm adv}) \sim \frac{h_c V_g^2}{c_i} , \end{eqnarray*}$$

(26)

where in the second expression we have used the approximate version of |$\dot{M}$| from equation (23) and have taken r_adv ≈ r₀, consistent with h_c ≲ 1. Equating equations (25) and (26) then yields h_c ∼ (c_i/V_g)(κ/r₀c_i)^1/2. For κ → ∞ this gives h_c > >1, inconsistent with the local approximation used here; the solution should be p_c∝1/r, i.e. h_c ≃ 1, so that

$$\begin{eqnarray*} h_c \sim \min \left (1,\frac{c_i}{V_g}\sqrt{\frac{\kappa }{r_0 c_i}}\right). \end{eqnarray*}$$

(27)

As an example, if κ ∼ 10²⁹ cm² s|$^{-1} \simeq 10 \, (r_0 c_i/30 \, {\rm km \, \, s}^{-1}\, {\rm kpc})$|⁠, |$c_i \sim 10 \, \, {\rm km \, \, s}^{-1}$|⁠, and |$V_g \sim 100 \, \, {\rm km \, \, s}^{-1}$|⁠, h_c ∼ 1/3. In our analytic scalings that follow in this section, we primarily normalize h_c to a value of 1/4 motivated by these fiducial parameters, but in several plots we will use equation (27). In Section 3, we also compare equation (27) to our numerical solutions and find good agreement.

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

$$\begin{eqnarray*} \dot{M}_w &\simeq & \frac{4 \pi r_0^2 \, \rho _0\, c_i}{h_c}\, \left(\frac{ c_{c,0}^2}{2{\rm V^2_{g,eff}}}\right)\left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/\left(2V_g^2\right)}\nonumber\\ &\simeq & 0.06 \, \, {\rm M_\odot \, yr^{-1}}\left(\frac{r_0}{1 \, \, {\rm kpc }}\right)^2 \left(\frac{1/4}{h_c}\right) \left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/\left(2V_g^2\right)}\nonumber\\ && \times \left(\frac{n_0}{1 \, {\rm cm^{-3}}}\right) \left(\frac{c_i \ c^2_{c,0}}{[10 \, \, {\rm km \, \, s}^{-1}]^3}\right) \left(\frac{{\rm V_{g,eff}}}{100 \, \, {\rm km \, \, s}^{-1}}\right)^{-2} , \end{eqnarray*}$$

(28)

where n₀ = ρ₀/m_p. For reference, if we scale for parameters appropriate to a galaxy like the Milky Way (with V_g ≃ 150 km s⁻¹, c_i = c_{c, 0} ≃ 10 km s⁻¹, n₀ ≃ 1 cm⁻³, and r₀ ∼ 5 kpc), equation (28) yields |$\dot{M}_w \simeq 1 \, {\rm M_\odot \, yr^{-1}}$|⁠, comparable to the star formation rate. Equation (28) also predicts |$\dot{M}_w \propto p_{c,0} c_i/h_c$|⁠, which scales ∝p_{c, 0}c_i or |$\propto p_{c,0} c_i^{1/2}$|⁠, depending on which regime of equation (27) is appropriate (these scalings assume r_s ≃ r₀ for simplicity). To the (uncertain) extent that the CR pressure is comparable in different phases of the ISM, the outflow is thus likely to be dominated by the warmer ISM phases, though only by a factor of a few.

Fig. 1 shows the mass-loss rate (equation (28) with h_c from equation (27) taking κ = 10r₀c_i) as a function of the two key dimensionless parameters in the problem, namely the strength of gravity (V_g/c_i = (ξ/2)^1/2; see equation 15) and the base CR sound speed (⁠|$p_{c,0}/\rho _0 c_i^2=c_{c,0}^2/c_i^2\propto A/\xi$|⁠); see equation 16). The mass-loss rate is given in units of

$$\begin{eqnarray*} \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\nonumber\\ &&\times \left(\frac{n_0}{1 \, {\rm cm^{-3}}}\right) \left(\frac{c_i}{10 \, \, {\rm km \, \, s}^{-1}}\right). \end{eqnarray*}$$

(29)

To express the mass-loss rate in CR driven winds in a more intuitive form we take V_g ≫ c_i and use the simplified expression for the mass-loss rate in equation (23). We then write the CR energy injection rate at the base of the outflow as

$$\begin{eqnarray*} \dot{E}_c = \epsilon _c \dot{M}_* c^2, \end{eqnarray*}$$

(30)

where |$\dot{M}_*$| is the star formation rate and ϵ_c ≡ 10^−6.3ϵ_{c, −6.3} is related to the fraction of SNe energy that goes into CRs: for 10⁵¹ ergs per SNe and 1 SNe per 100 M_⊙ of stars formed, ϵ_c = 10^−6.3 if |$10 {{\ \rm per\ cent}}$| of the SNe energy goes into primary CRs. Equation (23) can then be approximately rewritten as (taking r_s ≃ r₀ to simplify equation 21)

$$\begin{eqnarray*} \frac{\dot{M}_w}{\dot{M}_*} &\simeq & 0.8 \ \epsilon _{c,-6.3} \left(\frac{c_i r_0}{\kappa }\right) \left(\frac{100 \, \, {\rm km \, \, s}^{-1}}{V_{\rm g,\, eff}}\right)^{2}\nonumber \\ &\simeq & 0.08 \ \epsilon _{c,-6.3} \, \left(\frac{c_i r_0}{{30\, \rm km\, \, s^{-1} \, kpc}}\right) \nonumber \\ &&\times \, \left(\frac{10^{29}\, {\rm cm^2\, s^{-1}}}{\kappa }\right) \left(\frac{100 \, \, {\rm km \, \, s}^{-1}}{V_{\rm g,\, eff}}\right)^{2}. \end{eqnarray*}$$

(31)

Per equation (7) and the associated discussion, |$\kappa \sim \, 10 c_i r_0$| is plausible; we have scaled to representative values in the second line. This expression shows that significant wind mass loading relative to the global star formation rate is in principle possible, and that it should grow strongly with decreasing |$V_{\rm g,\, eff}$|⁠. However, if κ ∼ 10²⁹ cm² s⁻¹ is appropriate, the mass-loading factors in dwarf galaxies are relatively modest compared to the values |$\dot{M}_w \gg \dot{M}_*$| needed to reconcile the stellar and dark matter mass functions. We return in Section 4 to an observational calibration of the diffusion coefficients in other star-forming galaxies.

$Analytic mass-loss rates for CR-driven galactic winds in the limit of rapid diffusion (equation (28) with hc from equation 27), as a function of the strength of gravity relative to the gas sound speed in the disc (Vg/ci) and the base CR pressure (${\rm p_{c,0}/\rho _0 c_i^2}$). The mass-loss rates are normalized by equation (29). Labelled values of Vg/ci on the colour bar are logarithmically distributed and correspond to the curves on the plot.$

Figure 1.

Analytic mass-loss rates for CR-driven galactic winds in the limit of rapid diffusion (equation (28) with h_c from equation 27), as a function of the strength of gravity relative to the gas sound speed in the disc (V_g/c_i) and the base CR pressure (⁠|${\rm p_{c,0}/\rho _0 c_i^2}$|⁠). The mass-loss rates are normalized by equation (29). Labelled values of V_g/c_i on the colour bar are logarithmically distributed and correspond to the curves on the plot.

Open in new tab Download slide

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

$$\begin{eqnarray*} \dot{\Sigma }_* \simeq \frac{\sqrt{8} \pi G \Sigma _g^2 \phi }{v_*}, \end{eqnarray*}$$

(32)

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

$$\begin{eqnarray*} \frac{\dot{M}_{\rm w}}{\dot{M}_*} \simeq 8\left(\frac{1/4}{h_c}\right) \left(\frac{c_i v_*}{3 \times 10^4 \, {\rm km^2 \, s^{-2}}}\right) \left(\frac{100 \, \, {\rm km \, \, s}^{-1}}{V_g}\right)^{2} \frac{p_{c,0}}{\pi G \Sigma _g^2 \phi } \nonumber\\ . \end{eqnarray*}$$

(33)

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

2.1 Wind energetics, momentum flux, and velocity

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

$$\begin{eqnarray*} \dot{E}_{w} = \dot{M}_w \left(\frac{1}{2} v^2 + c_i^2 \ln \rho + 4 c_c^2 + \Phi \right) + \dot{E}_{\rm c} = {\rm constant} , \end{eqnarray*}$$

(34)

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

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

$$\begin{eqnarray*} v_\infty &\simeq& 2 {\rm V_{g,eff}}\left(\frac{3 \kappa }{r_0 c_i}\right)^{1/2} \left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{-c_i^2/(4V_g^2)} \nonumber\\ &\simeq& 10^3 \, \, {\rm km \, \, s}^{-1}\left(\frac{V_g}{100 \, {\rm km \, s^{-1}}}\right) \left(\frac{\kappa }{10^{29} {\rm cm^2 \, s^{-1}}}\frac{\rm 30 \, kpc \, {\rm km \, s^{-1}}}{r_0 c_i}\right)^{1/2}\nonumber\\ ,\end{eqnarray*}$$

(35)

where the second expression assumes V_g ≫ c_i. Recall that we require κ ≳ r₀c_i given our assumptions used in deriving the approximate wind solutions here. In this case, equation (35) shows that v_∞ ≳ V_{g, eff}, with v_∞ ∼ 1 − 10V_g plausible. We show below (Fig. 7) that equation (35) agrees very well with our numerical simulations.

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

$$\begin{eqnarray*} \dot{p}_w=\dot{M}_w v_\infty \simeq \frac{4\pi r_0^2 c_i p_{c,0}}{h_c V_{\rm g}}\left(\frac{3\kappa }{r_0c_i}\right)^{1/2}, \end{eqnarray*}$$

(36)

where we assume V_g ≫ c_i to simplify the expressions and highlight the key scalings. Assuming no pionic losses in the host galaxy, we can write this expression in terms of the star formation rate using equations (24) and (30) as

$$\begin{eqnarray*} \dot{p}_w=\dot{M}_w v_\infty \simeq \frac{\epsilon _c\dot{M}_* c^2}{V_{\rm g}}\left(\frac{r_0 c_i}{3\kappa }\right)^{1/2}. \end{eqnarray*}$$

(37)

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

$$\begin{eqnarray*} \frac{\dot{p}_w}{\dot{p}_{*}}&\simeq &\frac{\epsilon _c}{\epsilon _{*}}\frac{c}{V_{\rm g,\, eff}}\left(\frac{r_0 c_i}{3\kappa }\right)^{1/2} \nonumber \\ &\simeq &0.5\, \frac{\epsilon _{c,-6.3}}{\epsilon _{*,\, -3.3}}\left(\frac{100\, \rm km \, \, s^{-1}}{V_{\rm g}}\right) \nonumber \\ &&\times \left(\frac{r_0\, c_i}{30 \, \rm km \, \, s^{-1}\, \, kpc}\frac{\rm 10^{29}\, cm^2\, \, s^{-1}}{\kappa }\right)^{1/2}. \end{eqnarray*}$$

(38)

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

2.2 Maximum mass-loss rate

There is a strict upper limit to the mass-loss rate associated with CR-driven winds that is set by energy conservation. This is the regime in which gravity significantly modifies the energetics of the outflow and cannot be neglected as we did in deriving equation (35). If the asymptotic speed of the wind vanishes then |$\dot{E}_{c}(r_0) + \dot{M}_{\rm max} \Phi (r_0) \simeq 0$|⁠, i.e. |$\dot{M}_{\rm max} \simeq 2 \dot{E}_{c}(r_0)/v^2_{\rm esc}(r_0)$| where v_esc(r₀) is the escape speed from the base of the wind.³ This regime is analogous to photon-tired winds in stellar wind theory (Owocki & Gayley 1997). Using |$\dot{E}_c = \epsilon _c \dot{M}_* c^2$| (equation 30), we can write

$$\begin{eqnarray*} \frac{\dot{M}_{\rm max}}{\dot{M}_*} \simeq \, 9 \, \epsilon _{c,-6.3} \, \left(\frac{100 \, \, {\rm km \, \, s}^{-1}}{v_{\rm esc}(r_0)}\right)^2, \end{eqnarray*}$$

(39)

which has an ‘energy-driven’ velocity scaling (Murray, Quataert & Thompson 2005). Alternatively, using |$\dot{E}_{c}(r_0) \simeq 12 \pi r_0 \kappa p_{c,0}/h_c$| we find

$$\begin{eqnarray*} \dot{M}_{\rm max}&\simeq &\frac{24\pi r_0\kappa p_{c,\, 0}}{v_{\rm esc}^2(r_0)h_c}\simeq \frac{24\pi r_0\kappa \rho _0}{h_c}\left(\frac{c_{c,\, 0}}{v_{\rm esc}(r_0)}\right)^2 \nonumber \\ &\simeq &23\, {\rm M_\odot \, \, yr^{-1}}\left(\frac{r_0}{\rm kpc}\right)\left(\frac{\kappa }{10^{29}\rm cm^2\, \, s^{-1}}\right)\left(\frac{n_0}{\rm cm^{-3}}\right) \nonumber \\ &&\times \left(\frac{1/4}{h_c}\right) \left(\frac{c_{c,\, 0}}{10\, {\rm km\, \, s^{-1}}}\right)^2\left(\frac{100\, {\rm km\, \, s^{-1}}}{v_{\rm esc}(r_0)}\right)^2 \end{eqnarray*}$$

(40)

The maximum cosmic ray-driven wind momentum will be realized for mass-loss rates a bit below |$\dot{M}_{\rm max}$| when v_∞ ∼ v_esc(r₀) (at |$\dot{M}_{\rm max}$|⁠, v_∞ ≡ 0 and so |$\dot{p}_w = 0$|⁠). A rough estimate is

$$\begin{eqnarray*} \frac{\dot{p}_{w,\, {\rm max}}}{\dot{p}_*}\sim 5 \frac{\epsilon _{c,-6.3}}{\epsilon _{*,\, -3.3}} \left(\frac{100\, \rm km \, \, s^{-1}}{v_{\rm esc}(r_0)}\right). \end{eqnarray*}$$

(41)

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

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

$$\begin{eqnarray*} \dot{M}_w \simeq \dot{M}_{\rm max}\, \left(\frac{r_0 c_i}{3 \kappa }\right), \end{eqnarray*}$$

(42)

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

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

3 NUMERICAL SIMULATIONS

3.1 Equations

We solve the time-dependent cosmic ray hydrodynamic equations based on the two moment approach as developed by Jiang & Oh (2018) in 1D spherical polar coordinates.⁴ This algorithm has been implemented in the magnetohydrodynamic code ATHENA++ (Stone et al. 2020). As in Section 2, we use an isothermal equation of state for the gas with isothermal sound speed c_i and take the gravitational potential to be |$\Phi =2V_g^2\ln r$|⁠.

Neglecting CR streaming, the full set of equations for gas density ρ, flow velocity v, comic ray energy density E_c, and flux F_c in 1D spherical polar coordinates are

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v\right)&=&0,\nonumber \\ \frac{\partial \left(\rho v\right)}{\partial t} +\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2 \rho v^2\right)&=&-\rho \frac{\partial \Phi }{\partial r}-c_i^2\frac{\partial \rho }{\partial r}\nonumber\\ &&+\,\,\sigma _c\left[F_c-v(E_c+p_c)\right],\nonumber \\ \frac{\partial E_c}{\partial t}+\frac{1}{r^2}\frac{\partial \left(r^2F_c\right)}{\partial r}&=& - v \, \sigma _c [F_c-v(E_c+p_c)],\nonumber \\ \frac{1}{V_m^2}\frac{\partial F_c}{\partial t}+\frac{\partial p_c}{\partial r}&=&-\sigma _c\left[F_c-v(E_c+p_c)\right]. \end{eqnarray*}$$

(43)

Here the cosmic ray pressure is p_c = E_c/3. The reduced speed of light is V_m, which is chosen to be much larger than v in the whole simulation box. For CR transport by spatially independent diffusion, σ_c is a constant and is related to the diffusion coefficient used elsewhere in this paper by κ ≡ 1/(3σ_c).⁵

Substituting the right-hand-side of the fourth of equations (43) into the 3rd term on the right-hand-side of the 2nd of equations (43) produces the usual ∂p_c/∂r cosmic ray pressure gradient in the momentum equation, along with another term related to the time variation of the CR flux. Likewise, the term on the right-hand-side of the 3rd of equations (43) becomes the usual v∂p_c/∂r term related to CR pdV work, again with another term related to the time variation of the CR flux. The steady state versions of equations (43) are thus equivalent to the steady state equations (4)–(6) solved in Section 2.

3.2 Initial and boundary conditions

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

$$\begin{eqnarray*} \rho (r)&=&\rho _0 \left(r/r_0\right)^{-V_g^2/c_i^2},\nonumber \\ E_c(r)&=&3p_{c,0}\left(r/r_0\right)^{-V_g^2/c_i^2}. \end{eqnarray*}$$

(44)

We apply floor values of 10⁻⁴ for ρ/ρ₀ and p_c(r)/p_{c, 0} in the initial condition. The flow velocity and cosmic ray flux are set to be zero in the whole simulation domain initially.

For the bottom boundary condition at r₀, the cosmic ray energy density E_c and gas density ρ are fixed to be the initial values. The flow velocity in the ghost zones is set such that mass flux ρvr² is continuous from the last active zone to the ghost zones. We keep the gradient of the cosmic ray flux continuous across the bottom boundary. For the boundary condition at the top of the simulation domain, we keep the gradient of ρ, E_c, and F_c continuous across the boundary. The flow velocity at the top boundary is set by requiring the mass flux ρvr² to be continuous.

3.3 Overview of simulations

Table 1 summarizes our suite of simulations. The key physical parameters of each simulation are V_g/c_i and |$p_{c,0}/\rho _0 c_i^2$|⁠, as in the analytics of Section 2, as well as the diffusion coefficient κ (in units of r₀c_i). The key numerical parameters are the resolution, reduced speed of light, and box size.

Table 1.

Open in new tab

Parameters and properties of our numerical simulations. Numerical resolution is δr/r = 5.25 × 10⁻⁴ unless otherwise noted. Columns are diffusion coefficient κ, velocity V_g of the isothermal gravitational potential, base CR pressure p_{c, 0}, outer radius of domain r_out, reduced speed of light V_m, simulation mass-loss rate |$\dot{M}_{\rm sim}$| in units of |$\dot{M}_0$| (equation 29), simulation mass-loss rate in units of |$\dot{M}_{\rm max}$| (the maximum mass-loss rate allowed by energy conservation; equation 47), the dimensionless base CR scale-height h_c (equation 10), and the kinetic energy and CR enthalpy fluxes at the top of the domain relative to the diffusive CR flux at the base |$(\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out}))/\dot{E}_c(r_0)$|⁠. When the latter is ≃1, most of the CR energy flux at the base of the wind goes into kinetic energy at large radii, while when it is ≲1, most of the CR energy goes into work against gravity, and |$\dot{M}_w \simeq \dot{M}_{\rm max}$| (equations 40 and 47).

κ	V_g	p_{c, 0}	r_out	V_m	\|$\dot{M}_{\rm sim}$\|	\|$\dot{M}_{\rm sim}$\|	v(r_out)	h_c	\|$\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out})$\|
(r₀c_i)	(c_i)	(⁠\|$\rho _0 c_i^2$\|⁠)	(r₀)	(c_i)	(⁠\|$\dot{M}_0$\|⁠)	(⁠\|$\dot{M}_{\rm max}$\|⁠)	(V_g)	–	\|$(\dot{E}_c[r_0])$\|
20^a	200	400	7.6	60 000	0.21	0.032	13.8	0.023	0.95
10^b	100	100	7.6	30 000	0.16	0.067	9.5	0.03	0.94
3.3^c	33	11	7.6	10 000	0.09	0.19	5	0.053	0.8
30	10	1	5	3000	0.016	0.018	13.6	0.32	0.93
10	10	1	5	3000	0.022	0.054	8.2	0.23	0.93
3.3	10	1	5	3000	0.03	0.15	4.6	0.16	0.83
2	10	1	5	3000	0.03	0.42	3.3	0.13	0.57
1	10	1	15	3000	0.03	0.75	0.5	0.14	0.23
0.33	10	1	15	3000	0.007	0.78	0.4	0.21	0.16
0.11	10	1	25	3000	0.0019	0.85	0.43	0.23	0.08
10	10	0.3	5	3000	0.007	0.058	8.2	0.23	0.99
10	10	0.1	5	3000	0.0027	0.056	8.1	0.24	0.95
10	10	3	15	3000	0.065	0.089	9.4	0.23	0.9
10	6	1	5	3000	0.043	0.055	7.5	0.33	0.89
3.3	6	1	15	3000	0.055	0.26	4.9	0.24	0.67
10	3	1	15	3000	0.136	0.11	7.6	0.5	0.9
3.3	3	1	15	3000	0.156	0.3	4.2	0.4	0.74

κ	V_g	p_{c, 0}	r_out	V_m	\|$\dot{M}_{\rm sim}$\|	\|$\dot{M}_{\rm sim}$\|	v(r_out)	h_c	\|$\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out})$\|
(r₀c_i)	(c_i)	(⁠\|$\rho _0 c_i^2$\|⁠)	(r₀)	(c_i)	(⁠\|$\dot{M}_0$\|⁠)	(⁠\|$\dot{M}_{\rm max}$\|⁠)	(V_g)	–	\|$(\dot{E}_c[r_0])$\|
20^a	200	400	7.6	60 000	0.21	0.032	13.8	0.023	0.95
10^b	100	100	7.6	30 000	0.16	0.067	9.5	0.03	0.94
3.3^c	33	11	7.6	10 000	0.09	0.19	5	0.053	0.8
30	10	1	5	3000	0.016	0.018	13.6	0.32	0.93
10	10	1	5	3000	0.022	0.054	8.2	0.23	0.93
3.3	10	1	5	3000	0.03	0.15	4.6	0.16	0.83
2	10	1	5	3000	0.03	0.42	3.3	0.13	0.57
1	10	1	15	3000	0.03	0.75	0.5	0.14	0.23
0.33	10	1	15	3000	0.007	0.78	0.4	0.21	0.16
0.11	10	1	25	3000	0.0019	0.85	0.43	0.23	0.08
10	10	0.3	5	3000	0.007	0.058	8.2	0.23	0.99
10	10	0.1	5	3000	0.0027	0.056	8.1	0.24	0.95
10	10	3	15	3000	0.065	0.089	9.4	0.23	0.9
10	6	1	5	3000	0.043	0.055	7.5	0.33	0.89
3.3	6	1	15	3000	0.055	0.26	4.9	0.24	0.67
10	3	1	15	3000	0.136	0.11	7.6	0.5	0.9
3.3	3	1	15	3000	0.156	0.3	4.2	0.4	0.74

Notes. ^aδr/r = 8.25 × 10⁻⁶ for 1 ≤ r ≤ 1.14; c_i is 1/20 of κ = 1, V_g = 10 sim

^bδr/r = 2.5 × 10⁻⁵ for 1 ≤ r ≤ 1.11; c_i is 1/10 of κ = 1, V_g = 10 sim

^cδr/r = 6.6 × 10⁻⁵ for 1 ≤ r ≤ 1.11; c_i is 1/3 of κ = 1, V_g = 10 sim

Table 1.

Open in new tab

Parameters and properties of our numerical simulations. Numerical resolution is δr/r = 5.25 × 10⁻⁴ unless otherwise noted. Columns are diffusion coefficient κ, velocity V_g of the isothermal gravitational potential, base CR pressure p_{c, 0}, outer radius of domain r_out, reduced speed of light V_m, simulation mass-loss rate |$\dot{M}_{\rm sim}$| in units of |$\dot{M}_0$| (equation 29), simulation mass-loss rate in units of |$\dot{M}_{\rm max}$| (the maximum mass-loss rate allowed by energy conservation; equation 47), the dimensionless base CR scale-height h_c (equation 10), and the kinetic energy and CR enthalpy fluxes at the top of the domain relative to the diffusive CR flux at the base |$(\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out}))/\dot{E}_c(r_0)$|⁠. When the latter is ≃1, most of the CR energy flux at the base of the wind goes into kinetic energy at large radii, while when it is ≲1, most of the CR energy goes into work against gravity, and |$\dot{M}_w \simeq \dot{M}_{\rm max}$| (equations 40 and 47).

κ	V_g	p_{c, 0}	r_out	V_m	\|$\dot{M}_{\rm sim}$\|	\|$\dot{M}_{\rm sim}$\|	v(r_out)	h_c	\|$\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out})$\|
(r₀c_i)	(c_i)	(⁠\|$\rho _0 c_i^2$\|⁠)	(r₀)	(c_i)	(⁠\|$\dot{M}_0$\|⁠)	(⁠\|$\dot{M}_{\rm max}$\|⁠)	(V_g)	–	\|$(\dot{E}_c[r_0])$\|
20^a	200	400	7.6	60 000	0.21	0.032	13.8	0.023	0.95
10^b	100	100	7.6	30 000	0.16	0.067	9.5	0.03	0.94
3.3^c	33	11	7.6	10 000	0.09	0.19	5	0.053	0.8
30	10	1	5	3000	0.016	0.018	13.6	0.32	0.93
10	10	1	5	3000	0.022	0.054	8.2	0.23	0.93
3.3	10	1	5	3000	0.03	0.15	4.6	0.16	0.83
2	10	1	5	3000	0.03	0.42	3.3	0.13	0.57
1	10	1	15	3000	0.03	0.75	0.5	0.14	0.23
0.33	10	1	15	3000	0.007	0.78	0.4	0.21	0.16
0.11	10	1	25	3000	0.0019	0.85	0.43	0.23	0.08
10	10	0.3	5	3000	0.007	0.058	8.2	0.23	0.99
10	10	0.1	5	3000	0.0027	0.056	8.1	0.24	0.95
10	10	3	15	3000	0.065	0.089	9.4	0.23	0.9
10	6	1	5	3000	0.043	0.055	7.5	0.33	0.89
3.3	6	1	15	3000	0.055	0.26	4.9	0.24	0.67
10	3	1	15	3000	0.136	0.11	7.6	0.5	0.9
3.3	3	1	15	3000	0.156	0.3	4.2	0.4	0.74

κ	V_g	p_{c, 0}	r_out	V_m	\|$\dot{M}_{\rm sim}$\|	\|$\dot{M}_{\rm sim}$\|	v(r_out)	h_c	\|$\dot{E}_k(r_{\rm out}) + \dot{E}_{c,\rm adv}(r_{\rm out})$\|
(r₀c_i)	(c_i)	(⁠\|$\rho _0 c_i^2$\|⁠)	(r₀)	(c_i)	(⁠\|$\dot{M}_0$\|⁠)	(⁠\|$\dot{M}_{\rm max}$\|⁠)	(V_g)	–	\|$(\dot{E}_c[r_0])$\|
20^a	200	400	7.6	60 000	0.21	0.032	13.8	0.023	0.95
10^b	100	100	7.6	30 000	0.16	0.067	9.5	0.03	0.94
3.3^c	33	11	7.6	10 000	0.09	0.19	5	0.053	0.8
30	10	1	5	3000	0.016	0.018	13.6	0.32	0.93
10	10	1	5	3000	0.022	0.054	8.2	0.23	0.93
3.3	10	1	5	3000	0.03	0.15	4.6	0.16	0.83
2	10	1	5	3000	0.03	0.42	3.3	0.13	0.57
1	10	1	15	3000	0.03	0.75	0.5	0.14	0.23
0.33	10	1	15	3000	0.007	0.78	0.4	0.21	0.16
0.11	10	1	25	3000	0.0019	0.85	0.43	0.23	0.08
10	10	0.3	5	3000	0.007	0.058	8.2	0.23	0.99
10	10	0.1	5	3000	0.0027	0.056	8.1	0.24	0.95
10	10	3	15	3000	0.065	0.089	9.4	0.23	0.9
10	6	1	5	3000	0.043	0.055	7.5	0.33	0.89
3.3	6	1	15	3000	0.055	0.26	4.9	0.24	0.67
10	3	1	15	3000	0.136	0.11	7.6	0.5	0.9
3.3	3	1	15	3000	0.156	0.3	4.2	0.4	0.74

Notes. ^aδr/r = 8.25 × 10⁻⁶ for 1 ≤ r ≤ 1.14; c_i is 1/20 of κ = 1, V_g = 10 sim

^bδr/r = 2.5 × 10⁻⁵ for 1 ≤ r ≤ 1.11; c_i is 1/10 of κ = 1, V_g = 10 sim

^cδr/r = 6.6 × 10⁻⁵ for 1 ≤ r ≤ 1.11; c_i is 1/3 of κ = 1, V_g = 10 sim

The units for the results of our numerical simulations are summarized in Table 2: gas density is in units of the base density ρ₀, 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.

Table 2.

Open in new tab

Summary of units for numerical simulations. c_i is the gas isothermal sound speed (assumed to be a constant), ρ₀ is the gas density at the base of the wind (at radius r = r₀), and p_{c, 0} is the base CR pressure.

Quantity	Symbol	Units
Radial velocity	v	c_i
‘Isothermal’ CR sound speed	\|$c_c \equiv \left(\frac{p_c}{\rho }\right)^{1/2}$\|	c_i
Gravitational velocity	V_g	c_i
Density	ρ	ρ₀
CR pressure	p_c	p_{c, 0}
CR flux	F_c	c_ip_{c, 0}
CR diffusion coefficient	κ	c_ir₀

Quantity	Symbol	Units
Radial velocity	v	c_i
‘Isothermal’ CR sound speed	\|$c_c \equiv \left(\frac{p_c}{\rho }\right)^{1/2}$\|	c_i
Gravitational velocity	V_g	c_i
Density	ρ	ρ₀
CR pressure	p_c	p_{c, 0}
CR flux	F_c	c_ip_{c, 0}
CR diffusion coefficient	κ	c_ir₀

Table 2.

Open in new tab

Summary of units for numerical simulations. c_i is the gas isothermal sound speed (assumed to be a constant), ρ₀ is the gas density at the base of the wind (at radius r = r₀), and p_{c, 0} is the base CR pressure.

Quantity	Symbol	Units
Radial velocity	v	c_i
‘Isothermal’ CR sound speed	\|$c_c \equiv \left(\frac{p_c}{\rho }\right)^{1/2}$\|	c_i
Gravitational velocity	V_g	c_i
Density	ρ	ρ₀
CR pressure	p_c	p_{c, 0}
CR flux	F_c	c_ip_{c, 0}
CR diffusion coefficient	κ	c_ir₀

Quantity	Symbol	Units
Radial velocity	v	c_i
‘Isothermal’ CR sound speed	\|$c_c \equiv \left(\frac{p_c}{\rho }\right)^{1/2}$\|	c_i
Gravitational velocity	V_g	c_i
Density	ρ	ρ₀
CR pressure	p_c	p_{c, 0}
CR flux	F_c	c_ip_{c, 0}
CR diffusion coefficient	κ	c_ir₀

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

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

3.4 Numerical solutions for κ ≳ r₀c_i

Fig. 2 shows the resulting density, CR pressure, fluid velocity, and CR sound speed profiles for V_g/c_i = 10, varying κ from 0.33 to 30 (in units of r₀c_i; Table 2). The vertical dotted line in Fig 2 (as well as Figs 3 and 4 below) is the critical point predicted in equation (21), namely r_s ≃ 1.047r₀.

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

Figure 2.

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

Open in new tab Download slide

Numerator and denominator (equations 45) of the steady state wind equation evaluated for our κ = 10, Vg = 10 simulation. N2 is the approximation to the numerator of the wind equation used in our analytics (equation 46). The vertical dotted line is the analytic estimate of the location of the critical point (equation 21). The numerical simulation reaches a steady state that passes through a critical point close to that predicted by our (approximate) analytics.

Figure 3.

Numerator and denominator (equations 45) of the steady state wind equation evaluated for our κ = 10, V_g = 10 simulation. N₂ is the approximation to the numerator of the wind equation used in our analytics (equation 46). The vertical dotted line is the analytic estimate of the location of the critical point (equation 21). The numerical simulation reaches a steady state that passes through a critical point close to that predicted by our (approximate) analytics.

Open in new tab Download slide

$Energetics of the wind evaluated for our κ = 10, Vg = 10 simulation. The total energy flux ($\dot{E}_w$; see equation 34) and mass flux ($\dot{M}_w$) are independent of radius. At the base of the wind, nearly all of the energy is carried by the diffusing CRs ($\dot{E}_c$). Advection of CR enthalpy by the motion of the gas ($\dot{E}_{c,\rm adv}$) is important at intermediate radii while at large radii most of the energy is in the kinetic energy of the outflowing gas ($\dot{E}_{k}$). For the logarithmic potential used in our simulations, we define $\dot{E}_{\rm grav}$ using $\Phi =2 V_g^2 \ln (r/r_0)$ so that $\dot{E}_{\rm grav} \gt 0$.$

Figure 4.

Energetics of the wind evaluated for our κ = 10, V_g = 10 simulation. The total energy flux (⁠|$\dot{E}_w$|⁠; see equation 34) and mass flux (⁠|$\dot{M}_w$|⁠) are independent of radius. At the base of the wind, nearly all of the energy is carried by the diffusing CRs (⁠|$\dot{E}_c$|⁠). Advection of CR enthalpy by the motion of the gas (⁠|$\dot{E}_{c,\rm adv}$|⁠) is important at intermediate radii while at large radii most of the energy is in the kinetic energy of the outflowing gas (⁠|$\dot{E}_{k}$|⁠). For the logarithmic potential used in our simulations, we define |$\dot{E}_{\rm grav}$| using |$\Phi =2 V_g^2 \ln (r/r_0)$| so that |$\dot{E}_{\rm grav} \gt 0$|⁠.

Open in new tab Download slide

The left-hand panel of Fig. 2 shows that for all of the solutions, the CR pressure falls off much more slowly than the density near the base. This is because the gas density scale-height is initially quite small, ∼r₀(c_i/V_g)² ∼ 10⁻²r₀ while the CR pressure scale-height is significantly larger due to CR diffusion ∼h_cr₀ ∼ 0.25r₀ (equation 27). For all of the solutions, the hydrostatic equilibrium approximation for the density profile (equation 14) is reasonable at very small radii. There is, however, a bifurcation in the density profiles as a function of κ. The solutions with κ = 0.33 and 1 are closer to hydrostatic over the entire radial range while for κ ≳ 1 the density falls off much more rapidly at large radii. The right-hand panel of Fig. 2 shows that this bifurcation in density profiles corresponds to a bifurcation in velocity profiles. The solutions with κ ≳ 3 drive strong winds that accelerate to a velocity larger than the escape speed by r ∼ 1.2r₀. By contrast, for κ = 1 the velocities are much lower and for κ = 0.33 the solution is subsonic for all of the radii shown in Fig. 2. We return to the lowest κ simulation in Section 3.5 and show that it does produce a wind, but one whose character is very different from the higher κ simulations focused on here.

Fig. 3 quantifies the extent to which the numerical solution satisfies the time steady critical point conditions, for the κ = 10 numerical simulation. Specifically, from equation (11) we define the denominator and numerator of the critical point equation to be

$$\begin{eqnarray*} D \equiv v^2 - c_i^2 \ \ \ {\rm and} \ \ \ N \equiv \frac{2 c_i^2}{r} - \frac{1}{\rho }\frac{\mathrm{ d}p_c}{\mathrm{ d}r} - \frac{2 V_g^2}{r} \end{eqnarray*}$$

(45)

and we further define

$$\begin{eqnarray*} N _2\equiv \frac{2 c_i^2}{r} + \frac{p_{c,0} r_0}{\rho r^2 h_c} - \frac{2 V_g^2}{r} \end{eqnarray*}$$

(46)

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

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

Fig. 4 shows that interior to the critical point, the energy flux is dominated by CR diffusion, as assumed in the analytic calculation in Section 2. At intermediate radii the CR advective flux dominates the energy transport while eventually at large radii, when the CR pressure has accelerated the gas to well above the CR sound speed and escape speeds, the gas kinetic energy flux dominates. Note that the gravitational contribution to the energy flux is small at all radii, consistent with |$\dot{M}_w \ll \dot{M}_{\rm max}$| and v > V_g. The qualitative features of Fig. 4 are the same for κ = 3−30. For κ = 1, however, most of the energy at the outer edge of the domain is in CR enthalpy rather than gas kinetic energy.⁶ For κ = 0.33, the energetics of the flow is yet more different from that in Fig. 4, as we discuss in Section 3.5.

Fig. 5 shows how the κ = 10 solutions depend on the depth of the gravitational potential V_g (in units of the isothermal sound speed c_i). The solutions are qualitatively very similar. The primary difference is that weaker gravity (lower V_g) implies a larger density scale-height at the base of the wind. Correspondingly, the wind accelerates significantly more slowly.

$Dependence of density and pressure profiles (left-hand panel) and velocity and CR sound speed profiles (right-hand panel) on the depth of the gravitational potential Vg (in units of ci, the isothermal gas sound speed). All solutions are for κ = 10 and equal gas and CR pressures at the base of the wind ($p_{c,0}=\rho _0 c_i^2$). The CR pressure profiles are relatively independent of the gravitational potential, being largely set by rapid CR diffusion. By contrast, the density falls off more gradually, and the flow accelerates more gradually, for weaker gravity.$

Figure 5.

Dependence of density and pressure profiles (left-hand panel) and velocity and CR sound speed profiles (right-hand panel) on the depth of the gravitational potential V_g (in units of c_i, the isothermal gas sound speed). All solutions are for κ = 10 and equal gas and CR pressures at the base of the wind (⁠|$p_{c,0}=\rho _0 c_i^2$|⁠). The CR pressure profiles are relatively independent of the gravitational potential, being largely set by rapid CR diffusion. By contrast, the density falls off more gradually, and the flow accelerates more gradually, for weaker gravity.

Open in new tab Download slide

3.5 Low κ ≲ r₀c_i solutions

The analytic derivations in Section 2 primarily assumed κ ≳ r₀c_i. Our numerical simulations with κ ≳ r₀c_i are consistent with the properties of these analytic solutions, as we discuss further in the next section. Here we present numerical wind models with κ ≲ r₀c_i, which differ dramatically from their high κ counterparts. We stress that the higher κ solutions are likely the most astrophysically relevant, given estimates of diffusion coefficients from MW CR data (Section 2) and gamma-ray data from pion decay (Section 4).

Fig. 6 shows the kinematics (top) and energetics (bottom) of the wind for our κ = 0.33, V_g = 10 simulation, which differs significantly from the higher κ simulations. The velocity profile in Fig. 6 is shown over a larger range of radii than was plotted in Fig. 2. The flow accelerates very gradually reaching the sonic point only at r_s ∼ 8r₀,⁷ compared to r_s ≃ r₀ for the high κ solutions. The most striking property of the energetics is that |$\dot{E}_{\rm grav} \simeq \dot{E}_w$| (see equation 34) at large radii while at the base of the wind, nearly all of the energy is carried by CR diffusion (⁠|$\dot{E}_c$|⁠); recall that for our logarithmic potential |$\Phi = 2 V_g^2 \ln (r/r_0)$| and so |$\dot{E}_{\rm grav} \gt 0$|⁠. The fact that |$\dot{E}_{\rm grav}(r_{\rm out}) \simeq \dot{E}_c(r_0)$| implies that most of the energy supplied to CRs at the base goes into lifting (almost hydrostatically) the gas out to large radii. As a result, the gas kinetic energy is a minor contribution to the energy flux at large radii, in marked contrast to the high κ simulation shown in Fig. 4. The condition |$\dot{E}_{c} \simeq \dot{E}_{\rm grav}$| is precisely the condition for the CR-analogue of photon-tired winds; this is also the regime in which |$\dot{M}_w \sim \dot{M}_{\rm max}$| (equation 40). For our simulations, we can define the escape velocity from the base of the wind to be the velocity needed to reach the outer edge of the box, i.e, |$v_{\rm esc}(r_0) = 2 V_g \sqrt{\ln (r_{\rm out}/r_0)}$|⁠. Equation (40) then becomes

$$\begin{eqnarray*} \frac{\dot{M}_{\rm max}}{\dot{M}_0} \simeq \frac{3}{2 h_c \ln (r_{\rm out}/r_0)} \frac{\kappa }{r_0 c_i} \frac{c_{c,0}^2}{V_g^2} \end{eqnarray*}$$

(47)

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

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

Figure 6.

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

Open in new tab Download slide

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

3.6 Synthesis of analytics versus numerics

Fig. 7 compares our numerical and analytic solutions for the mass-loss rate, terminal velocity, and base CR scale-height. We focus on the κ > r₀c_i solutions for which we have the most detailed analytic estimates. For κ < r₀c_i, |$\dot{M}_{w} \simeq \dot{M}_{\rm max}$| (equation 47), but we do not have a prediction for the terminal velocity in this regime.

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

Figure 7.

Ratio of the analytically predicted mass-loss rate (equation 28 with h_c from equation 27), terminal velocity (equation 35), and base CR scale-height (equation 27) compared to the simulation results. The agreement is good for the full range of simulated parameters. Multiplying h_c in equation (27) by ≃0.75 would remove the small systematic offset in h_c and |$\dot{M}_w$|⁠.

Open in new tab Download slide

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

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

4 IMPLICATIONS OF GAMMA-RAY OBSERVATIONS FOR COSMIC-RAY WIND MODELS

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

4.1 Non-thermal emission from pion decay

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

We consider a simple one-zone model in which a galaxy is characterized by its size r₀, gas surface density Σ_g, and isothermal/turbulent velocity c_i. The cosmic ray scale-height is H_c while the gas scale-height is H_g. For CR diffusion we expect H_c > H_g, as is indeed the case in our simulations in Section 3. The gamma-ray luminosity is set by the ratio of the time-scale for pion losses |$t_\pi \simeq t_0 \, (n_\pi /n_{eff})$| (where t₀ = 5 × 10⁷ yr and n_π = 1 cm⁻³ are set by the cross-section for pion production and the effective density n_eff is defined below) to the time-scale for CRs to escape by diffusion, |$t_{\rm diff} \simeq H_c^2/\kappa$|⁠. The effective density n_eff is the density averaged over the region the CRs are diffusing through and so is given by ≃nH_g/H_c where n is the mid-plane density of the galaxy. For reasons that will become clear, we choose to express t_π in terms of gas surface density by using Σ_g ≃ 2ρH_g. The ratio of the diffusion time-scale to the pion loss time-scale is then

$$\begin{eqnarray*} \frac{t_{\rm diff}}{t_{\pi }} \simeq \frac{\Sigma _g H_c}{2 \kappa \rho _\pi t_0 }, \end{eqnarray*}$$

(48)

where ρ_π = n_πm_p ≃ 1.67 × 10⁻²⁴ cm⁻³. Equation (48) shows that for a fixed ratio of hadronic losses (t_π) to CR escape (t_diff) there is a degeneracy between the CR scale-height and diffusion coefficient, with κ∝H_c. This is consistent with the known degeneracy between CR diffusion coefficient and ‘halo size’ in the literature (e.g. fig. 3 of Trotta et al. 2011).

The gamma-ray emission from pion decay in a star-forming galaxy is given by |$L_\gamma \simeq 1/3 \dot{E}_\pi$| where |$\dot{E}_\pi \simeq \dot{E}_c \min (t_{\rm diff}/t_\pi , 1)$| is the rate of energy loss to pion decay and the factor of 1/3 quantifies the fraction of pion decay in neutral versus charged pions. Defining |$\dot{E}_c = E_{cr} \dot{M}_*/m_*$| with⁸ E_cr = 10⁵⁰E_{cr, 50} the energy per SN going into CRs and m_* ≃ 100M_⊙ the total stellar mass formed per core-collapse SNe, we find (Lacki et al. 2011)

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

(49)

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

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

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

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

(50)

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

$Empirically constrained cosmic ray properties in star-forming galaxies inferred from gamma-ray observations of pion decay by Fermi, as a function of galaxy star formation rate. All calculations assume gas isothermal sound speed $c_i=10 \, {\rm km \, \, s}^{-1}$; scalings to other values of ci are shown in Fig. 9 and discussed in the text. Star formation rates of example galaxies are indicated near the x-axis. Three examples are considered corresponding to nuclear starbursts (black lines; r0 = 0.3 kpc, $V_g=150 \, {\rm km \, \, s}^{-1}$), dwarf galaxies (red lines; r0 = 1 kpc, $V_g= 50\, {\rm km \, \, s}^{-1}$), and star-forming disc galaxies (blue lines; r0 = 3 kpc, $V_g=150 \, {\rm km \, \, s}^{-1}$). Left-hand panel: tdiff/tπ∝Lγ/L* is the ratio of the CR diffusion time to the pion loss time and sets the gamma-ray luminosity of the galaxy (equation 49). We infer tdiff/tπ for different galaxy star formation rates using the Fermi correlation between Lγ and star formation rate. We then calculate the CR diffusion coefficient κ and the dimensionless CR scale-height hc = Hc/r0 using equations (48), (32), and (27). Middle: The fractional contribution of CR pressure to the pressure in the galactic disc ($p_{HE}=\pi G \Sigma _g^2 \phi)$ depends on the assumed size of the star-forming disc r0, with nuclear starbursts having suppressed pc/pHE because of rapid pion losses and CR diffusion in the dense nuclear regions. The predicted mass-loss rate relative to the star formation rate $\dot{M}_w/\dot{M}_*$ is largest for the disc galaxy model. Right-hand panel: Asymptotic wind speed v∞ and momentum flux in the wind $\dot{p}_w$ relative to the stellar radiation field ($\dot{p}_* = L_*/c$).$

Figure 8.

Empirically constrained cosmic ray properties in star-forming galaxies inferred from gamma-ray observations of pion decay by Fermi, as a function of galaxy star formation rate. All calculations assume gas isothermal sound speed |$c_i=10 \, {\rm km \, \, s}^{-1}$|⁠; scalings to other values of c_i are shown in Fig. 9 and discussed in the text. Star formation rates of example galaxies are indicated near the x-axis. Three examples are considered corresponding to nuclear starbursts (black lines; r₀ = 0.3 kpc, |$V_g=150 \, {\rm km \, \, s}^{-1}$|⁠), dwarf galaxies (red lines; r₀ = 1 kpc, |$V_g= 50\, {\rm km \, \, s}^{-1}$|⁠), and star-forming disc galaxies (blue lines; r₀ = 3 kpc, |$V_g=150 \, {\rm km \, \, s}^{-1}$|⁠). Left-hand panel: t_diff/t_π∝L_γ/L_* is the ratio of the CR diffusion time to the pion loss time and sets the gamma-ray luminosity of the galaxy (equation 49). We infer t_diff/t_π for different galaxy star formation rates using the Fermi correlation between L_γ and star formation rate. We then calculate the CR diffusion coefficient κ and the dimensionless CR scale-height h_c = H_c/r₀ using equations (48), (32), and (27). Middle: The fractional contribution of CR pressure to the pressure in the galactic disc (⁠|$p_{HE}=\pi G \Sigma _g^2 \phi)$| depends on the assumed size of the star-forming disc r₀, with nuclear starbursts having suppressed p_c/p_HE because of rapid pion losses and CR diffusion in the dense nuclear regions. The predicted mass-loss rate relative to the star formation rate |$\dot{M}_w/\dot{M}_*$| is largest for the disc galaxy model. Right-hand panel: Asymptotic wind speed v_∞ and momentum flux in the wind |$\dot{p}_w$| relative to the stellar radiation field (⁠|$\dot{p}_* = L_*/c$|⁠).

Open in new tab Download slide

4.2 Models calibrated to gamma-ray data

Given t_diff/t_π as a function of galaxy star formation rate (purple line Fig. 8), we now use equation (48) to constrain the CR diffusion coefficient. To do so, however, we need an estimate of the CR scale-height and gas surface density as a function of star formation rate. We model h_c = H_c/r₀ using equation (27) which assumes that the CR scale-height is set by diffusion in a CR-driven galactic wind. We model the gas surface density Σ_g using equation (32) which yields |$\dot{M}_* \simeq \pi r_0^2 \sqrt{8} \pi G \Sigma _g^2 \phi /v_*$| where r₀ is the size of the galactic disc. The free parameters of our model are thus the size of the star-forming galactic disc r₀, the galaxy circular velocity set by V_g, the gas isothermal sound speed c_i, and the dimensionless stellar contribution to the gravitational potential ϕ (as well as several ‘microphysics’ parameters such as the CR energy per supernovae). Given choices for these parameters, as well as the observed L_γ − L_* correlation, we solve equations (27), (32), and (48) for κ, h_c, and Σ_g. In what follows, we initially assume |$c_i=10 \, {\rm km \, \, s}^{-1}$| and consider parameters approximating three different classes of galaxies. Results for these classes are given in Fig. 8:

Nuclear starbursts: black lines in Fig. 8; r₀ = 0.3 kpc, |$V_g=150 \, {\rm km \, \, s}^{-1}$|⁠, ϕ = 1. This model is meant to approximate the nuclear starbursts M82 and NGC 253 with |$\dot{M}_*\simeq 10$| M_⊙ yr⁻¹, as well as ultraluminous starbursts like Arp 220 with |$\dot{M_*} \sim 10^2$| M_⊙ yr⁻¹.
Dwarf galaxies: red lines in Fig. 8; r₀ = 1 kpc, |$V_g= 50\, {\rm km \, \, s}^{-1}$|⁠, ϕ = 5. This model is meant to approximate normal star-forming dwarf galaxies like the LMC or SMC and their potentially much more rapidly star-forming counterparts.
Star-forming disc galaxies: blue lines in Fig. 8; r₀ = 3 kpc, |$V_g=150 \, {\rm km \, \, s}^{-1}$|⁠, ϕ = 5. This model is meant to approximate the Milky Way, M31, and other star-forming spirals in the local (⁠|$\dot{M}_* \sim \, {\rm M_\odot \, yr^{-1}}$|⁠) and high-redshift (⁠|$\dot{M}_* \sim 10-100 \, {\rm M_\odot \, yr^{-1}}$|⁠) Universe.

We discuss below (Fig. 9 and equations 55–65) the scaling of our results to other gas sound speeds c_i.

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

Figure 9.

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

Open in new tab Download slide

The dotted lines in the left-hand panel of Fig. 8 show the CR diffusion coefficient in units of 10²⁹ cm² s⁻¹ required to explain the gamma-ray luminosities of star-forming galaxies, per the method described in the previous paragraph. Fig. 8 also shows the dimensionless CR scale-height h_c we derive (dashed lines). For the dwarf and starburst models in Fig. 8, h_c ≃ 1 and the diffusion coefficient is ∼2 × 10²⁹ cm² s⁻¹ ≫r₀c_i with only a modest factor of few variation with star formation rate. For the star forming disc model, however, h_c ∼ 0.1 and κ ∼ r₀c_i ∼ 10²⁸ cm² s⁻¹. The dimensionless scale-height for our ‘disc’ model is comparable to that inferred from synchrotron emission in nearby star-forming disc galaxies (e.g. Krause et al. 2018), though the latter depends on both the magnetic field and CR scale-heights. Our inferred values of h_c ∼ 0.1 and κ ∼ 10²⁸ cm² s⁻¹ for disc galaxies in Fig. 8 are, however, factors of ∼5 smaller than preferred in phenomenological models constrained by Milky Way CR data (e.g. Trotta et al. 2011). However, the estimated gamma-ray luminosity of the Milky Way ≃ 8 × 10³⁸ erg s⁻¹ (Strong et al. 2010) is a factor of ≃ 3 lower than what would be predicted by its infrared luminosity given the correlations from Ackermann et al. (2012) and Griffin et al. (2016) used here. h_c and κ would then be 3 and 10 times larger, respectively (see equations 55 and 61 below), in better agreement with detailed Milky-Way modelling. In addition, we show below that for h_c ≲ 1, h_c∝c_i and κ∝c_i (equations 55 and 61). Our disc galaxy model in Fig. 8 would thus have a larger scale-height and diffusion coefficient if we assumed that the CRs primarily coupled to the warm-hot phase of the ISM, as is quite plausible.

The key dimensionless number from Section 3 that determines the properties of CR-driven galactic winds is κ/(r₀c_i); winds accelerate rapidly and reach speeds significantly larger than V_g if κ/(r₀c_i) ≳ 1 (Fig. 2). For our three galaxy models in Fig. 8, we note that r₀c_i ≃ 10²⁷ cm² s⁻¹ (starburst; black lines), r₀c_i ≃ 3 × 10²⁷ cm² s⁻¹ (dwarf; red lines), and r₀c_i ≃ 10²⁸ cm² s⁻¹ (disc; blue lines). A key conclusion from Fig. 8 is that the gamma-ray data on star-forming galaxies are consistent with diffusion coefficients that are in the regime of κ ≳ r₀c_i, though only marginally so for our star-forming disc model. This does not prove that such diffusion coefficients are correct, at a minimum because it is possible that CR escape is not set by diffusion but rather by streaming or advection in winds driven by other mechanisms (e.g. supernovae), in which case our constraint on the CR diffusion coefficient using equation (49) would not apply. But the diffusion coefficients inferred in Fig. 8 are none the less a useful and instructive observational check on diffusive CR-driven galactic wind models.

Given the diffusion coefficients inferred in the left-hand panel of Fig. 8, we can estimate the CR energy density in the galactic disc – which sets the base conditions for the wind – as follows: if the star formation rate per unit area is |$\dot{\Sigma }_*$|⁠, the CR pressure in the galaxy is given by (e.g. Thompson & Lacki 2013)⁹

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

(51)

where min(t_diff, t_π) determines the effective loss/escape time for CRs in the galaxy. Using equation (32) and that fact that t_diff ≲ t_π even at the highest star formation rates in Fig. 8, we find:

$$\begin{eqnarray*} \frac{p_{c,0}}{\pi G \Sigma ^2_g \phi } \simeq \frac{\sqrt{2} E_{cr} h_c r_0}{3 \kappa p_*} . \end{eqnarray*}$$

(52)

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

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

$$\begin{eqnarray*} \frac{\dot{M}_w}{\dot{M}_*} \simeq \frac{E_{cr}/m_*}{{\rm V^2_{g,eff}}} \frac{c_i r_0}{3 \kappa } \left(\frac{4 h_c V_g^4}{c_{c,0}^2 c_i^2}\right)^{c_i^2/(2V_g^2)} , \end{eqnarray*}$$

(53)

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

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

The right-hand panel of Fig. 8 shows the terminal velocity (equation 35) and momentum flux |$\dot{p}_w=\dot{M}_w v_\infty$| (in units of the photon momentum flux |$\dot{p}_*=L/c$|⁠) for the same three galaxy models (compare with Lochhaas, Thompson & Schneider 2021). For the disc and dwarf models |$v_\infty \sim 1000 \, {\rm km \, \, s}^{-1}$| and |$\dot{p}_w \sim \dot{p}_*$| while for the starburst model |$v_\infty \sim 10^4 \, {\rm km \, \, s}^{-1}$| and |$\dot{p}_w \sim 0.1 \dot{p}_*$|⁠. The large velocities and correspondingly lower momentum fluxes for the starburst model are a consequence of κ ≫ r₀c_i needed to avoid overproducing the gamma-ray luminosities.

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

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

4.3 Analytic scalings

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

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

(54)

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

$$\begin{eqnarray*} h_c &\simeq & \min \Biggr [1, \frac{0.9 \, E_{cr,50} v_{*,8.5}^{1/2}}{\alpha _{0.07} \, \phi ^{1/2}} \left(\frac{c_i}{10 \, {\rm km \, \, s}^{-1}}\right) \nonumber\\ && \times \left(\frac{\rm kpc}{r_0}\right) \left(\frac{100 \, {\rm km \, \, s}^{-1}}{V_g}\right)^2 \Biggr ] , \end{eqnarray*}$$

(55)

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

There are two regimes depending on whether h_c < 1 or h_c = 1 in equation (55). For h_c ∼ 1, which is the regime appropriate for dwarf galaxies and nuclear starbursts in our models in Figs 8 and 9,

$$\begin{eqnarray*} \kappa \simeq 3 \times 10^{29} \, {\rm cm^2 \, s^{-1}} \, \left(\frac{h_c E_{cr,50} v_{*,8.5}^{1/2}}{\alpha _{0.07} \, \phi ^{1/2}}\right) , \end{eqnarray*}$$

(56)

$$\begin{eqnarray*} \frac{p_{c,0}}{\pi G \Sigma _g^2 \phi } \simeq 0.03 \, \alpha _{0.07} \, \phi ^{1/2} v_{*,8.5}^{-3/2} \, m_{*,2}^{-1} \, \left(\frac{r_0}{3 \, {\rm kpc }}\right) \end{eqnarray*}$$

(57)

$$\begin{eqnarray*} \frac{\dot{M}_w}{\dot{M}_*} & \simeq& 0.1 \, \alpha _{0.07} \, \phi ^{1/2} \, v_{*,8.5}^{-1/2} \, m_{*,2}^{-1} \nonumber\\ &&\times \left(\frac{100 \, {\rm km \, \, s}^{-1}}{{\rm V_{g,eff}}}\right)^2 \, \left(\frac{c_i r_0}{h_c 30 \, {\rm km \, \, s}^{-1}\, {\rm kpc }}\right), \end{eqnarray*}$$

(58)

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

$$\begin{eqnarray*} v_\infty & \simeq& 1600 \, \, {\rm km \, \, s}^{-1}\, \alpha _{0.07}^{-1/2} \, E_{cr,50}^{1/2} \, \phi ^{-1/4} \, v_{*,8.5}^{1/4} \nonumber\\ && \times \left(\frac{{\rm V_{g,eff}}}{100 \, {\rm km s^{-1}}}\right) \left(\frac{\rm 30 \, h_c \, kpc \, {\rm km \, s^{-1}}}{r_0 c_i}\right)^{1/2} \end{eqnarray*}$$

(59)

Equations (58) and (59) then imply

$$\begin{eqnarray*} \frac{\dot{p}_w}{\dot{p}_{*}} & \simeq& 1 \, \frac{\alpha _{0.07}^{1/2} \, E_{cr,50}^{1/2} \, \phi ^{1/4}}{m_{*,2} \, \epsilon _{*,-3.3}\, v_{*,8.5}^{1/4}} \nonumber\\ && \times \left(\frac{100 \, {\rm km \, s^{-1}}}{{\rm V_{g,eff}}}\right) \left(\frac{r_0 c_i}{\rm 30 \, h_c \, kpc \, {\rm km \, s^{-1}}}\right)^{1/2} \end{eqnarray*}$$

(60)

In the opposite regime of h_c ≲ 1, which is the regime of star-forming disc galaxies in Fig. 8, we find

$$\begin{eqnarray*} \kappa & \simeq& 10^{29} \, {\rm cm^2 \, s^{-1}} \left(\frac{ E_{cr,50}^2 v_{*,8.5}}{\alpha _{0.07}^2 \, \phi } \right) \, \left(\frac{c_i}{10 \, {\rm km \, \, s}^{-1}}\right) \nonumber\\ && \times \left(\frac{100 \, {\rm km \, \, s}^{-1}}{V_g}\right)^2 \, \left(\frac{3 \, \, {\rm kpc }}{r_0}\right), \end{eqnarray*}$$

(61)

$$\begin{eqnarray*} \frac{p_{c,0}}{\pi G \Sigma _g^2 \phi } \simeq 0.03 \, \alpha _{0.07} \, \phi ^{1/2} v_{*,8.5}^{-3/2} \, m_{*,2}^{-1} \, \left(\frac{r_0}{3 \, {\rm kpc }}\right) \end{eqnarray*}$$

(62)

$$\begin{eqnarray*} \frac{\dot{M}_w}{\dot{M}_*} \simeq 0.2 \, \alpha ^2_{0.07} \, \phi \, v_{*,8.5} \, m_{*,2}^{-1} \, \left(\frac{r_0}{3 \, {\rm kpc }}\right)^2, \end{eqnarray*}$$

(63)

$$\begin{eqnarray*} v_\infty \simeq 1000 \, \, {\rm km \, \, s}^{-1}\, \alpha _{0.07}^{-1} \, E_{cr,50} \, \phi ^{-1/2} \, v_{*,8.5}^{1/2} \, \left(\frac{3 \, {\rm kpc }}{r_0}\right) \end{eqnarray*}$$

(64)

and

$$\begin{eqnarray*} \frac{\dot{p}_w}{\dot{p}_{*}} \simeq 1.3 \frac{\alpha _{0.07} \, E_{cr,50} \, \phi ^{1/2} \, v_{*,8.5}^{3/2}}{m_{*,2} \, \epsilon _{*,-3.3}} \, \left(\frac{r_0}{3 \, {\rm kpc }}\right) \end{eqnarray*}$$

(65)

We reiterate that ϕ in equations (55)–(65) quantifies the contribution of stars and dark matter to the gravitational acceleration with ϕ ∼ 1 for gas-dominated systems and ϕ ∼ 3−10 for typical star-forming galaxies in the local Universe. Also note that equations (57) and (62) are identical, i.e. the same for h_c < 1 and h_c = 1. This is because p_c depends only on the ratio H_c/κ (equation 52) which is uniquely determined by our calibration of t_π/t_diff (equations 48 and 54).

Equations (55)–(65) do a good job of reproducing many key results shown in Fig. 8: (1) the characteristic diffusion coefficient ∼10²⁸⁻²⁹ cm² s⁻¹ relatively independent of galaxy properties required to reproduce the observed gamma-ray-star formation rate correlation, and (2) a base CR pressure p_{c, 0} and wind mass-loading |$\dot{M}_w/\dot{M}_*$| that are relatively independent of star formation rate, but a strong function of the size of the star-forming disc r₀ at fixed star formation rate (because diffusion is much more rapid for smaller r₀ and decreases p_{c, 0} and |$\dot{M}_w$|⁠)..

Equations (55)–(65) also elucidate how our gamma-ray inferences depend on the gas isothermal sound speed. Fig. 8 takes |$c_i=10 \, {\rm km \, \, s}^{-1}$|⁠. Our star-forming disc model in Fig. 8 is in the regime with h_c ≲ 1, so that equations (61)–(65) apply. In this regime κ∝c_i so that κ/r₀c_i ∼ 1 from Fig. 8 is in fact true relatively independent of c_i. Moreover, for h_c ≲ 1, |$\dot{M}_w$| is independent of c_i (equation 63). The results in the middle and right-hand panels of Fig. 8 are thus applicable over a wide range of c_i (this is the reason that we do not plot the ‘disc’ model in Fig. 9). Our conclusion that CR-driven winds constrained by gamma-ray observations of disc galaxies are dynamically important with |$\dot{M}_w \sim \dot{M}_*$| is thus robust to uncertainties in c_i.

In contrast to the disc model, the dwarf and starburst models in Fig. 8 are in the regime where h_c ∼ 1 so that equations (56)–(60) apply. In this case the inferred κ is independent of c_i (equation 56); however, even |$c_i \sim 100 \, {\rm km \, \, s}^{-1}$| still corresponds to κ/r₀c_i ≳ 1 so that the qualitative physics of the wind does not change. Equations (58), (59), and (60) predict that the inferred mass-loss rate, terminal velocity, and momentum flux given gamma-ray constraints scale |$\propto c_i, c_i^{-1/2}$| and |$c_i^{1/2}$|⁠, respectively. This is consistent with, though weaker than, the trends in Fig. 9. The stronger dependence on c_i in Fig. 9 is because for larger values of c_i the approximation |$\left(4 h_c V_g^4/[c_{c,0}^2 c_i^2]\right)^{c_i^2/(2V_g^2)} \sim 1$| used to derive equations (56)–(65) is no longer as accurate. Physically, this is because for larger values of c_i and smaller c_{c, 0} the sonic point r_s is at somewhat larger radii ∼few ×r₀, which increases |$\dot{M}_w$| (equation 28) and decreases v_∞ (equation 35).

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

5 SUMMARY AND DISCUSSION

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

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

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

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

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

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

$Density, velocity, cosmic ray pressure profiles, and CR sound speed ($c_c = \sqrt{p_c/\rho }$) for our Vg = 200 simulation (see Table 1). For this plot, because of the very low base gas sound speed, we have normalized the velocity, CR sound speed, and and CR pressure using Vg, Vg, and $\rho V_g^2$, respectively. Note the onset of an instability and strong fluctuations at r ∼ 1.07 (the resolution decreases at r ≃ 1.14 due to a change in mesh refinement, which likely is responsible for suppressing the short wavelength fluctuations exterior to that radius). Despite the large density fluctuations, the mass-loss rate and terminal velocity in the simulation are well-described by our steady state analytic solutions.$

Figure A1.

Density, velocity, cosmic ray pressure profiles, and CR sound speed (⁠|$c_c = \sqrt{p_c/\rho }$|⁠) for our V_g = 200 simulation (see Table 1). For this plot, because of the very low base gas sound speed, we have normalized the velocity, CR sound speed, and and CR pressure using V_g, V_g, and |$\rho V_g^2$|⁠, respectively. Note the onset of an instability and strong fluctuations at r ∼ 1.07 (the resolution decreases at r ≃ 1.14 due to a change in mesh refinement, which likely is responsible for suppressing the short wavelength fluctuations exterior to that radius). Despite the large density fluctuations, the mass-loss rate and terminal velocity in the simulation are well-described by our steady state analytic solutions.

Open in new tab Download slide

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

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

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

One of the uncertain parameters in applying our results to observations is the appropriate isothermal gas sound speed. This depends on the phase of the ISM that the cosmic-rays most effectively couple to. Fig. 8 assumes |$c_i = 10 \, {\rm km \, \, s}^{-1}$|⁠, which is an appropriate mass-averaged value in the Milky Way. For typical star-forming disc galaxy parameters, we find that the properties of the winds using gamma-ray constrained diffusion coefficients are weakly dependent on c_i (also derived analytically in equations 63–65). Our conclusion that cosmic-rays are a significant source of winds in normal disc galaxies is thus reasonably robust to the uncertainty of the phase of the ISM that primarily determines CR transport. We note, however, that the CR diffusion coefficient we infer for disc galaxies using gamma-ray data scales linearly with the assumed isothermal sound speed. We find the best agreement with phenomenological models (⁠|$\kappa \sim 3 -10 \, 10^{28}$| cm² s⁻¹; e.g. Trotta et al. 2011; Evoli et al. 2020) for |$c_i \sim 100 \, {\rm km \, \, s}^{-1}$| (see Section 4.3) while |$c_i \sim 10 \, {\rm km \, \, s}^{-1}$| yields a significantly smaller diffusion coefficient (Fig. 8). This is consistent with the hypothesis that CR transport is set primarily by the hot volume-filling phase of the ISM. We also find that for star-forming disc galaxies, κ ∼ r₀c_i so that the winds are on the border between the rapidly accelerating high-speed solutions we find for κ ≳ r₀c_i and the slower more hydrostatic solutions obtained for κ ≲ r₀c_i.

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

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

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

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

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

ACKNOWLEDGEMENTS

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

DATA AVAILABILITY

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

Footnotes

1

The two-moment CR model we solve numerically in Section 3 evolves F_c as an independent variable and the flux reduces to equation (3) only when time variations are sufficiently slow (see equation 43).

2

Even given this restrictive assumption, the diffusion term in equation (3) in general depends on local plasma conditions. It can also depend on the cosmic ray energy density itself, in which case the ‘diffusion’ term in equation (3) is not even truly diffusive (e.g. Skilling 1971; Wiener, Oh & Guo 2013). We do not consider this complication in this work.

3

v_esc(r₀) is formally not defined for the ln (r) potential focused on here, but it is ∼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).

4

The two-moment approach enables a much more accurate treatment of CR streaming, which is not essential for this paper but will be in later publications in this series.

5

Jiang & Oh (2018) follow the conventions of the radiation transfer literature in which the diffusive flux is |$F_c = -\sigma _c^{-1} \nabla p_c$|⁠. In the rest of this paper, we follow the conventions of the CR literature and define the diffusive flux as F_c = −κ∇E_c (equation 2). This accounts for the factor of 3 relating κ and |$\sigma _c^{-1}$|⁠.

6

This implies that v(r_out) = 0.5V_g = 5c_i in Table 1 is a lower limit. If all of the CR enthalpy is converted to kinetic energy, the final velocity would be ≃ 1.8V_g, close to the analytic estimate in equation (35).

7

Note that there are small oscillations in the solution near the sonic point; these remain at late times even once the solution settles into a steady state. We do not have a definitive explanation for these oscillations but they are not present in any of our κ ≳ r₀c_i solutions (the majority in Table 1).

8

Note that ϵ_c defined in equation (30) is given by ϵ_c = E_cr/m_*c². In this section, it is convenient to separate out ϵ_c into a part that depends on the initial mass function (m_*) and a part that depends on the CR energy supplied per SN (E_cr).

9

Note that this expression for p_{c, 0} is larger than its spherical counterpart used in Section 2 (equation 9) by a factor of 2. The difference is the surface area of a sphere (⁠|$4 \pi r_0^2$|⁠) versus that of the bottom and top of the disc (⁠|$2 \times \pi r_0^2$|⁠). We use the disc version in this section but our results are not sensitive to factor of 2 changes in p_{c, 0}.

10

http://www.astropy.org

REFERENCES

Abdo

A. A.

et al. ,

2010a

,

A&A

,

523

,

L2

10.1051/0004-6361/201015759

Crossref

Search ADS

Abdo

A. A.

et al. ,

2010b

,

ApJ

,

709

,

L152

10.1088/2041-8205/709/2/L152

Crossref

Search ADS

Ackermann

M.

et al. ,

2012

,

ApJ

,

755

,

164

Crossref

Search ADS

Amato

E.

,

Blasi

P.

,

2018

,

Adv. Space Res.

,

62

,

2731

10.1016/j.asr.2017.04.019

Crossref

Search ADS

Andrews

B. H.

,

Weinberg

D. H.

,

Schönrich

R.

,

Johnson

J. A.

,

2017

,

ApJ

,

835

,

224

10.3847/1538-4357/835/2/224

Crossref

Search ADS

Astropy Collaboration

,

2013

,

A&A

,

558

,

A33

10.1051/0004-6361/201322068

Crossref

Search ADS

Astropy Collaboration

,

2018

,

AJ

,

156

,

123

10.3847/1538-3881/aabc4f

Crossref

Search ADS

Barcos-Muñoz

L.

et al. ,

2018

,

ApJ

,

853

,

L28

10.3847/2041-8213/aaa28d

Crossref

Search ADS

Bell

E. F.

,

2003

,

ApJ

,

586

,

794

10.1086/367829

Crossref

Search ADS

Booth

C. M.

,

Agertz

O.

,

Kravtsov

A. V.

,

Gnedin

N. Y.

,

2013

,

ApJ

,

777

,

L16

10.1088/2041-8205/777/1/L16

Crossref

Search ADS

Boulares

A.

,

Cox

D. P.

,

1990

,

ApJ

,

365

,

544

10.1086/169509

Crossref

Search ADS

Breitschwerdt

D.

,

McKenzie

J. F.

,

Voelk

H. J.

,

1991

,

A&A

,

245

,

79

Buckman

B. J.

,

Linden

T.

,

Thompson

T. A.

,

2020

,

MNRAS

,

494

,

2679

10.1093/mnras/staa875

Crossref

Search ADS

Chan

T. K.

,

Kereš

D.

,

Hopkins

P. F.

,

Quataert

E.

,

Su

K. Y.

,

Hayward

C. C.

,

Faucher-Giguère

C. A.

,

2019

,

MNRAS

,

488

,

3716

10.1093/mnras/stz1895

Crossref

Search ADS

Crocker

R. M.

,

Krumholz

M. R.

,

Thompson

T. A.

,

2021

,

MNRAS

,

503

,

2651

10.1093/mnras/stab502

Crossref

Search ADS

Drury

L. O.

,

Falle

S. A. E. G.

,

1986

,

MNRAS

,

223

,

353

10.1093/mnras/223.2.353

Crossref

Search ADS

Everett

J. E.

,

Zweibel

E. G.

,

Benjamin

R. A.

,

McCammon

D.

,

Rocks

L.

,

Gallagher III

J. S.

,

2008

,

ApJ

,

674

,

258

10.1086/524766

Crossref

Search ADS

Evoli

C.

,

Morlino

G.

,

Blasi

P.

,

Aloisio

R.

,

2020

,

Phys. Rev. D

,

101

,

023013

10.1103/PhysRevD.101.023013

Crossref

Search ADS

Griffin

R. D.

,

Dai

X.

,

Thompson

T. A.

,

2016

,

ApJ

,

823

,

L17

10.3847/2041-8205/823/1/L17

Crossref

Search ADS

Hopkins

P. F.

et al. ,

2020

,

MNRAS

,

492

,

3465

10.1093/mnras/stz3321

Crossref

Search ADS

Hopkins

P. F.

,

Chan

T. K.

,

Ji

S.

,

Hummels

C. B.

,

Kereš

D.

,

Quataert

E.

,

Faucher-Giguère

C.-A.

,

2021a

,

MNRAS

,

501

,

3640

10.1093/mnras/staa3690

Crossref

Search ADS

Hopkins

P. F.

,

Squire

J.

,

Chan

T. K.

,

Quataert

E.

,

Ji

S.

,

Kereš

D.

,

Faucher-Giguère

C.-A.

,

2021b

,

MNRAS

,

501

,

4184

10.1093/mnras/staa3691

Crossref

Search ADS

Hunter

J. D.

,

2007

,

Comput. Sci. Eng.

,

9

,

90

10.1109/MCSE.2007.55

Crossref

Search ADS

Ipavich

F. M.

,

1975

,

ApJ

,

196

,

107

10.1086/153397

Crossref

Search ADS

Ji

S.

et al. ,

2020

,

MNRAS

,

496

,

4221

10.1093/mnras/staa1849

Crossref

Search ADS

Jiang

Y.-F.

,

Oh

S. P.

,

2018

,

ApJ

,

854

,

5

10.3847/1538-4357/aaa6ce

Crossref

Search ADS

Krause

M.

et al. ,

2018

,

A&A

,

611

,

A72

10.1051/0004-6361/201731991

Crossref

Search ADS

Kulsrud

R.

,

Pearce

W. P.

,

1969

,

ApJ

,

156

,

445

10.1086/149981

Crossref

Search ADS

Lacki

B. C.

,

Thompson

T. A.

,

2013

,

ApJ

,

762

,

29

10.1088/0004-637X/762/1/29

Crossref

Search ADS

Lacki

B. C.

,

Thompson

T. A.

,

Quataert

E.

,

2010

,

ApJ

,

717

,

1

10.1088/0004-637X/717/1/1

Crossref

Search ADS

Lacki

B. C.

,

Thompson

T. A.

,

Quataert

E.

,

Loeb

A.

,

Waxman

E.

,

2011

,

ApJ

,

734

,

107

10.1088/0004-637X/734/2/107

Crossref

Search ADS

Linden

T.

,

2017

,

Phys. Rev. D

,

96

,

083001

10.1103/PhysRevD.96.083001

Crossref

Search ADS

Linden

T.

,

Profumo

S.

,

Anderson

B.

,

2010

,

Phys. Rev. D

,

82

,

063529

10.1103/PhysRevD.82.063529

Crossref

Search ADS

Lochhaas

C.

,

Thompson

T. A.

,

Schneider

E. E.

,

2021

,

MNRAS

,

504

,

3412

10.1093/mnras/stab1101

Crossref

Search ADS

Mao

S. A.

,

Ostriker

E. C.

,

2018

,

ApJ

,

854

,

89

10.3847/1538-4357/aaa88e

Crossref

Search ADS

Muratov

A. L.

,

Kereš

D.

,

Faucher-Giguère

C.-A.

,

Hopkins

P. F.

,

Quataert

E.

,

Murray

N.

,

2015

,

MNRAS

,

454

,

2691

10.1093/mnras/stv2126

Crossref

Search ADS

Murray

N.

,

Quataert

E.

,

Thompson

T. A.

,

2005

,

ApJ

,

618

,

569

10.1086/426067

Crossref

Search ADS

Ostriker

E. C.

,

Shetty

R.

,

2011

,

ApJ

,

731

,

41

10.1088/0004-637X/731/1/41

Crossref

Search ADS

Owocki

S. P.

,

Gayley

K. G.

,

1997

, in

Nota

A.

,

Lamers

H.

, eds,

ASP Conf. Ser. Vol. 120, The Physics of Stellar Winds Near the Eddington Limit

.

Astron. Soc. Pac

,

San Francisco

, p.

121

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Owocki

S. P.

,

Townsend

R. H. D.

,

Quataert

E.

,

2017

,

MNRAS

,

472

,

3749

10.1093/mnras/stx2251

Crossref

Search ADS

Pakmor

R.

,

Pfrommer

C.

,

Simpson

C. M.

,

Springel

V.

,

2016

,

ApJ

,

824

,

L30

10.3847/2041-8205/824/2/L30

Crossref

Search ADS

Pavlidou

V.

,

Fields

B. D.

,

2001

,

ApJ

,

558

,

63

10.1086/322455

Crossref

Search ADS

Pohl

M.

,

1994

,

A&A

,

287

,

453

Recchia

S.

,

Blasi

P.

,

Morlino

G.

,

2016

,

MNRAS

,

462

,

4227

10.1093/mnras/stw1966

Crossref

Search ADS

Ruszkowski

M.

,

Yang

H. Y. K.

,

Zweibel

E.

,

2017

,

ApJ

,

834

,

208

10.3847/1538-4357/834/2/208

Crossref

Search ADS

Salem

M.

,

Bryan

G. L.

,

2014

,

MNRAS

,

437

,

3312

10.1093/mnras/stt2121

Crossref

Search ADS

Simpson

C. M.

,

Pakmor

R.

,

Marinacci

F.

,

Pfrommer

C.

,

Springel

V.

,

Glover

S. C. O.

,

Clark

P. C.

,

Smith

R. J.

,

2016

,

ApJ

,

827

,

L29

10.3847/2041-8205/827/2/L29

Crossref

Search ADS

Skilling

J.

,

1971

,

ApJ

,

170

,

265

10.1086/151210

Crossref

Search ADS

Socrates

A.

,

Davis

S. W.

,

Ramirez-Ruiz

E.

,

2008

,

ApJ

,

687

,

202

10.1086/590046

Crossref

Search ADS

Somerville

R. S.

,

Davé

R.

,

2015

,

ARA&A

,

53

,

51

10.1146/annurev-astro-082812-140951

Crossref

Search ADS

Stone

J. M.

,

Tomida

K.

,

White

C. J.

,

Felker

K. G.

,

2020

,

ApJS

,

249

,

4

10.3847/1538-4365/ab929b

Crossref

Search ADS

Strong

A. W.

,

Porter

T. A.

,

Digel

S. W.

,

Jóhannesson

G.

,

Martin

P.

,

Moskalenko

I. V.

,

Murphy

E. J.

,

Orlando

E.

,

2010

,

ApJ

,

722

,

L58

10.1088/2041-8205/722/1/L58

Crossref

Search ADS

Thompson

T. A.

,

Lacki

B. C.

,

2013

,

Astrophys. Space Sci. Proc. Vol. 34, The FIR-Radio Correlation in Rapidly Star-Forming Galaxies: The Spectral Index Problem and Proton Calorimetry

.

Springer-Verlag

,

Berlin

, p.

283

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Thompson

T. A.

,

Quataert

E.

,

Murray

N.

,

2005

,

ApJ

,

630

,

167

10.1086/431923

Crossref

Search ADS

Thompson

T. A.

,

Quataert

E.

,

Waxman

E.

,

2007

,

ApJ

,

654

,

219

10.1086/509068

Crossref

Search ADS

Torres

D. F.

,

2004

,

ApJ

,

617

,

966

10.1086/425415

Crossref

Search ADS

Trotta

R.

,

Jóhannesson

G.

,

Moskalenko

I. V.

,

Porter

T. A.

,

Ruiz de Austri

R.

,

Strong

A. W.

,

2011

,

ApJ

,

729

,

106

10.1088/0004-637X/729/2/106

Crossref

Search ADS

Uhlig

M.

,

Pfrommer

C.

,

Sharma

M.

,

Nath

B. B.

,

Enßlin

T. A.

,

Springel

V.

,

2012

,

MNRAS

,

423

,

2374

10.1111/j.1365-2966.2012.21045.x

Crossref

Search ADS

Veilleux

S.

,

Maiolino

R.

,

Bolatto

A. D.

,

Aalto

S.

,

2020

,

A&AR

,

28

,

2

10.1007/s00159-019-0121-9

Crossref

Search ADS

Wiener

J.

,

Oh

S. P.

,

Guo

F.

,

2013

,

MNRAS

,

434

,

2209

10.1093/mnras/stt1163

Crossref

Search ADS

Wiener

J.

,

Pfrommer

C.

,

Oh

S. P.

,

2017

,

MNRAS

,

467

,

906

10.1093/mnras/stx127

Crossref

Search ADS

Yan

H.

,

Lazarian

A.

,

2002

,

Phys. Rev. Lab.

,

89

,

281102

10.1103/PhysRevLett.89.281102

Crossref

Search ADS

Yoast-Hull

T. M.

,

Everett

J. E.

,

Gallagher

J. S. I.

,

Zweibel

E. G.

,

2013

,

ApJ

,

768

,

53

10.1088/0004-637X/768/1/53

Crossref

Search ADS

Yoast-Hull

T. M.

,

Gallagher

J. S. I.

,

Zweibel

E. G.

,

Everett

J. E.

,

2014

,

ApJ

,

780

,

137

10.1088/0004-637X/780/2/137

Crossref

Search ADS

APPENDIX A: LINEAR STABILITY

In this Appendix, we study the linear stability of the CR magnetohydrodynamic equations. Since the simulations are 1D, we restrict ourselves to 1D perturbations and also consider a local Cartesian approximation instead of the global spherical geometry used in Section 3. Physically, this system of equations admits longitudinal sound waves, in which both gas and CR pressure are the restoring force, as well as gas and CR entropy modes. In what follows, we show that ignoring background gradients, the two-moment CR system is linearly stable in the presence of cosmic ray diffusion. There is, however, an instability driven by a background cosmic ray pressure gradient that is present in the one-moment CR system (Drury & Falle 1986), i.e. the instability does not rely on the finite speed of light. We show, however, that this instability grows too slowly to be dynamically important in galactic winds, consistent with the laminar numerical solutions we found in Section 3. The only exception to this is if the gas sound speed is very low, as we show in Fig. A1.

A1 Instabilities of the two-moment system

We assume here that perturbations are ∝exp (− iωt + ikr) and that kH ≫ 1 (where H is a characteristic length-scale in the equilibrium state) so that a WKB analysis is appropriate. For now, we neglect the background gradients in the problem.

The key frequencies in the problem are the isothermal gas sound wave frequency

$$\begin{eqnarray*} \omega _g = k c_i, \end{eqnarray*}$$

(A1)

the adiabatic CR sound wave frequency

$$\begin{eqnarray*} \omega _c = k c_{\rm eff} \equiv k \sqrt{4 p_c/3 \rho }, \end{eqnarray*}$$

(A2)

the CR diffusion frequency associated with the assumed constant diffusion coefficient κ

$$\begin{eqnarray*} \omega _d = k^2 \kappa \end{eqnarray*}$$

(A3)

and a characteristic frequency in the problem due to the finite speed of light, which we define as

$$\begin{eqnarray*} \omega _M = \frac{v_M^2}{\kappa } . \end{eqnarray*}$$

(A4)

Note that in the simulations described in Section 3, κ ∼ 1 − 30c_ir so that ω_d/ω_g ∼ kr(κ/c_ir) ≫ 1. The same inequality holds for ω_d/ω_c. By contrast, |$\omega _d/\omega _M \sim \kappa ^2 k^2/v_M^2 \sim (k \ell)^2 \lt 1$| is required for the fluid approximation to the CR dynamics to be valid (where ℓ here is the CR mean free path).

Working in the WKB limit, the 1D linear dispersion relation for equations (43) is given by

$$\begin{eqnarray*} 0 &= & \frac{3 \omega ^4}{\omega _M} + i \omega ^3\left(1 + 3 \frac{\omega _c^2}{\omega _d \omega _M}\right) - \omega ^2 \left(\omega _d + 3 \frac{\omega _g^2}{\omega _M}\right) \nonumber\\ &&- \,\, i \omega \left(\omega _g^2 + \omega _c^2\right) + \omega _g^2 \omega _d \end{eqnarray*}$$

(A5)

In the rapid diffusion limit of ω_d ≫ ω_g, ω_c and ω_M ≫ ω_d, the 4 solutions to equation (A5) are all stable:

$$\begin{eqnarray*} \omega &\simeq&-i \, \frac{\omega _M}{3} \nonumber\\ \omega &\simeq& -i \, {\omega _d}\nonumber\\ \omega &\simeq& \pm \omega _g -i \, \frac{\omega _c^2}{2 \omega _d} \end{eqnarray*}$$

(A6)

The first two solutions in equation (A6) are strongly damped entropy modes. The last is a weakly damped gas sound wave. Physically, the latter wave arises because in the limit ω_d → ∞, CR pressure gradients are completely wiped out by diffusion and the only restoring force for a sound wave is the gas pressure. At finite ω_d, there is a small residual CR pressure gradient, the diffusion of which leads to damping of the associated sound wave.

We reiterate that the rapid diffusion ordering used to derive equation (A6) is the appropriate one for our simulations in Section 3. The absence of any growing modes in equation (A6) is consistent with the numerical solutions which find laminar wind solutions.

In the limit of slow CR diffusion, ω_d ≪ ω_g, ω_c ≪ ω_M the solutions of equation (A5) are also damped, namely

$$\begin{eqnarray*} \omega &\simeq& -i \, \frac{\omega _M}{3} \nonumber\\ \omega &\simeq& -i \, {\omega _d} \frac{\omega _g^2}{\omega _g^2 + \omega _c^2} \nonumber\\ \omega &\simeq& \pm \sqrt{\omega _g^2 + \omega _c^2} -i \, \frac{\omega _d}{2} \frac{\omega _c^2}{\omega _g^2 + \omega _c^2} \end{eqnarray*}$$

(A7)

A2 Instabilities of the one-moment CR system with background gradients

Instabilities of the one-moment CR system for a homogeneous background can be derived using the results in Section A1 by taking v_M → ∞. The sound and entropy modes are both stable in this limit. Including background gradients in the calculation, however, leads to an instability of sound waves that was discussed by Drury & Falle (1986). We briefly summarize a derivation of this instability for completeness and then discuss its relevance to our galactic wind simulations. The Drury & Falle (1986) instability is present in the one-moment CR system and so we restrict our analysis to this limit for ease of algebra.

We consider an isothermal gas plus CR system that satisfies the following conservation laws

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t} + \frac{\mathrm{ d} \rho v}{\mathrm{ d} z} = 0 \end{eqnarray*}$$

(A8)

$$\begin{eqnarray*} \rho \frac{\partial v}{\partial t} + \rho v \frac{\mathrm{ d} v}{\mathrm{ d} z} = - c_i^2 \frac{\mathrm{ d} \rho }{\mathrm{ d} z} - \frac{\mathrm{ d} p_c}{\mathrm{ d} z} - \rho g \end{eqnarray*}$$

(A9)

We linearize equations (A8) and (A9). To start we assume that all perturbations, labelled by δ, are ∝exp (− iωt) but we do not Fourier transform in z. We do the latter only at the end of the calculation to ensure that all background gradient terms are properly kept. The linearly perturbed equations are then

$$\begin{eqnarray*} i \omega \delta \rho = \frac{\mathrm{ d} (\rho \delta v)}{\mathrm{ d} z} \end{eqnarray*}$$

(A10)

$$\begin{eqnarray*} -i \omega \rho \delta v = -c_i^2 \frac{\mathrm{ d} \delta \rho }{\mathrm{ d} z} - \frac{\mathrm{ d} \delta p_c}{\mathrm{ d} z} - \delta \rho g \end{eqnarray*}$$

(A11)

Equations (A10) and (A11) can be combined to yield

$$\begin{eqnarray*} \omega ^2 \delta \rho = -c_i^2 \frac{\mathrm{ d}^2 \delta \rho }{\mathrm{ d}^2 z} - \frac{\mathrm{ d}^2 \delta p_c}{\mathrm{ d}^2 z} - g \frac{\mathrm{ d} \delta \rho }{\mathrm{ d} z} \end{eqnarray*}$$

(A12)

In the limit of rapid CR diffusion, the linearized CR energy equation with diffusion (equations 2 and 3) simply becomes |$\kappa \, \mathrm{ d}^2 \delta p_c/\mathrm{ d}^2 z \simeq 0$|⁠. Substituting this into equation (A12), assuming perturbations ∝exp [ikz − z/(2H)], where H is the density scale-height, and using hydrostatic equilibrium in the background yields

$$\begin{eqnarray*} \omega = |k c_i| + i \frac{ c_{\rm c}^2}{2 c_i} \frac{k}{|k|} \frac{\mathrm{ d} \ln p_c}{\mathrm{ d}z} \end{eqnarray*}$$

(A13)

to O(1/H) (⁠|$c_{\rm c}^2 = p_c/\rho$| as in the main text). Equation (A13) is equivalent to the dispersion relation in Drury & Falle (1986) in the limit of rapid CR diffusion. Drury & Falle (1986) further show that the rapid diffusion approximation leading to equation (A13) only applies if κ ≳ 4/3|dln p_c/dr|⁻¹c_i; otherwise the system is stable. Our lowest κ simulations in Table 1 with κ/r₀c_i = 0.33, 0.11 are stable at most radii per this condition; otherwise, the rapid diffusion approximation is a good one in our simulations.

The number of e-foldings for the Drury & Falle (1986) instability can be estimated as A(r) ≃ Im(ω)H_ρ/c_i where H_ρ is the density scale-height on which the background structure changes and the flow accelerates. Using equations (A13) and (27), we find

$$\begin{eqnarray*} A(r) \sim \sqrt{\frac{r_0 c_i}{\kappa }}\frac{c_c^2}{c_i V_g} \end{eqnarray*}$$

(A14)

near the base of the outflow where the instability derivation is applicable. For our fiducial simulation with κ ∼ 10r₀c_i and c_c ≃ c_i ≃ 0.1V_g near the base, we find A ∼ 0.03, i.e. very little growth of the instability. This is consistent with our laminar numerical simulations. Fundamentally, the reason for this is that the CR pressure gradient that drives the Drury & Falle (1986) instability is very shallow in galactic winds driven by CR diffusion, with a CR pressure scale-height much larger than the density scale-height in the subsonic portion of the wind at small radii where equation (A13) applies (see Fig. 2). The large CR pressure scale-height in the present context means that the growth of the Drury & Falle (1986) instability is slow and is the key reason why nearly all of our simulations do not show any sign of this linear instability.

From equation (A14), the Drury & Falle (1986) instability is most likely to grow when c_i is small and/or V_g is large (gravity is strong), both of which decrease the CR scale-height (equation 27). Indeed, we find that our simulation with the smallest value of the gas isothermal sound speed does show evidence of an instability. In this case (the first row in Table 1), we predict A(r) ≃ 0.5 near the base of the wind, and a somewhat larger value at the sonic point where c_c is larger. Fig. A1 shows that there is indeed evidence of an instability that sets in at r ∼ 1.07 in this simulation. This may be a manifestation of the Drury & Falle (1986) instability. However, the instability in Fig. A1 sets in at radii well exterior to the sonic point and even exterior to where the flow speed equals the CR sound speed. We suspect that these are fluctuations generated by the Drury & Falle (1986) instability at small radii and advected out to large radii where they become non-linear due to conservation of wave action. Despite the large density fluctuations, however, the wind mass-loss rate and terminal velocity in this simulations are still well-described by the analytic solution in Section 2. We note that the resolution of the simulation in Fig. A1 decreases at r ≃ 1.14 due to a change in mesh refinement, which likely is responsible for suppressing the short wavelength fluctuations exterior to that radius.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
November 2021	3
December 2021	9
January 2022	29
February 2022	10
March 2022	10
April 2022	12
May 2022	3
June 2022	10
July 2022	7
August 2022	11
September 2022	14
October 2022	8
November 2022	23
December 2022	6
January 2023	13
February 2023	29
March 2023	10
April 2023	22
May 2023	10
June 2023	14
July 2023	20
August 2023	17
September 2023	7
October 2023	11
November 2023	16
December 2023	30
January 2024	13
February 2024	9
March 2024	21
April 2024	23
May 2024	26
June 2024	26
July 2024	20
August 2024	20
September 2024	27
October 2024	25
November 2024	15
December 2024	14
January 2025	28
February 2025	23
March 2025	22
April 2025	22
May 2025	5

Article Contents

The physics of galactic winds driven by cosmic rays I: Diffusion

ABSTRACT

1 INTRODUCTION

2 ANALYTIC APPROXIMATIONS FOR COSMIC RAY-DRIVEN WINDS WITH DIFFUSION

2.1 Wind energetics, momentum flux, and velocity

2.2 Maximum mass-loss rate

3 NUMERICAL SIMULATIONS

3.1 Equations

3.2 Initial and boundary conditions

3.3 Overview of simulations

3.4 Numerical solutions for κ ≳ r₀c_i

3.5 Low κ ≲ r₀c_i solutions

3.6 Synthesis of analytics versus numerics

4 IMPLICATIONS OF GAMMA-RAY OBSERVATIONS FOR COSMIC-RAY WIND MODELS

4.1 Non-thermal emission from pion decay

4.2 Models calibrated to gamma-ray data

4.3 Analytic scalings

5 SUMMARY AND DISCUSSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: LINEAR STABILITY

A1 Instabilities of the two-moment system

A2 Instabilities of the one-moment CR system with background gradients

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

The physics of galactic winds driven by cosmic rays I: Diffusion Free

ABSTRACT

1 INTRODUCTION

2 ANALYTIC APPROXIMATIONS FOR COSMIC RAY-DRIVEN WINDS WITH DIFFUSION

2.1 Wind energetics, momentum flux, and velocity

2.2 Maximum mass-loss rate

3 NUMERICAL SIMULATIONS

3.1 Equations

3.2 Initial and boundary conditions

3.3 Overview of simulations

3.4 Numerical solutions for κ ≳ r0ci

3.5 Low κ ≲ r0ci solutions

3.6 Synthesis of analytics versus numerics

4 IMPLICATIONS OF GAMMA-RAY OBSERVATIONS FOR COSMIC-RAY WIND MODELS

4.1 Non-thermal emission from pion decay

4.2 Models calibrated to gamma-ray data

4.3 Analytic scalings

5 SUMMARY AND DISCUSSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: LINEAR STABILITY

A1 Instabilities of the two-moment system

A2 Instabilities of the one-moment CR system with background gradients

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

The physics of galactic winds driven by cosmic rays I: Diffusion

3.4 Numerical solutions for κ ≳ r₀c_i

3.5 Low κ ≲ r₀c_i solutions