Quantifying Supernovae-driven Multiphase Galactic Outflows

The simulations are performed using the Eulerian hydrodynamical code Enzo (Bryan et al. 2014). We set up a rectangular box with z-dimension of 5 kpc (−2.5 ≤ z ≤ 2.5 kpc). The midplane of the disk is located at z = 0. The horizontal cross section, i.e., x–y plane, is a square. The length of the horizontal dimension, l_x, varies with Σ_gas, as listed in Table 1. The idea is that we adopt higher resolutions for larger Σ_gas, while keeping the corresponding l_x smaller to gain computational efficiency (but sufficiently large to include many SNRs). The grid is refined near the midplane, with two refinement levels. Each refinement increases the resolution by a factor of two. The first level is within $| z| \leqslant 1\,\mathrm{kpc}$, and the second is $| z| \leqslant 0.5\,\mathrm{kpc}$. The finest resolution for each run is so chosen that the cooling radius of an SNR R_cool is resolved by approximately 12 cells for the initial midplane density ρ_mid. (For the definition of R_cool, see Equation (1) of Li et al. 2015.) Kim & Ostriker (2015) have shown that resolving R_cool by 10 cells is necessary to well-capture the evolution of an SNR in the Sedov–Taylor phase. Once the ISM becomes multiphase, SNe exploding in the dense region could be under-resolved. But, as we discuss in Section 4.6, this is likely a minor issue.

Table 1.  Model Description

Run	Σgas	ρmid	${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ a	PEH	Σ*	hSN,cc	vΔϕb	Tminc	Resd	lx	Rinj	tsime
	(M⊙ pc−2)	(cm−3)	(M⊙ kpc−2 yr−1)	(erg s−1)	(M⊙)	(pc)	(km s−1)	(K)	(pc)	(pc)	(pc)	(Myr)
Σ1-KS	1	0.011	1.26E-4	2.8E-28	0.5	150	52	300	12.5	1200	80.0	300
Σ1-3KS	⋯	⋯	3.78E-4	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	240
Σ1-10KS	⋯	⋯	1.26E-3	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	130
Σ10-KS	10	0.822	6.31E-3	1.4E-26	35	⋯	89	⋯	2.0	350	12.0	160
Σ10-KS-4g	⋯	1.78	⋯	⋯	180	⋯	178	⋯	⋯	⋯	⋯	64
Σ55-KS-h75	55	8.2	0.158	3.5E-25	35	75	117	⋯	0.75	150	4.5	40
Σ55-KS(-h150)	⋯	⋯	⋯	⋯	⋯	150	⋯	⋯	⋯	⋯	⋯	⋯
Σ55-KS-h300	⋯	⋯	⋯	⋯	⋯	300	⋯	⋯	⋯	⋯	⋯	⋯
Σ55-KS-h450	⋯	⋯	⋯	⋯	⋯	450	⋯	⋯	⋯	⋯	⋯	⋯
Σ55-KS-noPEH	⋯	⋯	⋯	0	⋯	150	⋯	⋯	⋯	⋯	⋯	⋯
Σ55-KS-5PEH	⋯	⋯	⋯	1.75E-24	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
Σ55-KS-1e4K	⋯	⋯	⋯	0	⋯	⋯	⋯	104	⋯	⋯	⋯	⋯
Σ150-KS	150	26.3	0.708	1.6E-24	⋯	⋯	160	300	0.60	122	3.0	15

Notes.

^aSN rate is ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$/m₀, where m₀ = 150 M_⊙. ^b ${v}_{{\rm{\Delta }}\phi }^{2}\equiv 2{\rm{\Delta }}\phi ({z}_{\max })=2{\int }_{z\,=\,0}^{{z}_{\max }}{g}_{\mathrm{tot}}{dz}$. ^cTemperature cut-off of the cooling curve. ^dResolution at z ≤ 500 pc. ^eSimulation time.

Download table as:  ASCIITypeset image

The boundary conditions are periodic for the x- and y-directions, and outflowing for z. We use the finite-volume piece-wise parabolic method (Colella & Woodward 1984) as the hydro-solver. We use the cooling curve as in Rosen & Bregman (1995), for the temperature range of 300–10⁹ K. PEH is time-independent and uniform across the box. The rate of PEH scales linearly with the star formation surface density ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$; for the solar neighborhood model Σ10-KS (see below), we adopt a PEH rate of 1.4 × 10⁻²⁶ erg s⁻¹ per H atom (Draine 2011). We explore the variations of the PEH that deviate from the fiducial settings in Section 4.2.

We have four fiducial runs: Σ1-KS, Σ10-KS, Σ55-KS, Σ150-KS. The number after Σ indicates the gas surface density in units of M_⊙ pc⁻². The SFRs associated with those runs are along the KS relation. Figure 1 shows the (Σ_gas, ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$) adopted in our simulations, indicated by blue triangles. They are plotted on top of Figure 15 of Bigiel et al. (2008), which shows the observed correlations of Σ_gas and ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ for nearby galaxies on sub-kpc scales. Table 1 summarizes the setups of the simulations. Σ10-KS is the model for the solar neighborhood. Variations of fiducial runs are described in Section 4. The gravitational field, initial gas distribution and the model of SN feedback are detailed in the next two subsections. For the fiducial run Σ10-KS, we have carried out a resolution convergence check, where we lower the spatial resolution by a factor of 2. The results agree very well, including the ISM properties, volume fraction of difference gas phases, and outflow fluxes.

Zoom In Zoom Out Reset image size

The gravitational field ("g-field" hereafter) has two components: a baryonic disk and a DM halo. The disk is modeled as self-gravitating with an iso-thermal velocity dispersion, so its g-field has the form g = 2πGΣ_* tanh (z/z_*), where z_* is the scale height of the stellar disk: ${z}_{* }\equiv {\sigma }_{* }^{2}/(\pi G{{\rm{\Sigma }}}_{* })$, in which σ_* and Σ_* are the velocity dispersion and the surface density of stars, respectively. The height z_* = 300 pc is observed for the solar neighborhood (Gilmore & Reid 1983; Binney & Tremaine 2008); we keep it unchanged for all our runs. Since we do not include self-gravity of the gas, we multiply the stellar gravitational field by a factor of 1/f_*, where f_* = Σ_*/(Σ_* + Σ_gas). The g-field from the baryonic disk is

$$\begin{eqnarray}&&{g}_{\mathrm{disk}}(z)=\displaystyle \frac{1}{{f}_{* }}2\pi G{{\rm{\Sigma }}}_{* }\,\tanh \left(\displaystyle \frac{z}{{z}_{* }}\right),\end{eqnarray} \tag{ 1 }$$

The g-field from the DM halo is modeled as an NFW profile projected to the z-direction,

$$\begin{eqnarray}&&{g}_{\mathrm{DM}}(z)=\displaystyle \frac{{{GM}}_{\mathrm{DM}}(r)\,z}{{r}^{3}},\end{eqnarray} \tag{ 2 }$$

where ${M}_{\mathrm{DM}}(r)$ is the enclosed mass of an NFW halo within radius r, so M_DM(r) = 4πρ_DM ${R}_{s}^{3}\,\{\mathrm{ln}(1+r/{R}_{s})-r/(r+{R}_{s})\}$, R_s = R_vir/c, ${\rho }_{\mathrm{DM}}=200\,{\bar{\rho }}_{\mathrm{DM}}\,c{(1+c)}^{2}$; ${\bar{\rho }}_{\mathrm{DM}}$ is the mean cosmic DM density at redshift 0, and c is the concentration parameter. For the Milky Way (MW) case, we take R_vir = 200 kpc, and c = 12 (Navarro et al. 1997). Note that r and z are related through ${r}^{2}={z}^{2}+{R}_{d}^{2}$, in which R_d is the distance from the location of the ISM patch we simulate to the galactic center. The total gravitational field is therefore

$$\begin{eqnarray}&&{g}_{\mathrm{tot}}={g}_{\mathrm{disk}}+{g}_{\mathrm{DM}}.\end{eqnarray} \tag{ 3 }$$

Table 1 lists Σ_* for all of our runs. We keep R_d = 8 kpc for all simulations, except for Σ10-KS-4g which has R_d = 3 kpc (see Section 4.3 for details). In Table 1 we also include v_Δϕ (see the table footnote for its definition) as an indicator for the total potential well for each simulation box.

Note that in the literature the adopted g-field can vary by a factor of a few even when the same "solar neighborhood" is claimed. In Figure 2 we compare our value to a few others. Walch et al. (2015) do not include the DM halo potential, so they have a smaller g; at z = 5 kpc, their g-field is about 1/3 of our value. Joung et al. (2009) uses the observed g-field in the solar neighborhood from Kuijken & Gilmore (1989), and extrapolates it into the halo. This works for z ≲ 1–2 kpc, but above that a simple extrapolation is likely too large. Δϕ for each curve shows the gravitational potential ${\rm{\Delta }}\phi ({z}_{\max })={\int }_{z=0}^{{z}_{\max }}{g}_{\mathrm{tot}}{dz}$. The numerical values of Δϕ are not negligible compared to the kinetic energies of the outflows, which are typically 100–500 km s⁻¹ (see Section 3.2). Consequently, the gravitational field is dynamically important. Indeed, the value of g-field turns out to be important for the loading efficiencies of the outflows (Section 3.3). Any meaningful comparisons between different works should take into account the difference in g-fields.

Zoom In Zoom Out Reset image size

Initially the gas has a uniform temperature T₀ = 10⁴ K. We set up the gas initial density to be in hydrostatic equilibrium in the g-field g_disk(z), i.e.,

$$\begin{eqnarray}&&\rho ={\rho }_{\mathrm{mid}}\,{\left\{\sech \left(\displaystyle \frac{z}{{z}_{* }}\right)\right\}}^{2\alpha },\end{eqnarray} \tag{ 4 }$$

where $\alpha =\gamma {\sigma }_{* }^{2}/({f}_{* }{c}_{s,0}^{2})$, γ = 5/3 is the adiabatic index of the gas and c_s,0 is the sound speed for T₀. The power-law decay at large z can result in very low density, so we set up a density floor of 3 × 10⁻²⁸ g cm⁻³. Due to the density floor and an extra gravitational field from the DM, the gas is not in perfect hydrostatic equilibrium, but in practice this has little consequence, since the outflowing gas will soon dominate the space above the midplane.

The SN surface density is related to ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ by assuming one SN explosion per m₀ = 150 M_⊙ star formation. There are some uncertainties associated with m₀; different works have adopted m₀ = 100–200 M_⊙. The distribution of SNe over time is uniform. SNe are randomly located horizontally; in the z-direction, the distribution is stratified. We distinguish two components of SNe: Type Ia and core-collapse SNe. Type Ia constitutes 10% of SN occurrence and core-collapse the rest. Type Ia SNe have an exponential distribution in the z-direction, with a scale height of 325 pc, similar to the old stellar disk (Freeman 1987). Core-collapse SNe have a Gaussian distribution vertically with a scale height h_SN,cc = 150 pc. We note that, due to runaway OB stars, core-collapse SNe may explode outside of the dense gas layer. We test the sensitivity of the results on h_SN,cc, as described in Section 4.1.

Each SN is implemented as injecting E_SN = 10⁵¹ erg energy, m_SN = 10 M_⊙ mass, and m_0,met metals (metals are modeled as "color tracers" that passively follow the mass flux, in arbitrary units), evenly distributed in a sphere. The energy added is 100% thermal. The injection radius R_inj varies for Σ_gas, and is chosen to be 0.45–0.50 of the cooling radius for the initial midplane density ρ_mid. Kim & Ostriker (2015) argued that R_inj/R_cool ≤ 1/3 is the robust criterion to capture the evolution of an SNR in the Sedov–Taylor phase. Our choice is slightly larger than that.

Figure 3 shows slices of the four fiducial runs in the x–z plane. Note that the physical scales of the slices are different from each other. The horizontal lengths are same as the simulation boxes; the vertical dimensions are shown partially. The actual dimensions each slice represents are indicated at the bottom of the temperature slices.

Zoom In Zoom Out Reset image size

Most of the dense gas stays near the midplane. The medium has multiple phases—a cold phase at a few hundred K, a warm phase at around 10⁴ K, and a hot phase at T ≳ 10⁶ K. At the boundary between the hot and warm/cold phase, gas with intermediate temperature, 10^5–6 K, is also seen. For all four runs, the hot gas volume fraction is about 60%–80% for the midplane. Hot gas occupies more volume in the halo for higher Σ_gas. Multiphase outflows are being launched from the midplane for all four runs. Cool clouds in the halo are clearly being stripped by the hot, faster gas. The hot phase appears hotter for higher Σ_gas. These qualitative results agree with previous works (McKee & Ostriker 1977; Joung et al. 2009; Creasey et al. 2013).

Figure 4 shows the phase diagram for the run Σ10-KS at t = 100 Myr. The color coding indicates the fractional mass in each (density, pressure) bin. The three phases, hot, warm, and cold, are clearly seen. Each of the three phases has some spread in the density distribution. But the majority of them are in pressure equilibrium, with P/k_B ∼ 5 × 10³ cm⁻³ K. This is in good agreement with the observations near the solar neighborhood (Cox 2005). Some mass, which lies in between the two diagonal lines that indicate the standard "warm" and "cold" phases, is out of thermal equilibrium (Heiles & Troland 2003).

Zoom In Zoom Out Reset image size

Figure 5 shows the density, weighted by volume, of different phases. Hotter phases have progressively smaller densities. The gas density for each phase near the plane is higher than that in the outflows. The warm-hot phase with 10⁴ < T < 3 × 10⁵ K has a slightly lower density than the warm phase. As shown in Figure 3, the warm-hot phase is mostly at the interface between warm clouds and hot gas. Even the coldest phase is seen at large z, even though the volume fraction can be very small. The densities for each phase agree with the observations of the local ISM (see, e.g., Draine 2011)

Zoom In Zoom Out Reset image size

The velocity structure of the gas determines how far the gas can travel in a gravitational potential. Moreover, the velocity distribution can be observed from the profiles of emission/absorption lines. Since different gas phases have drastically different velocities, and observationally they are detected through different line tracers, we hereby show the velocity distribution for each gas phase separately.

Figure 6 shows the z-velocity distribution for the run Σ10-KS at t = 160 Myr. The y-axis indicates the fractional mass in each velocity bin. Each curve is normalized to unity. Hotter phases have larger velocities, agreeing with the general observational trend (e.g., Heckman et al. 2001; Rupke et al. 2002), and other simulation works (Creasey et al. 2013; Girichidis et al. 2016b). The hottest phase has the broadest range of velocities, up to ≳600 km s⁻¹. A fraction of the warm phase can reach >100 km s⁻¹. The velocities of the cold phase remain small at ≲50 km s⁻¹.

Zoom In Zoom Out Reset image size

Our simulations only capture the gas evolution that is relatively close to the midplane, i.e., $| z| \leqslant 2.5\,\mathrm{kpc}$. One way to relate the "local" outflows to their large-scale evolution is to estimate how far the gas can travel in a given potential. First, let us consider a ballistic evolution. For a parcel of gas with a velocity v₀ at the bottom of a potential well, the furthest distance it can reach, R, is simply determined by 1/2 ${v}_{0}^{2}={\rm{\Delta }}\phi (R)$. We use the function R(v₀) to describe such a relation. Figure 7 shows R(v₀) for the MW, for a single stream line that is perpendicular to the disk and goes through the center of the disk. The potential of the DM halo is the same as described in Section 2, and the disk is modeled as a 2D razor-thin disk with a mass M_D = 5 × 10⁵ M_⊙ and a radius of R_D = 9.5 kpc. The mass distribution within the disk is uniform. Thus, along the stream line mentioned above, the g-field from the disk has a simple analytic form: ${g}_{{\rm{D}}}=2{{GM}}_{D}(1-z/\sqrt{{z}^{2}+{R}_{{\rm{D}}}^{2}})/{R}_{{\rm{D}}}^{2}$. From Figure 7, gas with v₀ ≳ 620 km s⁻¹ can escape from the DM halo; gas with v₀ = 300 km s⁻¹ can travel to R ∼ 50 kpc, and so on.

Zoom In Zoom Out Reset image size

We now discuss what should be used as v₀. The naive answer, i.e., the bulk velocity projected to the direction of g, may only give a lower bound. For a compressible fluid, it is likely that the gas motion is not ballistic, but thermal energy can later convert to bulk motions. According to the Bernoulli principle, the Bernoulli constant ${ \mathcal B }\equiv {v}_{z}^{2}/2+\gamma /(\gamma -1)P/\rho +\phi$, remains unchanged along a stream line in a steady-state flow (for constant γ). We thus define a modified "Bernoulli velocity" ${v}_{\widetilde{{ \mathcal B }}}\equiv \sqrt{2}{\widetilde{{ \mathcal B }}}^{1/2}$, where $\widetilde{{ \mathcal B }}\equiv { \mathcal B }-\phi$. So for a parcel of gas with a bulk velocity v_z and a Bernoulli velocity ${v}_{\widetilde{{ \mathcal B }}}$, the approximate range of radii it can reach is roughly $R({v}_{0}={v}_{z})\sim R({v}_{0}={v}_{\widetilde{{ \mathcal B }}})$. Note we only aim at a very rough estimate, ignoring cooling, interaction among different gas phases, etc., and assuming γ = 5/3.

Figure 8 shows the mass-wighted v_z and ${v}_{\widetilde{{ \mathcal B }}}$ for different phases for the fiducial run Σ10-KS. The y-axis on the right shows R corresponding to the velocities on the left axis. Only gas at $| z| \geqslant 1\,\mathrm{kpc}$ is included. The data are averaged over the last 20% of the simulation time. The error bars indicate time variations. The hot gas is affected more by each SN explosion, thus its properties vary stronger with time. Both v_z and ${v}_{\widetilde{{ \mathcal B }}}$ increase with gas temperature. The hottest phase, given its large v_z (∼150 km s⁻¹) and ${v}_{\widetilde{{ \mathcal B }}}$ (∼370 km s⁻¹), would travel much further into the halo, to about 30–70 kpc. Since the majority of the hot gas would not escape from the DM halo, large-scale fountain flows would form. The small velocities of the cool phase imply that they would fall back at below 10 kpc, unless being accelerated significantly. The velocity of the warm-hot phase, with T = 10⁴–3 × 10⁵K, is much closer to the warm phase than the hot. The ratio ${v}_{\widetilde{{ \mathcal B }}}$/v_z is largest for the hot phase, meaning that a significant fraction of the energy is thermal, which may convert to the bulk motion at large radii. For cooler phases, in contrast, most energy is kinetic. Note that the fiducial run is for the solar neighborhood, and is not representative of the MW disk in general. We discuss the model for the MW-average, Σ10-KS-4g, in Section 4.3. Hot flows are much faster there than the solar neighborhood, and can thus have a much broader impact on the CGM (see Section 4.3 for details).

Zoom In Zoom Out Reset image size

Figure 9 shows the mass-weighted v_z and ${v}_{\widetilde{{ \mathcal B }}}$ for the hot outflows, as a function of ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ for the four fiducial runs. Again, the error bars show time variation. Both v_z and ${v}_{\widetilde{{ \mathcal B }}}$ increase with SFR. The velocities for the Σ1-KS run are 60–200 km s⁻¹, and rise to 600–900 km s⁻¹ for Σ150-KS. The large velocities imply that the hot outflows can travel far, and even escape from the halo potential. This suggests that hot outflows play a critical role in regulating the CGM and even the IGM.

Zoom In Zoom Out Reset image size

We find that ${v}_{z}\propto {\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}^{0.24}$, and ${v}_{\widetilde{{ \mathcal B }}}\propto {\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}^{0.16}$. Observationally, while there is little direct constraint on velocities for hot gas, for the warm/cool phases, Martin (2005) and Weiner et al. (2009) found $v\propto {\mathrm{SFR}}^{0.3\mbox{--}0.35}$ for galactic-scale outflows. Our findings seem to indicate that the velocity dependence on SFR for the hot gas is weaker than for the cooler phases.

In this section, we discuss the loading capability of the outflows. We define outflows to be at $| z| \geqslant 1\,\mathrm{kpc}$ and with outgoing z-velocity. We find in our simulations that  the outflow fluxes show little variation with z at $| z| \geqslant 1\,\mathrm{kpc}$.

The mass loading factor η_m is defined as the ratio between the outflowing mass flux and ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$, that is,

$$\begin{eqnarray}&&{\eta }_{{\rm{m}}}\equiv \displaystyle \frac{\langle \rho {v}_{z}\rangle }{{\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}}.\end{eqnarray} \tag{ 5 }$$

The outflow flux includes both sides of the plane, and "$\langle ...\rangle$" denotes averaging over space ($1\leqslant | z| \leqslant 2.5\,\mathrm{kpc}$) and time (last 40% of t_sim).

The energy loading factor η_E is the ratio between the z-component of the energy flux and the energy production rate by SNe, that is,

$$\begin{eqnarray}&&{\eta }_{{\rm{E}}}\equiv \displaystyle \frac{\langle ({e}_{{\rm{k}}}+{e}_{\mathrm{th}}){v}_{z}\rangle }{{\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}{E}_{\mathrm{SN}}/{m}_{0}},\end{eqnarray} \tag{ 6 }$$

where e_k and e_th are the kinetic and thermal energy per unit volume.

The metal loading factor η_met is the ratio of the z-component of the metal flux to the metal production rate by SNe:

$$\begin{eqnarray}&&{\eta }_{\mathrm{met}}\equiv \displaystyle \frac{\langle {\rho }_{\mathrm{met}}{v}_{z}\rangle }{{\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}\ {m}_{0,\mathrm{met}}/{m}_{0}},\end{eqnarray} \tag{ 7 }$$

where ρ_met is the density of metals, and ${m}_{0,\mathrm{met}}$ is the mass of metals each SN produces (in arbitrary units, see Section 2.3). Note that we assume the metals are solely produced by SNe, and the ISM is otherwise pristine.

Figure 10 summarizes the loading factors and the volume fraction of each gas phase in the outflows as functions of Σ_gas, for the four fiducial runs. For the loading factors we show the total loading, which includes all gas phases, as well as that of the hot outflows only. The error bars indicate the standard deviation of the time variation.

Zoom In Zoom Out Reset image size

We find that η_m decreases monotonically with increasing Σ_gas. The largest mass loading is about 6 for Σ1-KS. For the solar neighborhood case, i.e., Σ10-KS, our η_m is around 2–3. For our highest density case Σ150-KS, η_m is only 0.2. The fraction of the mass loading contributed from the hot gas is about 1/3, except for Σ150-KS, where most of the mass flux is hot. The warm phase dominates the outflowing mass flux except when ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ is very high.

We fit our mass loading factor by a simple power-law function of Σ_gas:

$$\begin{eqnarray}&&{\eta }_{{\rm{m}}}=7.4{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{1{{M}}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{{\alpha }_{\mathrm{ml}}},{\alpha }_{\mathrm{ml}}=-0.61\pm 0.03.\end{eqnarray} \tag{ 8 }$$

For the hot gas, the mass loading factor is

$$\begin{eqnarray}&&{\eta }_{{\rm{m}},{\rm{h}}}=2.1{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{1{{M}}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{{\alpha }_{\mathrm{ml},{\rm{h}}}},{\alpha }_{\mathrm{ml},{\rm{h}}}=-0.61\pm 0.03.\end{eqnarray} \tag{ 9 }$$

Creasey et al. (2013) have found a sharper decline, with α ≈−1.1. Our results agree with theirs for Σ_gas ≲ 10 M_⊙ pc⁻², but there are relatively large discrepancies at higher densities; see Section 4.5 for a discussion. The X-ray emission from the halo of edge-on galaxies suggests a decreasing mass loading of hot gas for higher SFR (Zhang et al. 2014; Bustard et al. 2016; Wang et al. 2016), consistent with our results.

The energy loading factor η_E shows surprisingly little dependence on Σ_gas. Despite a factor of 150 span of Σ_gas, η_E stays at about 10%–30%. This means a significant fraction of SN energy goes into the outflows. There is no obvious trend of η_E as a function of Σ_gas. The hot gas contains the majority, >90%, of the outflow energy.

The metal loading factor η_met shows somewhat larger variation than η_E, although again we do not find an apparent dependence on Σ_gas. Overall, a quite large fraction of metals goes into the outflows, about 40%–90%. Hot outflows carry 35%–60% of the metals produced by SNe. While the warm/cool phase may fall back to the disk later, the hot gas has the potential to travel much further, even to escape the halo (see Section 3.2), and metals will be carried along. The mass–metallicity relation of galaxies implies that a significant fraction of all metals ever produced are no longer in galaxies (Tremonti et al. 2004; Erb et al. 2006). Our numbers agree with this general picture.

The volume in outflows is progressively occupied by the hot gas as Σ_gas becomes larger. The cold phase with T < 10³ K has a negligible volume fraction, thus we omit it in the plot. For Σ1-KS, the volume is equally shared by the warm and hot phases; for Σ55-KS, more than 90% of volume is hot; for Σ150-KS, the outflows are completely dominated by the hot phase.

Where SNe explode is critical for feedback efficiency. An SN exploding in a dense medium quickly radiates away its energy, and has little impact on the large-scale ISM, let alone contributing to driving winds (Girichidis et al. 2016b). On the other hand, if an SN explodes in an environment dominated by tenuous gas, then the cooling is much less efficient, and a significant fraction of energy can be preserved (e.g., Hennebelle & Iffrig 2014; Simpson et al. 2014; Gatto et al. 2015; Li et al. 2015; Walch et al. 2015). One key factor to determine where SNe explode is the fact that a significant fraction of OB stars are "runaways," that is, they have high velocities. A simple calculation shows that OB runaways can migrate a few tens to a few hundreds of parsecs before exploding as SNe (Li et al. 2015). This greatly facilitates SN feedback by allowing some of them to release their energy outside the dense SF regions.

In principle, the locations of core-collapse SNe depend on the velocities of OB stars, their lifetimes, the external gravitational field, close encounters with other stars, etc. One can also infer the SNe explosion sites from the spatial distribution and the velocities of pulsars (Narayan & Ostriker 1990). In this paper, we do not aim to model the location of SNe from first principles, but simply explore how sensitively the outflow properties depend on the vertical distribution of SNe. For the fiducial runs we have the ${h}_{\mathrm{SN},\mathrm{cc}}=150\,\mathrm{pc}$. Now we experiment with h_SN,cc. We take the run of Σ55-KS and change h_SN,cc to 75, 300, and 450 pc. In Table 1, they are identified by the names Σ55-KS-h75, Σ55-KS-h300, and so on.

Figure 11 shows the mass, energy, and metal loading factors for different h_SN,cc. Interestingly, η_m depends on h_SN,cc in a different way from η_E and η_met. As h_SN,cc increases, η_m first increases, reaches a peak at h_SN,cc = 150 pc, and then declines. The loading factors η_E and η_met increase monotonically with h_SN,cc, and reach a plateau at h_SN,cc ≳ 150–300 pc. For h_SN,cc = 75 pc, most SNe are buried in the midplane, and radiate their energy there, so the feedback is least efficient. As h_SN,cc becomes larger, more SNe explode in the low-density disk–halo interface, resulting in a more effective energy and metal loading. The non-monotonic dependence on h_SN,cc of η_m can be understood in this way: when h_SN,cc is too small, most energy radiates away in the disk, so the energy insufficiency is the limiting factor for the mass loading; when h_SN,cc is too large, the outflowing mass is simply SNe ejecta, with little ISM involved. Consequently, the maximum η_m occurs in between those two extremes. We find that η_m achieves unity when h_SN,cc = 150 pc, while η_m ≲ 0.2 for other cases. For h_SN,cc ≳ 300 pc, 40%–50% of energy and 80% of metals produced by SNe end up in the outflows.

Zoom In Zoom Out Reset image size

In the absence of SN explosions, PEH maintains a two-phase warm/cold ISM (Draine 1978; Wolfire et al. 1995). The value of the PEH rate Γ_PEH determines the relative amount of mass in the two phases and the pressure of the ISM (Wolfire et al. 2003). With SNe, Γ_PEH is an important factor in determining whether the ISM is in a thermal runaway state or not (Li et al. 2015). A higher Γ_PEH keeps more gas in the warm phase and increases the ISM pressure, thus limiting the size of the hot bubbles of SNRs. As a result, SNRs may not effectively overlap, and each SNR loses the majority of its energy at the cooling stage. Therefore little energy is left to drive an outflow. We thus expect that Γ_PEH is important in determining the outflow properties.

In this section we study the effect of different values of Γ_PEH. We note that Γ_PEH depends on many factors: the far-UV background, dust abundance, work function of the dust grains, ionization fraction of the gas, etc. (Draine 2011). A star-forming region has a very complex structure with strong and time-varying radiation background with both ionizing and non-ionizing photons. The radiation background also varies in space, and is much more intense around OB stars. The exact condition is thus hard to determine. For simplicity, we keep Γ_PEH constant in time and uniform in space for each simulation, but just change Γ_PEH to explore its effect. Note that some previous works adopt a cooling curve with a cut-off at 10⁴ K, which prohibits the formation of the cold phase. This is similar to the effect of a very high Γ_PEH. We include a discussion of the cooling curve cut-off as well.

We compare four runs that have Σ_gas = 55 M_⊙ pc⁻². The set-ups are identical (including the SN rate) except Γ_PEH:

• (a)  
  ${{\rm{\Gamma }}}_{\mathrm{PEH}}=0$ ("Σ55-KS-noPEH");
• (b)  
  Γ_PEH = 3.5 × 10⁻²⁵ erg s⁻¹ (fiducial, "Σ55-KS");
• (c)  
  Γ_PEH = 1.75 × 10⁻²⁴ erg s⁻¹ ("Σ55-KS-5PEH");
• (d)  
  cooling curve has a cut-off at ${T}_{\min }={10}^{4}\,{\rm{K}}$ ("Σ55-KS-1e4K").

In Figure 12 we show the slices for the fiducial run and Σ55-KS-1e4K. Adopting T_min = 10⁴ K results in a much larger scale height of gas (defined as enclosing 80% of the mass in the box), h_gas ∼ 200 pc, in contrast to h_gas ∼ 10 pc for the fiducial case. Assuming hydrostatic equilibrium, i.e., gravity is balanced by the thermal and turbulence pressure,

$$\begin{eqnarray}&&{h}_{\mathrm{gas}}\sim 140\mathrm{pc}\,\left(\displaystyle \frac{1+{{ \mathcal M }}^{2}}{2}\right)\left(\displaystyle \frac{{T}_{\mathrm{gas}}}{{10}^{4}\,{\rm{K}}}\right)\left(\displaystyle \frac{5\times {10}^{-9}\,\mathrm{cm}\,{{\rm{s}}}^{-2}}{g}\right),\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal M }$ is the local Mach number of the gas, which is on the order of unity. Note that in our simulations, SNe have a Gaussian distribution with h_SN,cc = 150 pc. So for the fiducial case, once the multiphase medium is formed, most SNe explode outside of the gas layer, whereas for Σ55-KS-1e4K, most SNe explode within the gas layer. For the latter, since the ISM is not in a thermal runaway state, most of the energy released from SNe is radiated away. Therefore, a cooling curve with T_min = 10⁴ K gives much smaller energy, mass, and metal loading. It is true that we are adopting a temperature cut of 300 K, and the actual temperature of the cold phase can be even lower. But our temperature cut is low enough to allow the ISM at midplane to undergo a thermal runaway—as mentioned in Section 3.1, all fiducial runs have a volume fraction of hot gas of 60%–80% at midplane. We thus believe that our results do not suffer from a qualitatively erroneous cooling loss, while a temperature cut at 10⁴ K may do so.

Zoom In Zoom Out Reset image size

Figure 13 compares the loading factors of all four runs in this section. The simulation Σ55-KS-1e4K gives an energy loading two orders of magnitude smaller than the fiducial run; η_m is smaller by a factor of 10, and η_met by a factor of 30. This indicates that if the formation of the cold phase is prohibited, the power of SN feedback is severely underestimated. Comparing the three runs with different Γ_PEH, we find that when Γ_PEH is higher, the energy loading is smaller, as expected, since h_gas is increasingly larger for stronger PEH. The mass and metal loadings do not show significant variation. This is likely due to the following two effects counteracting each other: a smaller η_E means that less energy is available to drive the mass out, while a larger h_gas is favorable to loading more gas, as discussed in Section 4.1.

Zoom In Zoom Out Reset image size

To explore the effect of external gravitational field on the outflows, we take the MW as an example. The fiducial run Σ10-KS uses the gravitational field in the solar neighborhood, which has Σ_gas = 10 M_⊙ pc⁻², Σ_* = 35 M_⊙ pc⁻², and a displacement R_d = 8 kpc from the center of the DM halo. This g-field is smaller than the inner part of the MW disk. We set up a higher gravity run Σ10-KS-4g, which uses a g-field more typical for the inner MW disk, with Σ_gas = 10 M_⊙ pc⁻², Σ_* = 180 M_⊙ pc⁻² and R_d = 3 kpc. The g-field is approximately four times that of the solar neighborhood. Figure 14 shows the different components of the two g-fields.

Zoom In Zoom Out Reset image size

Naively, one would think a larger gravity would make the feedback less effective, as gravity drags the outflows toward the disk. While this is generally true for regions away from the launching region, the situation near the disk is more complex. We plot the ratios of the loading factors between Σ10-KS-4g and Σ10-KS in Figure 15. The comparison is done for the time interval t = 40–64 Myr, and the error bar shows the standard deviation of time variation. While η_m is indeed smaller by a factor of 3–4 in the higher gravity case, the energy loading η_E is, nevertheless, a factor of 3–5 larger. The metal loading η_met is also larger by a factor of 1.5.

Zoom In Zoom Out Reset image size

How can we understand this? It turns out that the dominant impact of larger g-field here is to reduce h_gas. As we also show in Figure 15, a factor of 4 increase in gravity results in roughly the same factor of decrease in h_gas. This is expected for gas in hydrostatic equilibrium (Equation (10)). Since we keep h_SN,cc the same for the two runs, a smaller h_gas exposes more SNe in the low-density halo. As discussed in Section 4.1, more SNe above the gas layer can lead to a smaller η_m but larger η_E and η_met. Since less mass is heated by more energy, ${v}_{\widetilde{{ \mathcal B }}}$ of the outflowing gas is much larger in Σ10-KS-4g, by a factor of 3 than in the solar neighborhood. The values of v_z and ${v}_{\widetilde{{ \mathcal B }}}$ are about 175 km s⁻¹ and 980 km s⁻¹, respectively; the latter is even larger than the escape velocity of the MW halo ∼620 km s⁻¹. Thus the outflows for the MW-average are much more vigorous, which can broadly impact the CGM and even the IGM (see further discussion in Section 5.1).

For our fiducial runs, we assume that ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$ scales with Σ_gas as in the KS relation. Although the KS relation is well-established at scales ≳kpc, variation appears at smaller scales (Heiderman et al. 2010). In particular, star formation tends to occur in groups, and the OB stars are clustered in space and time. The sizes of our simulation boxes are in the sub-kpc regime, so it will be interesting and relevant to discuss the variation on the SN rate. In this section, we discuss the effect of enhanced SN rates on the outflows. We are interested in how the energy, mass, and metal loading efficiencies depend on the SN rates. Since the interaction of blast waves is highly nonlinear, it is non-trivial to predict whether the impact of multiple SNRs would be a simple add-up, to reinforce, or to cancel each other out.

We take the run Σ1-KS, and increase the SN rate by 3× and 10×. These runs are listed in Table 1 as "Σ1–3KS" and "Σ1–10KS." Figure 16 shows the mass, energy, and metal fluxes of the outflows relative to the fiducial run. The mass flux scales with the SN rate in a sub-linear manner. A factor of 3 and 10 increase in the SN rate only results in, on average, a factor of 1.5 and 3 enhancement in the mass flux, respectively. The energy and metal fluxes, on the other hand, show a roughly linear correlation with the SN rate. This means that the mass loading is less efficient when we increase the SN rate, while the energy and metal loading factors remain roughly constant for different SN rates.

Zoom In Zoom Out Reset image size

We caution that, even for the fiducial run, which has the lowest SN rate, most of the midplane is in a hot-dominated multiphase state. The sub-linear dependence of the mass flux and roughly linear dependence of the energy and metal flux are likely to be the features in this regime. If, for example, one starts with an SN rate sufficiently small so that the SNRs in the disk would not overlap, then the enhancement of the SN rate would lead to a transition from a steady-state ISM to forming outflows. As a result, the dependence of all fluxes on the SN rate would be super-linear. But since we are, in this paper, interested in the regime where outflows are generated, we do not explicitly explore the parameter space that leads to a steady-state ISM. Indeed, for all four fiducial runs, which are along the KS relation, the ISM at the midplane is in a thermal runaway state and outflows are being launched. Therefore, the scaling relations as shown in Figure 16 should hold for other Σ_gas cases as well.

Girichidis et al. (2016b) find that clustering of some SNe does not affect the mass outflow rate. We find a very mild increase in mass flux, although with large fluctuations. Within error bars our results are consistent with each other. Overall, once the ISM is hot-dominated, clustering of SNe does not help with the loading factors, and may even be negative for loading mass.

Girichidis et al. (2016b) have found that for the solar neighborhood, SNe can blow away most of the gas in the midplane, and drive outflows with a mass loading up to 10. This is higher than our value, which is around 2–3. Compared to our model, their SN scale height is smaller, 50 pc, which causes an "explosive" thermal-runaway at the midplane. Initially there is no leak of those hot gases, whose high-pressure propels the neutral gas layer up, like the formation of a super-bubble. Later on, the Rayleigh–Taylor instability will develop, the shell will fragment, and hot gas leaks to form winds (Mac Low & Ferrara 1998). For our case, h_SN,cc is larger, meaning that SNe are more spread-out in the z-direction. We do see warm shells of gas being driven out initially, but that does not involve too much mass, and would later either go beyond the box, or fragment and fall back. Our loading factors are calculated after those initial transient stages. Additionally, a weaker g-field (see Section 2.2) may also partially account for their relatively large mass loading.

Creasey et al. (2013) study SNe-driven outflows covering a broad parameter space of Σ_gas and g-field. The g-field can be expressed using the gas fraction, f_g ≡ Σ_gas/(Σ_gas + Σ_*). We here conduct a one-on-one comparison between our fiducial runs and their models. We note that their boxes are smaller in the vertical direction, $| z| \leqslant 0.5\,\mathrm{kpc}$, and their outflow fluxes are measured at the outer boundaries; whereas ours are averaged over $1\leqslant | z| \leqslant 2.5\,\mathrm{kpc}$. Our fiducial runs that overlap with their models are Σ10-KS, Σ55-KS, and Σ150-KS, which corresponds to (Σ_gas, f_g) = (10, 0.15), (55, 0.5), and (150, 0.7), respectively. We convert our g-field from the DM halo to an equivalent surface density of ≈25 M_⊙ pc⁻². We interpolate their data if there is no direct comparison.

For the solar neighborhood, their mass loading factor is about unity, and energy loading ("thermalization factor" in their terminology) is around 0.1. These are slightly smaller, by a factor of 1.5–2, compared to our simulation. For the higher Σ_gas cases, however, the discrepancies are larger. For Σ55-KS and Σ150-KS, our η_m are 0.9 ± 0.3 and 0.2 ± 0.05, whereas theirs are 0.2 and 0.03, smaller by a factor of 4–7 than our values; for η_E, our values are around 0.2 for both cases, whereas theirs are around 0.05, smaller by a factor of 4. In a follow-up paper, Creasey et al. (2015) measure the metal loading efficiency of the outflows for some runs in Creasey et al. (2013). Our model parameters only overlap with theirs for the solar neighborhood, in which we have η_met = 0.65 ± 0.2 and they have a smaller ≈0.2. Note also that we both assume a KS relation to relate Σ_gas and ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$, but we convert the SFR to SN rate by assuming m₀ = 150 M_⊙ (definition of m₀ in Section 2.3), whereas they have m₀ = 100 M_⊙. This means the difference is even larger by a factor of 1.5. We attribute the discrepancies mainly to their adoption of a cooling cut-off at 10⁴ K. As a result, the neutral gas layer in their runs is thicker, and most SNe lose their energy there, thus the loading factors are smaller (see Section 4.2 for detailed discussions).

Under the impact of many SN blast waves, the ISM becomes multiphase. The cooler, denser phase settles down near the midplane, whereas the hotter phases escape and form outflows. There are two regimes of the media: (i) a warm/cool-dominated ISM where, if an SN explodes within, it would lose most energy by radiative cooling, and (ii) a hot-dominated ISM, where SN shock waves would propagate much faster and further, while the cooling is inefficient (Cowie et al. 1981; Li et al. 2015). The fraction of SNe that explode in a hot-dominated ISM, f_SN,h, is key in determining the efficiency of the loading efficiency of energy and metals. A stronger external g-field or a weaker PEH leads to a smaller h_gas, leaving more SNe exploding in the low-density medium, thus a more powerful loading of energy and metals; a larger h_SN,cc has a similar effect. Simply enhancing the SN rate, without changing h_SN,cc or h_gas, yields unchanged energy and metal loading. The mass loading factor, on the other hand, has a more complex dependence on f_SN,h: a f_SN,h that is either too small or too large would result in a small mass loading factor.

We briefly discuss the possible impact of the physics that we do not include in this work.

Under-resolved SNe. The resolution and R_inj for each run are fixed, and are chosen based on R_cool for the initial ρ_mid. Later, when the ISM becomes multiphase, SNe exploding in the tenuous/hot phase continue to be well resolved. For all fiducial runs, the hot gas volume fraction at the midplane is 60%–80%. Since SNe are randomly located, a similar fraction of SNe would explode in the hot phase. The remaining SNe which explode in the denser phase are likely under-resolved, but we argue that these SNe are unimportant in driving large-scale outflows. Taking Σ10-KS for example, the cold phase has a density of about 10 cm⁻³, corresponding to a R_cool ∼ 7 pc. So even when the evolution of the SNR is resolved spatially and temporally, it will lose the majority of its energy at R_cool. This means that these SNe only have a very localized impact, and contribute little to large-scale outflows.

Magnetic fields. In the solar neighborhood, the magnetic pressure is overall similar to the thermal pressure of gas, and even larger for dense phases (Heiles & Crutcher 2005). This extra pressure, if included, would be likely to enhance h_gas, and also make SN bubbles smaller (Slavin & Cox 1993). The net effect on the loading factors would be similar to having a stronger PEH—the energy and metal loadings would be smaller, while the change of the mass loading is not certain. The magnetic field on the dynamics of the diffuse ISM is rather mild for the solar neighborhood (de Avillez & Breitschwerdt 2005; Hill et al. 2012; Walch et al. 2015), except possibly for providing support for the vertical distribution of the gas (Cox 2005). Magnetic fields alone are unlikely to play an active role in driving the outflows.

Self-gravity. We did not incorporate self-gravity in our simulations. For the solar neighborhood, only about 1% of the volume is in a self-gravitating state (Draine 2011). For larger gas surface density cases, however, self-gravity is probably more important. Including self-gravity would make the cold phase smaller in volume, thus may facilitate feedback for those SNe that explode in the inter-cloud space (e.g., Girichidis et al. 2016b; Kim & Ostriker 2016). But we note that (i) the external gravitational field in our simulations does include that from the gas, so the effect of self-gravity is not fully missed, at least in the vertical direction. (ii) As mentioned in Section 3.1, essentially in all runs, the volume fraction of hot gas at the midplane is 60%–80%, thus we believe shrinking the volume fraction of cold gas by some percentage would not lead to a qualitative change in the loading factors of outflows. (iii) Self-gravity should be included with caution; that is, one should also resolve the counterbalancing force—feedback acting below the Jeans scale of the dense clumps, which is challenging given the resolution of the current ISM simulations; otherwise, most gas would collapse into a few small clumps, which is also not realistic.

Cosmic rays (CRs). SNe are considered to be the main acceleration sites for CRs. Around 5%–15% of SNe energy may go into CRs (Hillas 2005). Recent simulations show CRs are promising candidates to drive galactic scale outflows, with a mass loading around 0.5 (Uhlig et al. 2012; Booth et al. 2013; Hanasz et al. 2013; Salem & Bryan 2014). Recent high-resolution, local simulations indicate that for the solar neighborhood, both CRs and SN thermal feedback can drive an outflow with mass loading around unity (Girichidis et al. 2016a; Simpson et al. 2016). CR-driven outflows are cooler, slower, and smoother than thermally driven ones. Peters et al. (2015) point out that compared to thermal feedback, CR-driven outflows give too little hot gas to match the soft X-ray background. Many questions remain open regarding how CRs propagate in, and interact with, a multiphase ISM. Further work is needed to determine whether thermal feedback or CRs are the dominant force launching outflows.

Radiation pressure. Radiation pressure can propel gas out if it is optically thick, especially for infrared light. The photon energy from the stars can then be effectively utilized to drive outflows (Murray et al. 2005; Hopkins et al. 2012). This can be the case for the extremely dense and dusty SF regions, such as in the ultra-luminous infrared galaxies (ULIRGs). We do not actually explore those extreme cases (Figure 1). For the Σ_gas adopted in our models, the optical depth for the infrared is much less than unity, therefore the radiation pressure is not important.

We also do not include preprocessing of the ISM by other stellar feedbacks, such as stellar winds, ionizing photons, etc. In general these feedbacks alone do not contribute directly to launching outflows (unless the radiation is highly trapped), but they may make the ISM inhomogeneous. Therefore it is possible for SNe to explode in a less dense environment, thus facilitating SN feedback to some extent (e.g., Gatto et al. 2016; Peters et al. 2017).

Hot, X-ray-emitting coronae have been observed around the MW (Snowden et al. 1998) and other star-forming disk galaxies, up to a few tens of kpc away from the galaxy (Strickland et al. 2004; Tüllmann et al. 2006; Anderson & Bregman 2011; Dai et al. 2012). Hot gas around the MW is also detected through O vii and O viii absorption lines (Fang et al. 2006; Bregman & Lloyd-Davies 2007; Gupta et al. 2012).The COS–Halos survey finds the warm-hot metal-enriched gas out to ≳100 kpc for MW-like galaxies, as traced by the O vi absorption lines (Tumlinson et al. 2011). Our simulations indicate that the hot phase in the outflows has the largest v_z and ${v}_{\widetilde{{ \mathcal B }}}$, and can potentially travel distances comparable to the size of DM halos. Hot outflows carry 10%–50% of energy and 30%–50% of metals produced by SNe (Section 3.3). Thus, the hot, metal-enriched CGM may at least be partly due to these SNe-driven outflows. Those outflows should have an important dynamic impact on the CGM (Hopkins et al. 2012).

Recently, Faerman et al. (2016) create a two-phase phenomenological model of the halo gas for MW-like galaxies, which simultaneously fits the absorption features of O vi, O vii, and O viii. The two phases in their model have volume-weighted medium temperatures of 3 × 10⁵ K and 1.8 × 10⁶ K. Nearly 90% of the mass is in the hotter phase. The total masses of the two phases are comparable to that of the baryons in the DM halo but not in the galaxy. Our hot phase in Σ10-KS and Σ10-KS-4g has a temperature comparable to that in their model. We note, though, that the cooling luminosity of the CGM in their model is about twice the SNe heating rate from the MW, and given that η_E ≈ 0.45 for the MW-average, the total discrepancy is about a factor of 4.

Warm/cool outflows are observed through interstellar absorption lines, which have velocities of about several hundred km s⁻¹ (Steidel et al. 1996; Shapley et al. 2003; Martin 2005; Weiner et al. 2009; Heckman et al. 2015). Mass loading factors around a few are frequently reported. In our simulations, the cool phase in general does not have such high velocities, and the mass loading is usually smaller than unity, especially for high Σ_gas. We here discuss the possible reasons for this apparent discrepancy.

We first note that observations usually focus on the star-bursting regime, such as (U)LIRGs and Lyman break galaxies, which we do not cover in our simulations (Figure 1). As mentioned before, radiation pressure may account for the potential high mass loading in those extreme environments, provided that the infrared light can be trapped by dust. Note that those bursting SF systems are very rare; even at higher redshifts where they are more common, they only account for ≲10% of cosmic SF density (Rodighiero et al. 2011).

Second, the observed mass loading factors have large uncertainties, due to the unconstrained metallicity, ionization fraction, geometry, etc. This can lead to an error bar as large as the mass loading factor itself. Recent work by Chisholm et al. (2016), which claims a better constraints on the above quantities, suggests a small η_m ≈ 0.1. This is, in fact, not inconsistent with our value for the highest Σ_gas.

On the simulation side, we note that even though we adopt high enough resolution to make sure the Sedov–Taylor phase of SNRs is well-resolved, the interaction between different gas phases may still not be sufficiently resolved, once the ISM becomes multiphase. To what extent this may affect the mass loading is not clear.

Another apparent discrepancy is that we find that the cooler phases dominate less of the outflow mass flux for higher ${\dot{{\rm{\Sigma }}}}_{\mathrm{SFR}}$. It is possible, however, that some cool clouds may form on the way due to thermal instability (Field 1965; Thompson et al. 2016), which we do not capture because of the relatively small box size. But again, how much cool mass would form under this mechanism is not clear.

Overall, there remains much uncertainty about the observational interpretation and theoretical modeling of cool outflows.

Stellar feedback is one of the key ingredients for galaxy formation in cosmological simulations and semi-analytic models (Somerville & Davé 2015). The general consensus is that strong feedback/outflows are needed to reproduce various observed scaling relations of galaxies (Springel & Hernquist 2003; Stinson et al. 2006; Hopkins et al. 2012; Hummels & Bryan 2012). Due to the limited resolution in cosmological simulations, however, a multiphase ISM and individual SNRs cannot be resolved. Ad hoc models are widely used in the field, with free parameters fine-tuned to match observations. Different groups adopt very different recipes, usually invoking shutting off cooling or hydrodynamics. To evaluate the real impact of feedback, however, a physically grounded model is necessary. While formulating such a model is beyond the scope of this paper, we compare the loading factors found in our work to current cosmological simulations, and briefly discuss the implications.

In cosmological simulations, a significant fraction of SNe energy is used to generate outflows, roughly 30%–50% (e.g., Oppenheimer & Davé 2008). Our η_E varies from 15% to 50%, broadly consistent with that fraction. To match the observed mass–metallicity relation of galaxies, a significant fraction of metals has to be driven out of the galaxies. Our η_met is about 40%–80%, in general agreement with what observations require and cosmological simulations adopt. The main difference is the mass loading of the outflows or, in other words, the amount of mass that those energy and metals couple to. Our η_m is in general smaller than what is adopted in cosmological simulations, by roughly an order of magnitude. A smaller mass loading is not surprising for the simulations where the multiphase ISM is resolved: SN blast waves tend to vent through the least-obstructed channel, thus preferentially heating and accelerating the tenuous phase (Cowie et al. 1981; Li et al. 2015). Cosmological simulations cannot resolve the multiphase ISM/outflows generally. Mixing the fast/tenuous and the slow/dense phase due to insufficient resolution would result in slower, cooler, and mass-loaded outflows. Therefore current cosmological simulations may not accurately predict the impacts of galactic outflows on the CGM and galaxy formation. From our simulations, hot outflows, though not carrying a great amount of mass, may be able to suppress the inflows given their vigor (Section 3.2), therefore restricting galaxy masses. The implications of tenuous yet vigorous hot outflows to galaxy formation and the CGM are not clear (although see recent attempts by Davé et al. 2016 & Fielding et al. 2017).