Galactic Winds and Fountains

Numerical Simulations of Multiphase Winds and Fountains from Star-Forming Galactic Disks: I. Solar Neighborhood Tigress Model

[ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA [ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA,

Gas blown away from galactic disks by supernova (SN) feedback plays a key role in galaxy evolution. We investigate outflows utilizing the solar neighborhood model of our high-resolution, local galactic disk simulation suite, TIGRESS. In our numerical implementation, star formation and SN feedback are self-consistently treated and well resolved in the multiphase, turbulent, magnetized interstellar medium. Bursts of star formation produce spatially and temporally correlated SNe that drive strong outflows, consisting of hot () winds and warm () fountains. The hot gas at distance from the midplane has mass and energy fluxes nearly constant with . The hot flow escapes our local Cartesian box barely affected by gravity, and is expected to accelerate up to terminal velocity of . The mean mass and energy loading factors of the hot wind are 0.1 and 0.02, respectively. For warm gas, the mean outward mass flux through is comparable to the mean star formation rate, but only a small fraction of this gas is at velocity . Thus, the warm outflows eventually fall back as inflows. The warm fountain flows are created by expanding hot superbubbles at ; at larger neither ram pressure acceleration nor cooling transfers significant momentum or energy flux from the hot wind to the warm outflow. The velocity distribution at launching near better represents warm outflows than a single mass loading factor, potentially enabling development of subgrid models for warm galactic winds in arbitrary large-scale galactic potentials.

galaxies: ISM – galaxies: star formation – magnetohydrodynamics (MHD) – methods: numerical

0000-0003-2896-3725]Chang-Goo Kim 0000-0002-0509-9113]Eve C. Ostriker

1 Introduction

Galactic scale gas outflows (or winds) are ubiquitous in star forming galaxies (see Veilleux et al., 2005; Heckman & Thompson, 2017, for reviews) and believed to be essential to distribution of the gas and metals in galaxies and the circumgalactic/intergalactic medium (CGM/IGM) and hence to regulating cosmic star formation history (see Somerville & Davé, 2015; Naab & Ostriker, 2017, for reviews). Theoretical models of the stellar mass-halo mass relation constructed by abundance matching of observational stellar mass functions to simulated halo mass functions (e.g., Moster et al., 2013; Behroozi et al., 2013; Rodríguez-Puebla et al., 2017) indicate that galaxies (or dark matter halos) are very inefficient at converting gas into stars. At low redshift, at most of the available baryonic mass has been converted into stars at halo mass of , while the ratio of stellar mass to halo mass declines steeply toward both lower and higher masses. Recent cosmological hydrodynamic simulations of large volumes of the Universe require strong outflows driven by both star formation and active galactic nuclei feedback to explain low baryonic abundance in galaxies compared to the cosmic fraction (e.g., Springel & Hernquist, 2003; Vogelsberger et al., 2014; Schaye et al., 2015).

Direct inclusion of feedback processes in large-volume cosmological galaxy formation simulations is still not feasible in practice. For star formation feedback by supernovae (SNe), implementation via simple thermal energy dumps suffers “overcooling,” with energy radiated away without preventing in situ star formation or driving winds (e.g., Katz, 1992); this is because resolving the cooling radii of SN remnants requires much higher resolution (Kim & Ostriker, 2015a) than is practicable in large-volume simulations. In cosmological zoom-in simulations, more careful treatments of SNe allowing for a “momentum prescription” at low resolution can solve at least some aspects of the overcooling problem (e.g., Kimm & Cen, 2014; Kimm et al., 2015; Hopkins et al., 2014, 2017), especially for dwarfs. However, given the constraints of computational expense, treating unresolved physics with parameterized models is unavoidable in many situations, including in simulations of galaxy groups/clusters, and in the large boxes needed for fully-sampled statistics of cosmic galaxy populations (e.g., Schaye et al., 2015; Pillepich et al., 2017).

When star formation feedback physics is not directly simulated, galactic winds are not an outcome but an input that is part of the “subgrid” treatment (Somerville & Davé, 2015). Currently, however, subgrid models of wind driving by stellar feedback often either adopt highly simplified scaling prescriptions for wind mass loss rates (relative to the star formation rate) and velocities (relative to the halo potential depth), or else are calibrated using empirical results from a limited set of galaxies (and hence are not fully predictive). Better theoretical models are clearly needed. Towards this end, the first step is to provide a physical understanding and detailed characterization of outflowing gas (including both winds and fountains) in galaxies, informed and calibrated based on high-resolution three-dimensional numerical magnetohydrodynamic (MHD) simulations. To fully capture the interaction between stellar feedback and a realistic multiphase interstellar medium (ISM), it is crucial to self-consistently include the gravitational collapse that produces star clusters and to resolve the local injection of energy from individual massive stars.

In classical theoretical models of galactic winds motivated by observed starburst galaxies (e.g. Chevalier & Clegg, 1985), a steady, adiabatic flow is produced by a central energy source. In this approach, hot, overpressured gas flows are characterized by “mass loading” and “energy loading” factors, respectively defined by the ratios of mass and energy outflow rates to star formation rates and energy injection rates at the wind base. Simple spherical wind models can also be constructed that allow for radiative cooling, such that the temperature precipitously drops at some radius in certain parameter regimes (e.g. Wang, 1995; Bustard et al., 2016; Thompson et al., 2016).

Observed galactic outflows are multiphase in nature. Systematic observations reveal prevalent multiphase structure of galactic winds with cold molecular (e.g., Weiß et al., 1999; Leroy et al., 2015), neutral (e.g., Heckman et al., 2000; Martin, 2005; Rupke et al., 2005; Chen et al., 2010; Contursi et al., 2013), ionized (e.g., Pettini et al., 2001; Shapley et al., 2003; Steidel et al., 2010; Erb et al., 2012; Heckman et al., 2015; Chisholm et al., 2017), and hot gas phases (e.g., Strickland & Stevens, 2000; Strickland & Heckman, 2007; López-Cobá et al., 2017). For the best studied example, local starburst M82, Leroy et al. (2015, see also ) have shown a clear signature of decreasing outward mass fluxes in molecular and neutral gas as a function of the distance from the disk midplane, implying a fountain (Shapiro & Field, 1976; Bregman, 1980) rather than a wind for the cooler gas. In one conceptual framework, warm and cold gas in outflows results when a hot medium accelerated by its own pressure gradient cools radiatively; an alternative concept is that overdense warm and cold ISM clouds are “entrained” by a high-velocity, low-density hot wind. More realistically, both effects can in principle occur, and in general there is a complex interaction between the multiple phases that are present.

The mass and energy loading factors are key quantities that characterize winds and quantify their significance in controlling baryonic mass cycles of galaxies. Measuring these loading factors has been of intense observational interest, but uncertainties are still large. In particular, the reported mass loading factor ranges widely from 0.01 to 10 (Veilleux et al., 2005). Depending on the assumed geometry, metallicity, and ionization state, the mass outflow rate can easily be reduced by a factor of 10 (e.g., Chisholm et al., 2016b, 2017). In addition, uncertain deprojection may result in an overestimate the velocity, incorrectly leading to interpretation of a low-temperature outflow as a wind rather than a fountain. If gas is not really escaping, the outward mass flux will be a decreasing function of distance from the wind launching region, and mass fluxes estimated at small radii would exceed the true losses from a galaxy.

Predicting wind loading factors theoretically requires modeling the interaction between SN remnants (including from clustered SNe) and the ISM. Expansion of superbubbles driven by multiple SNe has been studied by idealized analytic models and simulations (e.g., McCray & Kafatos, 1987; Mac Low & McCray, 1988; Mac Low et al., 1989; Mac Low & Ferrara, 1999; Kim et al., 2017). While these idealized models provide essential physical insight and quantitative estimates, firm theoretical measurements of mass loading in multiphase winds from galactic disks require ISM models with realistic spatio-temporal distribution of SNe and vertical stratification. A number of local stratified-disk simulations, with increasingly high resolution, have modeled the multiphase ISM with SN feedback (e.g., Korpi et al., 1999; de Avillez, 2000; de Avillez & Breitschwerdt, 2004; Joung & Mac Low, 2006; Hill et al., 2012; Gent et al., 2013; Walch et al., 2015; Girichidis et al., 2016b; Li et al., 2017), albeit with SN rates and locations imposed rather responding self-consistently to star formation. Some recent numerical work has focused specifically on the outflow properties driven by SN feedback (e.g., Creasey et al., 2013; Martizzi et al., 2016; Fielding et al., 2017), although with a cooling cutoff at that does not allow for the full range of ISM phases.

Very recently, it has become possible to evolve the turbulent, magnetized, multiphase ISM in local galactic disks with cooling and heating, galactic differential rotation, and self-gravity, including fully self-consistent resolved star formation and SN feedback over durations of several 100 Myr (Kim & Ostriker, 2017, Paper I hereafter). A few other recent simulations have also included self-gravity to model SN rates and positions self-consistently with star formation (Hennebelle & Iffrig, 2014; Gatto et al., 2017; Peters et al., 2017; Iffrig & Hennebelle, 2017), but given their relatively short simulation duration (), they have not reached a statistically quasi-steady state and wind properties may be subject to transient effects from the simulation start-up.

In this paper, we analyze our fiducial model from the TIGRESS (Three-phase ISM in Galaxies Resolving Evolution with Star formation and SN feedback) simulation suite introduced in Paper I, in order to provide more comprehensive understanding of multiphase gas outflows in the realistic ISM. Our analysis here mainly focuses on characterizing differences among outflows of different thermal phases. In a subsequent paper, we will analyze models with different galactic conditions to investigate systematic scaling relations of wind properties (e.g., Muratov et al., 2015; Heckman et al., 2015; Heckman & Borthakur, 2016; Chisholm et al., 2017).

In Section 2, we review equations for steady, adiabatic flows and summarize key physical quantities to be measured from the simulation. We then present an overview of gas flows in the simulation, including overall mass fluxes and vertical profiles of each gas phase; this demonstrates the necessity of a phase-by-phase analysis. In Section 3, we analyze the hot gas component, showing that it is consistent with a wind having well-defined mass flux and specific energy (or Bernoulli parameter) that are approximately constant as a function of distance from the midplane. Section 4 presents an analysis of the warm gas component, showing characteristics of a fountain flow that has decreasing mass and energy fluxes as a function of distance from the midplane. Section 5 provides mass and energy loading factors of each phase, comparing these to previous work and to observations. Section 6 summarizes our main conclusions.

2 Overview of Gas Flows in TIGRESS

2.1 Outflow terminology, vertical profiles, and classical adiabatic winds

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

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

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


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

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


consisting of specific kinetic energy, specific enthalpy , and gravitational potential terms. Note that here, for simplicity, we consider adiabatic, unmagnetized gas, but Appendix A presents the full equations for general MHD flows in a local shearing box, and shows that horizontal- and time-averaging yields a set of simple steady-state 1D flow equations, which can be applied to our simulations. In this paper, the reference point is at the midplane.

For flows in a local Cartesian box, like ours, with energy sources near the midplane, any time-averaged steady winds that may exist are launched vertically along the axis. Through Gauss’s Law, volume integration of Equations (1) and (2) gives mass and total energy fluxes through surface area at as




Here, is the distance from the disk midplane, and total mass and energy injection rates within are and . The quantities and stand for mass and total energy fluxes averaged over horizontal area at height .111 For the energy flux, we hereafter use subscripts KE, TE, GE, and ME to denote different energy components. These are respectively (kinetic), (thermal, with enthalpy), (gravitational), and (magnetic; see Section A and Equation (A10)). The subscript E will denote the total energy term, the sum of all four components. The above relations assume periodic boundary conditions in both of the horizontal directions; energy terms associated with the background shear from integrals over faces perpendicular to are discussed in Appendix A.

Since SN explosions are usually concentrated within a thin layer near the midplane, and are expected to approach constant values for , where H is the disk scale height. Hence, if is constant as in the local Cartesian coordinates, for a steady wind the (areal) mass and energy fluxes and would also approach constants at above the source region. For general geometries (e.g., Chevalier & Clegg 1985 for spherical coordinates), however, where the area that encloses a given set of streamlines varies as a function of distance (e.g., for the spherical case), the fluxes would also vary (e.g., ).

To the extent that horizontal correlations (or anticorrelations) among variables at a given may be neglected, , and if we may assume symmetry across the midplane this yields . This implies that beyond the source region where and have reached their final values, the Bernoulli parameter becomes constant for a steady (or time-averaged) flow. In the more general case with an increasing function of distance, because the area perpendicular to streamlines is the same for both mass and energy flux, and would vary but would still be conserved along streamlines (beyond the source region) irrespective of geometry. Therefore, for steady adiabatic winds (or equivalently for any time-averaged, adiabatic portion in a more general outflow), the Bernoulli parameter is a key quantity that enables robust extrapolation of flow evolution to large distances.222This holds true for MHD flows if the Poynting flux is negligible, as in our simulations; see Figure 5. As applied to the present problem, this suggests that evaluation of with our local disk simulations should provide predictions for properties of the hot portion of the wind at large distance (outside our simulation domain) where wind streamlines open up (becoming more radial than vertical). This motivates the need to quantify mass and energy fluxes in the launching region – just above the source region – in our simulations.

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

2.2 TIGRESS simulation model and analysis

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

Our simulations yield realistic, fully self-consistent three-phase ISM models with self-regulated star formation.333Of course, the absence of global geometry means that we are unable to follow effects of strong noncircular flows in the disk, transport of gas from one radius to another in a fountain (which would also require significant angular momentum exchange), or the transition of hot winds through a sonic point. Nevertheless, the high resolution afforded by our local scope is extremely valuable for limiting artificial mixing, which is essential for understanding key characteristics of multiphase flows. Paper I presented results from a fiducial model with parameters similar to those of the Solar neighborhood. There, we showed that after a quasi-steady state is reached. When stars form, massive stars enhance heating in the warm and cold ISM, and the SN rate increases, driving turbulence throughout the ISM. Both feedback processes disperse dense clouds and puff the gas disk up, temporarily shutting star formation off. With the corresponding reduction in star formation feedback, the gas can settle back to the midplane and collect into dense clouds which then form a new generation of stars. In Paper I, we evaluated several basic ISM and wind properties, and demonstrated their convergence as a function of the numerical spatial resolution.

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

Figure 1: Sample snapshots illustrating (left) outflow-dominated, and (right) inflow-dominated periods, at and , respectively. Gas density is shown in color scale volume rendering and the velocity field is shown with vectors; vector colors (rather than length) indicate velocity magnitudes. A fast-moving, low-density (dark blue) outflow is evident in the left snapshot, while moderate-density (green) gas that was previously blown out to large distances is falling back toward the midplane in the right. Velocity fields are turbulent near the midplane (), but ordered in either outflowing or inflowing directions at large distances. For visual clarity, only the upper half of the simulation box ( to ) is shown, with a full horizontal crossection .

Figure 2: Sample slices through showing temperature and vertical velocity at (left) and (right), representing outflow- and inflow-dominated periods. During outflow-dominated periods (left), the hot component fills most of the volume at high and flows outward at high velocity, while the warm component is confined in small cloudlets. During inflow-dominated periods (right), the warm component occupies most of the volume and falls toward the midplane. As in the volume renderings of Fig. 1, only the upper half of the simulation box ( to ) is shown.
Figure 3: Time evolution of (a) mass and (b) energy fluxes driven by bursts in the star formation rate. Mass and energy fluxes per unit area at distances 1, 2, and 3 kpc from the midplane are measured (Equations 4 and 5). Each burst in (shown in grey in upper panel) leads to outflows in both mass and energy near the midplane. The net mass flux decreases at larger because most of the outflow near the midplane is warm “fountain” gas with low velocity that turns around (resulting in periods of negative mass flux at ). Energy fluxes are substantial at all heights, due to the dominance of hot gas that escapes as a wind.

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

Although star formation and hence SN feedback are impulsive rather than continuous, the system approaches a quasi-equilibrium state. This state is a limit cycle mediated by the feedback loop, in which epochs of cooling and collapse alternate with epochs of heating and expansion.444We note, however, that a quasi-steady state is not guaranteed for all galactic conditions (e.g., Torrey et al., 2017), and may only hold in parts of the parameter space in which the vertical oscillation time (which controls collapse and star formation) is sufficiently long compared to the stellar evolution timescale (which controls feedback and expansion). Given that a quasi-steady state exists in the present simulation, horizontal- and temporal averages can be constructed to characterize this mean state. Since different thermal phases coexist at all heights, to understand the outflows and inflows of mass and energy it is further necessary to separately construct horizontal averages of each thermal phase.

Figure 4: Vertical distributions of mass, vertical momentum, and energy densities, averaged horizontally and over upper () and lower () sides of the disk and over time . Colored lines for separate thermal components show profiles of (a) hydrogen number density , (b) vertical momentum density (which is the same as the mass flux), and (c) total energy density (excluding gravity) , each as a function of distance from the midplane. Color-coded shaded regions represent one-sigma temporal fluctuations. In order to properly visualize both the magnitude and sign of the momentum density, we use a linear scale for and a log scale for . The warm and hot phases respectively dominate mass and energy densities above the disk scale height (indicated by vertical dotted line), and the hot component also has the largest vertical momentum density (net outward mass flux).

In Figure 4 we plot horizontally- and temporally-averaged profiles of mass, momentum, and energy distributions for thermally separated gas phases; these profiles average over and also average over upper () and lower () sides. Profiles show hydrogen number density, , outward vertical momentum density, , and total energy density (excluding gravity)


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

The four thermal phases plotted in Figure 4 are cold (), warm (), ionized (), and hot ().555 Note that we omit the “unstable” phase () defined in Paper I and merge it into the “cold” phase since (1) these phases are not of primary interest in this paper since we are focusing on gas above the disk scale height, and (2) the sum of these two phases numerically converges better than the individual phases. Above the warm/cold layer (, where ), the mass density is dominated by the warm component and the energy density is dominated by the hot component. As the individual terms in the energy density are proportional to corresponding terms in the momentum flux (, ), the hot medium also dominates the momentum flux away from the midplane. The hot medium is the largest contributor to the time-averaged vertical momentum density, which is the same as the time-averaged net mass flux. Although the mean value of vertical momentum density (or net mass flux) of the warm medium is effectively zero, there are large temporal fluctuations (indicated by the green shaded region) at small and intermediate because the warm gas contributes significantly to both outgoing and incoming mass fluxes at different times (as in Figure 3).

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

3 Hot Winds

In this section, we focus on the hot component, defined by , representing gas that has been shock heated by SN blastwaves. The hot medium fills most of the volume above the disk scale height, and cooling in tenuous hot gas is inefficient. With pressure gradients that accelerate it outward (cf. Figure 4(c)) and a source from SN shocks propagating into the surrounding warm and cold medium near the midplane, the horizontally and temporally averaged hot medium naturally fits the criteria for a (quasi-)steady, adiabatic wind. We therefore might expect that the mass flux and energy flux (and thus the Bernoulli parameter) of the time-averaged hot component would be (approximately) conserved as the gas flows outward.

Figure 5 explicitly shows vertical profiles of (a) the mass flux, and (b) the specific energy, , based on time averages of each horizontally averaged flux. To distinguish between outflows and inflows, for the mass flux in addition to net flux we separately show the flux of outflowing () and inflowing () gas.

Figure 5: Time averaged () vertical profiles of hot gas (a) mass fluxes and (b) specific energies. (a) The net mass flux is shown as a blue line with one-sigma temporal fluctuations as the blue shaded region. The hot gas does not show any significant inflows (yellow line), implying that any decrement of the mass flux is due to a phase transition to cooler phases and not direct inflows of the hot gas. (b) Total Bernoulli parameter as well as individual components (i.e. kinetic, thermal, gravitational, and magnetic terms; see Equations (3) and (A11)). Mean values are shown as colored lines, and one-sigma temporal fluctuations are shown as shaded regions. The specific enthalpy and kinetic energy are much larger than the gravitational potential. The magnetic term plays a minor role. In both (a) and (b), the vertical dotted line indicates the disk scale height ().
Figure 6: Mass flux and specific energy profiles for the hot wind as in Figure 5, except at a single instant .

The mass flux of the hot gas shows net outflow (blue line in Figure 5(a)), with negligible inflow flux at all heights (yellow line in Figure 5(a)). As the SNe that create the hot medium are concentrated near the midplane, the hot gas mass flux within the warm/cold layer (; below the vertical dotted line in Figure 5) first increases with and then decreases as hot SNRs have a maximum size before the onset of cooling. Above the point where the hot gas in the interior of SN remnants (or superbubbles) breaks out of the warm/cold layer the hot gas mass flux is nearly constant with . The slight decrease of the mass flux at large is caused by non-zero cooling (radiative and mixing of hot gas with cooler phases).

Figure 5(b) shows near constancy of the time-averaged total specific energy of the hot gas above , as expected for a steady, adiabatic wind (given the negligible Poynting flux). We also calculate individual components of the (averaged) energy flux, and divide them by the (averaged) mass flux. This provides the components of the mass-flux-weighted specific energy , consisting of kinetic, thermal (enthalpy), gravitational, and magnetic terms (see Equations 3 and A10). Once the hot gas breaks out of the warm/cold layer, the gas flow approximately preserves mass and energy fluxes because except for limited cooling there are no sources (or sinks) for the hot gas mass and energy. We have checked that the individual cooling and heating terms, including Reynolds and Maxwell stresses arise from the shearing box (see Equation (A7)), are indeed small compared to the SN energy injection.

From , the enthalpy of the hot gas implies a temperature in the range (or sound speed in the range ), and the kinetic energy of the hot gas translates to a velocity of . The Poynting flux contribution to is negligible, corresponding to an Alfvén speed of .

In SN-driven hot winds, the enthalpy (specific heat) term dominates over other components of the Bernoulli parameter, including the gravitational potential, at heights less than the disk radius. At large distance where streamlines open up in angle, one could expect hot galactic winds with constant to accelerate past a sonic point to an asymptotic velocity of as they adiabatically cool, similar to classical Parker stellar wind solutions.666Note that for galactic potentials, is generally computed relative to the midplane, whereas for point mass potentials the reference point is at infinite distance; the zero point in is irrelevant as long as it is consistent in and . Global simulations indeed show the expected behavior at large distance (e.g. Fielding et al., 2017).

In general, hot winds are accelerated by pressure gradients at the same time as enthalpy is converted to kinetic energy, and a sonic transition in a steady wind is only possible if the crossectional area increases as the flow moves outward (e.g. Shu, 1992). In the Cartesian geometry of the present simulations, streamlines cannot open up and there is no associated adiabatic cooling, limiting the pressure and density gradients and therefore the acceleration of the flow. However, the constancy of both the Bernoulli parameter and mass flux with for the hot medium in our simulation suggests that it properly represents the near-disk regions for a generalized galactic disk wind, in which streamlines emerge from the disk vertically (with when ) and would open up ( increasing with distance) when . The mass and energy fluxes carried by the hot wind are controlled by the interaction between SN shocks and the warm-cold medium that creates the hot ISM well within the disk, in processes that are unlikely to be affected by large-scale global galactic and CGM properties and geometry. Therefore, the Bernoulli parameter we calculate should be a robust estimator of asymptotic hot wind speed irrespective of the constraints of our Cartesian box.

As shown in Figure 3 (and demonstrated by the shaded region in Figure 5(b)), star formation and hence outflows are bursty, resulting in large temporal fluctuations. Based on analysis of for the hot component within short time ranges, the maximum asymptotic velocity777Here, this maximum wind velocity is calculated using at a height , although in reality the wind velocity would continue to decrease slowly with distance due to the logarithmic increase of with at large distances in dark matter halo potentials. of the hot wind could reach up to , while the mean asymptotic wind speed would be . Defining the speed required to escape to as , in our simulation domain , so the hot wind easily escapes, i.e. even at large for the hot component. A hot wind launched with the local conditions of our simulation would also be able to propagate far into the halo for the Milky Way, where the escape velocities are and using MWPotential2014 in galpy (Bovy, 2015). More generally, the far-field velocity for a hot wind with given local launching conditions can be estimated based on and the large-scale galactic potential.

Finally, we note that while we have mostly based the discussion above on temporal averages of horizontally-averaged profiles, the properties of instantaneous profiles are quite similar. To illustrate this, Figure 6 shows the mass flux and specific energy of the hot component at Myr, the “outflow” snapshot shown in Fig. (2). Except for local fluctuations, the instantaneous mass flux and specific energy profiles are overall very similar to the corresponding time-averaged profiles. In particular, the specific energy profiles in Fig. 5 and 6 are quantitatively almost the same, while the mass flux is higher in Fig. 6 than in Fig. 5 because the latter includes non-outburst epochs as well as outburst epochs in the temporal averaging.

4 Warm Fountains

Over the duration of the simulation, the mean net mass flux in warm gas out of the simulation domain is , about 28 percent of the net mass flux in hot gas. As we shall show, the outward mass flux of the warm medium secularly decreases with , such that if our box were taller we would expect the mean net mass flux out of the simulation domain to be even smaller. With negligible time-averaged net outward mass flux at large distance, the warm medium at in our simulation is not a true galactic wind. Even though the inflows and outflows of warm gas over long timescales are essentially balanced, time variations in the warm gas flux are quite large and include both positive and negative values (green shaded region of Figure 4(b)). The fluctuating behavior of the warm medium can be contrasted with the much smaller temporal fluctuations of the mass flux in the hot gas (red shaded region of Figure 4(b)). Evidently, the warm gas does not produce an escaping wind like that in the hot gas but a fluctuating fountain that at any time consists of both outflow and inflow. Figure 3 shows that alternating inflow and outflow dominance in the warm gas is reflected in the alternating signs of the mass flux for the whole medium at . In addition, the decrease in the magnitude of the total (phase-integrated) mass flux with increasing reflects the secular decrease in the net mass flux of the warm fountain (difference between inward and outward fluxes) with distance.

Figure 7: Vertical profiles of warm gas mass fluxes (top) and specific energies (bottom) averaged over outflow (left; ) and inflow (right; ) periods. Inflows and outflows are comparable with each other, implying little net mass and energy outflows (or inflows) from the simulation domain associated with the warm medium. Note that, for visualization purpose, we take absolute values of the inflow mass flux and net flux during the inflow period; these have negative signs by definition. For the warm medium, the kinetic term exceeds both the thermal and magnetic terms in the specific energy, but is always lower than the gravitational term above . This explains why the warm medium creates a fountain rather than a wind. The vertical dotted lines in all panels indicates the gas scale height ().

To quantify the characteristics of warm fountain flows, we take time averages selectively for outflow () and inflow () periods. Outflows of warm gas occur when many correlated SNe from a star formation burst lead to a superbubble expanding into the warm and cold layer, while inflows occur when the disk is in a quiescent state with reduced star formation after the cold medium has been dispersed by a previous burst. Figure 7 plots time averaged mass flux (top) and specific energy (bottom) for outflow (left) and inflow (right) periods. Although one flow dominates the other during each period, the opposite flows always exist at all heights and are more significant compared to the case for hot gas (Figure 5). For the warm medium, the kinetic term in the specific energy exceeds the enthalpy, but remains below the gravitational potential term at . This explains why most of the warm gas outflow turns around at and falls back toward the midplane. During the outflow period, the mean velocity of the warm gas is in the range for .

Overall, the warm medium occupies more volume above during the inflow period than the outflow period. This is because the hot gas is mainly generated during the outflow period by shocking the warm gas, and the high-pressure hot gas confines the warm gas into small cloudlets. When the disk becomes quiescent, the warm gas expands into previous hot wind channels (see Figure 2). About a factor of 5 to 10 more volume is occupied by the warm gas in the inflow period than the outflow period for and slabs.

Figure 8: Outward vertical velocity () probability distributions of the warm medium. Distributions of volume (top), mass (middle), and mass flux (bottom) as a function of velocity are calculated for thickness slabs above and below the midplane () centered at . In the bottom row, the mass flux within each velocity bin is in units of , and negative mass fluxes are shown with dotted lines. Warm gas with both signs of is present during both “outflow” and “inflow” periods, but the corresponding mean velocity (for ) changes sign. Acceleration of the warm gas to high velocities, especially during outflow periods, is evident in the difference of profile shapes between the midplane () and higher latitude, with an exponential tail at high velocity developing by . However, an overall deceleration with height in the increasing gravitational potential is also evident in comparison of profiles at increasing . The magenta dashed lines in the left column are fits to the fast-moving gas at given by Equations (7) and (8).

Although in the current simulations the warm gas is almost entirely confined by the galactic gravitational potential, this would not necessarily be true if the potential were shallower, as in dwarf galaxies. Direct simulations for different galactic conditions, including shallower potentials, are underway using the same TIGRESS framework. However, we can also use our current simulation to provide information on what might be expected by quantifying the fraction of fast-moving warm gas. Figure 8 shows outward velocity probability distribution functions (PDFs) weighted by volume (top), mass (middle), and mass flux (bottom) within slabs of thickness during outflow (left) and inflow (right) periods. We show results at several distances averaged over both sides of the disk (). Note that the volume-weighted PDF is normalized by slab volume (), and the mass-weighted PDF is normalized by the mean mass within the same volume () ( is the mass in the whole domain at a given time), while the mass flux PDF is in physical units of . During both outflow and inflow periods, Figure 8 shows that the warm gas velocity has a broad distribution with both outward and inward velocities.

During the outflow period, Figures 8(a) and (c) show that the volume and mass of high velocity () warm gas increases from the midplane to , where the specific kinetic energy is larger than the gravitational potential (see Figure 7(c)). The increase in mass of high-velocity, high-altitude warm gas between the midplane and is due to acceleration of the warm medium pushed by expanding superbubbles; this includes warm gas that was shock-heated to the hot phase (and accelerated to high velocity) and subsequently cooled back down. Figures 8(c) and (e) show that as increases the peaks of the distributions of mass and mass flux move to higher velocity and the overall outflowing gas fraction decreases. This general trend represents dropout/turnaround of warm gas fluid elements with low (and decelerating) velocities that are unable to climb to large in the gravitational potential.

In principle, acceleration of warm clouds driven by hot-gas ram pressure, cooling of fast hot gas, and dropout of low-velocity warm fluid elements could all contribute to the gradual increase of warm-medium specific kinetic energy at shown in Figure 7(c). Figure 8(c) shows, however, that overall the mass of high-velocity warm gas is decreasing with increasing . This has the important implications that in our simulation (1) warm clouds at large are not significantly accelerated by ram pressure of the hot, high-velocity gas that is flowing out around them, and (2) relatively little hot gas is converted to the warm phase through cooling at large . Rather, the warm medium is primarily accelerated via direct energy input from SNe at , and at higher altitudes warm fluid elements slow and turn around according to the competition between the gravitational potential and the kinetic energy they initially acquired at small . Figure 4(c) and Figure 5(b) are also telling in this regard: the total energy density (and individual components) of the hot medium declines very slowly for , while the energy density of the warm medium declines steeply; since momentum flux terms are proportional to energy density terms, this indicates that there is no significant transfer of momentum from the hot to the warm gas.

During the inflow period, the majority of the warm gas is falling. Since SN feedback is never completely turned off, however, some warm gas is still accelerated outward. The fraction of (outgoing) fast-moving warm gas () is reduced by a factor of 5 to 10 in inflow compared to outflow periods. Combined with the total outflowing mass fluxes of and at during outflow and inflow periods (see Figure 7(a,b)), respectively, the mass fluxes of the fast-moving warm gas are about and (see Figure 8(e) and (f)). This can be compared to a mean mass flux of hot gas at the same height of and during outflow and inflow periods, respectively. Although the “fast” warm outflow has comparable mass flux at to that of the hot medium during outflow periods, even at the warm medium yields very little mass escaping from the large-scale potential in the present simulation, while the hot medium mostly escapes ( and of hot and warm gas have respectively escaped during the time interval of ). This emphasizes the importance of measuring not just mass fluxes, but mass fluxes and specific energies in comparison to galactic escape speeds.

At large distances the large-scale gravitational potential strongly affects the warm-gas velocity distribution by enforcing dropout of lower-velocity material, but closer to the midplane this is less of an issue. At around the disk scale height, the gravitational potential is small compared to the specific kinetic energy of the gas, and global geometric effects are not important yet. We therefore consider the velocity distribution at as representative of the launching conditions for a warm outflow, which would apply relatively independently of the global galaxy (e.g. in a dwarf as well as a large galaxy for given local conditions). Here, we find the PDFs during outflow periods of the fast-moving warm gas () at are well fitted by a single exponential function,


for or , where the normalization factors for volume and mass weighted PDFs are and , respectively, and the characteristic velocities are and . Using the mass PDFs, the mass-flux PDF is given by


where the mean density, , is given by , corresponding to hydrogen number density .

By quantifying the mass flux PDF in warm gas for models with different local conditions, it should be possible to develop a comprehensive quantitative characterization of warm wind launching by star formation feedback. These local results could then be used to make global predictions. For example, integration of Equation (8) for velocities for would yield a mass flux . If this holds in general, it means that measurement of the launching-region warm mass flux for given local disk conditions could be used to predict the flux that actually escapes into the halo for a galaxy with arbitrary halo velocity . This statistical characterization can be combined with measurements of mass flux and Bernoulli parameter for the hot medium to develop subgrid models of multiphase wind driving for implementation in galaxy formation simulations.

5 Mass and Energy Loading

5.1 Simulation Results

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


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


where is the total energy per SN, and is the total mass of new stars per SN (see Paper I). To compute , , and we select only zones with outflowing gas, i.e. with . The areal star formation rate averaged over is ; we use this average in computing all loading factors. Note that a rolling mean of with 100Myr time window gives 50% variation with respect to the mean. Because there is generally a temporal offset between star formation bursts and winds (see Figure 3), the instantaneous ratio of mass or energy fluxes to the star formation rate can significantly over- or under-estimate the true physical loading. The ratio of time-averaged wind fluxes to the time-averaged star formation rate is more meaningful.

Figure 9: Mass and energy loading factors in each thermal phase. For both cold and warm gas, both mass and energy loading factors stiffly decrease as a function of . Mass loading also drops off with for the ionized phase. Only the hot gas is a true wind, with well-defined mass and energy loading factors of and , respectively.

We decompose the gas into four thermal components and present the loading factors of each thermal phase as a function of (Figure 9). The energy flux is always dominated by the hot gas, with the energy loading factor of above . SN feedback causes large outgoing mass fluxes of warm and cold gas within the disk scale height , but the majority of the warm and cold medium has low velocity and cannot travel far from the midplane; this mass flux at small is best thought of as the “upwelling” of turbulent motions within the disk driven by expanding SNRs and superbubbles.888Allowing for the work done by shock-heated hot gas, both isolated and spatially correlated SNe inject a mean spherical momentum/SN to the ISM of (see Kim & Ostriker, 2015a; Kim et al., 2017, and references therein), an order of magnitude greater than the momenta of the initial SN ejecta. Most of this momentum goes into maintaining quasi-equilibrium force balance with gravity in the bulk of the ISM, rather than driving a wind. The mass loading factor of the ionized gas also decreases significantly with . By the time the flow reaches the edge of the simulation domain, is at least a factor larger than all the other components. Therefore, as we have concluded in Sections 3 and 4, only the hot gas forms a genuine“galactic wind,” with a mass loading factor of that is nearly constant above . We note that with the above definitions of and , the Bernoulli parameter for the hot wind is given by .

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

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

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

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

5.2 Comparison with Observations

The low mass loading factor of the hot gas in our simulations is comparable to (or slightly smaller than) the observed mass loading factor estimated in the best studied local starburst M82. Using Chandra X-ray observations, Strickland & Heckman (2009) constrained the “central” mass loading factor of the hot gas to about , with the “central” energy loading factor about . Here, they constrained quantities using a large number of hydrodynamical models to explain the observed diffuse and hard X-rays, which come from the central 500pc region. 999However, we note that in their modeling the hot wind freely expands into a very tenuous medium rather than expanding into a dense ISM. By comparison, it is evident e.g. in Fig. 8 of Kim et al. (2017) that before superbubble breakout from the warm/cold ISM, the hot gas has very high velocity but mostly remains subsonic. This suggests that lack of a dense surrounding medium in a hydrodynamic wind comparison model might lead to a Mach number higher than would be realistic, and hence for density and temperature constrained by emission properties the mass loss rate could be overestimated by a factor of a few. It is difficult to make a direct comparison, but if we consider the state of the hot medium within the energy-injecting layer () in our simulations, we have and . Although there are no systematic observational studies of mass loading factors of the hot gas, is suggested for a wide range of star formation rates from dwarf starbursts to ultraluminous infrared galaxies, utilizing the Chevalier & Clegg (1985) wind model with observational constraints of the X-ray luminosity and star formation rates (Zhang et al., 2014).

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

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

5.3 Wind Driving Simulation Context

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

Girichidis et al. (2016b) ran a set of simulations with solar neighborhood conditions, focusing on the effect of the SN placement and of SN clustering. The authors emphasize that due to the short duration of their simulations (), definitive conclusions cannot be drawn regarding wind driving. However, their measurements of outflow in the plane can be compared to our fountain flow measurements. For a range of different SN feedback treatments, they report mass loading factors . Although this is larger than what we find ( at ), their simulations do show a decline in outflow over time, and our measurements are at . As pointed out by the authors, the majority of the outflowing mass is relatively dense () and moves slowly (), similar to the properties of our warm fountain at low . Their fraction of high-velocity gas () is only , which correspondingly reduces the mass loading factor of gas that is certain to escape to large distance.

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

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

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

Li et al. (2017) conducted simulations using fixed SN rates, but decomposed by thermal phase in reporting mass and energy (as well as metal) loading factors. Their measurements are at , with their “warm” component including most of what we consider “ionized,” and their “hot” extending down to , slightly lower than our “hot” definition. For their solar neighborhood model, the hot and total mass loading factors are about and , which are larger than ours by a factor of 3 to 8 (see also their discussion in comparison to work of Creasey et al. 2013, 2015; Girichidis et al. 2016b). The energy loading factor is also about an order of magnitude larger in their model.

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

Finally, we remark that a fine enough grid to spatially resolve both low filling-factor warm gas and high filling-factor hot gas in the wind launching region is crucial for proper physical characterization of galactic winds. If warm and hot gas are artificially mixed, e.g. if the flow in AMR and semi-Lagrangian simulations moves from a higher resolution region near the midplane to a lower resolution region at high latitude, the result can be unphysical in ways that would compromise the implications for real galactic systems. For example, consider mixing of warm and hot flows that have total horizontally-averaged mass, momentum, and energy vertical fluxes of , , and . If we further assume a steady state and neglect gravity and magnetic fields, the outgoing fluxes of the mixed gas must be the same as the sum of the horizontally-averaged incoming fluxes of the warm and hot gas. That is, , and similarly and are the sum of “warm” and “hot” terms based on horizontal averages over multiple zones.

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


Typically, (near the disk) will be dominated by the warm medium contribution (see Fig. 9a), will be dominated by the hot medium contribution (see Fig. 4c), and will be dominated by the hot medium contribution (see Fig. 9b).

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

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

Figure 10: The hypothetical post-mixing velocity (black) if artificial numerical mixing were to happen at a distance from the midplane. We use horizontally-averaged outgoing fluxes to calculate using Equation (11), and take a time average over the outflow period. Also shown are mean values of the warm-gas and hot-gas outflowing velocities. The nominal value of is similar to the warm medium velocity at small , and gets closer to the hot medium velocity at large .

6 Summary

Gas flows blown out of galactic midplane regions by energetic stellar feedback, most notably type II SNe, will either escape as a wind or turn around as a fountain, depending on the specific kinetic and thermal energy of the gas compared to the gravitational potential. Due to the multiphase structure of the ISM, the simplest steady, adiabatic solutions (e.g., Chevalier & Clegg, 1985) are not directly applicable, although aspects of these solutions are informative when considering a phase-decomposed analysis. More generally, the question of “how much mass and energy are carried out by each phase of outflowing gas?” depends strongly on when and where star formation and SNe occur and how feedback transfers energy to the ISM. Given the complexity involved, the properties of galactic winds and fountain flows created by the star-forming ISM must be investigated via fully self-consistent numerical simulations.

The TIGRESS implementation described in Paper I provides comprehensive, self-consistent simulations of the multiphase star-forming ISM that can be used for high-resolution investigations of outflows in a wide variety of environments. In this paper, we analyze the solar neighborhood model of Paper I to characterize vertical gas flows phase-by-phase. Our principal finding is that the outflowing gas consists of a superposition of a hot wind and a warm fountain flow. We summarize key conclusions from our analysis of each component below.

6.1 Hot wind

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

2. The hot gas mass loading factor, defined as the ratio of outflowing mass flux to the star formation rate per unit area, is , decreasing by from to . The measured value of is consistent with estimates from idealized numerical experiments of superbubble expansion in a warm-cold ISM.

3. Above , the hot outflow contains a tiny fraction of the energy originally injected in the ISM by SNe. Another tiny fraction is converted to kinetic energy of the warm and cold medium accelerated by superbubble expansion, but most of the original SN energy is radiated away within the disk scale height. The ratio of outflowing hot-gas energy flux to the SN energy injection rate per unit area is , decreasing by from to . For , the mean temperature of the hot gas is , and its mean velocity is . Magnetic energy in the outflow is negligible, with mean Alfvén speed of .

4. In the present simulation, the energy of the hot medium is carried mainly as heat because the outflow is constrained to remain vertical in our local Cartesian box. Allowing for the opening of streamlines at larger distance (beyond our simulation domain), much of the enthalpy would be converted to kinetic energy. The hot medium is barely affected by gravity within the simulation domain, and with an asymptotic maximum wind speed of would easily escape to very large distances. More generally, local simulations with varying galactic conditions can provide measures of and the Bernoulli parameter that can be used to provide predictions for hot wind properties at large distances in a wide range of galaxies.

6.2 Warm Fountain

1. The warm (including both neutral and ionized) medium, which dominates the total gas mass in the simulation, acquires energy from SN feedback, mostly in kinetic form. However, the typical outgoing velocity of the warm medium above is only , which is not enough to overcome the large-scale gravitational potential in the present simulation. As a result, most of the warm gas that is blown out of the midplane eventually turns around and falls back, forming a fountain.

2. Star formation bursts lead to a succession of outflow-dominated and inflow-dominated periods of the fountain flow. During both outflow-dominated and inflow-dominated periods, gas with both signs of velocity (i.e. outflows and inflows) is present on both sides of the disk. Considering just the outflowing gas, owing to deceleration and turnaround the outgoing warm-medium mass loading factor is a decreasing function of height . The mean warm-gas mass flux in outflowing “fountain” gas at is , but this drops steeply to at .

3. The value of the warm-medium energy loading factor drops from to over to . Because drops slightly less than , the mean specific energy of the warm medium increases with . The corresponding (volume-weighted) mean velocities of the warm gas at and are and . However, it is important to note that the increase of mean velocity in the warm medium with is primarily due to dropout of low-velocity fluid elements, and does not reflect acceleration with height. Detailed distributions show a secular decrease in the mass and mass flux of high-velocity warm gas with increasing .

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

6.3 Caveats and prospects

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

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

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

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

We are grateful to the referee for helpful report, and to Miao Li and Drummond Fielding for fruitful discussions. This work was supported by grants AST-1312006 from the National Science Foundation, and NNX14AB49G and NNX17AG26G from NASA. Resources supporting this work were provided in part by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and in part by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center. \softwareThis work made use of the Athena MHD code (Stone et al., 2008; Stone & Gardiner, 2009). This work also made use of analysis and visualization softwares including yt (Turk et al., 2011), astropy (Astropy Collaboration et al., 2013), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), IPython (Perez & Granger, 2007), pandas (McKinney, 2010), and VisIt (Childs et al., 2012).


  • Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  • Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
  • Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
  • Bregman (1980) Bregman, J. N. 1980, ApJ, 236, 577
  • Bustard et al. (2016) Bustard, C., Zweibel, E. G., & D’Onghia, E. 2016, ApJ, 819, 29
  • Chen et al. (2010) Chen, Y.-M., Tremonti, C. A., Heckman, T. M., et al. 2010, AJ, 140, 445
  • Chevalier & Clegg (1985) Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44
  • Childs et al. (2012) Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance Visualization–Enabling Extreme-Scale Scientific Insight, 357–372
  • Chisholm & Matsushita (2016) Chisholm, J., & Matsushita, S. 2016, ApJ, 830, 72
  • Chisholm et al. (2017) Chisholm, J., Tremonti, C. A., Leitherer, C., & Chen, Y. 2017, MNRAS, 469, 4831
  • Chisholm et al. (2016a) Chisholm, J., Tremonti, C. A., Leitherer, C., Chen, Y., & Wofford, A. 2016a, MNRAS, 457, 3133
  • Chisholm et al. (2016b) Chisholm, J., Tremonti Christy, A., Leitherer, C., & Chen, Y. 2016b, MNRAS, 463, 541
  • Contursi et al. (2013) Contursi, A., Poglitsch, A., Grácia Carpio, J., et al. 2013, A&A, 549, A118
  • Creasey et al. (2013) Creasey, P., Theuns, T., & Bower, R. G. 2013, MNRAS, 429, 1922
  • Creasey et al. (2015) —. 2015, MNRAS, 446, 2125
  • de Avillez (2000) de Avillez, M. A. 2000, MNRAS, 315, 479
  • de Avillez & Breitschwerdt (2004) de Avillez, M. A., & Breitschwerdt, D. 2004, A&A, 425, 899
  • Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, MNRAS, 414, 3501
  • Erb et al. (2012) Erb, D. K., Quider, A. M., Henry, A. L., & Martin, C. L. 2012, ApJ, 759, 26
  • Fielding et al. (2017) Fielding, D., Quataert, E., Martizzi, D., & Faucher-Giguère, C.-A. 2017, MNRAS, 470, L39
  • Gatto et al. (2015) Gatto, A., Walch, S., Low, M.-M. M., et al. 2015, MNRAS, 449, 1057
  • Gatto et al. (2017) Gatto, A., Walch, S., Naab, T., et al. 2017, MNRAS, 466, 1903
  • Gent et al. (2013) Gent, F. A., Shukurov, A., Fletcher, A., Sarson, G. R., & Mantere, M. J. 2013, MNRAS, 432, 1396
  • Girichidis et al. (2016a) Girichidis, P., Naab, T., Walch, S., et al. 2016a, ApJ, 816, L19
  • Girichidis et al. (2016b) Girichidis, P., Walch, S., Naab, T., et al. 2016b, MNRAS, 456, 3432
  • Hanasz et al. (2013) Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJ, 777, L38
  • Hawley et al. (1996) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690
  • Heckman et al. (2015) Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147
  • Heckman & Borthakur (2016) Heckman, T. M., & Borthakur, S. 2016, ApJ, 822, 9
  • Heckman et al. (2000) Heckman, T. M., Lehnert, M. D., Strickland, D. K., & Armus, L. 2000, ApJS, 129, 493
  • Heckman & Thompson (2017) Heckman, T. M., & Thompson, T. A. 2017, ArXiv e-prints, arXiv:1701.09062
  • Hennebelle & Iffrig (2014) Hennebelle, P., & Iffrig, O. 2014, A&A, 570, A81
  • Hill et al. (2012) Hill, A. S., Joung, M. R., Mac Low, M.-M., et al. 2012, ApJ, 750, 104
  • Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581
  • Hopkins et al. (2012) Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522
  • Hopkins et al. (2017) Hopkins, P. F., Wetzel, A., Keres, D., et al. 2017, ArXiv e-prints, arXiv:1702.06148
  • Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  • Iffrig & Hennebelle (2017) Iffrig, O., & Hennebelle, P. 2017, A&A, 604, A70
  • Joung & Mac Low (2006) Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266
  • Katz (1992) Katz, N. 1992, ApJ, 391, 502
  • Kim & Ostriker (2015a) Kim, C.-G., & Ostriker, E. C. 2015a, ApJ, 802, 99
  • Kim & Ostriker (2015b) —. 2015b, ApJ, 815, 67
  • Kim & Ostriker (2017) —. 2017, ApJ, 846, 133
  • Kim et al. (2013) Kim, C.-G., Ostriker, E. C., & Kim, W.-T. 2013, ApJ, 776, 1
  • Kim et al. (2017) Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834, 25
  • Kimm & Cen (2014) Kimm, T., & Cen, R. 2014, ApJ, 788, 121
  • Kimm et al. (2015) Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900
  • Korpi et al. (1999) Korpi, M. J., Brandenburg, A., Shukurov, A., Tuominen, I., & Nordlund, Å. 1999, ApJ, 514, L99
  • Leroy et al. (2015) Leroy, A. K., Walter, F., Martini, P., et al. 2015, ApJ, 814, 83
  • Li et al. (2017) Li, M., Bryan, G. L., & Ostriker, J. P. 2017, ApJ, 841, 101
  • Li et al. (2015) Li, M., Ostriker, J. P., Cen, R., Bryan, G. L., & Naab, T. 2015, ApJ, 814, 4
  • López-Cobá et al. (2017) López-Cobá, C., Sánchez, S. F., Moiseev, A. V., et al. 2017, MNRAS, 467, 4951
  • Mac Low & Ferrara (1999) Mac Low, M.-M., & Ferrara, A. 1999, ApJ, 513, 142
  • Mac Low & McCray (1988) Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776
  • Mac Low et al. (1989) Mac Low, M.-M., McCray, R., & Norman, M. L. 1989, ApJ, 337, 141
  • Martin (2005) Martin, C. L. 2005, ApJ, 621, 227
  • Martizzi et al. (2016) Martizzi, D., Fielding, D., Faucher-Giguère, C.-A., & Quataert, E. 2016, MNRAS, 459, 2311
  • McCray & Kafatos (1987) McCray, R., & Kafatos, M. 1987, ApJ, 317, 190
  • McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  • Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121
  • Muratov et al. (2015) Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., et al. 2015, MNRAS, 454, 2691
  • Naab & Ostriker (2017) Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59
  • Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9, 21
  • Peters et al. (2017) Peters, T., Naab, T., Walch, S., et al. 2017, MNRAS, 466, 3293
  • Pettini et al. (2001) Pettini, M., Shapley, A. E., Steidel, C. C., et al. 2001, ApJ, 554, 981
  • Pillepich et al. (2017) Pillepich, A., Nelson, D., Hernquist, L., et al. 2017, ArXiv e-prints, arXiv:1707.03406
  • Rodríguez-Puebla et al. (2017) Rodríguez-Puebla, A., Primack, J. R., Avila-Reese, V., & Faber, S. M. 2017, MNRAS, 470, 651
  • Rupke et al. (2005) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115
  • Ruszkowski et al. (2017) Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, ApJ, 834, 208
  • Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521
  • Shapiro & Field (1976) Shapiro, P. R., & Field, G. B. 1976, ApJ, 205, 762
  • Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65
  • Shu (1992) Shu, F. H. 1992, The physics of astrophysics. Volume II: Gas dynamics.
  • Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51
  • Springel & Hernquist (2003) Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289
  • Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289
  • Stone & Gardiner (2009) Stone, J. M., & Gardiner, T. 2009, New A, 14, 139
  • Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142
  • Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  • Strickland & Heckman (2007) Strickland, D. K., & Heckman, T. M. 2007, ApJ, 658, 258
  • Strickland & Heckman (2009) —. 2009, ApJ, 697, 2030
  • Strickland & Stevens (2000) Strickland, D. K., & Stevens, I. R. 2000, MNRAS, 314, 511
  • Thompson et al. (2016) Thompson, T. A., Quataert, E., Zhang, D., & Weinberg, D. H. 2016, MNRAS, 455, 1830
  • Torrey et al. (2017) Torrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2017, MNRAS, 467, 2301
  • Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
  • van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  • Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769
  • Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177
  • Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238
  • Wang (1995) Wang, B. 1995, ApJ, 444, 590
  • Weiß et al. (1999) Weiß, A., Walter, F., Neininger, N., & Klein, U. 1999, A&A, 345, L23
  • Zhang et al. (2014) Zhang, D., Thompson, T. A., Murray, N., & Quataert, E. 2014, ApJ, 784, 93

Appendix A MHD equations in a Shearing-Box

The ideal MHD equations of mass, momentum, and total energy conservation in frame rotating at are




where the total energy density is


The total gravitational potential may include self-gravity and external-gravity, as well as the tidal potential (see Paper I). Note that the zero point of is at the center of the simulation domain.

With the help of shearing-periodic boundary conditions in the horizontal directions, the horizontally-averaged equations are then given by




where the angle brackets denote a horizontal average, . Here, specific energy of the gas, or the Bernoulli parameter, is defined by


where , and the Poynting vector is defined by


with its vertical component of


The right hand side of Equation A7 is the sum of energy source terms due to the Reynolds and Maxwell stresses (on the -surfaces of the box, from shearing boundary conditions – e.g., Hawley et al. 1996; Stone & Gardiner 2010) and the net cooling.

Although the evolution of our simulation is highly dynamic and time-dependent, the fluctuations in the horizontally-averaged gas properties are with respect to a well-defined equilibrium state (see Figure 4), after an early transient period. We thus take time averages to consider a quasi-steady equilibrium state, dropping the time derivative terms in Equations (A5)-(A7) to analyze characteristics of mean gas flows. We explicitly measured each energy source term for separate thermal phases and confirmed that for the hot gas the energy source terms are negligible, and for the warm gas the cooling dominates the energy source terms. The hot gas can thus be treated nearly adiabatically. Above the region where SN energy is injected, the hot gas mass and energy fluxes remain (roughly) constant, and hence the ratio of the two, , is also (roughly) constant. Here


is the specific magnetic energy carried by the magnetic field, which is shown to be small compared to (Figure 5). Thus, outside of the energy injection region where the hot component is created, its mean vertical mass flux and its Bernoulli parameter are expected to be roughly constant, independent of distance relative to the midplane.

We note that more generally, if we consider the time-averaged state of a flow (i.e. ) that is weakly magnetized and has negligible cooling, Equations (A1) and (A3) reduce to and . That is, the Bernoulli parameter is conserved along streamlines, regardless of the geometry of the flow.

Appendix B Numerical Convergence Tests

Due to complex and nonlinear interactions between SN feedback and the surrounding medium, satisfaction of simple physical conditions (e.g., resolving the Sedov stage of a single SN blast wave; Kim & Ostriker 2015a) will not guarantee numerical convergence of the simulated system as a whole. It is therefore important to check robustness of our main results for varying numerical resolution.

In addition, it has been suggested that the SN placement can significantly impact the resulting ISM structure and feedback efficiency (e.g., Gatto et al., 2015; Li et al., 2015). In TIGRESS, unlike in most other simulations, the SN locations are not set based on a pre-defined vertical distribution, but are self-determined based on the locations where star clusters form and migrate over time. A fraction of the SNe in our simulation are associated with runaways, and we adopt a distribution for the ejection velocity from clusters based on a binary population synthesis model (Eldridge et al., 2011). This distribution, together with the fraction of binary runaways, is in fact not very certain. It is therefore valuable to test how our results might depend on the placement of SNe by comparing different test runs.

b.1 Numerical Resolution

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

Figure 11: Mass- and mass flux-weighted velocity distributions of the warm medium measured at for different numerical resolutions during the outflow period (same as Figure 8(c) and (e)). The magenta dashed lines in each panel are fits given by Equations (7) and (8).
Figure 12: Mass and energy loading factors of the warm (green) and hot (red) phases (same as Figure 9), for different numerical resolutions.

b.2 SN Placement

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

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

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

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

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

Figure 13: Mass- and mass flux-weighted velocity distributions of the warm medium measured at for different SN prescriptions during the outflow period (same as Figure 8(c) and (e)).
Figure 14: Mass and energy loading factors of the warm (green) and hot (red) phases for different SN prescriptions (same as Figure 9).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description