Numerical Simulations of Multiphase Winds and Fountains from Star-forming Galactic Disks. I. Solar Neighborhood TIGRESS Model

In galactic disks, star formation takes place in dense gas near the midplane, within the scale height of the ISM. Prodigious energy is injected by SNe within this thin layer, and the high-entropy, overpressured gas expands outward. Strong shocks heat and accelerate both dense, cold cloudlets and the warm, diffuse intercloud medium surrounding individual SNRs and superbubbles, with most of this gas cooling back to its original temperature relatively quickly (e.g., Kim & Ostriker 2015a; Kim et al. 2017). Depending on the level of remaining specific energy with respect to the gravitational potential, outflows of cooled, SN-accelerated warm (or cold) gas may either keep moving out of the disk, or may turn around at some height. While most of the energy deposited by SNe is radiated away, some of the hot gas created in strong SN shocks is at low enough density that it has very limited cooling. This accelerates away from the midplane towards higher-latitude, lower-pressure regions, achieving high enough velocity that it can escape from the galactic potential well. In this paper, the term "galactic winds" refers to outflows launched with high enough energy to escape the galactic gravitational potential, while the term "galactic fountains" refers to outflows launched with insufficient energy that eventually fall back.

In our simulations (and in real galaxies), the motions of gas are three-dimensional. Within any given patch of the ISM in the disk, followed on its orbit about the galactic center, horizontal averaging can be used to define a mean density, velocity, and other flow properties as a function of height z. In general, the instantaneous mean velocity will have both horizontal (radial-azimuthal) and vertical components. The horizontally averaged flow quantities may be further averaged over time (with window comparable to an orbit time, so that epicyclic motions are averaged away). If the accretion rate is low, the resulting temporally averaged velocity at any height will be dominated by vertical motion. Thus, horizontal and temporal averaging of the gas yields an effectively one-dimensional profile as a function of height, consisting of average values $\langle \rho (z)\rangle$, $\langle {v}_{z}(z)\rangle$, $\langle P(z)\rangle$, and so on. If there is a net outflow, $\mathrm{sign}({v}_{z})=\mathrm{sign}(z)$, and if there is a net inflow, the sign is reversed.

Similar to horizontal- and time-averaged flows, classical SN-driven wind solutions are one-dimensional, and the simplest solutions are also adiabatic. For steady-state one-dimensional adiabatic gas flows, the equations of mass and total energy conservation, including source terms, can be written as

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})={\dot{\rho }}_{\mathrm{inj}}(z),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{ \mathcal B })={\dot{e}}_{\mathrm{inj}}(z),\end{eqnarray} \tag{ 2 }$$

where ${\dot{\rho }}_{\mathrm{inj}}$ and ${\dot{e}}_{\mathrm{inj}}$ are the volumetric mass and energy injection rates, respectively, arising from SN feedback. In Equations (1) and (2), z represents the vertical direction for a flow perpendicular to the plane of a galactic disk; an approximately spherical galactic center flow (e.g., Chevalier & Clegg 1985) would instead have ${\dot{\rho }}_{\mathrm{inj}}(r)$ and ${\dot{e}}_{\mathrm{inj}}(r)$. Note that ${\dot{\rho }}_{\mathrm{inj}}$ does not represent SN ejecta itself. Rather, shock heating of ambient ISM gas near SNe increases the mass injection rate above that of the initial SN ejecta, while cooling tends to reduce the rate. Allowing for this shock heating and cooling, ${\dot{\rho }}_{\mathrm{inj}}$ is simply the mean local rate at which hot material is added the steady flow. Similarly, ${\dot{e}}_{\mathrm{inj}}$ represents the mean local energy input rate to the flow, which is bounded above by the initial energy carried by SN ejecta.

In Equation (2), the total specific energy (the Bernoulli parameter) is defined by

$$\begin{eqnarray}&&{ \mathcal B }\equiv \displaystyle \frac{{v}^{2}}{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }+{\rm{\Phi }},\end{eqnarray} \tag{ 3 }$$

consisting of specific kinetic energy, specific enthalpy $h\,=[\gamma /(\gamma -1)]P/\rho$, and gravitational potential terms. Note that here, for simplicity, we consider adiabatic, unmagnetized gas, but Appendix A presents the full equations for general MHD flows in a local shearing box, and shows that horizontal- and time-averaging yields a set of simple steady-state 1D flow equations, which can be applied to our simulations. In this paper, the ${\rm{\Phi }}=0$ reference point is at the midplane.

For flows in a local Cartesian box, like ours, with energy sources near the midplane, any time-averaged steady winds that may exist are launched vertically along the $\hat{z}$ axis. Through Gauss's law, volume integration of Equations (1) and (2) gives mass and total energy fluxes through surface area $A={L}_{x}{L}_{y}$ at $z=\pm d$ as

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{\mathrm{wind}}(d)\equiv {F}_{M}(z=d)-{F}_{M}(z=-d)={\dot{M}}_{\mathrm{inj}}(d)/A\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{{ \mathcal F }}_{E}(d)\equiv {F}_{E}(z=d)-{F}_{E}(z=-d)={\dot{E}}_{\mathrm{inj}}(d)/A.\end{eqnarray} \tag{ 5 }$$

Here, d is the distance from the disk midplane, and total mass and energy injection rates within $| z| \lt d$ are ${\dot{M}}_{\mathrm{inj}}(d)\,={\int }_{-d}^{d}{\dot{\rho }}_{\mathrm{inj}}{dV}$ and ${\dot{E}}_{\mathrm{inj}}(d)={\int }_{-d}^{d}{\dot{e}}_{\mathrm{inj}}{dV}$. The quantities ${F}_{M}(z)\,=\langle \rho {v}_{z}\rangle$ and ${F}_{E}(z)=\langle \rho {v}_{z}{ \mathcal B }\rangle$ stand for mass and total energy fluxes averaged over a horizontal area at height z.³ The above relations assume periodic boundary conditions in both of the horizontal directions; energy terms associated with the background shear from integrals over faces perpendicular to $\hat{x}$ are discussed in Appendix A.

Since SN explosions are usually concentrated within a thin layer near the midplane, ${\dot{M}}_{\mathrm{inj}}(d)$ and ${\dot{E}}_{\mathrm{inj}}(d)$ are expected to approach constant values for $d\gg H$, where H is the disk scale height. Hence, if A is constant, as in the local Cartesian coordinates, for a steady wind the (areal) mass and energy fluxes ${\dot{{\rm{\Sigma }}}}_{\mathrm{wind}}$ and ${{ \mathcal F }}_{E}$ would also approach constants at $d\gg H$ above the source region. For general geometries (e.g., Chevalier & Clegg 1985 for spherical coordinates), however, where the area that encloses a given set of streamlines varies as a function of distance (e.g., A ∝ r² for the spherical case), the fluxes would also vary (e.g., ∝r⁻²).

To the extent that horizontal correlations (or anticorrelations) among variables at a given z may be neglected, ${F}_{E}={F}_{M}{ \mathcal B }$, and if we may assume symmetry across the midplane this yields ${ \mathcal B }={{ \mathcal F }}_{E}/{\dot{{\rm{\Sigma }}}}_{\mathrm{wind}}={\dot{E}}_{\mathrm{inj}}/{\dot{M}}_{\mathrm{inj}}$. This implies that beyond the source region where ${\dot{E}}_{\mathrm{inj}}$ and ${\dot{M}}_{\mathrm{inj}}$ have reached their final values, the Bernoulli parameter becomes constant for a steady (or time-averaged) flow. In the more general case with A an increasing function of distance, because the area perpendicular to streamlines is the same for both mass and energy flux, F_E and F_M would vary ∝A⁻¹, but ${ \mathcal B }$ would still be conserved along streamlines (beyond the source region), irrespective of geometry. Therefore, for steady adiabatic winds (or equivalently for any time-averaged, adiabatic portion in a more general outflow), the Bernoulli parameter is a key quantity that enables robust extrapolation of flow evolution to large distances.⁴ As applied to the present problem, this suggests that evaluation of ${ \mathcal B }$ with our local disk simulations should provide predictions for properties of the hot portion of the wind at large distance (outside our simulation domain) where wind streamlines open up (becoming more radial than vertical). This motivates the need to quantify mass and energy fluxes in the launching region—just above the source region—in our simulations.

The simple analysis above provides helpful intuition for gas flows driven by localized energy sources, but the real ISM—and our simulations—is not a single phase, adiabatic gas. In fact, material in the ISM spans a wide range of density and temperature. SN shocks are responsible for generating the hot gas phase ($T\sim {10}^{6}-{10}^{7}\,{\rm{K}}$), which interacts with surrounding warm ($T\sim {10}^{4}\,{\rm{K}}$) and cold ($T\sim {10}^{2}\,{\rm{K}}$) phases. Considering each thermal component individually, radiative heating and cooling as well as mass and energy transfers between components would act as source and sink terms in the conservation equations of each phase. Because SN ejecta strongly interact with the surrounding ISM in extremely complex ways to heat and accelerate gas, some of which may be able to escape from a galaxy, it is not at all obvious how one would estimate ${\dot{M}}_{\mathrm{inj}}$ and ${\dot{E}}_{\mathrm{inj}}$ for individual thermal components of a multiphase outflow. Moreover, star formation and hence SN events are very bursty, and this burstiness may affect yields. Clearly, the total and individual-phase ${\dot{M}}_{\mathrm{inj}}$ and ${\dot{E}}_{\mathrm{inj}}$ are only quantifiable with self-consistent numerical simulations that capture the full physics of the ISM. Nevertheless, while simulations are essential for obtaining the detailed properties of realistic multiphase outflows, we can still expect certain aspects of classical wind solutions to hold when suitably applied. In particular, as we shall show, the spacetime-averaged hot portion of the wind, when considered separately from other phases, shares many similarities with adiabatic one-dimensional winds.

In Paper I, we presented a novel framework for multi-physics numerical simulations of the star-forming ISM implemented in the Athena MHD code (Stone et al. 2008; Stone & Gardiner 2009). We solve the ideal MHD equations in a local, shearing box, representing a small patch of a differentially rotating galactic disk. This treatment allows us to achieve uniformly high spatial resolution compared to what is possible in a global simulation of an entire galaxy or galaxies (e.g., Hopkins et al. 2012; Muratov et al. 2015), which is crucial to resolve both star formation and SN feedback as well as all thermal phases of the ISM both near and far from the midplane. We include gaseous and (young) stellar self-gravity, a fixed external gravitational potential to represent the old stellar disk and dark matter halo, galactic differential rotation, and optically thin cooling and grain photoelectric heating. We utilize sink particles to follow formation of and accretion onto star clusters in dense, cold gas. Massive young stars in these star clusters feed energy back to the ISM, by emitting far-ultraviolet radiation (FUV) and exploding as SNe. The former heats the diffuse warm and cold ISM, while the latter creates hot ISM gas, drives turbulence, and induces outflows.

Our simulations yield realistic, fully self-consistent three-phase ISM models with self-regulated star formation.⁵ Paper I presented results from a fiducial model with parameters similar to those of the Solar neighborhood. There, we showed that after $t\sim 100\,\mathrm{Myr}$, a quasi-steady state is reached. When stars form, massive stars enhance heating in the warm and cold ISM, and the SN rate increases, driving turbulence throughout the ISM. Both feedback processes disperse dense clouds and puff the gas disk up, temporarily shutting star formation off. With the corresponding reduction in star formation feedback, the gas can settle back to the midplane and collect into dense clouds, which then form a new generation of stars. In Paper I, we evaluated several basic ISM and wind properties, and demonstrated their convergence as a function of the numerical spatial resolution (see also Appendix B.1).

In this paper, we analyze detailed properties of galactic winds and fountains for a vertically extended version of the fiducial model of Paper I. The simulation domain covers $1\,\mathrm{kpc}\times 1\,\mathrm{kpc}$ horizontally and $-4.5\,\mathrm{kpc}\lt z\lt 4.5\,\mathrm{kpc}$ vertically, at a uniform resolution ${\rm{\Delta }}x=4\,\mathrm{pc}$. Representative snapshots displaying a volume rendering of density and velocity vectors during outflow- and inflow-dominated periods are shown in Figure 1. Figure 2 displays slices of temperature and vertical velocity through the y = 0 plane for the same snapshots shown in Figure 1. The outflows and inflows seen in Figures 1 and 2 are part of an overall cycle that repeats, representing the response to large amplitude temporal fluctuations in star formation rates (see the gray line in Figure 3(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To characterize vertical gas flows, we first construct horizontally averaged quantities. We then calculate the net mass and energy fluxes that pass through horizontal planes (both upper and lower sides) at d = 1, 2, and $3\,\mathrm{kpc}$ (cf. Equations (4) and (5)). Figure 3 shows time evolution of (a) mass flux, along with the areal star formation rate (${{\rm{\Sigma }}}_{\mathrm{SFR}}$), and (b) energy flux. Every star formation burst is followed by a burst of energy injection, and this burstiness is reflected in large temporal variations in the mass and energy fluxes. The fluxes can become negative, meaning that the mass and/or energy of gas flowing inward exceeds that of the gas flowing outward, at a given height. As distance d from the midplane increases, the net mass flux significantly decreases, while the net energy flux even through d = 2 and $3\,\mathrm{kpc}$ remains large. While net negative mass fluxes (implying fallback) occur at $d=1\,\mathrm{kpc}$ after each burst of mass outflow, net energy fluxes almost always remain positive. The differing behavior of mass and energy fluxes is a signature of multiphase flows.

Although star formation and hence SN feedback are impulsive rather than continuous, the system approaches a quasi-equilibrium state. This state is a limit cycle mediated by the feedback loop, in which epochs of cooling and collapse alternate with epochs of heating and expansion.⁶ Given that a quasi-steady state exists in the present simulation, horizontal- and temporal averages can be constructed to characterize this mean state. Since different thermal phases coexist at all heights, to understand the outflows and inflows of mass and energy, it is further necessary to separately construct horizontal averages of each thermal phase.

In Figure 4 we plot horizontally and temporally averaged profiles of mass, momentum, and energy distributions for thermally separated gas phases; these profiles average over $t=250\mbox{--}500\,\mathrm{Myr}$ and also average over upper ($z\gt 0$) and lower ($z\lt 0$) sides. Profiles show hydrogen number density, ${n}_{{\rm{H}}}\equiv \rho /({\mu }_{{\rm{H}}}{m}_{{\rm{H}}})$, outward vertical momentum density, $\rho {v}_{\mathrm{out}}\equiv \rho {v}_{z}\mathrm{sign}(z)$, and total energy density (excluding gravity):

$$\begin{eqnarray}&&{ \mathcal E }\equiv \displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{P}{\gamma -1}+\displaystyle \frac{{B}^{2}}{8\pi }.\end{eqnarray} \tag{ 6 }$$

Note that this energy density differs from the gas density multiplied by the Bernoulli parameter of Equation (3), as it includes just thermal energy (rather than enthalpy) and also includes magnetic energy density.

Zoom In Zoom Out Reset image size

The four thermal phases plotted in Figure 4 are cold ($T\lt 5050\,{\rm{K}}$), warm ($5050\,{\rm{K}}\lt T\lt 2\times {10}^{4}\,{\rm{K}}$), ionized ($2\times {10}^{4}\,{\rm{K}}\lt T\lt 5\times {10}^{5}\,{\rm{K}}$), and hot ($T\gt 5\times {10}^{5}\,{\rm{K}}$).⁷ Above the warm/cold layer ($d\gt H$, where $H\approx 400\,\mathrm{pc}$), the mass density is dominated by the warm component and the energy density is dominated by the hot component. As the individual terms in the energy density are proportional to corresponding terms in the momentum flux ($\rho {v}^{2}$, P), the hot medium also dominates the momentum flux away from the midplane. The hot medium is the largest contributor to the time-averaged vertical momentum density, which is the same as the time-averaged net mass flux. Although the mean value of vertical momentum density (or net mass flux) of the warm medium is effectively zero, there are large temporal fluctuations (indicated by the green shaded region) at small and intermediate d because the warm gas contributes significantly to both outgoing and incoming mass fluxes at different times (as in Figure 3).

The phase-separated momentum and energy density profiles in Figure 4 (and corresponding profiles of mass and momentum flux) reflect essential differences of gas flow dynamics between the warm and hot phases, which we separately analyze in the following sections.

In this section, we discuss key quantities of gas outflows driven by SNe, the mass and energy loading factors. The mass loading factor, β, through the surfaces at d (including both sides of the disk plane) is conventionally defined by the ratio of "outgoing" mass flux to the star formation rate as

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{{\dot{{\rm{\Sigma }}}}_{\mathrm{wind}}^{+}}{{{\rm{\Sigma }}}_{\mathrm{SFR}}},\end{eqnarray} \tag{ 9 }$$

while the energy loading factor, α, is defined by the ratio of "outgoing" kinetic + thermal (enthalpy) energy flux to the energy production rate of SNe,

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{{{ \mathcal F }}_{\mathrm{KE}}^{+}+{{ \mathcal F }}_{\mathrm{TE}}^{+}}{{E}_{\mathrm{SN}}{{\rm{\Sigma }}}_{\mathrm{SFR}}/{m}_{* }},\end{eqnarray} \tag{ 10 }$$

where ${E}_{\mathrm{SN}}={10}^{51}\,\mathrm{erg}$ is the total energy per SN, and ${m}_{* }=95.5\,{M}_{\odot }$ is the total mass of new stars per SN (see Paper I). To compute ${\dot{{\rm{\Sigma }}}}_{\mathrm{wind}}^{+}$, ${{ \mathcal F }}_{\mathrm{KE}}^{+}$, and ${{ \mathcal F }}_{\mathrm{TE}}^{+}$, we select only zones with outflowing gas (i.e., with $\mathrm{sign}({v}_{z})=\mathrm{sign}(z)$). The areal star formation rate averaged over $t=250\mbox{--}500\,\mathrm{Myr}$ is ${{\rm{\Sigma }}}_{\mathrm{SFR}}=0.006\,{M}_{\odot }\,{\mathrm{kpc}}^{-2}\,{\mathrm{yr}}^{-1};$ we use this average ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ in computing all loading factors. Note that a rolling mean of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ with 100 Myr time window gives 50% variation with respect to the mean. Because there is generally a temporal offset between star formation bursts and winds (see Figure 3), the instantaneous ratio of mass or energy fluxes to the star formation rate can significantly over- or under-estimate the true physical loading. The ratio of time-averaged wind fluxes to the time-averaged star formation rate is more meaningful.

We decompose the gas into four thermal components and present the loading factors of each thermal phase as a function of d (Figure 9). The energy flux is always dominated by the hot gas, with the energy loading factor of ${\alpha }_{h}\sim 0.01\mbox{--}0.05$ above $d\gt 1\,\mathrm{kpc}$ (as before, the "h" subscript denotes the hot phase). SN feedback causes large outgoing mass fluxes of warm and cold gas within the disk scale height $d\lt H$, but the majority of the warm and cold medium has low velocity and cannot travel far from the midplane; this mass flux at small d is best thought of as the "upwelling" of turbulent motions within the disk driven by expanding SNRs and superbubbles.¹⁰ The mass loading factor of the ionized gas also decreases significantly with d. By the time the flow reaches the edge of the simulation domain, ${\beta }_{h}$ is at least a factor ∼3 larger than all the other components. Therefore, as we have concluded in Sections 3 and 4, only the hot gas forms a genuine "galactic wind," with a mass loading factor of ${\beta }_{h}\sim 0.1$ that is nearly constant above $d\gt 1\,\mathrm{kpc}$. We note that with the above definitions of α and β, the Bernoulli parameter for the hot wind is given by ${ \mathcal B }={\rm{\Phi }}+({\alpha }_{h}/{\beta }_{h}){E}_{\mathrm{SN}}/{m}_{* }$.

Zoom In Zoom Out Reset image size

The vertical dependence of the mass loading factors of non-hot phases shown in Figure 9(a) implies that it is important to provide careful distinctions when reporting on the outflows measured in numerical simulations, as "wind" mass loading can be greatly overestimated if this is not done. In particular, the steep decrease of the warm-gas mass loading factor ${\beta }_{w}$ with d implies that if this were measured at $d\sim 1\,\mathrm{kpc}$, it would overestimate the value at the edge of our box ($d\sim 4\,\mathrm{kpc}$) by a factor ∼30, and the true value at larger distance would be even smaller (see also Martizzi et al. 2016, who reached similar conclusions about total β). Thus, a measurement of ${\beta }_{w}$ from a simulation with small vertical domain cannot by itself provide a prediction for warm gas mass loss in a galactic wind. However, even with a limited vertical domain, it is possible to discriminate between fountain flow and wind by combining measurements of the warm gas mass flux and its vertical velocity (as in Section 4). Similar considerations apply to observations of gas at $T\sim {10}^{4}\,{\rm{K}}$ at high latitudes in edge-on galaxies, but in that case, uncertain projection effects make this even more problematic: without an unambiguous measurement of velocity (which is subject to assumed wind geometry), it is impossible to distinguish between fountain flow and wind from observations of the emission measure.

More generally, it is essential to decompose outflows in simulations into separate thermal phases to distinguish winds from fountain flows. The phase-integrated β cannot be taken as a true wind mass loading unless the measurement is made at very large distance in a global simulation. Nevertheless, individual measurements of mass fluxes in different phases, when combined with analysis of the components of their individual Bernoulli parameters, can be used to distinguish fountains from winds even within a limited vertical domain.

The mass loading factor of the hot gas obtained here is consistent with simple estimates based on idealized experiments of superbubbles driven by multiple SNe in the warm-cold ISM (Kim et al. 2017). The shock from an individual SN sweeps up $\sim {10}^{3}\,{M}_{\odot }$ before it cools and forms a shell (Kim & Ostriker 2015a). For superbubbles created by multiple SNe, the maximum mass in hot gas per SN is ${\hat{M}}_{h}\sim 400\mbox{--}2000\,{M}_{\odot }$ prior to shell formation, but subsequently this drops to ${\hat{M}}_{h}\,\sim 10\mbox{--}100\,{M}_{\odot }$ (lower at higher ambient density). When star formation rates are self-regulated (e.g., Kim et al. 2013; Kim & Ostriker 2015b), the mean interval between SNe within projected area $\pi {H}^{2}$ in the disk of a galaxy is always $\sim 0.3\,\mathrm{Myr}$, and breakout of superbubbles is expected to occur after shell formation (see Section 5 of Kim et al. 2017). For a SN interval $\geqslant 0.1\,\mathrm{Myr}$, Figure 11 of Kim et al. (2017) shows that by the time superbubbles reach a radius of $2H$, ${\hat{M}}_{h}\sim 10\mbox{--}30\,{M}_{\odot }$, depending only weakly on ambient density. This corresponds to ${\beta }_{h}={\hat{M}}_{h}/{m}_{* }\sim 0.1\mbox{--}0.3$. The same idealized superbubble simulations show a hot-gas energy loading factor of a few percent when superbubbles expand beyond $\sim H$, since most of the energy has already been transferred to acceleration and heating of ambient gas and lost via radiative cooling at the time of shell formation.

Kim et al. (2017) argued that ${\beta }_{h}\lt 1$ is expected quite generally unless the temporal and spatial correlations of SNe are extremely enhanced compared to their mean values, requiring more than a factor of 40 elevation compared to the average conditions in self-regulated galactic disks where ${{\rm{\Sigma }}}_{\mathrm{SFR}}\geqslant {10}^{-3}\,{M}_{\odot }\,{\mathrm{kpc}}^{-2}\,{\mathrm{yr}}^{-1}$. Although our simulation does exhibit large temporal fluctuations (see Figure 3), the peak upward fluctuation compared to the mean star formation rate is only a factor of 5. Further systematic investigations for different galactic conditions will be needed to confirm whether the predictions for low ${\beta }_{h}$ apply universally in star-forming galaxies. If local or global conditions make star formation inherently extremely bursty, then ${\beta }_{h}$ may be higher.

The low mass loading factor of the hot gas in our simulations is comparable to (or slightly smaller than) the observed mass loading factor estimated in the best studied local starburst M82. Using Chandra X-ray observations, Strickland & Heckman (2009) constrained the "central" mass loading factor of the hot gas to about ${\beta }_{h}\sim 0.3\mbox{--}1$, with the "central" energy loading factor about ${\alpha }_{h}\sim 0.3\mbox{--}1$. Here, they constrained quantities using a large number of hydrodynamical models to explain the observed diffuse and hard X-rays, which come from the central 500pc region.¹¹ It is difficult to make a direct comparison, but if we consider the state of the hot medium within the energy-injecting layer ($d\lt H$) in our simulations, we have ${\beta }_{h}\sim 0.5$ and ${\alpha }_{h}\sim 0.1\mbox{--}0.2$. Although there are no systematic observational studies of mass loading factors of the hot gas, ${\beta }_{h}\lt 1$ is suggested for a wide range of star formation rates from dwarf starbursts to ultraluminous infrared galaxies, utilizing the Chevalier & Clegg (1985) wind model with observational constraints of the X-ray luminosity and star formation rates (Zhang et al. 2014).

For the solar neighborhood conditions investigated in the present simulation, the warm gas accelerated by the SNe cannot reach velocities fast enough to escape the gravitational potential. Typically star formation bursts launch warm outflows with velocity up to $\sim 100\,\mathrm{km}\ {{\rm{s}}}^{-1}$, but with most gas at lower velocity (see Figure 8), which can climb to $\sim 1\,\mathrm{kpc}$ but no further (compare the kinetic and gravitational curves in Figure 7(c)). Without additional energy and momentum input at high-$| z|$, SN feedback alone cannot drive warm winds in the Milky Way-type galaxies. Above $d=1\,\mathrm{kpc}$, Figure 9(a) shows that ${\beta }_{w}\lt 1$. Even though relatively fast warm gas could escape from dwarf galaxies with a shallower potential well than the Milky Way, Equation (7) implies that ${\beta }_{w}\sim 1/3$ for ${v}_{\mathrm{out}}\gt 50\,\mathrm{km}\ {{\rm{s}}}^{-1}$ at most. This suggests that achieving $\beta \gt 1$ in dwarfs, as appears necessary to explain current-day galaxy-halo relationships and cosmic history, would require an additional acceleration mechanism such as interaction with an outflowing cosmic ray fluid (e.g., Hanasz et al. 2013; Girichidis et al. 2016a; Ruszkowski et al. 2017). In cosmic-ray-driven winds, SN feedback may still be crucial for pushing warm gas out to large d where cosmic ray pressure gradients are sufficient to produce efficient acceleration (Mao & Ostriker 2018).

Observational studies of warm phase outflows (possibly including the ionized phase according to our definition) indicate a wide range of the mass loading factor, ${\beta }_{w+i}\,\sim 0.1\mbox{--}10$ (e.g., Heckman et al. 2015; Chisholm et al. 2017), with systematically decreasing trends for increasing galaxy mass, circular velocity, and total SFR. As the level of wind mass loading may vary with both local conditions (including ${{\rm{\Sigma }}}_{\mathrm{gas}}$ and ${{\rm{\Sigma }}}_{* }$), as well as the global gravitational potential, it will be quite interesting to measure outflow properties in simulations under widely varying galactic conditions, for comparison with current empirical scaling trends and future observations.

Recently, a number of other research groups have performed simulations with similar local Cartesian box setups to study galactic outflows driven by SNe (e.g., Creasey et al. 2013, 2015; Girichidis et al. 2016a, 2016b; Martizzi et al. 2016; Gatto et al. 2017; Li et al. 2017). Most of these simulations have adopted fixed SN rates, while varying the SN placement (e.g., random versus in high density regions). In contrast, in our simulation, SN rates and locations are self-consistent with star formation, which we believe is crucial in creating realistic multiphase outflows.

Girichidis et al. (2016b) ran a set of simulations with solar neighborhood conditions, focusing on the effect of the SN placement and of SN clustering. The authors emphasize that due to the short duration of their simulations ($100\,\mathrm{Myr}$), definitive conclusions cannot be drawn regarding wind driving. However, their measurements of outflow in the $d=1\,\mathrm{kpc}$ plane can be compared to our fountain flow measurements. For a range of different SN feedback treatments, they report mass loading factors $\sim 5\mbox{--}10$. Although this is larger than what we find (${\beta }_{w}\sim 1$ at $d\sim 1\,\mathrm{kpc}$), their simulations do show a decline in outflow over time, and our measurements are at $t\gt 250\,\mathrm{Myr}$. As pointed out by the authors, the majority of the outflowing mass is relatively dense (${n}_{{\rm{H}}}\sim 0.1$) and moves slowly (${v}_{\mathrm{out}}\sim 20\mbox{--}40\,\mathrm{km}\ {{\rm{s}}}^{-1}$), similar to the properties of our warm fountain at low d. Their fraction of high-velocity gas (${v}_{\mathrm{out}}\gt 500\,\mathrm{km}\ {{\rm{s}}}^{-1}$) is only $2\mbox{--}8\times {10}^{-5}$, which correspondingly reduces the mass loading factor of gas that is certain to escape to large distance.

In simulations including cosmic rays, Girichidis et al. (2016a) found slowly moving (${v}_{\mathrm{out}}\sim 10\mbox{--}50\,\mathrm{km}\ {{\rm{s}}}^{-1}$) warm outflows with a mass loading factor near unity. The SN-only comparison model of Girichidis et al. (2016a) had an order of magnitude lower mass loading, with just a fast, low-density component passing through the $d=1\,\mathrm{kpc}$ surface where the mass flux is measured. While this comparison suggests that cosmic rays may be crucial to accelerating warm outflows, a concern is that the role of SNe may not have been properly captured in Girichidis et al. (2016a). Due to the high computational cost of including cosmic rays, relatively low spatial resolution ${\rm{\Delta }}x\sim 16\,\mathrm{pc}$ was adopted for this pair of comparison simulations. The reported mass outflow rate for SN-only models appears to differ significantly between the high resolution (${\rm{\Delta }}x\sim 4\,\mathrm{pc}$) simulation of Girichidis et al. (2016b), with mass flux at $d=1\,\mathrm{kpc}$ of $4\times {10}^{-2}\,{M}_{\odot }\,{\mathrm{kpc}}^{-2}\,{\mathrm{yr}}^{-1}$, and the low-resolution (${\rm{\Delta }}x\sim 16\,\mathrm{pc}$) simulation of Girichidis et al. (2016a), with mass flux of ${10}^{-3}\,{M}_{\odot }\,{\mathrm{kpc}}^{-2}\,{\mathrm{yr}}^{-1}$. In our own resolution study for a similar parameter regime (see Paper I), we found that numerical convergence of SN-driven ISM properties is not guaranteed at ${\rm{\Delta }}x\sim 16\,\mathrm{pc}$, and the low-resolution outcomes are quite sensitive to the exact prescription for SN feedback. This suggests that higher resolution simulations will be needed to assess the role of interactions between the cosmic ray fluid and gas in driving galactic winds.

The Gatto et al. (2017) models are most similar to our simulations in that they self-consistently model star formation and SN feedback, rather than adopting a fixed SN rate. The main conclusion they draw is that the outflow rate strongly depends on the volume filling factor of hot gas. Above a 50% volume filling factor, their measured mass loading factors at $d=1\,\mathrm{kpc}$ are $\sim 1\mbox{--}100$ (but note that this may be exaggerated since loading is evaluated instantaneously rather than based on temporal averages; see the discussion after Equation (10) in Section 5.1). While the simulation durations of Gatto et al. (2017) are $\lt 100\,\mathrm{Myr}$, such that star formation and wind mass-loss rates are likely both subject to "start-up" transient effects, we agree with the conclusion that significant driving of fountain flows (which dominate at $d=1\,\mathrm{kpc}$ based on our work) is associated with prominent superbubbles near the disk midplane. In Paper I, we found the hot gas fills $\sim 20 \% \mbox{--}60 \%$ of the volume at $| z| \lt H$.

Martizzi et al. (2016) performed simulations for solar neighborhood conditions, as well as environments with higher gas surface density and SN rates. For high-velocity gas (${v}_{\mathrm{out}}\gt 300\,\mathrm{km}\ {{\rm{s}}}^{-1}$) that could potentially escape, they reported mass loading factors of 0.02–0.005 (lower for higher gas surface density models). They also noted the absence of wind acceleration from subsonic to supersonic velocities as a limitation of local Cartesian box simulations. We have argued in Section 3 that provided the hot-gas Bernoulli parameter and mass flux have both approached constant values at large d in a Cartesian simulation, they can be combined to make predictions for wind properties at large distance. In this case, the asymptotic mass flux at high velocity would be larger than the near-disk value as high-enthalpy gas is further accelerated when streamlines open, but the total mass flux of hot gas would change little. However, it is not possible to compare hot-gas mass fluxes, as Martizzi et al. (2016) did not separate by phase.

Li et al. (2017) conducted simulations using fixed SN rates, but decomposed by thermal phase in reporting mass and energy (as well as metal) loading factors. Their measurements are at $d=1\mbox{--}2.5\,\mathrm{kpc}$, with their "warm" component including most of what we consider "ionized," and their "hot" extending down to $T=3\times {10}^{5}\,{\rm{K}}$, slightly lower than our "hot" definition. For their solar neighborhood model, the hot and total mass loading factors are about 0.8 and 2–3, which are larger than ours by a factor of 3 to 8 (see also their discussion in comparison to work of Creasey et al. 2013, 2015; Girichidis et al. 2016b). The energy loading factor is also about an order of magnitude larger in their model.

The reason for this large difference in mass and energy loading with respect to our findings is mainly due to the difference in the vertical scale height of SNe (relative to the gas scale height), as also pointed out in Li et al. (2017). By placing SNe randomly with a fixed SN scale height of 250pc, we find that we are able to reproduce their results (see Appendix B.2). When we adopt a fixed SN scale height of 250pc, the majority of SN explosions occur outside of the main gas layer. Each SN remnant can then expand into the hot, rarefied disk atmosphere with little interaction with warm gas. Most of the injected energy is carried outward before cooling. The resulting ISM properties are also quite different: the gas scale height is smaller ($H\sim 150\,\mathrm{pc}$), star formation rates are higher, and the hot and warm/cold phases are almost completely segregated (single phase outflow). In Appendix B.2, we also test models without runaway O stars, and with random SN placement with a smaller scale height of 50pc. Both alternatives result in loading factors and warm gas velocity distribution very similar to the fiducial model. Within the standard TIGRESS framework, the spatial distribution of SNe is subject to the adopted prescription for runaways, but is otherwise self-consistently determined with respect to the gas by the distribution of star formation sites. Better theoretical and observational constraints for runaway OB stars would lead to more accurate modeling of the SN distribution, which in principle could change the wind mass loading. However, our tests suggest that a very large proportion of high-velocity runaways would be needed to significantly increase the wind mass loading, while the corresponding star formation rate and ISM phase segregation in that case might not be consistent with observations.

Finally, we remark that a fine enough grid to spatially resolve both low filling-factor warm gas and high filling-factor hot gas in the wind launching region is crucial for proper physical characterization of galactic winds. If warm and hot gas are artificially mixed (e.g., if the flow in AMR and semi-Lagrangian simulations moves from a higher resolution region near the midplane to a lower resolution region at high latitude), the result can be unphysical in ways that would compromise the implications for real galactic systems. For example, consider mixing of warm and hot flows that have total horizontally averaged mass, momentum, and energy vertical fluxes of F_ρ, ${F}_{\rho v}$, and F_E. If we further assume a steady state and neglect gravity and magnetic fields, the outgoing fluxes of the mixed gas must be the same as the sum of the horizontally averaged incoming fluxes of the warm and hot gas. That is, ${F}_{\rho }={\rho }_{\mathrm{mix}}{v}_{\mathrm{mix}}\,=\langle {\rho }_{w}{v}_{w}\rangle +\langle {\rho }_{h}{v}_{h}\rangle$, and similarly ${F}_{\rho v}={\rho }_{\mathrm{mix}}{v}_{\mathrm{mix}}^{2}+{P}_{\mathrm{mix}}$ and ${F}_{E}={\rho }_{\mathrm{mix}}{v}_{\mathrm{mix}}[{v}_{\mathrm{mix}}^{2}/2+\gamma {P}_{\mathrm{mix}}/((\gamma -1){\rho }_{\mathrm{mix}})]$ are the sum of "warm" and "hot" terms based on horizontal averages over multiple zones.

The post-mixing mean velocity will depend on the total fluxes as

$$\begin{eqnarray}&&{v}_{\mathrm{mix}}=\displaystyle \frac{\gamma }{\gamma +1}\displaystyle \frac{{F}_{\rho v}}{{F}_{\rho }}\left[1-{\left(1-\displaystyle \frac{2({\gamma }^{2}-1)}{{\gamma }^{2}}\displaystyle \frac{{F}_{E}{F}_{\rho }}{{F}_{\rho v}^{2}}\right)}^{1/2}\right].\end{eqnarray} \tag{ 11 }$$

Typically, F_ρ (near the disk) will be dominated by the warm medium contribution (see Figure 9(a)), ${F}_{\rho v}$ will be dominated by the hot medium contribution (see Figure 4(c)), and F_E will be dominated by the hot medium contribution (see Figure 9(b)).

Figure 10 shows this hypothetical post-mixing velocity as calculated from our simulation. This is based on horizontal averages of the outgoing (i.e., zones with $\mathrm{sign}({v}_{z})=\mathrm{sign}(z)$) fluxes, temporally averaged over the outflow period and over both sides of the disk. Also shown, for comparison, are the mean outflow velocities of the warm and the hot gas. The implication of the comparison shown in Figure 10 is that if numerical mixing were to occur at $d\gt 1\,\mathrm{kpc}$, the result would be ${v}_{\mathrm{mix}}\sim 100\,\mathrm{km}\ {{\rm{s}}}^{-1}$, which is intermediate between the mean hot and mean warm outflow velocities. Depending on the galaxy/halo global properties, this could have two (opposite) unphysical consequences. On the one hand, if ${v}_{\mathrm{mix}}$ is greater than the galaxy escape speed but the original v_w is less than the galaxy escape speed, then the artificially mixed flow would be able to escape with a much larger mass loading factor than would be realistic (i.e., $\sim {\beta }_{w}$ rather than ${\beta }_{h}$). On the other hand, if ${v}_{\mathrm{mix}}$ is smaller than the galaxy escape speed, then the artificially mixed flow would not be able to escape at all, whereas in reality the hot wind should escape with mass loading factor ${\beta }_{h}$, and potentially loaded with more than its share of metals.

Zoom In Zoom Out Reset image size

To avoid these unphysical consequences, it is necessary to separately resolve the multiphase gas even above the dense midplane region.

1. The hot gas component, which is created by SN shocks, behaves very similarly to expectations for a steady, adiabatic flow as it expands away from the midplane. For the horizontally and temporally averaged hot component, the mass flux and the Bernoulli parameter are close to constant as a function of the distance from the midplane d above $d\gt 1\,\mathrm{kpc}$, indicating that there is relatively little mass added or subtracted by heating or cooling from/to another phase, and little energy added or subtracted either through shocks, radiation, pressure work on other phases, or mixing.

2. The hot gas mass loading factor, defined as the ratio of outflowing mass flux to the star formation rate per unit area, is ${\beta }_{h}\sim 0.1$, decreasing by 40% from d = 1 to $4\,\mathrm{kpc}$. The measured value of ${\beta }_{h}$ is consistent with estimates from idealized numerical experiments of superbubble expansion in a warm-cold ISM.

3. Above $d=1\,\mathrm{kpc}$, the hot outflow contains a tiny fraction of the energy originally injected in the ISM by SNe. Another tiny fraction is converted to kinetic energy of the warm and cold medium accelerated by superbubble expansion, but most of the original SN energy is radiated away within the disk scale height. The ratio of outflowing hot-gas energy flux to the SN energy injection rate per unit area is ${\alpha }_{h}\sim 0.02$, decreasing by 65% from d = 1 to $4\,\mathrm{kpc}$. For $d=2\mbox{--}4\,\mathrm{kpc}$, the mean temperature of the hot gas is $\sim 1.5\times {10}^{6}\,{\rm{K}}$, and its mean velocity is $\sim 200\,\mathrm{km}\ {{\rm{s}}}^{-1}$. Magnetic energy in the outflow is negligible, with mean Alfvén speed of $\sim 30\,\mathrm{km}\ {{\rm{s}}}^{-1}$.

4. In the present simulation, the energy of the hot medium is carried mainly as heat because the outflow is constrained to remain vertical in our local Cartesian box. Allowing for the opening of streamlines at larger distance (beyond our simulation domain), much of the enthalpy would be converted to kinetic energy. The hot medium is barely affected by gravity within the simulation domain, and with an asymptotic maximum wind speed of $\sim \sqrt{2({ \mathcal B }-{\rm{\Phi }})}\sim 350\,\mathrm{km}\ {{\rm{s}}}^{-1}$ would easily escape to very large distances. More generally, local simulations with varying galactic conditions can provide measures of ${\beta }_{h}$ and the Bernoulli parameter ${ \mathcal B }$ that can be used to provide predictions for hot wind properties at large distances in a wide range of galaxies.

1. The warm (including both neutral and ionized) medium, which dominates the total gas mass in the simulation, acquires energy from SN feedback, mostly in kinetic form. However, the typical outgoing velocity of the warm medium above $d=1\,\mathrm{kpc}$ is only ${v}_{\mathrm{out}}\sim 60\,\mathrm{km}\ {{\rm{s}}}^{-1}$, which is not enough to overcome the large-scale gravitational potential in the present simulation. As a result, most of the warm gas that is blown out of the midplane eventually turns around and falls back, forming a fountain.

2. Star formation bursts lead to a succession of outflow-dominated and inflow-dominated periods of the fountain flow. During both outflow-dominated and inflow-dominated periods, gas with both signs of velocity (i.e., outflows and inflows) is present on both sides of the disk. Considering just the outflowing gas, owing to deceleration and turnaround, the outgoing warm-medium mass loading factor ${\beta }_{w}$ is a decreasing function of height d. The mean warm-gas mass flux in outflowing "fountain" gas at $d\sim 1\,\mathrm{kpc}$ is ${\beta }_{w}\sim 1$, but this drops steeply to ${\beta }_{w}\sim 0.03$ at $d=4\,\mathrm{kpc}$.

3. The value of the warm-medium energy loading factor ${\alpha }_{w}$ drops from ∼0.002 to 10⁻⁴ over d = 1 to $4\,\mathrm{kpc}$. Because ${\alpha }_{w}$ drops slightly less than ${\beta }_{w}$, the mean specific energy of the warm medium increases with d. The corresponding (volume-weighted) mean velocities of the warm gas at d = 1 and $4\,\mathrm{kpc}$ are ∼60 and $80\,\mathrm{km}\ {{\rm{s}}}^{-1}$. However, it is important to note that the increase of mean velocity in the warm medium with d is primarily due to dropout of low-velocity fluid elements, and does not reflect acceleration with height. Detailed distributions show a secular decrease in the mass and mass flux of high-velocity warm gas with increasing d.

4. A promising way to characterize warm outflows is via the velocity distribution where they are launched at $d\sim 2H$. Here, we find that the high-velocity warm gas has an exponential distribution. For a given PDF in velocity, the portion of the warm gas mass flux that is able to escape as a wind will depend on the halo potential depth. For the large-scale galactic potential in the present simulation, very little warm gas escapes, but in a dwarf galaxy, the same distribution could lead to a wind with ${\beta }_{w}\sim 1/3$.

Our simulations have two main caveats for direct comparison with observations. First, we use a local Cartesian box to achieve high resolution. The uniformly high resolution is crucial for distinguishing different phases and limiting numerical mixing. However, the local Cartesian box prevents us from following the hot outflow's evolution under global geometry with a realistic galactic potential (see Martizzi et al. 2016). Without the opening of streamlines, the hot medium cannot accelerate through a sonic point to reach its asymptotic velocity; we therefore cannot follow this process directly in our simulations. However, proper decomposition of the gas phases allows us to provide well-defined mass and energy loading factors for the hot gas, which can be robustly extrapolated to obtain predictions for asymptotic wind properties at large scales. The hot-gas loading factors at large d are slightly lower than the observational constraints deduced from M82 (Strickland & Heckman 2009). This is likely because in M82 the strong starburst has successfully cleared much of the cooler gas away from the midplane; we find that when SNe explode above the denser phases, the mass and energy loading of hot winds increases. As discussed above, to the extent that the warm phase is ballistic above a scale height, its measured velocity distributions could be used to extrapolate its properties to large distance. However, it is clear that the gravitational potential even within our local box affects the warm-medium properties at large d, so it is important to treat conclusions regarding warm-medium outflows cautiously.

Second, in the present simulations we do not include photoionization. This affects the properties of the warm and ionized phases, for which the photoionization heating can be important. Inclusion of photoionization is necessary for direct comparison with observations, where the line diagnostics are sometimes better explained by photoionization rather than a shock model (e.g., Chisholm et al. 2016a). As a first step toward this, it may be sufficient to compute ionization in post-processing, as photoionization is unlikely to be important to the dynamics of high-velocity warm gas, even though it may dominate heating. Most observations of highly ionized absorption lines, however, are based on large apertures, so direct theoretical comparisons would also require simulations with global geometry.

We are able to run our simulations over an extended period, long after initial transients that may affect quantitative results reported by others for outflows in simulations with similar physics and resolution to the TIGRESS implementation, but shorter durations. Other recent high-resolution simulations do not have self-consistent star formation and feedback, which may affect outflow loading because this is sensitive to the spatio-temporal correlation of supernovae with ISM gas of different phases. In the future, as more groups run high-resolution simulations with self-consistent star formation and SN feedback for an extended duration, it will be informative to compare phase-separated results for mass and energy loading, Bernoulli parameters, and velocity distributions, for both fountain flows and winds.

Application of the TIGRESS implementation to other galactic environments is currently underway, and promises to be quite interesting. Varying the basic model input parameters (especially gas and stellar surface density and metallicity) will enable predictions for multiphase outflow properties in a wide range of galaxies, and will provide detailed information needed to build subgrid models for winds in cosmological simulations of galaxy formation. By extending TIGRESS and other self-consistent multiphase ISM/star formation numerical implementations to include additional feedback (especially radiation and cosmic rays), understanding the complex physics behind galactic winds and fountains is within reach.

We analyze the same set of simulations with different numerical resolutions of ${\rm{\Delta }}x=4,8$, and 16 pc. Note that we showed in Paper I that ${\rm{\Delta }}x=16\,\mathrm{pc}$ is a marginal condition for convergence in star formation rates and the ISM properties. Here, we further show detailed velocity distributions of the warm medium (Figure 11) and mass and energy loading factors of the warm and hot medium (Figure 12) at varying resolution. These properties, which are key characteristics of hot winds and warm fountains, are evidently converged.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to compare the effects of SN placement, we consider, in addition to the fiducial model, three models with different prescriptions: (1) no-runaway—all SNe are in star cluster particles without runaways; (2) random-250—all SNe are randomly located horizontally and consistent with an exponential vertical profile with scale height of ${z}_{\mathrm{SN}}=250\,\mathrm{pc}$ vertically, but with the rate determined by star cluster particles; (3) random-50—the same as random-250, but ${z}_{\mathrm{SN}}=50\,\mathrm{pc}$.

By comparing the fiducial model with the no-runaway model for velocity distributions of the warm medium (Figure 13) and mass and energy loading factors of the warm and hot medium (Figure 14), we conclude that presence or absence of runaways with our adopted prescription does not significantly alter the results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With randomly located SNe, the results can change substantially, depending on the adopted scale height. For the random-250 model, the hot gas mass and energy loading factors are about 0.8 and 0.2, respectively. In this case, most SNe explode above the gaseous scale height without interacting with the warm and cold medium. The hot SN remnants simply expand outward in the low-density, hot atmosphere. A negligible warm fountain is created. In contrast, for the random-50 model, the majority of SN events happen within the warm/cold gas layer. The hot gas created in SN remnants strongly interacts with the surrounding warm and cold medium, and hot gas is lost in the process. The resulting hot gas mass and energy loading factors in the wind for the random-50 model are similar to the fiducial model. This trend with varying SN scale height was also shown in Li et al. (2017).

In short, the mass and energy loading factors of the hot wind are sensitive to the vertical distribution of SNe relative to the gas (mainly the warm medium). In our fiducial model, we assume 2/3 of SNe explode in binaries and eject runaway companions. This gives a runaway SNe fraction of 1/3. Since the ejection velocity distribution is an exponential function with characteristic velocity $50\,\mathrm{km}\ {{\rm{s}}}^{-1}$, the fraction of SNe that explode above $200\,\mathrm{pc}$ is 15%. This fraction increases to 44% in the random-250 model, while the no-runaway and random-50 models give fractions of 0.4% and 7%, respectively.

It is noteworthy that the random-250 model fails to regulate star formation rates at the same level as the fiducial model (the SFR is about 1.5 times higher for the random-250 model). In the random-50 model, an asymmetric vertical distribution of gas develops at later times (more gas is in the upper half of the simulation domain). Since the SNe distribution is set relative to the initial gas distribution rather than the instantaneous gas distribution, the asymmetric gas structure persists and allows more efficient hot and warm outflows in the lower half of the domain. This results in slightly higher mass and energy loading factors of the warm fountain for the random-50 model relative to the fiducial model.