How Supernovae Launch Galactic Winds

TL;DR

This study uses high-resolution hydrodynamic simulations to investigate how supernovae drive galactic winds, revealing the importance of supernova remnant radii, clustering, and analytic modeling in understanding wind properties.

Contribution

The paper introduces a self-consistent simulation approach resolving supernova remnant radii and presents an analytic model explaining wind launching based on SNR breakout conditions.

Findings

Supernovae launch highly supersonic galactic winds consistent with analytic expectations.

Wind energy loading converges with resolution, but mass loading decreases at higher resolutions.

Clustering supernovae significantly enhances wind power.

Abstract

We use idealized three-dimensional hydrodynamic simulations of global galactic discs to study the launching of galactic winds by supernovae (SNe). The simulations resolve the cooling radii of the majority of supernova remnants (SNRs) and thus self-consistently capture how SNe drive galactic winds. We find that SNe launch highly supersonic winds with properties that agree reasonably well with expectations from analytic models. The energy loading ($\eta_E = \dot{E}_{\rm wind} / \dot{E}_{\rm SN}$) of the winds in our simulations are well converged with spatial resolution while the wind mass loading ($\eta_M = \dot{M}_{\rm wind} / \dot{M}_\star$) decreases with resolution at the resolutions we achieve. We present a simple analytic model based on the concept that SNRs with cooling radii greater than the local scale height breakout of the disc and power the wind. This model successfully explains the dependence (or lack thereof) of $\eta_E$ (and by extension $\eta_M$) on the gas surface density, star formation efficiency, disc radius, and the clustering of SNe. The winds in the majority of our simulations are weaker than expected in reality, likely due to the fact that we seed SNe preferentially at density peaks. Clustering SNe in time and space substantially increases the wind power.

Figures

Tables

Table 1. Table 1: Full range of simulation parameters. Not all combinations are discussed in the text. The fiducial models are in bold. M ∘ subscript 𝑀 M_{\circ} is the mass of the external gravitational potential. R d subscript 𝑅 𝑑 R_{d} is the characteristic radius of the gravitational potential and the scale length of the exponential gas disc. Σ gas subscript Σ gas \Sigma_{\rm gas} is central gas surface density. The parameter f cl subscript 𝑓 cl f_{\rm cl} is the clustering factor of SNe: the total energy per SNR is equal to f cl × 10 51 subscript 𝑓 cl superscript 10 51 f_{\rm cl}\times 10^{51} ergs while the SN rate ∝ f cl − 1 proportional-to absent superscript subscript 𝑓 cl 1 \propto f_{\rm cl}^{-1} . The efficiency of star formation is encapsulated in f ⋆ = t ⋆ / t dyn subscript 𝑓 ⋆ subscript 𝑡 ⋆ subscript 𝑡 dyn f_{\star}=t_{\star}/t_{\rm dyn} , such that a lower f ⋆ subscript 𝑓 ⋆ f_{\star} yields more star formation.

$M_{\circ }$$\left[M_{\odot }\right]$	$R_d$[pc]	${\Sigma }_{\mathrm{gas}}$$\left[M_{\odot }{\mathrm{pc}}^{-2}\right]$	$f_{\mathrm{cl}}$	$f_{\star }/100$
$1.0\times {10}^9$	100	10, 30, 50	1	0.3, 1, 3
$2.8\times {10}^9$	300	10, 30, 100, 300	1, 3, 10, 30	0.3, 1, 3
$9.2\times {10}^9$	1000	10, 30	1	1

Equations

$$P(\mathrm{SN}\mathrm{in}\mathrm{cell})=\frac{M_{\mathrm{cell}}}{100{\mathrm{M}}_{\odot }}\frac{d{t}_{\mathrm{hydro}}}{t_{\star }}\propto \frac{n_{\mathrm{cell}}}{t_{\star }}=\frac{n_{\mathrm{cell}}}{f_{\star }{t}_{\mathrm{dyn}}}.$$

$${\eta }_{\mathrm{E}}=\frac{{\dot{E}}_{\mathrm{wind}}}{{\dot{E}}_{\mathrm{SN}}}\sim \frac{n_{\mathrm{crit}}}{n_{\mathrm{mid}}}\propto h^{-3/2}\text{leavevmode}{f}_{\mathrm{cl}}^{5/7}\text{leavevmode}{\Sigma }_{\mathrm{gas}}^{-1}.$$

$${\eta }_{\mathrm{E}}\propto R_d^{-3/2}{\left(\frac{\delta v}{v_{\mathrm{circ}}}\right)}^{-3/2}{f}_{\mathrm{cl}}^{5/7}\text{leavevmode}{\Sigma }_{\mathrm{gas}}^{-1}.$$

Full text

††pagerange: How Supernovae Launch Galactic Winds–References††pubyear: 2017

1 Introduction

Galactic winds help limit the efficiency with which galaxies turn gas into stars by expelling material from the ISM and by halting gas inflow into galaxies (e.g., Dekel & Silk, 1986; Springel & Hernquist, 2003; Faucher-Giguère et al., 2011). They are also responsible for enriching and heating the CGM (e.g., Aguirre et al., 2001; Oppenheimer & Davé, 2006; Hummels et al., 2013; Fielding et al., 2017). As a result galactic winds are at the heart of many of the cornerstone relationships of modern astronomy such as the stellar mass function, stellar mass to halo mass relation, and the mass-metallicity relation. Many processes are capable of launching galactic winds, but in star forming galaxies energy deposition by SNe is often thought to be a key driver.

Constraints on the nature of star formation powered galactic winds come from numerous sources. First, extensive observations have directly measured the energy and mass loading of galactic winds in different environments (e.g., Heckman et al., 1990; Veilleux et al., 2005; Chisholm et al., 2017). Second, analytic considerations predict the density, temperature, and velocity profiles of a galactic wind for a given energy and mass loading (e.g., Chevalier & Clegg 1985 hereafter CC85, Thompson et al. 2016). Third, cosmological simulations demonstrate that winds with a particular range of efficiencies are required to reproduce many observations (e.g., Finlator & Davé, 2008; Somerville & Davé, 2015; Muratov et al., 2015). Given all we know about galactic winds there is nonetheless a surprising degree of disagreement on if/how SNe are capable of launching winds that meet all the necessary constraints.

Numerical simulations of isolated galaxies inform how SNe drive winds—commonly using local stratified box simulations (e.g., Joung & Mac Low, 2006; Creasey et al., 2013; Girichidis et al., 2016; Li et al., 2016; Kim & Ostriker, 2016). These stratified box simulations generally predict winds that are subsonic, which may be a result of the geometry, namely the lack of a well-defined escape speed and free-space for the wind to expand into (Martizzi et al. 2016 hereafter M16). To more faithfully address how SNe launch galactic winds we designed a new suite of simulations that adopts a global geometry, capturing an entire gaseous galactic disc while resolving most SNRs. Because we neglect self-gravity, molecular line cooling, and other important physics these are not the final word on the true energy and mass loading of SNe-driven winds. But they do significantly sharpen our understanding of the origin and properties of such winds.

2 Method

We ran a series of simulations designed to study the launching of galactic winds by SNe in a global geometry using the Eulerian hydrodynamics code athena (Stone et al., 2008). We evolve a gaseous galactic disc that is stratified by an external potential—representing the gravitational field from baryons and dark matter—in which intermittent, discrete SNe go off at a given rate. The setup of the numerical experiment is simple and provides a useful counterpoint to analogous experiments that differ essentially only in their use of local Cartesian simulation domains; we compare primarily to our earlier work (M16). The simulations evolve an ideal fluid in three dimensions with cooling and without self-gravity. The gas has solar metallicity everywhere and cooling proceeds assuming collisional ionization equilibrium. Cooling below $10^{4}$ K and photoelectric heating are not included. The gas is initialized in a rotating, $10^{4}$ K disc that is in radial centripetal balance and vertical hydrostatic equilibrium with the background gravitational potential that is given by a Hernquist profile, $\Phi(r)=-{GM_{\circ}}/({r+a})$, where $r$ is the spherical radius ($R$ represents the cylindrical radius) (Hernquist, 1990). The disc initially has an exponentially declining surface density profile characterized by a central surface density $\Sigma_{\rm gas}$ and a scale length $R_{d}$ that is set to match the characteristic radius of the gravitational potential $a$.

The simulation parameters are listed in Table 1. In all simulations the parameters are chosen so that $v_{\rm circ}(r=a=R_{d})=\sqrt{GM(r)/r}=100$ km s*-1*. For our fiducial simulations we adopt the relatively small disc scale length of $R_{d}=300$ pc, although we study the differences resulting from using a ${\sim}3\times$ smaller and larger disc. Our disc sizes and $v_{\rm circ}$ are similar to those in M82 and NGC 253, but our $\Sigma_{\rm gas}$ are not quite as high.

The major drawback of running global rather than local simulations is the computational expense of resolving the SNR injection when we are primarily interested in the large-scale structure and wind dynamics. This partially motivates the smaller disc scale-lengths considered here, so that SNR evolution can be resolved. To further this goal, we use five levels of static-mesh-refinement. We space our cubical $256^{3}$ grids logarithmically, enabling us to use a domain side length of $L_{\rm box}\gtrsim 40R_{d}$ and have high resolution (3 pc) in the disc mid-plane where most of the mass resides (the highest resolution region is $L_{\rm HR}\gtrsim 2R_{d}$ on a side).

We seed and inject SNRs in our simulations as in the “SC” model used by M16, so we refer the reader there for more details. In short, the probability of a SN being set off is proportional to the local gas density and inversely proportional to the local star formation timescale, which is chosen to be proportional to the dynamical time, so that

[TABLE]

This assumes that for every 100 ${\rm M}_{\odot}$ of stars that form there is one SN. Our fiducial choice of $f_{\star}$ is 100, which corresponds to a 1% star formation efficiency and results in star formation rate surface densities $\dot{\Sigma}_{\star}$ that are similar to observations.

SNRs are injected using the subgrid model developed by Martizzi et al. (2015), which accounts for subgrid cooling and injects both kinetic and thermal energy at a value calibrated to high resolution single SNR simulations. Additionally, 3 ${\rm M}_{\odot}$ of ejecta is added to the SNR per SN, so the ISM mass loading ($\dot{M}_{\rm ej}/\dot{M}_{\star}$) is 0.03. One of the primary aims of this study is to determine wind mass and energy loss rates when the SNRs’ cooling radii ($r_{\rm cool}$) are explicitly resolved, so we chose parameters to ensure that this occurs for our higher resolution simulations. In this limit our SNe injection model corresponds to $2.9\times 10^{50}$ ergs of kinetic energy and $7.1\times 10^{50}$ ergs of thermal energy. One new feature we added to the injection scheme relative to M16 is a somewhat crude model for the clustering of SNe in space and time (future work will expand this feature). We allow the injected energy per SN to be scaled up by an integer clustering factor $f_{\rm cl}$, which represents multiple SNe going off simultaneously (the SN rate is correspondingly reduced by $f_{\rm cl}$, so that the total injected energy by SNe is unchanged). The cooling radius and other radii in the Martizzi et al. (2015) fits for subgrid injection are scaled up by $f_{\rm cl}^{2/7}$ in accordance with analytic expectations (Cioffi et al., 1988).

3 Results

We begin our presentation of the simulation results with a qualitative description to ground the readers’ intuitions. In Fig. 1 we show density, temperature, and spherical radial velocity images from simulations with $R_{d}=300$ pc, $\Sigma_{\rm gas}=10\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$, $f_{\star}=100$ and both $f_{\rm cl}=1$ and $f_{\rm cl}=30$, after 300 Myr of evolution (${\sim}16\leavevmode\nobreak\ t_{\rm orb}$). These images are slices through the computational domain along the rotation axis of the disc. Clearly shown are the strong biconical outflows driven by the SN. Along the rotation axis of the disc densities are low, temperatures are high, and velocities reach upwards of 300 km s*-1*($\sim 500$ km s*-1in the clustered model). In the midplane of the disc densities remain high, the temperature remains at roughly $10^{4}$ K, and the gas is turbulent with a mass-weighted velocity of ${\sim}10$s km s-1*. These images show several supernova remnants breaking out of the disc. These breakout events carve out a region of the disc and dump thermal energy into the low density wind region, thereby powering the wind. The discrete breakouts lead to an inhomogeneous outflow composed of a series of hot, dense, and fast fronts of material that are trailed by gas which has expanded, cooled, and slowed down. This is reminiscent of what is observed in the M82 wind.

In Fig. 2 we show the time averaged radial profiles of $T$, $n$, $v_{r}$, $c_{s}$, $\eta_{\rm E}$, and $\eta_{\rm M}$ for the fiducial $\Sigma_{\rm gas}{=}10\,{\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$ simulation. The averages are volume-weighted and computed in a biconical region centered on the disc rotation axis with a half opening angle of $45^{\circ}$. Several features are immediately apparent. The outflow is supersonic. The mass and energy outflow rates are roughly independent of radius beyond a certain point, indicating that we have a steady state outflow. The density of the wind material falls off rapidly, and the temperature decreases with radius slower than expected for just adiabatic expansion, which predicts $T\propto r^{-4/3}$. We have omitted the profiles from simulations with other parameters because they are all sufficiently similar and show the same trends.

Comparing the profiles in Fig. 2 to those from local Cartesian box simulations demonstrates the simulation geometry’s effect on the wind structure. For example, Fig. 4 of M16 shows that the sound speed and velocity of the wind are roughly independent of height beyond the scale height of the disc and the flow is always subsonic. Additionally, Fig. 9 and B1 of M16 show that $\eta_{\rm M}$ decreases dramatically with distance from the disc whereas we find $\eta_{\rm M}$ to be roughly constant with radius. It is, therefore, the ratio of wind thermal to kinetic energy and the fraction of wind material that escapes that are primarily affected by the simulation geometry (important quantities for galactic winds!).

Standard theoretical arguments for the structure of a galactic wind of a given $\eta_{\rm E}$, $\eta_{\rm M}$, and $\dot{\Sigma}_{\star}$ assume spherical symmetry and a uniform injection of energy and mass (e.g., CC85, Thompson et al. 2016). Nevertheless, a comparison to the analytic work is instructive. A generic prediction of these models is that the wind will be supersonic beyond a sonic point that is approximately the radius of the star forming region (unless cooling is too strong). Our simulations agree with this prediction very well as can be seen in Fig. 2 where the sonic point is at ${\sim}R_{d}/3$. CC85 predict the asymptotic velocity to be $v_{\infty}\approx 10^{3}$ km s*-1*$(\eta_{\rm E}/\eta_{\rm M})^{1/2}$ and the temperature at the sonic point to be $T\approx 2\times 10^{7}$K $\eta_{\rm E}/\eta_{\rm M}$, which agrees strikingly well with our simulations. We indicate the best fit power law slope for the $T$, $n$, and $v_{r}$ profiles at large radii in Fig. 2 for comparison to observations and analytic models. Interestingly, $T$ falls off slower with distance than is expected for adiabatic expansion ($T\propto r^{-4/3}$), which could be due to cooling at small radii, additional heating beyond $R_{d}$ either by the rare (but effective) distant SNe, or by internal shocks. The biconical $n$ profile falls off roughly as $r^{-2}$ as expected for a freely expanding constant $\dot{M}$ wind. When averaging over a spherical region the $n$ profile falls off much more quickly as $r^{-\alpha}$ with $\alpha{\sim}4$-$5$ for different simulations. The steeper fall off of the spherically averaged profile indicates there may be some fall back of wind material on the side of the biconical outflows creating a fountain flow (see Fig. 1). This steeper slope is also consistent with the values inferred for local starburst galaxy M 82 (Leroy et al., 2015).

In Fig. 3 we show the time evolution of $\eta_{\rm E}$ and $\eta_{\rm M}$ for the fiducial $\Sigma_{\rm gas}{=}10$ and $100\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$ simulations at the highest and half the highest central spatial resolution. The energy and mass outflow rates of the wind are measured at 4 kpc. After the initial transient the outflow settles into a steady state (quite different from local Cartesian box simulations; e.g., M16 Fig. 9 and B1). We performed extensive resolution testing on our simulations and found that with even half the resolution $\eta_{\rm E}$ is converged. However, as we increased the resolution $\eta_{\rm M}$ continued to decrease. This is likely due to more mixing from unresolved shear layers with worse resolution, and more efficient venting at high resolution. The convergence of $\eta_{\rm E}$ does not depend sensitively on the degree to which we resolve the SNRs’ $r_{\rm cool}$. In our calculations, a resolved SNR has $r_{\rm cool}{>}r_{\rm inj}{=}2\Delta x$. In the fiducial $\Sigma_{\rm gas}=10\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$ simulation with $\Delta x=3,6,$ and 12 pc 97.3, 73.5, and 31.8 per cent of the SNR are resolved, respectively, and yet $\eta_{\rm E}$ is the same. In the fiducial $\Sigma_{\rm gas}=100\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$ simulation with $\Delta x=3,6,$ and 12 pc 49.0, 12.1, and 1.4 per cent of the SNRs are resolved, respectively. In these higher surface density simulations $\eta_{\rm E}$ is roughly the same for the two higher resolution cases, but at $\Delta x=12$ pc, $\eta_{\rm E}$ drops by a factor of ${\sim}2$. Overall, because $\eta_{\rm M}$ measured at large radii depends sensitively on the structure of the disc and the properties of the surrounding circumgalactic medium swept up by the wind (Sarkar et al., 2015) we consider $\eta_{\rm M}$ of secondary importance relative to $\eta_{\rm E}$. Nevertheless, the lack of $\eta_{\rm M}$ convergence is something that should be considered whenever simulations similar to ours are compared to observations.

A standard physical picture of models such as CC85 is that winds are launched when the volume filling fraction of hot gas in the disc is large enough so that the number of SNe per cooling time and cooling volume—known as the porosity $Q_{c}$ (McKee & Ostriker, 1977)—is greater than 1. Martizzi et al. (2015) found that $r_{\rm cool}\approx 20.8\leavevmode\nobreak\ \mathrm{pc}\leavevmode\nobreak\ {n_{H}}^{-2/5}\leavevmode\nobreak\ f_{\rm cl}^{2/7}$, and $t_{\rm cool}\approx 2.9\times 10^{4}$ yrs ${n_{H}}^{-0.54}$. Therefore, the porosity of the disc is $Q_{c}=(4\pi/3)r_{\rm cool}^{3}t_{\rm cool}\dot{n}_{SN}=6\times 10^{-5}(n_{H}/100\leavevmode\nobreak\ \mathrm{cm}^{-3})^{-4/5}(100/f_{\star})(10^{6}\leavevmode\nobreak\ \mathrm{yrs}/t_{\rm dyn})$. In all of the simulations carried out here $Q_{c}\ll 1$. As can be seen in Fig. 1, the SNe that contribute to the launching of the wind are the (rare) ones whose SNRs are able to breakout of the disc before radiating away their energy. Working under the assumption that SNRs that breakout satisfy $r_{\rm cool}\gtrsim h$, where $h$ is the local scale height of the gaseous disc, we now provide a simple argument for the expected scaling of the wind properties with the disc and injection properties. There is a critical hydrogen number density $n_{\mathrm{crit}}=(h/20.8\mathrm{\leavevmode\nobreak\ pc})^{-5/2}f_{\rm cl}^{5/7}\leavevmode\nobreak\ \mathrm{cm}^{-3}$ that satisfies $r_{\rm cool}=h$. Since, by design, the SN rate is proportional to the local density, the fraction of SNe that satisfy $n_{H}\leq n_{\mathrm{crit}}$ should be roughly equal to the ratio of $n_{\mathrm{crit}}$ to the midplane density $n_{\rm mid}\leavevmode\nobreak\ {\sim}\leavevmode\nobreak\ \Sigma_{\rm gas}/2hm_{p}$ (this need not be true for models with different SN seeding schemes or that account for additional physics). Neglecting further radiative loses post breakout, the fraction of the injected energy that goes into the wind should be equal to the same ratio, which yields

[TABLE]

The total energy injection rate is $\dot{E}_{\rm SN}{=}10^{51}\mathrm{\leavevmode\nobreak\ ergs\leavevmode\nobreak\ }\pi R_{d}^{2}\Sigma_{\rm gas}/(t_{\star}\,100{\rm M}_{\odot}),$ which follows from the definition of $t_{\star}$ and our model’s assumption that there are $10^{51}$ ergs released per SN and there is one SN per $100{\rm M}_{\odot}$ of stars formed. Combining this expression with equation (2) we find that $\dot{E}_{\rm wind}=\eta_{\rm E}\dot{E}_{\rm SN}\propto R_{d}^{2}\,f_{\rm cl}^{5/7}/(h^{3/2}\,t_{\star}).$ One particularly interesting feature of this expression is the lack of any dependence on $\Sigma_{\rm gas}$. Finally, the ratio of the scale height to the disc radius is approximately equal the ratio of the velocity dispersion to the circular velocity ${h}/{R_{d}}\sim{\delta v}/{v_{\rm circ}}$, which is $\sim 0.1$ in our model. With this assumption in hand we end up with the following expected scalings for the wind energetics

[TABLE]

The scalings above have been derived for $\eta_{\rm E}$ because the condition $r_{\rm cool}\gtrsim h$ is explicitly a statement about energetics, but similar reasoning can be applied to mass-loading. Assuming that each SNR that breaks out contributes to the wind all of the mass it swept up prior breakout, $M_{\rm swept}=(4\pi/3)r_{\rm cool}^{3}\rho\approx 1200\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ n_{\rm crit}^{-1/5}$, yields $\eta_{\rm M}\approx(n_{\rm crit}/n_{\rm mid})(M_{\rm swept}\,\dot{N}_{\rm SN}/\dot{M}_{\star})=12\,\eta_{\rm E}\,n_{\rm crit}^{-1/5}$. Therefore, we expect $\eta_{\rm M}$ to scale similarly to $\eta_{\rm E}$ and to be one to two orders of magnitude larger than $\eta_{\rm E}$.

The top panel of Fig. 4 demonstrates that, for all other properties being equal, $\eta_{\rm E}$ is inversely dependent on $\Sigma_{\rm gas}$ exactly as predicted by equation (3). Increasing the star formation efficiency by decreasing $f_{\star}$ leads to an increase in $\dot{E}_{\rm wind}$, but no appreciable change in $\eta_{\rm E}$ as expected. Likewise, equation (3) captures roughly the correct behavior for the scaling of $\eta_{\rm E}$ with $R_{d}$ seen in Fig. 4, which demonstrates that more compact systems launch more powerful winds.

The bottom panel of Fig. 4 shows how $\eta_{\rm E}$ and $\eta_{\rm M}$ scale with $f_{\rm cl}$ for $R_{d}=300$ pc discs with $\Sigma_{\rm gas}=10$ and $30\leavevmode\nobreak\ {\rm M}_{\odot}\leavevmode\nobreak\ {\rm pc}^{-2}$. In these simulations $\eta_{\rm E}$ increases roughly linearly with $f_{\rm cl}$—somewhat more strongly than predicted in equation (2). Regardless of the exact scaling, the strong $f_{\rm cl}$ dependence of $\eta_{\rm E}$ may be critical for understanding the launching of real galactic winds that are powerful enough to match observations and satisfy the requirements from cosmological simulations. This is because stars—massive stars in particular—are expected to form in clusters and it is likely that their SNRs will be nested or overlap rather than being spatially and temporally separated as in our $f_{\rm cl}=1$ simulations. Although $\eta_{\rm M}$ also increases with $f_{\rm cl}$, the measured scaling is less robust due to the lack of $\dot{M}_{\rm wind}$ convergence. Nevertheless, it is worth noting that for nearly all cases $\eta_{\rm M}$ surpasses 0.03—the value corresponding to all wind material coming from SNe ejecta.

4 Discussion & Conclusion

Using idealized global galactic disc simulations we have quantified the properties of galactic winds driven solely by SNe for a range of disc and star formation properties. We have focused on small discs (${\sim}0.1{-}1$ kpc in size) in order to ensure that the cooling radii of most of the SNRs in our simulations can be resolved. Our simulations roughly reproduce the supersonic wind structure expected from analytic models (e.g., CC85, ). Previous simulations that attempted to study galactic winds launched by SNe in a stratified medium often adopted local Cartesian domains (with periodic and outflow boundary conditions in the disc plane and perpendicular to it, respectively) and found subsonic outflows with outflow rates that depend on box height (e.g., Martizzi et al., 2016; Girichidis et al., 2016; Kim & Ostriker, 2016). The more physical global geometry we adopt allows the winds to adiabatically expand causing them accelerate to supersonic velocities, and the gravitational potential with a well-defined escape velocity leads to outflows with radially constant $\dot{M}_{\rm wind}$ and $\dot{E}_{\rm wind}$ (Fig. 3). Other numerical models for studying galactic winds inject a fixed $\dot{E}$ uniformly in a given volume (Strickland & Heckman, 2009; Sarkar et al., 2015). Our calculations compliment these by addressing the key question of how discrete SNe collectively drive a wind.

In analytic galactic wind models the mass ($\eta_{\rm M}$) and energy ($\eta_{\rm E}$) loading of the wind are free parameters. Our simulations determine these wind properties as a function of the underlying disc structure (e.g., $\Sigma_{\rm gas}$ and $R_{d}$) and the SN seeding model (e.g., degree of clustering $f_{\rm cl}$). In our simulations the winds are driven by SNe that go off in low density regions where the cooling radius $r_{\rm cool}$ is larger than the local scale height $h$; this enables SNRs to drive the wind without radiative losses sapping their energy. In general only a small fraction of the SNe satisfy this constraint because at the disc midplane $r_{\rm cool}\ll h$. We present a simple model based on this concept (see eq. 3) that predicts, among other things, that $\dot{E}_{\rm wind}$ should be independent of $\Sigma_{\rm gas}$, increase with the degree of SNe clustering $f_{\rm cl}$ and the star formation efficiency $f_{\star}$, and decrease with increasing disc size. This simple model successfully explains many of the trends we find in our simulations (Fig. 4). Although this analytic model and the numerical scaling of wind efficiency with disc and SNe parameters are likely to be somewhat modified with different SNe seeding schemes, we expect that the general trends found here are likely to be more robust—in particular the simple criterion that $r_{\rm cool}\gtrsim h$ for the SNe that drive the wind.

The mass and energy loading of the galactic winds driven by SNe we find are likely lower than suggested by observations. This may be due to the fact that our SNe are set off preferentially at density peaks and that the ISM is relatively homogeneous. A more realistic (or a spatially random) SNe seeding scheme that separates the SN locations from density peaks and/or clustering the SNe would increase the efficiency of the outflows (e.g., Sharma et al., 2014; Girichidis et al., 2016; Kim & Ostriker, 2016; Gentry et al., 2017). Indeed, in our calculations, we implemented a simple model of clustering in which each SNR’s energy is increased by a factor of $f_{\rm cl}$, and the SN rate decreased by the same factor, leaving the total energy input rate the same. The resulting galactic wind energy loss rate increases roughly linearly with $f_{\rm cl}$ (Fig. 4). The wind power might well be further increased if additional physics were included such as molecular line cooling (Li et al., 2016) and stronger ISM turbulence possibly enhanced by self-gravity and/or galactic inflows (Sur et al., 2016). These would result in larger density inhomogeneities, causing more of the ISM volume to be filled with low density gas, and thereby allowing more SNe to go off in low density regions and breakout of the disc.

Acknowledgements

We thank Prateek Sharma for useful conversations during the development of this work. This work was supported in part by NASA ATP grant 12-APT12-0183 and a Simons Investigator award from the Simons Foundation to EQ. DF was supported by the NSF GRFP under Grant # DGE 1106400. CAFG was supported by NSF through grants AST-1412836 and AST-1517491, and by NASA through grant NNX15AB22G. DM was supported by the Swiss National Science Foundation as an Advanced Postdoc Mobility Fellow until Nov. 2016; grant number P300P2_161062.

This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at UC, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor of Research, and Office of the CIO). In addition, this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575, via grant number TG-AST160020.