The Physical Nature of Starburst-driven Galactic Outflows

The primary difference between this work and earlier simulations in the CGOLS suite is the implementation of the stellar feedback. In Paper II, we described a method of clustered feedback where mass and energy were deposited in eight spherical regions within the disk, each with R = 150 pc. The values of the mass and energy injection were defined arbitrarily such that the volume-filling gas in the resulting hot outflow either would have high enough density to cool efficiently down to 10⁴ K at a relatively small radius (the "high mass-loading state") or would not cool efficiently (the "low mass-loading state"), as described in Thompson et al. (2016). We defined these mass and energy injection rates M_inj and E_inj in relation to the star formation rate: mass loading is quantified by β,

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{\dot{M}}_{\mathrm{inj}}}{{\dot{M}}_{\mathrm{SFR}}},\end{eqnarray} \tag{ 3 }$$

and energy loading is quantified by α, such that

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{inj}}=3\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,\alpha \left[\displaystyle \frac{{\dot{M}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right],\end{eqnarray} \tag{ 4 }$$

where we have assumed that there is one supernova per 100 M_⊙ of star formation, and each supernova generates 10⁵¹ erg of energy. In each state, we held the values of α and β steady for many millions of years.

In this work, in addition to using a significantly smaller cluster radius, R_cl = 30 pc, we relax these assumptions about the values of α and β. Now, we set the mass and energy injection rates in the clusters using physically motivated values determined by running a separate "superbubble" simulation. This simulation will be described in detail in a future paper (E. E. Schneider & E. C. Ostriker 2020, in preparation), but in brief, we use a Starburst99 single-burst stellar population synthesis model (Leitherer et al. 1999) to inject mass and energy into a box containing clouds that represent a well-resolved multiphase ISM, and we then track the interaction between phases as the resulting bubble driven by the cluster feedback propagates through the box (see Kim et al. 2017, for an example of such simulations). The ISM gas that is swept up as the bubble expands before breaking out of the disk can therefore contribute to the mass loading at early times.

In order to compare with our earlier, more idealized models, we use very large clusters for the feedback in this simulation—each cluster is assumed to host 10⁷ ${M}_{\odot }$ of stars—and each cluster is turned "on" for 10 Myr, as in our previous work. For the first 10⁵ yr that the cluster is on, the mass and energy injection rates are defined by the equations

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{inj}}=1.2\times \exp \left(-{\left[\displaystyle \frac{{\left(t-0.05\right)}^{2}}{(2\times {0.028}^{2})}\right]}^{2}\right)\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{inj}}=4\times {10}^{41}\displaystyle \frac{{t}^{2}}{({t}^{2}+0.035t)}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 6 }$$

where t is measured in Myr; thereafter, $\dot{M}=0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and $\dot{E}=3\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. These mass and energy injection rates are plotted in Figure 1.

Zoom In Zoom Out Reset image size

In terms of an injected α_inj and β_inj for each cluster, this corresponds to a mass loading that peaks at β_inj = 1.2 as the cluster is breaking out of the disk, accounting for interactions with swept-up cold clouds in the multiphase ISM (not present in the CGOLS simulation), then lowers to a value of β_inj = 0.1, which is approximately the average for the mass return from supernovae and stellar winds at late times (Leitherer et al. 1999). With this model, approximately 10% of the total mass injection for each cluster occurs during the initial breakout phase. Meanwhile, the energy loading is a steadily increasing function that reaches α_inj = 1 after 10⁵ yr, by which time a cluster of this size will have broken out of the disk, and energy losses within the R_cl = 30 pc injection region should be very low (Fielding et al. 2018). Note that these values for α_inj and β_inj are defined in terms of an assumed "star formation rate" of 1 ${M}_{\odot }$ yr⁻¹ in each cluster, defined as the total mass divided by the time the cluster is turned on, ${\dot{M}}_{\mathrm{SFR}}={10}^{7}\,{M}_{\odot }/{10}^{7}\,\mathrm{yr}$. These functions for $\dot{M}$ and $\dot{E}$ are fits to the actual measured $\dot{M}$ and $\dot{E}$ values at a radius corresponding to a cluster in the CGOLS simulation (R_cl = 30 pc) as measured in a superbubble simulation with a similar average gas density (E. E. Schneider & E. C. Ostriker 2020, in preparation). We emphasize that these are the injected mass and energy rates within the clusters; the effective α and β measured in the wind can be very different as a result of interactions with gas in the disk.

Clusters are turned on every million years at a rate corresponding to 20 ${M}_{\odot }$ yr⁻¹ of star formation from 5 to 35 Myr and 5 ${M}_{\odot }$ yr⁻¹ for the remainder of the simulation, which runs to 70 Myr. The clusters are randomly distributed in radius and azimuthal angle throughout the central R = 1 kpc of the disk and up to z = 5 pc above or below the midplane. (We note that this 1 kpc radius for cluster distribution is significantly larger than that of the central starburst in M82; despite using it to model our initial conditions, we are not attempting to replicate that system in detail.) After being turned on, each cluster rotates with the disk at a speed set by the gravitational potential. At the end of its 10 Myr life span, each cluster is turned off. This cluster life span was chosen to match our previous simulations, and we will investigate the effects of a longer cluster lifetime in future work.

While this model represents a rather extreme mode of centralized star cluster feedback, it is not wholly unphysical. In any given starburst, the energetics of the wind will be driven by the most massive clusters. In the CGOLS simulations, we do not resolve the details of the multiphase ISM, and therefore any smaller clusters that would not break out of the disk may be neglected for the purposes of studying the outflow—their primary effect would be to inject momentum into the ISM. On the other hand, including a large number of more moderately sized (M_* ∼ 10⁵–10⁶ ${M}_{\odot }$), longer-lived clusters would have the effect of significantly increasing the surface area of interaction between clusters and disk gas and could therefore plausibly increase the mass-loading rate substantially while lowering the effective value of α in the outflow for a given star formation rate. We will test the effects of these assumptions in future work by examining a simulation with a distribution of cluster masses.

Like the previous CGOLS simulations, the simulation presented in this work was carried out using the Cholla hydrodynamics code (Schneider & Robertson 2015), using a PLM reconstruction scheme, an HLLC Riemann solver, and an unsplit predictor-corrector integration method (Stone & Gardiner 2009). We also employ a dual-energy method to track the internal energy of the gas, given the high Mach numbers attained by cool gas in the outflow (Bryan et al. 1995; Schneider & Robertson 2017). Unlike the previous simulations, in addition to evolving the conserved quantities of density, momentum, and total energy, in this simulation, we also evolve a passive scalar variable, s, which is advected with the fluid. The primary purpose of this variable is to trace where the gas in the outflow originated. Gas that is present at the start of the simulation, i.e., disk and halo gas, is initialized with a scalar value of zero. Gas that is injected in the cluster regions is given a scalar value of 1. Thus, at any later time, the fraction of mass in a cell that was originally injected via a cluster can be determined by the value of the scalar in that cell; that is, if a cell has s = 0.5, then half the mass in the cell was injected. Therefore, the scalar directly represents the fraction of gas in a given cell that was injected within a cluster as "hot." Because our cluster injection accounts for unresolved interactions between the superbubble and the ISM, this value is distinct from the fraction of the mass that represents pure supernova ejecta, which varies as a function of β_inj.

Unless otherwise noted, all of the variables in this subsection should be understood to have an implied "h" subscript, denoting that they refer to the hot phase of the outflow.

Mass. The mass continuity equation can be written

$$\begin{eqnarray}&&{\partial }_{t}(\rho )+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=\dot{\rho }.\end{eqnarray} \tag{ 7 }$$

Here the source term on the right-hand side captures changes in mass per unit time per unit volume and can include both mass gained in the hot phase from shocked ambient material and cool gas destruction and losses from hot gas cooling out or mixing with the cool phase. Under the assumptions of a time-steady flow and spherical symmetry, this equation becomes

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}{d}_{r}({r}^{2}\rho {v}_{r})=\dot{\rho }.\end{eqnarray} \tag{ 8 }$$

Neither of these assumptions have to be true in an outflow; we will evaluate their validity in our simulations later. By integrating from r = 0 to r = r over a cone, we can relate the radial profiles of density and velocity to the net mass per unit time flowing outward within the cone,

$$\begin{eqnarray}&&{\rm{\Omega }}{r}^{2}\rho {v}_{r}={\dot{M}}_{\mathrm{net}},\end{eqnarray} \tag{ 9 }$$

where all variables are assumed to be functions of r. Here Ω is the actual solid angle within the cone that is filled with hot gas, i.e., allowing for a filling factor <1.

For clarity, we denote mass flow rates measured at an arbitrary point in the flow as ${\dot{M}}_{\mathrm{net}}$ (and similar rates for other quantities) to indicate that this allows for both gains and losses since the injection at the base of the flow. We denote the rate of mass and energy injected in the cluster region as ${\dot{M}}_{\mathrm{inj}}$ and ${\dot{E}}_{\mathrm{inj}}$, respectively.

Energy. The energy continuity equation is

$$\begin{eqnarray}&&{\partial }_{t}\left(\displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{1}{\gamma -1}P\right)+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\left(\displaystyle \frac{1}{2}{{\boldsymbol{v}}}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)\right)={\dot{{ \mathcal E }}}_{\mathrm{tot}}.\end{eqnarray} \tag{ 10 }$$

The energy source term on the right-hand side accounts for any energy added, in addition to energy lost from cooling of the gas or mixing into the cooler phases. Assuming steady-state spherical symmetry and integrating over the cone, we get

$$\begin{eqnarray}&&{\rm{\Omega }}{r}^{2}\rho {v}_{r}\left(\displaystyle \frac{1}{2}{v}_{r}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)={\dot{E}}_{\mathrm{net}},\end{eqnarray} \tag{ 11 }$$

where again, all variables are assumed to be functions of radius, and the right-hand side is the net energy per unit time flowing through the cone, consisting of the initial value in the injection region plus any gains and minus any losses.

Scalar mass. As described in Section 2, all gas in the simulation has an associated scalar value, s, which is passively advected with the fluid. For material injected in a cluster, s = 1, while material that was originally in the disk was assigned s = 0. Thus, the scalar value within a given cell identifies how much of the mass in that cell was originally injected in the clusters (when β_inj = 0.1, this injected mass is the same as supernova ejecta). We can write a continuity equation for the scalar density, ρ_s, that is identical to the density continuity equation,

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{s})+{\rm{\nabla }}\cdot ({\rho }_{s}{\boldsymbol{v}})={\dot{\rho }}_{s}.\end{eqnarray} \tag{ 12 }$$

Following the same procedure as above relates the scalar mass flux to the scalar density and velocity profile:

$$\begin{eqnarray}&&{\rm{\Omega }}{r}^{2}{\rho }_{s}{v}_{r}={\dot{M}}_{{\rm{s}},\mathrm{net}}.\end{eqnarray} \tag{ 13 }$$

The right-hand side measures the total scalar mass per unit time flowing outward in the cone, which will allow for any initial material injected, minus losses.

An advantage of tracking this scalar mass is that it allows us to examine the relationships between total measured outflow rates at a given radius, ${\dot{M}}_{\mathrm{net}}$ and ${\dot{E}}_{\mathrm{net}}$, and the mass and energy that was injected by the clusters at the base of the flow. We note that for r approaching the source region, in the idealized perfectly spherical case, we would have ${\dot{M}}_{\mathrm{net}}={\dot{M}}_{{\rm{s}},\mathrm{net}}={\dot{M}}_{\mathrm{inj}}$. Taking the ratio of Equation (13) to Equation (9) gives

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{s}}{\rho }=\displaystyle \frac{{\dot{M}}_{{\rm{s}},\mathrm{net}}}{{\dot{M}}_{\mathrm{net}}},\end{eqnarray} \tag{ 14 }$$

where ρ_s/ρ ≡ s is the definition of the scalar. We therefore have

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{net}}=\displaystyle \frac{{\dot{M}}_{{\rm{s}},\mathrm{net}}}{s},\end{eqnarray} \tag{ 15 }$$

and in the case that none of the hot gas has cooled out, ${\dot{M}}_{{\rm{s}},\mathrm{net}}={\dot{M}}_{\mathrm{inj}}$. Thus, if s < 1 in the hot medium at some large radius and there have been negligible losses to cooling, ${\dot{M}}_{\mathrm{net}}$ will exceed the initial injected value due to mixing in of originally cool gas. Likewise, the ratio of Equation (11) to Equation (13) gives

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{net}}=\left(\displaystyle \frac{1}{2}{v}_{r}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)\displaystyle \frac{{\dot{M}}_{{\rm{s}},\mathrm{net}}}{s}.\end{eqnarray} \tag{ 16 }$$

Momentum. Momentum continuity states

$$\begin{eqnarray}&&{\partial }_{t}(\rho {\boldsymbol{v}})+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{\boldsymbol{v}}+{\boldsymbol{I}}P)=-\rho {\rm{\nabla }}{\rm{\Phi }}+\dot{q}.\end{eqnarray} \tag{ 17 }$$

While there is no (significant) momentum directly injected by the clusters, there are potentially momentum source and sink terms associated with mixing. For example, when hot gas is mixed into the cool medium (thereafter cooling), there will be a momentum sink term (negative $\dot{q}$) for the hot medium. If a small enough amount of cool gas is mixed into the hot such that it does not subsequently cool, this will be a positive $\dot{q}$ for the hot medium. The spherical and steady-state assumptions yield

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}({r}^{2}\rho {v}_{r}^{2})+{\partial }_{r}(P)=-\rho {\partial }_{r}({\rm{\Phi }})+\dot{q}.\end{eqnarray} \tag{ 18 }$$

There is no straightforward way to integrate the pressure term, so we will continue with the differential version of the momentum equation. A bit of algebra translates it into a more usable form,

$$\begin{eqnarray}&&{v}_{r}{\partial }_{r}({v}_{r})+\displaystyle \frac{{\partial }_{r}(P)}{\rho }=-{\partial }_{r}({\rm{\Phi }})-\displaystyle \frac{\dot{\rho }}{\rho }{v}_{r}+\displaystyle \frac{\dot{q}}{\rho },\end{eqnarray} \tag{ 19 }$$

from which we can assess the importance of each term in setting the velocity profile.

Entropy. Substituting Equation (9) into Equation (11) gives

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{v}_{r}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }=\displaystyle \frac{{\dot{E}}_{\mathrm{net}}}{{\dot{M}}_{\mathrm{net}}}.\end{eqnarray} \tag{ 20 }$$

Meanwhile, from Equation (19), we have

$$\begin{eqnarray*}&&{\partial }_{r}\left(\displaystyle \frac{1}{2}{v}_{r}^{2}\right)+\displaystyle \frac{{\partial }_{r}(P)}{\rho }=-{\partial }_{r}({\rm{\Phi }})-\displaystyle \frac{\dot{\rho }}{\rho }{v}_{r}+\displaystyle \frac{\dot{q}}{\rho }.\end{eqnarray*}$$

Combining the two gives a differential equation describing the entropy profile:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{r}\left(\displaystyle \frac{\mathrm{ln}[P{\rho }^{-\gamma }]}{\gamma -1}\right) & = & \displaystyle \frac{1}{P}\left[\rho {\partial }_{r}\left(\displaystyle \frac{{\dot{E}}_{\mathrm{net}}}{{\dot{M}}_{\mathrm{net}}}\right)+\rho {\partial }_{r}{\rm{\Phi }}\right.\\ & & +\ \left.\displaystyle \frac{{v}_{r}}{{r}^{2}}{\partial }_{r}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{net}}}{{\rm{\Omega }}}\right)-\dot{q}\right].\end{array}\end{eqnarray} \tag{ 21 }$$

In general, we do not expect the gravitational term on the right-hand side to affect the properties for the hot phase, so it may be omitted.

Effective α and β. We can use the above relationships to determine an "effective" α and β at each radius. (Recall that α is the ratio of energy in the outflow to the total energy injected by supernovae, and β is the ratio of mass flux in the outflow to star formation rate.)

We can define an effective mass loading in the wind at any radius by

$$\begin{eqnarray}&&{\beta }_{\mathrm{eff}}\equiv \displaystyle \frac{{\dot{M}}_{\mathrm{net}}}{{\dot{M}}_{\mathrm{SFR}}}={\beta }_{\mathrm{inj}}\displaystyle \frac{{\dot{M}}_{\mathrm{net}}}{{\dot{M}}_{\mathrm{inj}}}=\displaystyle \frac{{\beta }_{\mathrm{inj}}}{s}\displaystyle \frac{{\dot{M}}_{{\rm{s}},\mathrm{net}}}{{\dot{M}}_{\mathrm{inj}}}.\end{eqnarray} \tag{ 22 }$$

Similarly,

$$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}\equiv \displaystyle \frac{{\dot{E}}_{\mathrm{net}}}{{E}_{\mathrm{SN}}({\dot{M}}_{\mathrm{SFR}}/{m}_{* })}={\alpha }_{\mathrm{inj}}\displaystyle \frac{{\dot{E}}_{\mathrm{net}}}{{\dot{E}}_{\mathrm{inj}}},\end{eqnarray} \tag{ 23 }$$

where E_SN is the energy injected per supernova, m_* is the mass of stars formed per supernova, and ${\dot{M}}_{\mathrm{SFR}}$ is the total star formation rate. We can use Equation (20) to rewrite this as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}=\left(\displaystyle \frac{1}{2}{v}_{r}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)\displaystyle \frac{{m}_{* }}{{E}_{\mathrm{SN}}}\displaystyle \frac{{\dot{M}}_{{\rm{s}},\mathrm{net}}}{{\dot{M}}_{\mathrm{inj}}}\displaystyle \frac{{\beta }_{\mathrm{inj}}}{s}.\end{eqnarray} \tag{ 24 }$$

Hence, the scalar variable allows us to distinguish between the injected mass and energy rates set by our cluster prescription and the effective mass and energy loading actually measured in the simulation at any point in the wind.

In the absence of cooling, ${\dot{E}}_{\mathrm{net}}={\dot{E}}_{\mathrm{inj}}$, and there is no reduction of total scalar mass in the hot medium with distance, so that ${\dot{M}}_{{\rm{s}},\mathrm{net}}={\dot{M}}_{{\rm{s}},\mathrm{inj}}={\dot{M}}_{\mathrm{inj}}$. In this case, we have α_eff = α_inj and β_eff = β_inj/s. That is, the effective mass loading can be different from the injected value due to mixing of previously cool gas into the hot medium.

Completely analogous equations to those of Section 3.1 could be written for the cool gas. However, beyond the injection region, we know that for mass, scalar mass, and momentum, there are no "exogenous" sources or sinks. Thus, any losses from the hot must be gains for the cool, and vice versa. That is, ${\dot{\rho }}_{{\rm{c}}}=-{\dot{\rho }}_{{\rm{h}}}$, ${\dot{\rho }}_{s,{\rm{c}}}=-{\dot{\rho }}_{s,{\rm{h}}}$, and ${\dot{q}}_{{\rm{c}}}=-{\dot{q}}_{{\rm{h}}}$. For the energy, however, the total energy summed over the phases is not conserved due to radiation.

In order to accelerate the cool gas, momentum must be transferred to it from the hot phase. We assume that the cool gas is supersonic, and therefore we can neglect the pressure term in the "cool" version of Equation (18). The continuity equation for momentum then states

$$\begin{eqnarray*}&&{\left[\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}({r}^{2}\rho {v}_{r}^{2})+\rho {\partial }_{r}({\rm{\Phi }})\right]}_{{\rm{c}}}=-{\left[\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}({r}^{2}\rho {v}_{r}^{2})+{\partial }_{r}(P)\right]}_{{\rm{h}}},\end{eqnarray*}$$

where the source term on the right-hand side is the rate of momentum transferred per unit volume from the hot gas (subscript h) to the cool gas (subscript c). This is ${\dot{q}}_{{\rm{c}}}=-{\dot{q}}_{{\rm{h}}}$, and we have neglected the gravitational potential term for the hot medium. We can rewrite the right-hand side of this equation as

$$\begin{eqnarray*}&&=\,-{\left[{v}_{r}\displaystyle \frac{\dot{{\rho }_{s}}}{s}+\rho {v}_{r}{\partial }_{r}({v}_{r})+{\partial }_{r}(P)\right]}_{{\rm{h}}},\end{eqnarray*}$$

where we have used ${\dot{\rho }}_{{\rm{h}}}={\dot{\rho }}_{s,{\rm{h}}}/{s}_{{\rm{h}}}$. Conservation of the scalar variable combined with mass conservation for the cool (see Equation (8)) yields

$$\begin{eqnarray*}&&-{\dot{\rho }}_{s,{\rm{h}}}={\dot{\rho }}_{s,{\rm{c}}}={\left[\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}({r}^{2}{\rho }_{s}{v}_{r})\right]}_{{\rm{c}}}.\end{eqnarray*}$$

Thus, we can rewrite the momentum equation as

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{r}{\left({r}^{2}\rho {v}_{r}^{2}\right)}_{{\rm{c}}} & = & -{r}^{2}{\rho }_{{\rm{c}}}{\partial }_{r}({\rm{\Phi }})+\displaystyle \frac{{v}_{r,{\rm{h}}}}{{s}_{{\rm{h}}}}{\partial }_{r}{\left({r}^{2}{\rho }_{s}{v}_{r}\right)}_{{\rm{c}}}\\ & & -\ {r}^{2}{\left[\rho {v}_{r}{\partial }_{r}({v}_{r})+{\partial }_{r}(P)\right]}_{{\rm{h}}}.\end{array}\end{eqnarray} \tag{ 25 }$$

The left-hand side is the rate of increase of momentum flux in the cool gas as a function of r, while the first term on the right-hand side accounts for deceleration of the cool gas due to gravity. The second term on the right-hand side describes the rate of increase of the cool momentum flux transferred from the hot medium at a rate proportional to the rate of increase of scalar in the cool gas resulting from mixing in the hot gas. The last term on the right-hand side is the thermal pressure work from the hot on the cool. The second-to-last term on the right-hand side describes the deceleration of the hot gas; this is the term that represents the ram pressure work of the hot on the cool. If the hot gas is a supersonic wind, its deceleration and pressure gradients will be small.

If we assume that the mixing is the dominant source term for the cool gas momentum flux, and that the hot wind velocity and scalar have reached constant values, Equation (25) gives a linear relationship between the velocity of the cool gas and its scalar value,

$$\begin{eqnarray}&&{\rho }_{{\rm{c}}}{v}_{r,{\rm{c}}}^{2}=\displaystyle \frac{{v}_{r,{\rm{h}}}}{{s}_{{\rm{h}}}}{\rho }_{s,{\rm{c}}}{v}_{r,{\rm{c}}},\end{eqnarray} \tag{ 26 }$$

or

$$\begin{eqnarray}&&{v}_{r,{\rm{c}}}=\displaystyle \frac{{v}_{r,{\rm{h}}}}{{s}_{{\rm{h}}}}{s}_{{\rm{c}}},\end{eqnarray} \tag{ 27 }$$

which is normalized by the velocity of the hot gas.

As with previous CGOLS simulations, we endeavor to make the most of this (expensive) simulation by running in two different states. We begin the cluster feedback at 5 Myr in the high star formation rate state, with ${\dot{M}}_{\mathrm{SFR}}=20\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, equivalent to turning on two clusters every million years for 30 Myr. After that point, we turn on clusters at a rate corresponding to ${\dot{M}}_{\mathrm{SFR}}=5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, our low star formation rate state. These values were chosen to match our previous simulations. In the following, we primarily present results for two characteristic times: at 35 Myr, when the star formation rate has been high for its maximum time, and at 65 Myr, when it has been low for the same length of time.

Figure 2 shows a rendering of the density field at 35 Myr that highlights many general features of the outflow. At small radii, high-density clouds of disk gas are being driven out by the central clusters. At larger radii, the outflowing material is more diffuse, but there are still high-density filaments permeating the volume-filling low-density phase. These high-density filaments are clearly associated with the cool gas in the outflow, as can be seen in Figure 3, which shows off-axis density and temperature projections of the full volume at 35 Myr. Movies showing the density and temperature projections for the full time evolution can be found online.⁶

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the projections, slices through the box reveal interesting relationships between the density, velocity, temperature, and scalar values of gas in the outflow. In Figure 4, we show density, temperature, velocity, and scalar slices along the x−z plane during the high star formation rate state, while Figure 5 shows the same slices during the low star formation rate state. A few salient features of the outflow can be determined immediately upon inspection of these slices. First, at both times, the outflow is multiphase, characterized by a volume-filling hot phase at T > 10⁶ K punctuated by small, dense clouds of cool gas at T ∼ 10⁴ K. There is a clear correlation between the gas density, velocity, and temperature, with the lower-density channels corresponding to the hottest, fastest-moving gas. The outflow features are roughly biconical, though determining an opening angle is not straightforward and would likely depend on which snapshot was being examined. Regardless, the opening angle appears to be large, and the outflow does not show evidence of any sort of fountain, with cool gas raining back down onto the disk at larger radii. We note that this may be a result of our choice to position the clusters within the central R = 1 kpc of the disk, and that our limited horizontal volume prevents us from assessing whether a fountain flow would arise at large angles at larger radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast with our previous simulations, at no point in this simulation does the volume-filling hot phase become mass-loaded enough to undergo significant radiative cooling. While the simulations in Paper II employed an arbitrary β_inj = 0.4 for the hot phase, in this simulation, we set the mass loading for the injected material to a physically motivated value determined by massive star winds, supernova ejecta, and cluster breakout. This average β is close to 0.1 most of the time, and mixing with the disk gas at small radii is not efficient enough to increase the density of the majority of the hot gas enough for it to cool. This does not mean none of the hot gas cools, however. As we will demonstrate later, there is strong evidence that hot gas that does interact with cool clouds creates a mixed phase with intermediate temperatures, which can cool under the right circumstances.

Comparing Figures 4 and 5, we see that overall, the simulation during the higher star formation rate state is characterized by a higher-density, lower-velocity outflow, and the biconical outflow region is filled with more clouds of cool gas. At late times, when the star formation rate is lower, there is not enough interaction between the clusters and the disk gas to contribute significant mass loading to the outflow, as we will show directly when we examine the fluxes in Section 4.3. This lower mass loading leads to a lower-density hot phase and fewer cool clouds in the outflow.

Throughout the remainder of the results, we will examine the properties of the gas split by phase. We define three phases of interest: hot gas (T > 5 × 10⁵ K), intermediate-temperature gas (2 × 10⁴ K < T < 5 × 10⁵ K), and cool gas (T < 2 × 10⁴ K). Most of our analysis focuses on the hot and cool phases, as they are the longest-lived and represent the majority of the mass in the outflow (see Section 4.5).

In Figures 6 and 7, we show radial profiles for the hot and cool phases, respectively, for a number of physical parameters of interest: the number density n, passive scalar s, pressure P, temperature T, radial velocity v_r, sound speed c_s, Mach number ${ \mathcal M }$, and entropy K ∝ Pρ^−γ, where γ is the adiabatic index of the gas, taken to be 5/3 through this work. Note that in Figures 6 and 7, we plot the entropy as $\frac{K}{{k}_{b}}={{Tn}}^{-\tfrac{2}{3}}$. Each profile is measured within a cone with a half-opening angle Ω = 30° above and below the disk, as shown with dashed white lines in the density panels of Figures 4 and 5. Within the cone, we expect the properties of the outflow to be approximately spherically symmetric once we reach a radius equal to our cluster seeding radius within the disk, R = 1 kpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each panel, we show the median and mean values, as well as the first and third quartiles, to give an idea of the spread in the gas properties. We use density-weighted values for all quantities, unless otherwise noted. For example, the velocity average is calculated as $\langle {v}_{r}\rangle =\sum ({v}_{r}n)/(\langle n\rangle {N}_{\mathrm{cells}})$, with $\langle n\rangle =\sum (n)/{N}_{\mathrm{cells}}$ for all cells within the cone that meet the relevant temperature criterion. The choice of density-weighted versus volume-weighted averages has no effect on the cool gas profiles, but there are some differences for the hot phase. Thus, we also include the volume-averaged median quantities on the relevant panels in Figure 6. However, we note that the largest differences are at small radii, where we do not interpret our results as complete, because we are still within the gain region. The biggest differences between the mass- and volume-weighted profiles arise in the scalar variable, which tends to be higher in the hot gas if we use a volume-weighted average for the hot phase, and the velocities, which are also higher, on average. We further explore the relationship between the scalar and velocity in Section 4.6.

The hot phase. On average, the profiles for the hot phase follow relationships that are close to, but not exactly, what one would expect given adiabatic expansion of the gas. For example, Chevalier & Clegg (1985) calculated the radial solution for a hot outflow given mass and energy injection within a spherical region. In that model, density decreases as r⁻², while velocity quickly asymptotes to a value set by ${v}_{\mathrm{term}}=\sqrt{2{\dot{E}}_{\mathrm{inj}}/{\dot{M}}_{\mathrm{inj}}}$, as would be predicted by Equation (20) (assuming mass and energy are conserved after injection). With an adiabatic index γ = 5/3, adiabatic cooling implies that pressure decreases as r^−10/3 (as Equation (21) implies when source terms are zero), and thus temperature decreases as r^−4/3. We have plotted these relationships with black lines in Figure 6 to demonstrate their deviations, with the adiabatic profiles normalized to the measured values at 1 kpc. As the first panel shows, the density is falling off with radius but at a shallower rate than r⁻². This reflects the fact that mass is being added to the hot phase as a function of radius, as we will show in Section 4.3. Similarly, the pressure and temperature are decreasing with radius but at a slower rate than is implied by pure adiabatic expansion.

If we account for the mass addition to the hot phase by allowing a shallower density profile, we can accommodate changes in the pressure and temperature profiles as well. For example, we can scale the radial dependence of the density profile such that it is an arbitrary function of r that provides a good fit to the density profile. In the case of the plots in Figure 6, we use n ∝ r^f, with f = −0.05r − 1.08, and r measured in kiloparsecs. Then P ∝ r^fγ, and T ∝ r^f(γ−1). We have plotted these additional relationships as black dashed lines in Figure 6. While the slope of the density profile is now an arbitrary function, the slopes of the pressure and temperature profiles are set by the adiabatic physics and provide a much better fit to the data than the pure expansion wind model. In this framework, the hot phase of the wind can be understood entirely as adiabatic expansion with a mass source that depends on mixing in gas from cooler phases.

Although these new scalings do provide a much better fit to the density, pressure, and temperature profiles, they still underestimate the pressure and temperature. The entropy profiles displayed in the final panel of Figure 6 provide an explanation. While the volume-averaged entropy profile for the hot gas is quite flat, the density-weighted median entropy in the hot phase is an increasing function of radius. This increasing slope reflects the fact that low-entropy gas from the cool phase is being mixed into the hot gas, with the highest-density hot gas having experienced the most mixing. As the cool gas is added to the hot, its entropy rises, as would be expected in the case where the cool gas experiences a shock. This rising entropy profile then violates the assumption that the relationship between n, P, and T is adiabatic and, in particular, will result in a flatter temperature slope than predicted by adiabatic physics.

The expectation for the terminal value of the velocity is set by our choice of input $\dot{E}$ and $\dot{M}$. If we use our values of α_inj and β_inj, we find that ${v}_{\mathrm{term}}=\sqrt{2\dot{E}/\dot{M}}\approx 3000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, about a factor of 2 too high. However, if we measure the effective α and β in the wind at R = 1 kpc, we find that α_eff = 0.5, while β_eff = 0.2 (see Section 4.4), which leads to a terminal velocity of v_term ≈ 1500 km s⁻¹, approximately the measured velocity in the volume-averaged hot phase. The hot gas velocity as measured by the density-weighted median is slightly lower, which can naturally be explained by higher energy losses associated with higher-density gas (resulting in a lower α_eff).

Another direct indication of the mass transfer between phases is the radial profile of the scalar variable in Figures 6 and 7. By a radius of 1 kpc, the median density-weighted scalar value of the hot phase has already decreased from 1 to a value of 0.4, indicating significant mixing from the cool phase. Similarly, the median scalar value in the cool phase has risen from zero to 0.15. At first look, the relatively flat value of the density-averaged scalar in the hot phase combined with the rising scalar value in the cool phase might seem to be discrepant. However, the volume-averaged scalar in the hot phase tells the full story. Here we see a rapid decrease down to s = 0.5 at R < 2 kpc. At larger radii, s continues to decrease, although at a slower rate. This is consistent with rapid mixing between the hot and cool gas as the clusters initially accelerate cool material out of the disk, followed by more gradual mixing at larger radii as the volume density of cool gas goes down. The outflow rates presented in Section 4.3 bear out this explanation. However, we caution that this interpretation is dependent on the assumption that the outflow is steady-state. While we have purposely selected a time for this analysis where that appears to be the case, it is certainly not true for the entire simulation.

The cool phase. The profiles for the cool phase can be interpreted as a radially expanding isothermal gas. Density falls off approximately as n ∝ r⁻², and the temperature is a constant value set by our temperature floor. (In reality, we expect photoionization from the local ionizing sources within the galaxy to play this role at small radii and from the cosmic UV background to provide heating of the cool medium at larger radii.) The pressure also falls off approximately as P ∝ r⁻². However, the pressure in the cool phase is consistently lower than that in the hot phase, with larger differences at smaller radii. We highlight this by plotting the fit to the hot phase pressure from Figure 6 on the pressure panel in Figure 7.

Physically, the lower pressure of the cool phase can be understood if the cooling time for intermediate-temperature gas in individual cool clouds is shorter than the time it takes for them to equilibrate to the background wind pressure via pressure waves or shocks. In this scenario, cool gas in the simulation continuously interacts with the hot phase via mixing or shocks. If the interaction heats the cool gas to temperatures at the low end of the intermediate-temperature range, the relatively high density leads to very fast cooling. Because the interactions between the hot and cool phase are frequent, there is insufficient time for the cool gas to equilibrate to the background hot pressure before experiencing a subsequent interaction. In particular, given the average values of the profiles at 1 kpc, the sound-crossing time for cool gas at our resolution limit Δx = 5 pc is t_sc ≈ 5 × 10⁵ yr, which is an order of magnitude longer than the cooling time of gas at a density of n ∼ 1 cm⁻³ (see Section 4.5). The relevant timescale may instead be the slightly shorter cloud-crushing time, t_cc = (n_c/n_h)^1/2 (Δx/v_h) ≈ 10⁵ yr (Klein et al. 1994), but this is still longer than the cooling time. In either case, the cool gas would then experience compression from the hot phase, as it is out of equilibrium. The fact that the slopes of both the density and pressure for the cool phase are slightly shallower than r⁻² could be explained by this compression. We note that it is plausible that with higher resolution, the cool clouds may further "shatter" (McCourt et al. 2018), leading to faster equilibration with the pressure of the hot phase. This explanation is consistent with our previous work investigating cool clouds at much higher resolution (Schneider & Robertson 2017).

An important implication of this work is shown in the velocity panel of Figure 7. In this panel, we see that the velocities of the cool gas are a rising function of radius, with maximum velocities reaching 1000 km s⁻¹. In fact, almost all cool gas at large radius is moving with velocities of at least 400 km s⁻¹, which is approximately the escape velocity of the system. The scalar value of the cool gas is also increasing with radius; we will explore the origin of the cool gas acceleration in Section 4.6.

As discussed in Section 3, the radial profiles of the quantities shown in Section 4.2 are related to total mass and energy outflow rates. We have measured these outflow rates in the simulation within the same cone used to calculate the radial profiles. Each flux is averaged within a radial bin of width Δr = 0.125 kpc. So, to calculate the mass flux, for example, we sum the total $\dot{M}$ in a given phase within the cone between r and r + Δr and then divide by the bin width, $\dot{M}\,=\sum ({{Mv}}_{r})/{\rm{\Delta }}r$, where the sum is taken over individual cells. Our cone opening angle corresponds to a solid angle of Ω = 2 × 2π[1 − cos(30°)] ≈ 1.68, or just over 1/8 of the full spherical area. If the outflow were perfectly spherical, then we might expect these fluxes to be 1/8 of the total. However, when we measure over the full 4π, we find fluxes that are only a factor of ∼3 larger. Hence, we conclude that the outflow is preferentially biconical, with a modest opening angle.

We show the radial outflow rates in the biconical region for mass, scalar mass, momentum, and energy at 35 Myr in Figure 8 and 65 Myr in Figure 9. For the momentum outflow rate, we only include the kinetic component, that is, $\dot{p}=\sum ({{Mv}}_{r}^{2})/{\rm{\Delta }}r$. The energy outflow rate is calculated using the total energy and including both the kinetic and thermal components, $\dot{E}\,=\sum ({{Ev}}_{r})/{\rm{\Delta }}r$, where $E=\rho \left(\tfrac{1}{2}{v}^{2}+{e}_{\mathrm{th}}\right)$, and e_th is the specific internal energy.⁷ The outflow rates are split according to the temperature of the gas into hot (T > 5 × 10⁵ K), intermediate (2 × 10⁴ K < T < 5 × 10⁵ K), and cool (T < 2 × 10⁴ K) phases. Including the intermediate phase allows us to better understand the amount of mixing that is happening in the outflow, since much of the gas in this temperature range has a relatively short cooling time and is therefore a transient phase (see Section 4.5). We scale down the vertical axes by a factor of 4 in Figure 9 relative to Figure 8 to better compare given the lower "star formation rate" at late times: ${\dot{M}}_{\mathrm{SFR},\mathrm{early}}\,=20\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ versus ${\dot{M}}_{\mathrm{SFR},\mathrm{late}}=5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Accounting for that normalization, we see that the total outflow rates are fairly consistent between the early and late stages of the outflow, but the intermediate and cool fractions are much smaller at lower outflow rates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 8 and 9 also demonstrate some natural relationships for the outflow. Most of the energy and momentum are carried by the hot phase of the wind at all times, as would be expected for an energy-driven wind model. At small radii and early times, when the average density in the wind is high, there are comparable amounts of mass in both the hot and cool phases (particularly if we include the intermediate-temperature gas with the cool). There are times when the mass outflow rate in the cool gas is larger than that of the hot gas (not shown in these figures), but they are always within a factor of ∼2 in our simulations. That is to say, the mass outflow rate is never dominated by the cool gas. In each of the energy outflow rate panels, we have added a dotted line showing the total energy available given our input energy injection rate and assuming no losses (so, α_eff = 1). At both early and late times, we see that the total energy outflow rate is below this line by a factor of ∼2. It is interesting to note that the total momentum outflow rate in the cone, despite representing only a fraction of the total outflow, is nevertheless within a factor of 2 of L/c that would be expected for a freely expanding hot wind, which would be $\dot{p}\approx 4000\,{M}_{\odot }\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{yr}}^{-1}$.

We can also use these outflow rates to assess the validity of the time-steady assumption that went into our analytic model. At 35 Myr, the total outflow rates are all fairly flat as a function of radius outside of R = 2–3 kpc. Thus, we interpret the outflow at this time as an approximately time-steady flow. The flatness does depend on phase, however. While the mass outflow rate of the hot phase is increasing with radius, the cool phase is decreasing, and the intermediate-temperature gas remains roughly constant. This is a clear indication that mass is being transferred from the cool phase to the hot as a result of mixing between the phases. In addition, the slight decrease in the total mass outflow rate between R = 4 and 8 may indicate the presence of a fountain flow in the cool phase, as lower-velocity gas drops out at higher radii (see, e.g., Kim & Ostriker 2018).

Returning to the momentum outflow rates in Figure 8, we see that although the total momentum outflow rate is roughly constant outside of R ∼ 4 kpc, the momentum outflow rate in the cool phase is actually decreasing. This may seem counterintuitive, since the profiles in Figure 7 showed that the velocity of the cool gas is increasing as a function of radius, and we have attributed this increase to a transfer of momentum from the hot phase. In addition, the energy outflow rate in the hot phase stays roughly constant outside of R ∼ 2 kpc, which would seem to call into question an explanation in which the mixed intermediate-temperature gas radiates away its thermal energy and adds mass to the cool phase with an increased velocity (i.e., Gronke & Oh 2018). In that case, we might expect the hot phase to show a decreasing energy flux with radius, reflecting this lost energy.

To shed some light on the matter, we show in Figure 10 versions of the momentum and energy outflow rates normalized by the mass outflow rate in each phase. This corresponds to a velocity as a function of radius in the top panel and a specific total energy (${e}_{{\rm{T}}}=\tfrac{1}{2}{v}^{2}+{e}_{\mathrm{th}}$) as a function of radius in the bottom panel. With the total mass normalized out, we see that the gas in each phase behaves exactly as we would expect given the mixing model we have outlined. While the velocity in the hot phase stays roughly constant as a function of r, the velocity in the intermediate and cool phases is consistently rising. The specific energy of the cool and intermediate phases is also rising, reflecting the higher kinetic energy at larger radii. Meanwhile, the specific energy in the hot phase decreases slightly as a function of radius, reflecting the fact that low-energy cool gas is being mixed in with the hot, lowering the average specific energy of the hot gas as a function of radius. The energy loss from cooling by the hot phase is fractionally quite small, such that the outward decrease in the specific energy is mainly due to the outward mass increase of the hot, and thus we do not see a decrease in Figure 8.

Zoom In Zoom Out Reset image size

Using the profiles measured in Section 4.2 and the relationships defined in Section 3, we can calculate the effective mass and energy loading in the different phases of the outflow. Using Equation (22), we can straightforwardly calculate the effective mass loading in the hot outflow as the injected mass-loading factor (β_inj ∼ 0.1) divided by the passive scalar, s, and multiplied by the factor ${\dot{M}}_{s,\mathrm{net},{\rm{h}}}/{\dot{M}}_{\mathrm{inj}}$. Because we have measured the outflow rates in a cone, we must "correct" the denominator, ${\dot{M}}_{\mathrm{inj}}$, which represents the total injected mass rate (all hot), to the total mass injected in the cone, which is just $\dot{{M}_{s}}$ (all phases). This correction factor is plotted in Figure 11. Carrying out this calculation, we find β_eff,h = 0.17, approximately a factor of 6 below the star formation rate. We note, however, that β_eff,h ranges from 0.1 at small radii to 0.2 at 8 kpc as a result of the additional mass loading of the hot phase taking place as some of the cool gas is mixed in. We plot these values as a function of radius in Figure 12. Similarly, we can calculate the effective energy loading in the hot outflow using Equations (23) or (24), again corrected for the amount of mass injected into the cone. Using Equation (23), we find a roughly constant value of energy loading within the hot phase, α_eff = 0.4–0.5, implying that ∼50%–60% of the input energy has been radiated away in the outflow.⁸ While this implies a robust energy-driven wind, we caution against overinterpreting the exact measured values, since they may be sensitive to the cluster input scheme, in the sense that having only a few massive clusters minimizes the interaction of supernovae with the ISM and thus is likely to maximize the value of α_eff while potentially minimizing β_eff. We will explore this possibility in future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following the same procedure (applying Equations (22) and (24) corrected for the mass injected into the cone), we also calculate the effective mass and energy loading for the cool phase, as shown in Figure 12. The loading factors for the hot and cool phases demonstrate the same relationships seen in the profiles and outflow rates: at small radii, a significant fraction of the outflow is in the cool phase, but at larger radii, that gas has mostly been incorporated into the hot phase or dropped out of the flow. As a result, the mass-loading factor for the hot phase increases, while the cool phase decreases. As demonstrated by α_eff,c, at the maximum point around r = 4 kpc, the cool phase has coupled a few percent of the total energy available.

Given these mass- and energy-loading factors, we can now explain why this CGOLS simulation does not lead to the large-scale cooling of the hot wind seen in Paper II. Following Thompson et al. (2016), we estimated the cooling radius of the hot wind as

$$\begin{eqnarray}&&{r}_{\mathrm{cool}}\simeq 4\,\mathrm{kpc}\displaystyle \frac{{\alpha }^{2.13}}{{\beta }^{2.92}}{\mu }^{2.13}{R}_{0.3}^{1.79}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{4\pi }}{{\dot{M}}_{\mathrm{SFR},10}}\right)}^{0.789}.\end{eqnarray} \tag{ 28 }$$

In contrast with our previous work, where mass- and energy-loading values were assumptions of the model, we can now use the measured values of α and β in the hot phase to calculate the derived cooling radius for the wind. With α = 0.5 and β = 0.2, and using R = 1 kpc for the injection radius, we calculate ${r}_{\mathrm{cool}}=170\,\mathrm{kpc}\,{{\rm{\Omega }}}_{4\pi }^{0.789}$. While we have demonstrated that Ω may be significantly smaller than 4π, a conservative lower limit of 1/8 still gives r_cool = 33 kpc, well outside the bounds of our simulation box. Thompson et al. (2016) also estimated a minimum β below which cooling of the hot wind would not take place at any radius (see their Equation (7)), which with these same values would be β = 0.17. Given that Ω is likely underestimated in this case, the volume-filling hot wind in this simulation is therefore unlikely to cool. Nevertheless, we note that this cooling radius depends strongly on the exact values of α, β, and Ω. Using the same β and Ω but a slightly lower α = 0.1, r_cool ∼ 1 kpc, so we do not conclude that the simulation analyzed in this work rules out cooling of the volume-filling hot wind in other scenarios.

In addition to allowing us to calculate effective loading factors within the cone, the correction factor ${\dot{M}}_{s}/{\dot{M}}_{\mathrm{inj}}$ also gives us an estimate for what fraction of the total outflow is contained within the cone. As shown in Figure 11, this factor steadily increases at small radii, plateauing at around 0.35 at r ∼ 3–4 kpc. This means that roughly 1/3 of the total injected cluster mass is escaping within the biconical region, more than twice the rate that would be expected in a spherical outflow, implying that the disk is having a significant collimating effect. Of course, this analysis would not capture the impact of a higher rate of mass loading outside the biconical region or a lower energy loading. In fact, when we look at the mass, momentum, and energy outflow rates outside of the bicone, we do find higher values of β_eff, particularly for the cool phase at small radii. However, we also find that the total mass outflow rate falls from $\dot{M}=6\,$ to $5\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$ between r = 1 and 5 kpc within the full sphere, implying that the effects of a cool fountain are perhaps more significant at larger angles, making the large-angle contributions to the total mass outflow rate less important at larger radii.

In addition to looking at the properties of the outflow as a function of radius, we can use 2D histograms to better characterize the relationships between different physical properties. In this section, we explore how density, velocity, and temperature relate to one another.

In Figure 14, we show a phase diagram of the gas within the biconical region during the high star formation state. The bins are weighted by total mass to highlight the regions where the bulk of the mass resides. We additionally exclude the regions within 0.5 kpc of the midplane to avoid including pure disk gas or the interior of a cluster in our analysis. The histogram shows two clear peaks, demonstrating the primarily two-phase nature of the outflow. Integrating these regions, we find that there is slightly more total mass in the cool phase than the hot: 2.5 × 10⁷ and 1.7 × 10⁷ ${M}_{\odot }$, respectively. The ratio reverses during the low star formation state, however. The same integration at 65 Myr yields 3.1 × 10⁶ and 4.6 × 10⁶ ${M}_{\odot }$ for cool vs. hot. These total values are consistent with there being about four times more outflowing gas during the high star formation state (SFR = 20 ${M}_{\odot }$ yr⁻¹) than the low (SFR = 5 ${M}_{\odot }$ yr⁻¹), indicating a roughly constant mass-loading factor, as indeed is seen in Figure 13.

Zoom In Zoom Out Reset image size

We can estimate the cooling time of the gas as

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{{k}_{{\rm{B}}}T}{n{\rm{\Lambda }}(n,T)},\end{eqnarray} \tag{ 29 }$$

where Λ(n, T) is the density and temperature-dependent cooling function measured in erg s⁻¹ cm³. We show these cooling times plotted over the phase diagram in Figure 14. While this demonstrates that the hot gas in the wind has very long cooling times, and the temperature floor of T = 10⁴ K means that the cool gas in our simulations has an effectively infinite cooling time, there exists significant mass in gas at intermediate temperatures (2 × 10⁴ K < T < 5 × 10⁵ K) that is in an interesting regime. Given an approximate number density of n = 10⁻¹ cm⁻³ and a temperature at the peak of the cooling curve of T ≈ 10⁵ K, the approximate cooling time (assuming solar metallicity) is t_cool ∼ 10 kyr. For gas traveling at the hot gas speed of v ≥ 1000 km s⁻¹, this is comparable to the dynamical time of the wind,

$$\begin{eqnarray}&&{t}_{\mathrm{adv}}=10\,\mathrm{kyr}\left[\displaystyle \frac{r}{10\,\mathrm{kpc}}\right]{\left[\displaystyle \frac{v}{1000\mathrm{km}{{\rm{s}}}^{-1}}\right]}^{-1}.\end{eqnarray} \tag{ 30 }$$

We therefore conclude that a significant fraction of the intermediate-temperature gas (which tends to be moving at v < 1000 km s⁻¹; see Figure 15) can cool to T ∼ 10⁴ K while in the outflow. In particular, the gas that is most prone to cooling is the highest-density, lowest-temperature portion of the intermediate-temperature gas that is created by mixing and shocks. The flat nature of the ${\dot{M}}_{\mathrm{int}}$ flux seen in Figures 8 and 9 is an indication that this mixing is taking place at all radii, and the mixing rate is roughly balanced by the cooling rate. We expect this to hold as long as there is T ∼ 10⁴ K gas present at that radius.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the relationship between density and temperature, we can also investigate how velocity and temperature correlate in the outflowing gas. In Figure 15, we show a mass-weighted velocity temperature histogram. Again, the two peaks highlight that the bulk of the gas is either in the hot, fast-moving wind (v > 1000 km s⁻¹, T > 10⁶ K) or a cooler, slower phase (v < 800 km s⁻¹). By mass, the bulk of the cool gas has velocities less than 500 km s⁻¹, but there is a substantial tail with velocities all the way up to 1000 km s⁻¹, as was also demonstrated in the velocity profiles shown in Figure 7. This stands in contrast to our previous work using idealized cloud–wind simulations (Schneider & Robertson 2017), where only high-temperature gas attained high velocities. However, we note that these simulations are lower resolution than that work, and hence there is more numerical diffusion. We address this discrepancy in Section 5.

In Figure 16, we show a mass-weighted density velocity histogram, again at 35 Myr. In addition to the correlation between velocity and radius demonstrated in Figure 7, this figure shows that there is also a correlation between density and velocity. In other words, lower-density gas appears to be more effectively accelerated, which is certainly consistent with intuition. The sharp cutoff for high-velocity (v > 1000 km s⁻¹) gas at n ∼ 10⁻² cm⁻³ mirrors the cutoff in density for cool gas in Figure 14, possibly indicating that this is the maximum velocity that can be achieved by cool gas in a starburst-driven outflow. We investigate the physical mechanism by which this cool gas is accelerated in the following subsection.

Zoom In Zoom Out Reset image size

A goal of this work is to assess the potential contributions of both mixing and ram pressure to acceleration of the cool gas, as those are the two primary theories for cool gas acceleration via hydrodynamic processes. We have excluded other potential sources of acceleration in this work (e.g., magnetic drag, radiation pressure, and cosmic-ray pressure) in order to focus on the two processes most commonly invoked in the M82-like outflows we are studying here.

First, we consider the velocities we would expect for the cool gas if all of its momentum were gained via mixing with the hot phase of the wind. In this case, we would expect the momentum density of the cool gas, ρ_cv_c, to reflect the fraction of hot gas that has been mixed into the cool medium (and then cooled). Since the cool gas initially has zero "scalar mass," we can use the scalar variable to determine the fraction of material in any cool cell that was previously hot. From Equation (27), under the assumption of mixing-driven acceleration, we obtain a linear relationship between the cool gas velocity and its scalar value:

$$\begin{eqnarray}&&{v}_{{\rm{c}}}=\displaystyle \frac{{v}_{{\rm{h}}}}{{s}_{{\rm{h}}}}{s}_{{\rm{c}}}.\end{eqnarray} \tag{ 31 }$$

This relation is normalized by the velocity of the hot wind and its scalar value (this s_h does not have to be 1, because it reflects mixing of cool material into the hot gas at small radii).

To test the validity of this hypothesis, we plot in Figure 17 a mass-weighted histogram showing the radial velocity of the cool gas versus its scalar value for all gas with T < 2 × 10⁴ within the same cone used for the earlier profiles and fluxes (see also Melso et al. 2019). We also plot the linear relationship expected given Equation (31). Because v_h and s_h are radially dependent quantities, we use their average values outside of 1 kpc to normalize the slope of the line, which yields a slope v_h/s_ıh ≈ 2900 km s⁻¹.

Zoom In Zoom Out Reset image size

We see from Figure 17 that the slope of the line is a reasonably good fit to the locus of the cool gas, particularly at larger values of s_c, although the whole line is displaced upward slightly. The fact that very little gas sits above the line demonstrates that ram pressure cannot have a dominant effect on the acceleration of the cool gas. Ram pressure without mixing would transfer momentum from the hot gas to the cool by accelerating it without raising its scalar value, which would move gas straight up on the plot. On the other hand, several effects could push gas below the line, including gravity, nonradial acceleration, and a lower value of the normalization v_h/s_h. Gravity is important only for the larger column density clouds near the disk that are moving slowly; these also have the lowest value of s, and this is the part of the locus that lies the furthest from the line. Nonradial acceleration is expected given that the clusters are not distributed in a spherically symmetric way, though they are concentrated at the center of the disk. Within the distribution radius, R < 1 kpc, mixing can accelerate gas in all directions, leading to an initial shift below the linear relation in Figure 17 as scalar mass gets mixed in without producing radial acceleration. Lastly, the slope v_h/s_h is not necessarily constant; indeed, we see from Figure 6 that it should increase with radius. Because the histogram contains gas at all radii, we cannot directly represent this in the line plotted in Figure 17. That said, we see from Figure 7 that s_c is also increasing slightly with radius. In this case, we would expect the slope of the line to be shallower at low values of s_c and to steepen at higher values. This is consistent with the nature of the curve traced by the locus of the cool gas in Figure 17, which does appear flattest at small s and increases in slope at larger s.

Given that these results depend on the hydrodynamic mixing of gas with a large range of temperatures and densities, the question of whether they would hold at even higher resolutions is an important one. In order to assess convergence, we have run the same simulation at two additional resolutions, with Δx ≈ 10 and 20 pc. Because the input radius of the clusters themselves is R = 30 pc in this simulation, we expect that they are marginally resolved even in the lowest-resolution case. Any resolution dependence in the results is therefore likely to be a result of not capturing the mixing between the hot and cool gas in the disk and wind, as opposed to an unresolved feedback model.

We investigate the issue of convergence by comparing the mass, momentum, and energy outflow rates as a function of radius for each of the simulations. Figure 18 shows the resulting outflow rates for each simulation averaged between 33 and 37 Myr. In addition to averaging in time, we have additionally smoothed the radial fluxes using a third-order Savitzky–Golay filter before time-averaging them to damp some of the small-scale spatial variability and ease comparisons.

Zoom In Zoom Out Reset image size

A primary result demonstrated in Figure 18 is that the mass, momentum, and energy outflow rates vary by less than 50% for all phases at all radii between the three simulations. The biggest differences are in the cool gas, and the smallest are in the hot gas. The mass outflow rates do not show an obvious trend with resolution, though it is the case that the highest-resolution simulation has the highest cool mass outflow rate. However, there are some potential trends that become apparent in the momentum and energy panels. At intermediate radii, the momentum and energy outflow rates in the cool and intermediate phases are higher in the Δx = 5 pc simulation than at lower resolution. This indicates that momentum transfer from the hot to the cool gas is more efficient in the higher-resolution simulation, as would be expected if that transfer is primarily a result of mixing. Such a relationship would not be expected if the momentum transfer was primarily a result of ram pressure, since the lower-resolution simulations tend to have larger clouds with more surface area perpendicular to the wind. At higher resolution, the cool clouds break up more, leading to less (perpendicular) surface area relative to total column density and hence less efficient acceleration (see Schneider & Robertson 2017). The fact that the mass outflow rates for all phases are approximately equal at larger radii implies that although the details of the mixing process depend on the spatial resolution, the overall trend toward mass transfer from cool to hot gas as a function of radius holds.

We can see this relationship more clearly in Figure 19, where we again normalize the momentum and energy outflow rates by the corresponding mass outflow rate. Here the trend with resolution is obvious in both panels. As the resolution increases, the momentum in the hot phase decreases, while the momentum in the cooler phase increases. Similarly, the specific energy decreases faster in the hot phase at high resolution and increases faster in the cool phase, indicating that the momentum transfer is more efficient at higher resolution, because mixing is more efficient.

Zoom In Zoom Out Reset image size

In Section 4.6 we demonstrate using the scalar mass variable that acceleration of cool gas can be accounted for by mixing in some of the fast hot gas. Those results also imply an expected relationship between the metallicity of the gas and its velocity. In Figure 20, we show a velocity slice normalized by the associated scalar mass from the simulation at 35 Myr. While the correspondence between scalar mass and velocity is clear from Figure 4, Figure 20 shows this directly; when normalized by the scalar, most of the large-scale velocity structure in the outflow is wiped out.

Zoom In Zoom Out Reset image size

Because the injected gas in our model includes the directly deposited ejecta from supernovae, we expect the high scalar gas to be highly enriched in metals, while the disk gas has lower metallicity. Given a plausible metallicity for the disk and injected gas, it may thus be possible to back out a metallicity estimate for the cool gas observed in outflows based on its velocity, as follows. In the case that all acceleration is due to mixing with hot gas that originally had velocity v_h and scalar value s_h, it can be shown that

$$\begin{eqnarray}&&Z(v)=\displaystyle \frac{{Z}_{\mathrm{ej}}{s}_{{\rm{h}}}v+{Z}_{\mathrm{ISM}}{v}_{{\rm{h}}}}{{s}_{{\rm{h}}}v+{v}_{{\rm{h}}}},\end{eqnarray} \tag{ 32 }$$

where Z_ej = M_met,ej/M_ej is the metallicity of the supernova ejecta and Z_ISM is the original metallicity of the ISM. So, assuming that Z_ej and Z_ISM are known, we can use the roughly constant value of s_h = 0.4 from the simulation to estimate the metallicity of the hot gas:

$$\begin{eqnarray}&&{Z}_{{\rm{h}}}=\displaystyle \frac{{Z}_{\mathrm{ej}}{s}_{{\rm{h}}}+{Z}_{\mathrm{ISM}}}{{s}_{{\rm{h}}}+1}.\end{eqnarray} \tag{ 33 }$$

For the cool gas, one could use Equation (32) to obtain an estimate for its metallicity based on its observed velocity v and the observed or estimated values of v_h and s_h. Alternatively, in the case that Z_h is known from observations, Equation (33) can be used to estimate s_h in terms of Z_h, which can then be used in Equation (32) to estimate the metallicity of the cool gas.

Given the sophisticated treatment of the ISM in many current stellar feedback simulations (e.g., Gatto et al. 2017; Colling et al. 2018; Kannan et al. 2020; Kim & Ostriker 2018), the simplicity of the feedback model employed in this work may present a source of uncertainty. As mentioned in Section 2, we do not attempt to resolve low-temperature (T < 10⁴ K) gas in these simulations due to the additional computational expense, so we cannot model a full three-phase ISM or self-consistent star formation. Nevertheless, we expect our results for gas in the outflow to be reasonable, given that higher-resolution simulations that do take this lower-temperature gas into account find that the majority of outflowing gas is in the cool or hot phase (Kim & Ostriker 2018). It then remains to determine whether our simplified model of cluster feedback is sufficient to produce a realistic outflow.

As described in Section 2, the model employed in this work assumed that a small (<20) number of very massive (M_cl = 10⁷ ${M}_{\odot }$) clusters were responsible for driving the wind. Our rationale was that in systems like M82, a few of the biggest clusters tend to dominate the energy input in the wind. Additionally, only massive clusters are expected to break out from the disk (Mac Low & McCray 1988; Kim et al. 2017; Fielding et al. 2018), so while modeling smaller mass clusters is critical for accurately calculating momentum input in the ISM, we expect them to be less important in determining the physical characteristics of the outflow. Keeping the clusters centralized in the disk also allowed us to analyze these results in a spherical framework, considerably simplifying the analytic description. However, real star-forming galaxies have clusters with a range of masses, and they are often distributed more widely throughout the disk.

To test the extent to which our results depend on the simplifying assumptions made in our feedback model, we have run an additional CGOLS simulation with a more realistic cluster mass function ranging from M_cl = 10⁴ to  5 × 10⁶ ${M}_{\odot }$. We additionally allowed these clusters to be distributed throughout the disk, assigning their radial locations such that the integrated star formation rate surface density profile matched the gas density profile. Each cluster was turned "on" for a more realistic 40 Myr, and the mass and energy injection rates were set according to a Starburst99 single-burst population synthesis model. In a sense, this distributed model can be thought of as similar to the clustered model presented here but with a lower net star formation rate, given that a smaller fraction of the "star formation" is in clusters large enough to break out of the disk and contribute significant energy to driving the outflow.

While a full discussion of the results of that simulation is outside the scope of this paper, we note that the results are qualitatively very similar to those reported here. Both simulations produce multiphase outflows with similar profiles in mass, momentum, and energy outflow rate as a function of phase. The velocities in the hot gas are slightly lower (v ∼ 1000 km s⁻¹) in the distributed cluster simulation, which is expected given the smaller maximum cluster masses, and the net energy outflow rate, while still dominated by the hot gas, is lower by about a factor of 3, implying a lower α_eff. Nevertheless, cool gas is accelerated to similar velocities as those found here, and we observe a similar linear relationship between velocity and scalar mass in the cool gas. Thus, we conclude that the results presented in this work are reasonable, given that quantitative adjustments to our reported quantities would be expected for different star formation rates regardless.

Although the CGOLS simulations are unique in their resolution over the volume captured here, several recent numerical studies have addressed similar questions, and a comparison to their conclusions is warranted. These include both global wind simulations (Fielding et al. 2017; Vijayan et al. 2018) and higher-resolution simulations that focus on patches of the ISM (Kim et al. 2017; Li et al. 2017; Fielding et al. 2018; Kim & Ostriker 2018; Vasiliev et al. 2019; Vijayan et al. 2020). We are excluding from this discussion zoom-in simulations and cosmological models, as their subgrid prescriptions for supernova feedback and winds generally do not allow for a hot phase to form, making a comparison less relevant.

Generally speaking, the global disk simulations that have been carried out to date do not have self-consistent star formation and feedback (nor do the CGOLS simulations). Both Fielding et al. (2017) and Vijayan et al. (2018) prescribed a smooth galactic disk, similar to our work here, and then injected discrete supernovae (Fielding et al. 2017) or OB associations (Vijayan et al. 2018) in order to generate winds. Because both methods of feedback lacked the large, powerful clusters employed here, they both found less energetic winds. Fielding et al. (2017) found very low-energy winds with α ∼ 0.01, which they attribute to their single supernova injection model. This is consistent with our results, in which clustering reduces the energy losses and increases α. Both Fielding et al. (2018) and Vijayan et al. (2018) found values of the mass-loading factor, β, that are consistent with ours, with β ranging from a few percent to approximately 50% of the star formation rate and never exceeding it. While Fielding et al. (2017) did not have high enough resolution outside the disk to capture the multiphase nature of the winds, Vijayan et al. (2018) showed that their winds had a primarily two-phase structure, similar to those presented in this work. They additionally noted some acceleration of the cool gas in the wind to speeds of ∼300 km s⁻² but were unable to track the cool gas further due to the constraints of their box size and shorter run time of the simulation. The cool gas acceleration seen is consistent with that presented here. While Fielding et al. (2017) found maximum velocities of a few hundred km s⁻¹, Vijayan et al. (2018) found maximum velocities of 800 km s⁻¹ for the hot phase, in contrast to the maximum velocities presented here, which exceed 1500 km s⁻¹. This is also consistent with the increased energy of the winds resulting from our maximally clustered configuration versus the OB association clustering in Vijayan et al. (2018) and the discrete supernovae in Fielding et al. (2017).

We now turn our attention to an analysis of work that focuses on wind properties using simulations that do not capture a full galactic disk but rather simulate a smaller region (usually ∼1 kpc) at higher resolution. These studies generally have resolutions comparable to the CGOLS simulations, of order a few pc, and follow the winds out to distances of several kpc. Although most of the studies focus on Milky Way–like gas surface densities, Li et al. (2017) explored higher surface density environments more similar to our disk. In that work, the authors found low mass-loading rates (β ∼ 0.2) and moderate energy loading (α ∼ 0.3), which is consistent with our results. They also found that most of the energy stayed in the hot phase, leading to fast, hot outflows with velocities comparable to those presented here. However, little acceleration of the cooler phases was seen. The authors attributed to this to lower resolution outside the disk (they used adaptive mesh refinement to decrease the resolution above and below the midplane). This led to little mixing between the phases. While mostly focused on the physics of superbubble breakout, Fielding et al. (2018) also studied a higher surface density environment with explicitly added clusters that reach 10⁶ ${M}_{\odot }$. The results they found for their high surface density, high cluster mass simulations are very consistent with ours, with energy-loading factors between 10% and 50% that decrease slightly with radius and mass-loading factors of around 0.3 at their largest radii (0.5 kpc).

Although they focus on lower surface density environments representative of the solar neighborhood, we also compare our results to those of the TIGRESS model, presented in Kim & Ostriker (2018), and specifically the focus on wind properties in Vijayan et al. (2020), as these simulations represent the state of the art in physical realism, incorporating MHD, self-consistent star formation and stellar feedback, and shearing box periodic boundary conditions. The TIGRESS simulations were carried out in a 1 kpc × 1 kpc box that extended to 4 kpc on either side of the disk. In their analysis, the authors split the gas into the same temperature bins used here, allowing us to compare the results directly (though note that our "cool" gas is their "warm"). As in our simulations, the TIGRESS models show that the energy budget of the wind is dominated by the hot phase. At high velocities, the authors found that the acceleration of the cool gas is not consistent with a ballistic model, as the velocities are too high; hence, acceleration due to mixing with the hot phase must be at play. In our CGOLS simulations, the cool gas mass flux appears to increase out to a radius of ∼4 kpc, after which it decreases. Because this increase is present in all phases, we have attributed the increase primarily to our conical geometry, and as demonstrated in Figure 11, we are continuously capturing more of the total outflow out to approximately this radius. In TIGRESS, which features a planar geometry, the authors find a decreasing cool gas mass flux as a function of radius out to 4 kpc, which they attribute to a fountain flow. This is consistent with our findings at larger radii, where the total mass outflow rate decreases slightly. As in our simulations, this indicates that despite the transfer of momentum due to mixing, the total cool gas mass flux is not increasing with distance as a result of mixing with the hot phase.

At small radii, the TIGRESS simulation generally shows a larger loss of energy in the hot phase, with an effective energy loading of <0.1 by the time the hot gas reaches the scale height of 400 pc. By contrast, in our model, α_eff has decreased to 0.5 by 400 pc and only decreases by a small fraction at larger radii, indicating that our extreme clustering model leads to maximal hot gas energy loading. In the case where supernovae are less highly clustered, more remnants are expected to cool without breaking out of the disk, and more interaction between the hot bubbles and disk gas is expected in general, which appears to be the source of most of the energy loss in the wind. As a result of these increased energy losses, the TIGRESS models show an α of only a few percent at large radii, more than 10 times smaller than that predicted here. Between 1 and 4 kpc, however, the energy flux measured in the hot phase in TIGRESS decreases by less than a factor of 2, which is more consistent with our results at larger radii. Nevertheless (and despite a much lower gas surface density and star formation rate), they find a similar mass-loading factor for the hot medium, β_eff ∼ 0.1.

In recent years, a number of studies have focused particularly on the survival, acceleration, and mixing of cool gas moving through a hot background (see Section 1 for references). In particular, several studies, including Armillotta et al. (2016), Gritton et al. (2017), and Gronke & Oh (2018, 2020), predict the increasing mass of cool gas clouds as the cool gas mixes with the hot background and the resulting intermediate-temperature gas cools and condenses. While we observe evidence of substantial mixing in this CGOLS simulation, and the correspondence between the velocity of the cool gas and its scalar value indicates that the momentum is being transferred via this process, we do not see evidence of the total cool gas mass flux increasing with distance in the wind; rather, we find a decreasing cool gas outflow rate with a radius beyond a radius of 4 kpc. So, what might be the cause of this discrepancy?

First, we note that even with the unprecedented resolution we have employed for a global model, many of the cool gas clouds in our simulation are still unresolved relative to the criterion for mass growth outlined in Gronke & Oh (2018). There, the size required for a cool cloud to grow in mass is estimated as a function of the cloud size, overdensity relative to the background, background pressure, wind Mach number, and cloud temperature. Using the numbers relevant for cool clouds in our simulations as measured from the radial profiles in Figures 6 and 7, we find that at small radii (∼300 pc), clouds larger than roughly 0.1 pc are expected to grow in size and accrete material from the background wind. While we do not capture sizes this small, we do note several large clouds near the base of the outflow in, e.g., Figure 4, which should be resolved by at least eight cells per cloud radius, the resolution quoted in Gronke & Oh (2018) to capture cloud growth. At these radii, the cool gas mass outflow rate is indeed rising (see Figure 8). However, our simulation analysis at these radii is incomplete, since the spatial arrangement of our clusters means gas is still feeding into the cone at these radii. The fact that the total gas mass outflow rate continues to rise until R ∼ 3 kpc means that we cannot distinguish between cool cloud mass growth and increasing total mass outflow rates. However, at larger radii, R ∼ 6 kpc, $\dot{M}$ has plateaued, and this should not be a problem. Here we find that cool clouds with sizes larger than a few parsecs should be growing in mass according to the Gronke & Oh (2018) criterion. Again, we see clouds with radii larger than this that should be resolved in our simulations, but we find no evidence for increasing cool gas mass outflow rates at these radii; rather, we find the exact opposite (mass transfer from the cool phase to the hot).

Besides resolution, a primary difference between the clouds in our simulations and those in the simulations discussed above is the nature of the background wind. Whereas the winds are laminar and (in most cases) unchanging in the idealized simulations that see cloud mass growth, our winds are radially expanding and have substantial transverse velocity components that could contribute to disrupting cool clouds and preventing their growth. This is particularly true at small radii, as demonstrated in Figure 21, which shows Δv = v − v_r for a slice through the simulation at 35 Myr. As the figure demonstrates, the difference between the total and radial velocity can be as much as 1000 km s⁻¹, particularly at small radii. As a result, it is not clear whether our simulations actually contradict the results of these more idealized cloud–wind studies or the turbulent background conditions in winds imply that substantially larger cool cloud sizes are needed in order to see cool gas mass growth. In particular, it is entirely plausible that the process outlined in Gronke & Oh (2018) is happening in our winds and indeed is the reason for the cool gas acceleration, but the cool gas growth via mixing and cooling is balanced or overpowered by destructive processes owing to the turbulent structure of the background wind. In addition, Gronke & Oh (2020) also pointed out that the rate of cool gas mass growth is sensitive to the background pressure profile in the hot wind. While an adiabatically expanding hot wind leads to growth in their idealized simulations, a shallower pressure profile as measured in Figure 6 would tend to inhibit it. This reinforces the notion that global processes that are not captured in small box simulations may profoundly influence the survival of cool clouds in galactic winds.

Zoom In Zoom Out Reset image size