NIHAO XIX: How supernova feedback shapes the galaxy baryon cycle

TL;DR

This paper uses NIHAO simulations to analyze how supernova feedback influences the galaxy baryon cycle, regulating star formation, driving outflows, and affecting galaxy mass and baryon content across different scales.

Contribution

It provides a detailed characterization of supernova feedback effects on galaxy evolution, including outflow rates, gas recycling, and the impact on baryon content, based on high-resolution simulations.

Findings

SN feedback prevents HI condensation into H2, reducing star formation in dwarfs.

Most outflowing gas is cold and falls back, creating a fountain cycle that reduces galaxy mass.

Galactic winds divert inflowing gas, explaining low baryon content in ultradwarves.

Abstract

We have used the NIHAO simulations to explore how supernovae (SNe) affect star formation in galaxies. We find that SN feedback operates on all scales from the interstellar medium (ISM) to several virial radii. SNe regulate star formation by preventing condensation of HI into H$_2$ and by moving cold neutral gas to the hot HII phase. The first effect explains why the cold neutral gas in dwarf galaxies forms stars inefficiently. The second maintains the hot ISM of massive galaxies (HII vents out at lower masses). At $v_{\rm vir}\simeq 67{\rm\,km\,s}^{-1}$, the outflow rate follows the relation: $\dot{M}_{\rm out}=23\,(v_{\rm vir}/67{\rm\,km\,s}^{-1})^{-4.6}\,{\rm SFR}$. $20\%$ to $70\%$ of the gas expelled from galaxies escapes from the halo (ejective feedback) but outflows are dominated by cold swept-up gas, most of which falls back onto the galaxy on a $\sim 1\,$Gyr timescale. This `fountain feedback' reduces the masses of galaxies by a factor of two to four, since gas spends half to three quarter of its time in the fountain. Less than $10\%$ of the ejected gas mixes with the hot circumgalactic medium and this gas is usually not reaccreted. On scales as large as $6r_{\rm vir}$, galactic winds divert the incoming gas from cosmic filaments and prevent if from accreting onto galaxies (pre-emptive feedback). This process is the main reason for the low baryon content of ultradwarves.

Figures

Tables

Table 1. Table 1: Best-fit parameters for Eq. ( 7 ) at z = 0 𝑧 0 z=0 .

$\eta$	${\eta }_0$	$\alpha$
${\eta }_{\mathrm{cold}\mathrm{ISM}}$	7	-3.2$\pm$ 0.3
${\eta }_{\mathrm{gal}}$	4	-4.6$\pm$ 0.5
${\eta }_{\mathrm{gal}}^{\mathrm{cold}}$	2	-5.6$\pm$ 0.5
${\eta }_{\mathrm{halo}}$	2	-3.5$\pm$ 0.5

Table 2. Table 2: Virial masses, stellar masses, escape speeds, and mean radial velocities of the cold gas, the hot CGM and the hot ISM for the four NIHAO galaxies in Fig.  13 . All speeds are in km ​ s − 1 km superscript s 1 \rm km\,s^{-1} .

Name	$M_{\mathrm{vir}}$ [$M_{\odot }$]	$M_{\ast }$ [$M_{\odot }$]	$v_{\mathrm{esc}}$	${\overline{v}}_r^{\mathrm{cold}}$	${\overline{v}}_r^{\mathrm{hot}\mathrm{ISM}}$	${\overline{v}}_r^{\mathrm{hot}\mathrm{CGM}}$
g1.23e10	$1.05\times {10}^{10}$	$1.60\times {10}^6$	$79.4$	$-0.85$	not detected	$12.0$
g1.90e10	$2.70\times {10}^{10}$	$1.20\times {10}^7$	$114$	$-2.61$	$23.1$	$15.2$
g3.21e11	$3.61\times {10}^{11}$	$3.73\times {10}^9$	$230$	$6.67$	$57.0$	$13.6$
g1.12e12	$1.29\times {10}^{12}$	$7.95\times {10}^{10}$	$517$	$-41.0$	$32.4$	$-3.66$

Equations

$${\dot{\rho }}_{\star }={ϵ}_{\mathrm{sf}}\frac{\rho }{t_{\mathrm{dyn}}}$$

$$T=4200\mathrm{K}\times 10^{-0.25a(z)-0.86a^2(z)}{\left(\frac{\rho }{{\rho }_{\mathrm{b}}}\right)}^{\frac{2}{3}},$$

$$f_{\mathrm{ej}}^{\mathrm{gal}}=1-\frac{M_{\mathrm{gal}}}{M_{\mathrm{gal},\mathrm{acc}}},$$

$$f_{\mathrm{ej}}^{\mathrm{halo}}=1-\frac{M_{\mathrm{b}}}{M_{\mathrm{b},\mathrm{acc}}}.$$

$$t_{\mathrm{sf}}=\frac{M_{\mathrm{gas}}}{\mathrm{SFR}}$$

$$\eta =\frac{1}{\mathrm{SFR}}\frac{\Delta M_{\mathrm{out}}}{\Delta t},$$

$$\eta ={\eta }_0{\left(\frac{v_{\mathrm{vir}}}{100{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{\alpha },$$

$${\dot{M}}_{\mathrm{out},\max }=\frac{2{ϵ}_{\mathrm{SN}}{E}_{\mathrm{SN}}{\Psi }_{\mathrm{SN}}}{v_{\mathrm{esc}}^2}\mathrm{SFR},$$

$$v_{\mathrm{esc}}(r_{\mathrm{vir}})=\sqrt{\frac{2[c-\ln (1+c)]}{\ln (1+c)-c/(1+c)}}{v}_{\mathrm{vir}},$$

$$M_{\mathrm{gal},\mathrm{acc}}-M_{\mathrm{gal}}=\int {\dot{M}}_{\mathrm{out},\mathrm{cold}\mathrm{ISM}}\mathrm{d}t,$$

$$t_{\mathrm{reaccr}}=2\frac{r_{\mathrm{vir}}}{v_{\mathrm{vir}}}{\int }_0^{x_{\max }}\frac{\mathrm{d}x}{\sqrt{2\frac{c}{f(c)}\frac{\ln (1+x)}{x}}},$$

$$\frac{1}{2}{\left(\frac{v_0}{v_{\mathrm{vir}}}\right)}^2=\frac{c}{f(c)}\left[1-\frac{\ln (1+x_{\max })}{x_{\max }}\right]$$

$$t_{\mathrm{reaccr}}\simeq 1.1{\left(\frac{M_{\mathrm{vir}}}{10^{11}{M}_{\odot }}\right)}^{-0.18}Gyr.$$

$${\eta }_{\mathrm{gal}}=23{\left(v_{\mathrm{vir}}/67{\mathrm{km}\mathrm{s}}^{-1}\right)}^{-4.6}$$

$${\eta }_{\mathrm{gal}}^{\mathrm{cold}}=23{\left(v_{\mathrm{vir}}/67{\mathrm{km}\mathrm{s}}^{-1}\right)}^{-5.6}$$

$${\eta }_{\mathrm{gal}}={\eta }_{\mathrm{gal}}^{\mathrm{cold}}=23{\left(v_{\mathrm{vir}}/67{\mathrm{km}\mathrm{s}}^{-1}\right)}^{-2}$$

Full text

††pagerange: NIHAO XIX: How supernova feedback shapes the galaxy baryon cycle–References††pubyear: 2018

1 Introduction

Matching the abundances of galaxies and dark matter (DM) haloes (Leauthaud et al., 2012; Papastergis et al., 2012; Behroozi et al., 2013; Moster et al., 2013; Tollet et al., 2017), the Tully-Fisher relation (Geha et al., 2006; Pizagno et al., 2007; Springob et al., 2007) and the rotation curves of spiral galaxies (Persic et al., 1996) show that a halo of mass $M_{\rm vir}\sim 10^{11}\,M_{\odot}$ contains a galaxy with a typical stellar mass of $M_{\star}\sim 3\times 10^{8}\,M_{\odot}$. Many dwarf galaxies are rich in cold neutral gas due to inefficient star formation (e.g., Boselli et al., 2014; also see the simulations by Munshi et al., 2013 and Hopkins et al., 2014). Its presence can push the baryonic galaxy mass up to $M_{\rm gal}\sim 10^{9}\,M_{\odot}$ (Papastergis et al., 2012). However, even with gas, the galaxy mass is still an order of magnitude lower than the one expected from the universal baryon fraction $f_{\rm b}=0.16$.

The most common explanation is that the rest of the baryons have been blown out by supernova (SN) feedback (Dekel & Silk, 1986). However, a purely ejective scenario requires mass-loading factors of $\eta\sim 50$. For $1\,M_{\odot}$ that forms stars, $50\,M_{\odot}$ of gas need to be expelled, and the wind speed must be larger than the escape speed from the halo to prevent gas from coming back. These conditions are not easily met. Observed stellar outflows have typical mass-loading factors of order unity (Bouché et al., 2012; Cicone et al., 2014; Schroetter et al., 2015, 2016) if not smaller (Kacprzak et al., 2014; Pereira-Santaella et al., 2016; Falgarone et al., 2017).

Hydrodynamic simulations, too, have for a long time struggled to produce systems with mass-loading factors $\eta\gg 1$. Only in the present decade have they begun to form realistic disc galaxies and to reproduce the observed $M_{\star}-M_{\rm vir}$ relation (Guedes et al., 2011; Brook et al., 2012; Hopkins et al., 2014; Shen et al., 2014; Vogelsberger et al., 2014; Schaye et al., 2015; Crain et al., 2015; Wang et al., 2015).

Simulations by Brook et al. (2014), Nelson et al. (2015) and Mitchell et al. (2018) have shown that the combined action of heating by a photoionising UV background and outflows can prevent gas from falling into the shallow gravitational potential wells of haloes with $M_{\rm vir}\ll 10^{11}\,M_{\odot}$. However, galaxies with M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{11}\,M_{\odot} do accrete a baryon mass of order $f_{\rm b}M_{\rm vir}$ and expel from the halo less than half that value. Hence, massive spirals must expel a lot of gas that remains within the virial radius. This gas is potentially available for reaccretion, and Brook et al. (2014) find that almost $40\%$ of the stellar mass at $z=0$ is composed of reaccreted baryons.

Brook et al. (2014)’s $40\%$ figure is based on defining outflows as gas that flows out of a spherical surface of radius $r\sim 0.1r_{\rm vir}$. Nelson et al. (2015), Anglés-Alcázar et al. (2017) and Mitchell et al. (2018) adopted similar criteria111 Anglés-Alcázar et al. (2017) use two stellar half-mass radii, which contain $80\%$ of the total stellar mass on average. For most galaxies, this criterion gives a $r<0.1r_{\rm vir}$. This may explain why Anglés-Alcázar et al. (2017) find higher recycled fractions ($60\%$ to $90\%$) than both Brook et al. (2014) and Christensen et al. (2016).. However, gas may be heated by SNe, vent out of galactic discs, cool, and fall back onto galaxies without crossing this surface, so that the real contribution of ejected and reaccreted baryons may be even larger.

Christensen et al. (2016) added a density and temperature criteria to decide which gas particles belong to galaxies, and they found that between $20\%$ and $70\%$ of the ejected gas was later reaccreted. By tracking gas particles, they could also measure the average reaccretion timescale, which they found to be of the order of $1\,$Gyr. However, they did not follow the evolution of the ejected gas in density and temperature.

Our work is based on a Numerical Investigation of a Hundred Astrophysical Objects (NIHAO, Wang et al., 2015). Besides the size our sample, its novelty is that we track the motions of gas particles in both physical and thermodynamic phase space (that is, on a density-temperature diagram). Observed outflows have a multiphase structure with different speeds for gas that emits in CO, H$\alpha$ and X-rays (e.g., McKeith et al., 1995; Lehnert et al., 1999; Strickland & Heckman, 2009; Beirão et al., 2015; Leroy et al., 2015 for the winds in M82). Our approach allows us to separate the different thermodynamic phases of simulated outflows and to follow their individual reaccretion timescales. We shall show that this separation is important for a meaningful comparison with outflow rates measured from observations.

The primary goal of our study is to disentangle the relative contributions of the different mechanisms and pathways through which stellar feedback reduces the stellar masses of galaxies, since stellar feedback is not only a multiphase but also a multiscale phenomenon. Our study has also the secondary goal to improve the description of feedback in semianalytic models of galaxy formation.

The article begins by presenting the NIHAO simulations (Section 2) and by describing the procedure used to assign particles to the different thermodynamic phases present within haloes (Section 3). Then, it presents our results for the galaxy baryon cycle, which we separate into accretion (Section 4), star formation (Section 5), ejection (Sections 6 and 7) and recycling (Section 8). Section 9 summarises the conclusions of the article.

2 NIHAO simulations

The simulations used for this work are from Wang et al. (2015)’s Numerical Investigation of a Hundred Astrophysical Objects (NIHAO). They were run with cosmological parameters from the Planck Collaboration et al. (2014): Hubble constant $H_{0}=67.1\,\rm km\,s^{-1}\,\rm Mpc^{-1}$, matter density $\Omega_{\rm m}=0.3175$, dark-energy density $\Omega_{\rm\Lambda}=0.6824$, baryon density $\Omega_{\rm b}=0.049$, power-spectrum normalisation $\sigma_{\rm 8}=0.8344$, and power-spectrum slope $n=0.9624$.

From three cubic N-body volumes of side lengths $15$, $20$ and $60\,h^{-1}$Mpc (Dutton & Macciò, 2014), Wang et al. (2015) selected $\sim$ 90 haloes with masses ranging from $10^{9.5}\,M_{\odot}$ to $10^{12.3}\,M_{\odot}$ discarding all haloes with a companion of comparable mass within $3r_{\rm vir}$ at any time. The chosen haloes were resimulated at higher resolution with baryons.

At the highest level of refinement, the force softening length is $\simeq 50-80\,\rm pc$ for the gas particles and $\simeq 100-200\,\rm pc$ for the dark matter particles. The mass of the corresponding gas particles is of the order of $10^{3}\,\rm M_{\odot}$ (see Table 1 in Wang et al., 2015 for more details about the resolution of the simulations).

The initial conditions for the zoom simulations were computed using a modified version of the grafic2 code (Bertschinger, 2011; Penzo et al., 2014). The level of refinement was chosen to keep a roughly constant numerical resolution across the probed mass range. The hydrodynamic simulations were run with the gasoline2 smoothed-particle-hydrodynamics (SPH) code by Wadsley et al. (2004, 2017).

esf-gasoline2 contains a subgrid model for turbulent mixing of metals and energy (Wadsley et al., 2008). The calculation of heating and cooling includes the effects of a photoionising UV background on metal-line cooling (Shen et al., 2010). Ritchie & Thomas (2001)’s implementation of SPH allows a more accurate description of the multiphase structure of the gas and of the destruction of cold clumps. Limiting the timestep when it is much longer than the one for neighbouring particles (as proposed by Saitoh & Makino, 2009) improves accuracy in the description of strong shocks. Finally, the number of neighbours used for SPH has been increased from thirty-two to fifty.

Gas is available for star formation when it has temperature $T<T_{\rm th}=15000\,\rm K$ and number density $n>n_{\rm th}=10.3{\rm\,cm}^{-3}$. Gas that satisfies these conditions is converted into stars at the rate per unit volume:

[TABLE]

where $\epsilon_{\rm sf}=0.1$ is the efficiency of star formation, $\rho$ is the density of the gas available for star formation, and $t_{\rm dyn}$ is the dynamical time (Stinson et al., 2013 for technical details). Star formation is computed using timesteps of $0.8\,\rm Myr$.

The NIHAO simulations contains two types of stellar feedback: pre-SN feedback (radiative feedback, winds from massive OB stars), also referred to as early stellar feedback (ESF), and SN blastwaves.

Assuming a Chabrier (2003) initial mass function and the relation of Torres et al. (2010) between the mass sand the luminosity of an individual star, a stellar population releases a luminous energy of $2\times 10^{50}{\rm\,erg}\,M_{\odot}^{-1}$, most of which comes from short-lived massive stars. Pre-SN feedback is modelled assuming that a fraction $\epsilon_{\rm ESF}$ of this energy is deposited as thermal energy in the interstellar medium (ISM), while cooling continues to be allowed (Stinson et al., 2013). In the NIHAO simulations, $\epsilon_{\rm ESF}=0.13$ (Wang et al., 2015)222For comparison with previous work, the MaGICC simulations (Stinson et al., 2013) used a similar model with $T_{\rm th}=10^{4}\,\rm K$, $n_{\rm th}=9.3{\rm\,cm}^{-3}$, $\epsilon_{\rm sf}=0.017$ and $\epsilon_{\rm ESF}=0.1$..

SN feedback is implemented using the blastwave formalism described in Stinson et al. (2006) but without the superbubble feedback of Keller et al. (2014). The SN rate associated with each stellar particle is computed from its mass and age assuming a SN energy of $10^{51}\,$erg for all stars between $8\,M_{\odot}$ and $40\,M_{\odot}$. For a Chabrier (2003) initial stellar mass function, this assumption gives one SN every $120\,M_{\odot}$ of stars. This energy is distributed to the surrounding gas. To prevent its immediate radiation by the artificially dense surrounding gas (the gas would not be so dense if the multiphase structure of the ISM were properly resolved), radiative cooling is switched off for $\sim 30\,\rm Myr$ for some of the gas particles within the blastwave radius (Stinson et al., 2006 for details of how these particles are selected).

3 Phases and components

We classified gas particles into five categories: cold gas within the galaxy (the cold interstellar medium, ISM), cold gas outside the galaxy but within the halo (the cold circumgalactic medium, CGM), hot gas within the galaxy (the hot ISM), hot gas outside the galaxy but within the halo (the hot CGM), and gas outside the halo (the intergalactic medium, IGM). This classification keeps into account both a particle’s position in physical space (inside the galactic radius, between the galactic radius and the virial radius, outside the virial radius; Section 3.1) and its location in the thermodynamic phase-space (that is, on a density-temperature diagram; Section 3.2).

3.1 Virial radius and galactic radius

The first step of our analysis is to identify the centre, the virial quantities and the radius for each simulated galaxy.

We use the software Pynbody (Pontzen et al., 2013) to identify the main halo in each NIHAO simulation. Its centre is determined by the AHF halo finder (Knollmann & Knebe, 2009). The virial radius $r_{\rm vir}$ is the radius within which the mean density equals $\Delta_{\rm c}\rho_{\rm c}$, where $\rho_{\rm c}$ is the critical density of the Universe and $\Delta_{\rm c}$ is the critical overdensity contrast computed with the fitting formulae by Bryan & Norman (1998), which give $\Delta_{\rm c}\simeq 109$ at $z=0$.

To separate the cold ISM from the cold CGM, we define the galactic radius $r_{\rm g}$ as the radius $r_{\rm HI}$ within which the surface density $\Sigma_{\rm gas}$ of the cold gas ($T<15,000$) is larger than the threshold for star formation in Kennicutt (1998)’s law, $\Sigma_{\rm th}=9\,M_{\odot}{\rm pc}^{-2}$.

The galactic radius $r_{\rm g}$ differs from $r_{\rm HI}$ only when $r_{\rm HI}$ decreases ($r_{\rm g}$ is not allowed to decrease) or when $r_{\rm HI}$ is smaller than the half-mass radius of the stars, $r_{\rm half-mass}$, in which case we set $r_{\rm g}=r_{\rm half-mass}$, but that happens rarely and only in extremely gas poor objects. We also prevent the growth of $r_{\rm g}$ when $r_{\rm HI}$ increases by more than a factor of three from one output timestep to the next as this is invariably a sign of an ongoing merger. We impose these conditions to prevent the possibility that gas particles in the disc may be counted as expelled and reaccreted simply because of stochastic fluctuations in the value of $r_{\rm HI}$.

Figs. 1 to 3 are all for the same NIHAO galaxy g1.12e12, which exemplifies a typical massive spiral (it is the second most massive galaxy in the NIHAO series). We have given its identifier not only to be able to refer to it later on but also to enable comparison with other NIHAO publications. Fig. 1 illustrates the calculation of $r_{\rm g}$ graphically Fig. 2 shows the distribution of cold gas in g1.12e12 when the galaxy is viewed face-on; the red dashed circle corresponds to $r_{\rm HI}$ and shows that the outer edge of the Hi disc has been correctly identified. Fig. 3 shows the distribution of cold gas (contours) and hot gas (colour scale) edge-on; one can see that the hot gas is strongly concentrated within the galaxy’s inner region, and that the cold gas is in a disc that thickens at the centre due to entrainment of cold gas by the outflowing hot component.

3.2 Gas phases

Fig. 4a shows the density-temperature distribution for all gas particles in g1.12e12. Figs. 4b and 4c are identical to Fig. 4a except that they show only particles with $r<r_{\rm g}$ and $r_{\rm g}<r<r_{\rm vir}$, respectively.

The gas along the white dashed line $T\propto\rho^{2/3}$ in the bottom left part of Fig. 4a belongs to the IGM because it is still outside $r_{\rm vir}$ (hence its absence in Figs. 4b,c). This gas has such low density and temperature that it cannot radiate efficiently. Hence, its equation of state is virtually adiabatic.

When gas particles from the IGM reach the virial radius, they can take two paths (e.g., Kereš et al., 2005): they can streams onto the central galaxy in cold flows (cold accretion) or they can be shock-heated, form a hot low-density atmosphere (the hot CGM) and later cool onto the galaxy (Fig. 4a).

Galaxy g1.12e12 displays an atmosphere of shock-heated gas because it is a massive spiral with $M_{\rm vir}=1.29\times 10^{12}\,M_{\odot}$. We chose a massive galaxy as an example precisely so that we would be able to see all thermodynamic phases on a density-temperature diagram. The hot CGM is virtually absent in haloes with $M_{\rm vir}<10^{11}\,M_{\odot}$, where all the baryons are accreted cold. However, the rest of our considerations from Fig. 4 do extend to lower masses.

Comparing Figs. 4b and 4c shows that the hot high-density gas in the top right part of the density-temperature diagram and the hot low-density gas in top left part are different thermodynamic phases with different spatial distributions. The hot low-density gas is mainly at $r>r_{\rm g}$. The hot dense gas is confined to $r<r_{\rm g}$ and most of it is within the galaxy itself, although some of it is in a wind (Fig. 3). Hence, we refer to the hot low-density gas and the hot dense gas as the hot CGM and the hot ISM, respectively.333The hot ISM is identified on the basis of a purely thermodynamic criterion. Hence it includes both the genuine hot ISM and a hot dense wind component (the gas in the top right corner of Fig. 4c). We emphasize that the latter is artificial and exists only because Gasoline suppresses radiative cooling of gas particles heated by SN for $\sim 30\,$Myr (Stinson et al., 2006). Hence, its only use is to identify gas that has been recently heated by SNe. Separating cold gas accreting onto a galaxy from cold gas in the galaxy is less straightforward but Figs. 4b,c show that the star-forming gas is almost entirely within the galaxy.

Having identified the five phases that compose the gas in the NIHAO simulations (the IGM, the cold CGM, the cold ISM, the hot ISM and the hot CGM), we are now ready to set out precise criteria to assign a phase to each particle.

An obvious choice is to use a temperature criterion (the maximum temperature for star formation) to separate cold gas from hot gas. All the hot CGM (Fig. 4c) and the hot ISM (Fig. 4b) are above the lines AB and DE, corresponding to $T=15,000\,$K. The problem of using a temperature criterion alone is that this choice would attribute part of the cold gas to the hot CGM, since cold streams within haloes are in temperature range $T\sim 10,000-20,000\,$K.

To separate the cold filaments from the hot atmosphere, we observe that the adiabat BC corresponds to the entropy of IGM. Cooling pushes gas to the right of this adiabat. Shock-heating pushes it to the left. Hence, it is natural to use this adiabat at $T>15,000\,$K to separate the cold CGM from the hot CGM.

An adiabat is specified by its temperature at a given density. Let this density be the mean baryonic density of the Universe $\bar{\rho}_{\rm b}=\Omega_{\rm b}\rho_{\rm c}\simeq 6\,M_{\odot}{\rm\,kpc}^{-3}$. The adiabat that fits best the temperature-density relation around $\bar{\rho}_{\rm b}$ is:

[TABLE]

where $z$ is redshift and $a(z)=z/(1+z)$.

To separate cold gas from gas that is sufficiently hot to vent out of the galaxy, we use a pressure criterion. The line CD is an isobar with pressure equal to the pressure $p_{\rm ISM}\simeq 2\times 10^{-13}{\rm\,dyn\,cm}^{-2}$ of the cold dense star-forming ISM (i.e., the molecular gas). Gas with higher temperature and thus higher pressure will leave this phase but will remain within the galaxy as cold neutral gas if $p_{\rm ISM}$ is not exceeded by a large factor. Hence, gas with $T\leq 15,000\,$K is still considered part of the cold ISM even if it has pressure $p>p_{\rm ISM}$. The dotted-dashed line in Fig. 4a shows that this choice is equivalent to separating the cold ISM from the hot ISM at $p=100p_{\rm ISM}$.

We note that the ISM in the NIHAO simulations is very different from the one of McKee & Ostriker (1977), based on the assumptions of different phases in pressure equilibrium. The hot ISM is not in pressure equilibrium with the cold ISM. It is continually venting out of galaxies (and cooling back onto them, as we shall see in Section 9).

The distinction between cold gas in the galaxy and cold gas in the halo is the only one that cannot be performed based on the thermodynamic phase-space diagram alone and relies on our measurement of $r_{\rm g}$.

Particles at the borders between phases may complicate our analysis because small fluctuations in density and temperature may be misinterpreted as ejection and reaccretion. For this reason, a particle is moved from one phase to another only when the threshold value has been exceed by $10\%$. A particle initially outside the galaxy moves inside the galaxy when it moves to $r<0.9r_{\rm g}$. A particle within the galaxy is counted as ejected only when it moves to $r>1.1r_{\rm g}$. A particle crosses the lines AB or DE and moves to the hot CGM or ISM only when its temperature raises above $16,500\,K$. A particle crosses the lines AB or DE and moves to the cold phase only when its temperature drops below $13,500\,$K. The same $10\%$ elastic border is applied to the pressure $p_{\rm ISM}$ that separates the cold gas from the hot ISM (line CD). The only exception is the adiabat BC, to which we apply a tolerance of $5\%$ only, since many particles in the hot CGM have entropies just above the critical entropy that separates the hot CGM from the cold CGM.

4 Baryon content of galaxies and haloes

Having introduced the NIHAO simulations (Section 2) and our method to analyse them (Section 3), we are now ready to present our results. We start from the baryon fraction within galaxies and haloes (current section). Later we shall investigate the astrophysical processes that are responsible for the baryon fractions we find.

Fig. 5 shows the stellar mass $M_{\star}$ (filled blue circles), the baryonic mass $M_{\rm gal}$ within the galaxy (stars and cold gas, filled red circles), and the total baryonic mass $M_{\rm b}$ within $r_{\rm vir}$ (filled black circles) for the NIHAO galaxies at $z=0$. All masses are shown as a function of $M_{\rm vir}$ and in units of $f_{\rm b}M_{\rm vir}$, where $f_{\rm b}=0.16$ is the cosmic baryon fraction.

The filled blue circles (our $M_{\star}-M_{\rm vir}$ relation) are in reasonable agreement with the results from abundance matching (AM) by Papastergis et al. (2012; blue solid curve).The only small tension with the observations is that our ultradwarves tend to be of lower mass than in AM studies (but notice that our procedure to estimate $r_{\rm g}$ lowers $M_{\star}$ in gas-poor galaxies, and that passive ultradwarves fall into this category). Throughout this article, dwarves are galaxies with 3\times 10^{7}\,M_{\odot}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}M_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3\times 10^{9}\,M_{\odot} (corresponding to 3\times 10^{10}\,M_{\odot}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3\times 10^{11}\,M_{\odot}); ultradwarves and normal spirals are galaxies below and above this range of masses, respectively.

The filled red circles show the same relation for the baryonic galaxy mass $M_{\rm gal}$ (the mass of stars and the cold ISM). Comparing them to the AM relation by Papastergis et al. (2012; red shaded area), who measured $M_{\rm gal}=M_{\star}+1.4M_{\rm HI}$ from a sample with Hi data for each galaxy (the factor $1.4$ accounts for helium) shows that the agreement with the total galaxy mass is also satisfactory. The simulations’ ability to reproduce the stellar and baryonic masses of galaxies simultaneously for a same data set proves that gas fractions are reproduced reasonably well, at least for M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}3\times 10^{10}\,M_{\odot}.

The filled black circles show that the halo baryon content is $60-80\%$ of the cosmic baryon fraction for normal spirals with $M_{\rm vir}>3\times 10^{11}\,M_{\odot}$ but drops to $10-20\%$ of the cosmic baryon fraction for ultradwarves with $M_{\rm vir}<2\times 10^{10}\,M_{\odot}$. As already reported by Wang et al. (2017), these findings are consistent with those by Brook et al. (2014), Hopkins et al. (2014) and Mitchell et al. (2018), even though our baryon fractions tend to be slightly lower.

The open red circles replace $M_{\rm gal}$ with $M_{\rm gal,acc}$, the sum of the masses of all particles ever identified as belonging to the galaxy. $M_{\rm gal,acc}$ includes $M_{\rm gal}$, but it also includes the masses of all the particles that have been expelled from the galaxy. When $M_{\rm gal,acc}/M_{\rm gal}\sim 10$, that means that $90\%$ of the baryons accreted by the galaxy have been ejected from the galaxy.

The impact of the expulsion of matter on the baryon content of galaxies and haloes can be quantified more precisely by introducing the ejected fractions:

[TABLE]

[TABLE]

Fig. 6 (red curve) shows that $f_{\rm ej}^{\rm gal}$ decreases from about $95\%$ at $M_{\rm vir}\sim 4\times 10^{10}\,M_{\odot}$ to about $20\%$ at $M_{\rm vir}\sim 10^{12}\,M_{\odot}$. Low $f_{\rm ej}^{\rm gal}$ are also measured in ultradwarves where there has not been enough star formation to produce massive outflows.

In dwarves, the median ejected fraction from the halo is about half the ejected fraction from the galaxy (one third in the case of ultradwarves). The difference between $f_{\rm ej}^{\rm halo}$ and $f_{\rm ej}^{\rm gal}$ becomes smaller at higher masses. At $M_{\rm vir}>8\times 10^{11}\,M_{\odot}$, nearly all galaxies have $f_{\rm ej}^{\rm halo}>f_{\rm ej}^{\rm gal}$ (Fig. 6, black and red circles). This finding is only apparently paradoxical because most of the baryons ejected from haloes do not come from galaxies. They come from halo gas that has been pushed out by galactic winds or unbound by the thermal energy that SNe have deposited into it (Fig. 7). The NIHAO sample selects galaxies that are fairly isolated (Section 2) but some of the gas that leaves the halo without having ever been part of the galaxy may also come from satellites that have come in and out of the virial sphere.

One may wonder to what extend the baryons that are ejected from haloes but do not come from galaxies may be in fact winds coming from small satellites. We have verified that such winds contribute to the outflow rate from the halo but that our qualitative results do not change if we remove DM substructures and the associate baryons from our analysis.

The existence of haloes with $M_{\rm b}/(f_{\rm b}M_{\rm vir})\sim 0.06$ (Fig. 5) while f_{\rm ej}^{\rm halo}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.8 for all objects (Fig. 6) proves that ejective feedback is not the only mechanism that is responsible for the low baryon fractions of many NIHAO haloes.

To demonstrate that two different processes are at work, we have identified two haloes with $M_{\rm b}/(f_{\rm b}M_{\rm vir})\sim 0.1$. They correspond to the circles enclosed in squares in Fig. 5. The halo with $M_{\rm vir}\sim 1.5\times 10^{10}\,M_{\odot}$ has accreted $60-70\%$ of the universal baryon fraction but has lost $\sim 70\%$ of its baryons. The other halo (the one with $M_{\rm vir}\sim 8\times 10^{9}\,M_{\odot}$) has accreted less than $20\%$ of the universal baryon fraction but has retained most of the accreted baryons. Ejective feedback has played a major role in the first halo but not in the second one, where baryons never made it into the virial radius in the first place. The presence of many haloes with $M_{\rm b,acc}\ll f_{\rm b}M_{\rm vir}$ (Fig. 5, open circles) demonstrates that this is not an isolated case and that inability to accrete is a significant limiting factor for the growth of galaxies in haloes with $M_{\rm vir}<3\times 10^{10}\,M_{\odot}$ (Nelson et al., 2015; Wang et al., 2017; Mitchell et al., 2018).

To understand whether this inability is due to stellar feedback or heating by a photoionising UV background, we have resimulated two galaxies without stellar feedback (triangles in Fig. 5). Fig. 8 shows the flow lines and the temperature distribution of gas in the less massive of the two galaxies ($M_{\rm vir}=1.6\times 10^{10}\,M_{\odot}$) with and without stellar feedback (left and right panels, respectively). The right panel shows that reionisation feedback alone (the only feedback that is present when SN feedback and early stellar feedback are turned off) has limited effect at M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}2\times 10^{10}\,M_{\odot}. Hence, gas streams into the virial radius (dashed circle), accretes onto the galaxy (Fig. 8, right) and forms stars so efficiently that, in Fig. 5, the blue triangles ($M_{\star}$) hide the red triangles ($M_{\rm gal}$) almost entirely. The left panel of Fig. 8 shows that including stellar feedback changes the dynamics of the gas completely. Stellar-driven outflows terminate with a shock front at $r=6r_{\rm vir}$ (dashed circle) and prevent the infalling gas from coming in. The only gas that accretes onto the galaxy is warm ($10^{4}{\rm\,K}<T<10^{5}{\rm\,K}$) gas that cools within a sphere of radius $4r_{\rm vir}$. Fig. 9 shows that the hot outflowing gas drives the expansion of a shock front that propagates in the wind’s direction.

In conclusion, without stellar feedback, the NIHAO galaxies have $M_{\star}>0.4f_{\rm b}M_{\rm vir}$ down to $M_{\rm vir}\sim 2\times 10^{10}\,M_{\odot}$ (Fig. 5, blue triangles). SNe reduces $M_{\star}$ by two orders of magnitudes but not through ejective feedback alone. One order of magnitude is gained at a much lower energy cost by preventing gas from coming while it is still far away. This pre-emptive feedback is dominant mechanism in ultradwarves with M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{10}\,M_{\odot}.

5 Star formation

The efficiency of star formation $\epsilon_{\rm sf}=0.1$ adopted in the NIHAO simulations (Section 2) is calibrated on the local Kennicutt (1998) law and we have already seen that this assumption leads to integrated SFRs in reasonable agreement with the observed stellar masses of galaxies (Fig. 5). Here, however, we want to check the NIHAO galaxies’ star formation efficiencies directly by comparing the global star formation timescale

[TABLE]

with observations from the Herschel Reference Survey (HRS, Boselli et al., 2014).

In Fig. 10, the red circles show the $t_{\rm sf}-M_{\star}$ relation in the NIHAO simulations when $M_{\rm gas}$ (Eq. 5) is the total mass of the cold ISM. The red shaded area is a linear fit to the data for Hi-and-CO-detected late-type gas-rich galaxies (Boselli et al., 2014). In these observations, $M_{\rm gas}=1.3(M_{\rm HI}+M_{\rm H_{2}})$, where the factor $1.3$ accounts for the presence of helium. The NIHAO simulations are in reasonable agreement with the data and recover the observational scaling relation $t_{\rm sf}\propto M_{\star}^{-0.7}$.

For $M_{\star}<10^{9}\,M_{\odot}$, $t_{\rm sf}$ becomes longer than the age of the Universe. Hence, inefficient accretion and the expulsion of gas from gas galaxies and haloes are not the only mechanisms that slow down the conversion of baryons into stars. Inefficient star formation plays a role, too.

The fraction of the available gas that can be converted into stars in a time $t$ is $f=1-{\rm exp}(-t/t_{\rm sf})$. When $t_{\rm sf}$ is much shorter than the age of the Universe, $f$ is of order unity, but that is not true in galaxies with $M_{\star}\sim 10^{8}\,M_{\odot}$, where $t_{\rm sf}$ can take all values between $10$ and $200\,$Gyr (Fig. 10, red circles). Taking a median value $t_{\rm sf}\sim 40\,$Gyr and assuming $t\sim 10\,$Gyr limits the mass that can form stars to $\sim 20\,\%$ of the baryons available in the cold ISM.

The inefficiency of star formation at low $M_{\star}$ derives from the cold ISM’s inability to condense into molecular clouds: we have inspected the density-temperature diagram for low-mass galaxies and we have verified that dwarves with t_{\rm sf}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10\,Gyr contain very little gas with $T<1,000\,$K. We refer to this phenomenon with the name of anti-condensation feedback.

The same conclusion can also be reached from Fig. 10. The star-formation timescale for gas with $T<1,000$ is always short and decreases at low $M_{\star}$ (in agreement with Boselli et al., 2014). Hence, the molecular gas always forms stars effectively.

The proof that anti-condensation feedback is a genuine feedback effect is that nearly all the gas is converted into stars when stellar feedback is switched off (in Fig. 5, the blue triangles overlap with the red ones).

The data of Boselli et al. (2014) do not extend below $M_{\star}=10^{8}\,M_{\odot}$. In fact, they contain very few galaxies with $M_{\star}<10^{9}\,M_{\odot}$. However, the NIHAO simulations predict that the star-formation timescale for the Hi gas cannot keep rising like $t_{\rm sf}\propto M_{\star}^{-0.7}$ at low masses: dwarves galaxies with $M_{\star}<10^{8}\,M_{\odot}$ form so few stars that there are not enough SNe for effective feedback. Ultradwarves with $M_{\star}<10^{7}\,M_{\odot}$ are predicted to exhibit highly stochastic SFRs with intense starbursts followed by long periods of quiescence. Fig. 10 shows that in these systems $t_{\rm sf}$ can take all values between $0.3$ and $200\,$Gyr.

6 Mass loading

In Section 4, we have seen that $90\%$ of the baryons that accrete onto a dwarf galaxy ($M_{\star}\sim 10^{8}\,M_{\odot}$) are expelled from the galaxy and that this fraction goes down to $\sim 20\%$ for massive spirals such as Andromeda and the Milky Way. We have also seen that this is not the only reason why dwarf galaxies have very low $M_{\star}/M_{\rm vir}$ mass ratios. Inability to accrete (Section 4) and inefficient star formation (Section 5) play a role, too. We now want to take our analysis one step forward and establish a direct link between star formation and outflows. Our goal is to extract predictions that may be compared to observations and also to derive physically motivated prescriptions that may improve the description of feedback in semianalytic models of galaxy formation.

Let us start by introducing the mass-loading factor:

[TABLE]

where $\Delta t$ is the time between the consecutive output times $t$ and $t+\Delta t$, and where $\Delta M_{\rm out}$ is the mass of the gas particles ejected between $t$ and $t+\Delta t$.

In this article, we consider four different definitions of ejected particles:

particles that are in the cold ISM at $t$ but no longer belong to it at $t+\Delta t$; 2. 2.

particles that are in the cold ISM at $t$ but have moved to $r>r_{\rm g}$ at $t+\Delta t$; 3. 3.

particles that are in the cold ISM at $t$ but have moved or are in the process of moving to the cold CGM at $t+\Delta t$; 4. 4.

particles that are within the virial radius at $t$ but have moved to $r>r_{\rm vir}$ at $t+\Delta t$.

Each definition serves a different purpose. The first gives the rate at which gas ceases to be available for star formation. The second gives the rate at which gas flows out of the galaxy. The third is useful to compare with observations that trace the cold molecular component444Definition (iii) includes hot ISM’s particles in transit from the cold ISM to the cold CGM but few particles take this route. When they do, the time spent in the intermediate hot-ISM phase is very short. Therefore, particles identified with definition (iii) are effectively cold at all times.. The fourth gives the rate at which gas leaves the halo. We use $\eta_{\rm cold\,ISM}$, $\eta_{\rm gal}$, $\eta_{\rm gal}^{\rm cold}$ and $\eta_{\rm halo}$, respectively, to denote the mass-loading factors computed with these four definitions. Our definitions are such that, for each individual galaxy, $\eta_{\rm gal}^{\rm cold}\leq\eta_{\rm gal}\leq\eta_{\rm cold\,ISM}$ by construction.

Fig. 11 compares the mass-loading factors from the NIHAO simulations with those estimated by Cicone et al. (2014) for M82, NGC 253, NGC 2146, NGC 3256 and NGC 3628 (we took the sample by Cicone et al. and removed all objects in which they found evidence for an active galactic nucleus). The greatest difficulty of this comparison was to find reliable masses for such nearby objects. Their apparent sizes are so large that the SDSS and 2MASS miss a lot of the light from the outer regions. We therefore reverted to extinction-corrected B-, V- and R-band photometry from the Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al., 1991) and converted luminosities into masses using stellar mass-to-light ratios from pegase (Fioc & Rocca-Volmerange, 1997), which kept colours into account. The good agreement of $\eta_{\rm gal}$ (yellow circles) with the observations becomes excellent when the comparison is based on the cold component $\eta_{\rm gal}^{\rm cold}$ (filled blue circles) only, as it should be, since the points by Cicone et al. (2014) are based on CO data. We stress that no free parameter was tuned to obtain this result.

The point from Falgarone et al. (2017) refers to a different redshift range ($2<z<2.6$ instead of the local Universe) and carries considerable uncertainty555Falgarone et al. (2017) used ALMA data to measure mass-loading factors for six submillimetre galaxies with SFR$=300-13000\,M_{\odot}{\rm\,yr}^{-1}$ at $2<z<2.6$ and found values of $\sim 0.1$ for all of them (within the large uncertainties of their measurements). We have no information about the masses of these galaxies but we have estimated their order of magnitude from the SFR - $M_{\star}$ at $z\sim 2$ (Daddi et al., 2007; Tasca et al., 2015). $M_{\star}\sim 2\times 10^{11}\,M_{\odot}$ is the average stellar mass of a galaxy with SFR$=300\,M_{\odot}{\rm\,yr}^{-1}$ at $z=2$. This calculation is more likely to overestimate $M_{\star}$ than underestimate it.. We have added it nevertheless because our $\eta_{\rm gal}^{\rm cold}$ - $M_{\star}$ relation at $z=2.3$ is not significantly different from the one at $z=0$ and this point is our only constraint for $\eta_{\rm gal}^{\rm cold}$ at high masses. Falgarone et al. (2017)’s data point shows that the steep decrease of $\eta_{\rm gal}^{\rm cold}$ with $M_{\star}$ is real and that the NIHAO simulations reproduce the correct relation over two orders of magnitude in stellar mass.

The mass-loading factor $\eta_{\rm gal}^{\rm cold}$ depends on our definition of cold gas. The filled blue circles in Fig. 11 consider all the gas identified as belonging to the cold ISM in Fig. 4b. As the mass-loading factor from Cicone et al. (2014) and Falgarone et al. (2017) are inferred from Co and CH+ tracers, respectively one may wish to apply a more restrictive criterion and consider cold only the molecular gas ($T<1000\,\rm K$). When we do it the blue filled circles become the blue open circles, which are in good agreement with the observations at high masses but seems to underestimate mass-loading factor at lower masses. One should notice, however, that only massive NIHAO galaxies exhibit a dense molecular phase distinct from the atomic gas. In lower mass galaxies, the cold ISM does not contain enough particles for this phase to be properly resolved.

One should also be aware that Cicone et al. measured outflow rates by assuming that the outflowing clouds populate the spherical (or multi-conical) region affected by outflows uniformly. This assumption is reasonable if outflows are continually refilled with clouds ejected from the disc but it overestimates $\eta_{\rm gal}^{\rm cold}$ for stochastic outflows (which we find to be prevalent in dwarf galaxies but not in massive objects, see Section 7).

The comparison with observations has made us confident that the NIHAO simulations approximate the Universe reasonably well. With this confidence, we are now ready to study the dependence of $\eta$ on the virial velocity $v_{\rm vir}$ and the virial mass $M_{\rm vir}$ for each of our four definitions of the mass-loading factor. The red, yellow, blue and black symbols in Fig. 12a correspond to $\eta_{\rm cold\,ISM}$, $\eta_{\rm gal}$, $\eta_{\rm gal}^{\rm cold}$ and $\eta_{\rm halo}$, respectively, while the red, yellow, blue and black solid lines are log-log linear least-square fits of the form:

[TABLE]

to the symbols of matching colours.

Table LABEL:table_eta gives the best-fit parameters $\eta_{0}$ and $\alpha$ for each definition of $\eta$ when we stop the fits at $v_{\rm vir}=45{\rm\,km\,s}^{-1}$. Below this virial velocity, the scatter in $\eta$ becomes huge and the fits would be meaningless. This scatter is due to the stochastic SFR histories of dwarves, characterised by intermittent bursts of star formation interspersed between periods of quiescence (there can be no outflows without type II SNe and there can be no type II SNe without star formation).

The yellow line and the blue line in Fig. 12a have been computed fitting all the points at $v_{\rm vir}>45{\rm\,km\,s}^{-1}$ but we have stopped them at $v_{\rm vir}=67{\rm\,km\,s}^{-1}$ because it would be meaningless to extrapolate the fits where $\eta_{\rm gal}^{\rm cold}>\eta_{\rm gal}$ (the cold component cannot larger than the total outflow rate) or where $\eta_{\rm cold\,ISM}<\eta_{\rm gal}$ (including gas that is heated but remains within the galaxy cannot give lower outflow rates than excluding it). The steeper decline of the blue line compared to the yellow line and of the yellow line compared to the red line shows that at low masses most of the mass is the cold component, while in massive galaxies plenty of gas is heated by SNe but does not escape.

The cyan points with error bars in Fig. 12a compare the results of the NIHAO simulations to measurements from MgII absorption lines by Bouché et al. (2012, triangles), Kacprzak et al. (2014, square) and Schroetter et al. (2015, circles), which trace warm ionised gas with outflow speeds of a few hundred kilometres per seconds which are intermediate between those of cold molecular gas and hot X-ray gas. We also remark that the red line in Fig. 12b is consistent with the $\eta_{\rm c}$ - $M_{\rm vir}$ relation at $z=0$ in Anglés-Alcázar et al. (2017) and that the black line is consistent with the $\eta_{\rm loss}$ - $M_{\rm vir}$ relation at $z=0$ in the same article. Anglés-Alcázar et al. (2017) define $\eta_{\rm c}$ as the ratio of the total mass ejected in winds from the ISM (regardless of its fate after ejection) to the total stellar mass formed in situ from early times down to $z=0$, and $\eta_{\rm loss}$ as the same ratio counting only the particles that have been ejected and never reaccreted.

Besides the sporadic star formation of dwarves discussed above (also see Section 5), there is another physical reason why $\eta$ cannot not keep growing as $v_{\rm vir}^{\alpha}$ with $\alpha\ll-2$ at low $v_{\rm vir}$ (Cattaneo et al., 2017). Outflows require speeds greater than the escape speed $v_{\rm esc}$ and this sets an upper limit to outflow rate:

[TABLE]

(Silk, 2003), where $\Psi_{\rm SN}$ is the number of SNe per unit stellar, $E_{\rm SN}$ is the energy released by one SN, and $\epsilon_{\rm SN}$ is the efficiency with which the energy from SNe is converted into kinetic energy. As $v_{\rm esc}\propto v_{\rm vir}$, Eq.(8) implies $\eta_{\rm max}=\dot{M}_{\rm out,\,max}/{\rm SFR}\propto v_{\rm vir}^{-2}$ if we assume that $\epsilon_{\rm SN}$ is constant.

In reality, there is no reason why $\epsilon_{\rm SN}$ should be constant. On the contrary, there are physical reasons why $\epsilon_{\rm SN}$ should decrease with $v_{\rm vir}$. One of them is that massive haloes contain massive atmospheres of shock-heated gas (Dekel & Birnboim, 2006). The efficiency of SNe is limited by the radiative efficiency of the heated gas, which is higher if the gas is denser. SNe deposit their energy in the ISM, where even the dilute hot component is relatively dense. Only by expanding faster than it can cool can gas retain enough energy to produce an outflow. The hot CGM of massive galaxies confine the outflowing gas and slow down its expansion. Hence, the outflowing gas radiates a much larger fraction of its internal energy.

In Section 7, we shall present another physical reason why $\epsilon_{\rm SN}$ should decrease with $v_{\rm vir}$. However, the NIHAO simulations themselves show that this must be the case, since all fits in Table LABEL:table_eta decrease much faster than $v_{\rm vir}^{-2}$. So do semianalytic models of galaxy formation, since nearly all of them require $\alpha\ll-2$ to be consistent with the slope of the stellar mass function of galaxies at low masses (see, e.g., Guo et al., 2011, Lacey et al., 2016, and Cattaneo et al., 2017, but also Henriques et al., 2013; cyan lines in Fig. 12a666 Fig. 12a shows that the normalisation of $\eta$ varies wildly from one semianalytic model to another. GalICS 2.0 (Cattaneo et al., 2017) reproduces the stellar mass function of galaxies with a lower normalisation than other models because in GalICS 2.0 the ejected gas cannot be reaccreted.).

Still, $\epsilon_{\rm SN}$ cannot keep growing indefinitely at low $v_{\rm vir}$ because it cannot become larger than unity. Hence, there must be a virial velocity below which $\eta\leq\eta_{\rm max}\propto v_{\rm vir}^{-2}$ for all galaxies. Cattaneo et al. (2017) predicted this behaviour based on this consideration in the context of the GalICS 2.0 semianalytic model and argued that it could explain the upturn of the stellar mass function of galaxies at low masses, which is observed both in the local Universe (e.g., Yang et al., 2009; Baldry et al., 2012; Moustakas et al., 2013) and at higher redshifts (e.g., Tomczak et al., 2014).

The NIHAO simulations confirm this prediction. All galaxies with $v_{\rm vir}<45{\rm\,km\,s}^{-1}$ lie below the red dashed line $\eta_{\rm gal}\simeq 90(v_{\rm vir}/45{\rm\,km\,s}^{-1})^{-2}$ in Fig. 12a. This line corresponds to $\epsilon_{\rm SN}=1$ for $E_{\rm SN}\Psi_{\rm SN}=8\times 10^{48}{\rm\,erg\,M}_{\odot}^{-1}$ (Section 2) and $v_{\rm esc}(r_{\rm g})\simeq 2.1v_{\rm vir}$, which is a reasonable estimate for the escape speed from galaxies. Adding early stellar feedback increases the energy injected into the ISM by a factor of four. However, the fact that we find $\epsilon_{\rm SN}=1$ for a reasonable escape speed when we consider the contribution of SN feedback alone confirms the supposition that energy from early stellar feedback is dissipated much more efficiently and that its main role is to dilute the ISM in which SNe are about to explode.

The escape speed from the haloes is higher than the one from the galaxy. Hence, the outflow rate of the halo should be lower than the outflow rate from the galaxy by a factor $v_{\rm esc}^{2}(r_{\rm g})/v_{\rm esc}^{2}(r_{\rm vir})$ Assuming an NFW profile (Navarro et al., 1997), the escape speed to $r>r_{\rm vir}$ from the centre of the halo is:

[TABLE]

where $c$ is the concentration of the NFW profile. The fitting formulae by Dutton & Macciò (2014) give $c\simeq 14$ for a halo with $v_{\rm vir}=45{\rm\,km\,s}^{-1}$ in a cosmology with $\Delta_{\rm c}=109$ at $z=0$ (Section 3.1). For this concentration, Eq. (9) gives $v_{\rm esc}(r_{\rm vir})\simeq 3.6v_{\rm vir}$. Hence, we expect $\eta_{\rm halo}/\eta_{\rm gal}\sim(2.1/3.6)^{2}\simeq 0.3$ at low virial velocities. Fig. 12a shows that this very crude estimate is not far off from reproducing correctly the ratio of the black dashed line to the red dashed line.

Fig. 12a shows the relation between $\eta$ and $v_{\rm vir}$. As $M_{\rm vir}\propto v_{\rm vir}^{3}$, it is not surprising that we find a similar relation between $\eta_{\rm galaxy}$ and $M_{\rm vir}$ (Fig. 12b). However, the proportionality constant between $M_{\rm vir}$ and $v_{\rm vir}^{3}$ depends on redshift. Hence, the dependences of $\eta$ on $v_{\rm vir}$ and $M_{\rm vir}$ cannot both remain the same at $z=0$ and $z=2$. At least one has to evolve. The one that evolves less is the more fundamental one. Figs. 12c,d show that the $\eta$ - $v_{\rm vir}$ relation is more or less the same at at $z=2$ and $z=0$, whereas, for a same $M_{\rm vir}$, $\eta$ is lower at $z=2$, since the escape speed for a same virial mass is higher at high redshift. This finding confirms the theoretical expectation that the gas’s ability to escape is linked to the depth of the gravitational potential well more than it is to halo mass.

We remark that the picture presented in this section is in stark contrast with the one by Muratov et al. (2015, FIRE), who claimed that the logarithmic slope of $\eta(v_{\rm vir})$ is steeper for ultradwarves ($\alpha=-3.2$), shallower for dwarves ($\alpha=-1$), and then very steep again for massive spirals. Muratov et al. (2015) attribute the change of exponent from $\alpha=-3.2$ to $\alpha=-1$ around $v_{\rm vir}=60{\rm\,km\,s}^{-1}$ to a transition from energy-driven to moment-driven winds. Had the NIHAO simulations included momentum feedback (e.g., radiation pressure) and not only energy feedback (thermal feedback from SNe and earlier phases of stellar evolution), we might have discovered an intermediate range of $v_{\rm vir}$ where $\eta(v_{\rm vir})$ is shallower than what we find, although it is equally possible that scatter might have blurred any such effect (we have already remarked that Fig. 12b is consistent with a later FIRE publication by Anglés-Alcázar et al., 2017). We consider, however, that the exponent $\alpha=-3.2$ at $v_{\rm vir}<60{\rm\,km\,s}^{-1}$ is unphysical, and that it derives from defining as outflowing any particle with positive radial velocity, even if its speed is lower than the escape speed from the galaxy, whereas the behaviour that we find is much more consistent with the episodic nature of star formation and outflows in ultradwarves.

We conclude this section by comparing the lessons from mass-loading factors and ejected fractions. In Section 4, we saw that $f_{\rm ej}^{\rm gal}$ decreases with $M_{\rm vir}$ faster than $f_{\rm ej}^{\rm halo}$ (Fig. 6). Hence, we may have expected that $\eta_{\rm cold\,ISM}$ would decrease with $v_{\rm vir}$ faster than $\eta_{\rm halo}$ ($\eta_{\rm cold\,ISM}$ is the relevant quantity for this comparison because the calculation of $f_{\rm ej}^{\rm gal}$ is based on $M_{\rm gal}$, the total mass of stars and cold neutral gas within a galaxy). On the contrary, the $\alpha$ parameters for $\eta_{\rm cold\,ISM}$ and $\eta_{\rm halo}$ are very similar (Table LABEL:table_eta) and $\eta_{\rm halo}\sim\eta_{\rm galaxy}/3$ on average. If anything, it is $\eta_{\rm halo}$ that decreases slightly faster (compare the red line and the black line in Fig. 12a). This apparent paradox arises because $f_{\rm ej}^{\rm gal}$ is linked to the mass $M_{\rm gal,acc}-M_{\rm gal}$ of the baryons that are no longer in the galaxy at $z=0$ but have been inside it at some earlier epoch, while $\eta_{\rm cold\,ISM}$ is linked to the instantaneous rate $\dot{M}_{\rm out,\,cold\,ISM}$ at which gas leaves the cold ISM. If gas that leaves the cold ISM never came back, then:

[TABLE]

but we shall see that this is not the case (Section 8). Hence, the different behaviour of $f_{\rm ej}^{\rm gal}$ and $\eta_{\rm cold\,ISM}$ must be due to reaccretion: $f_{\rm ej}^{\rm gal}$ decreases more rapidly than $\eta_{\rm cold\,ISM}$ because, in massive galaxies, a lot of gas is ejected and reaccreted. This gas contributes to $\eta_{\rm cold\,ISM}$ but not to $f_{\rm ej}^{\rm gal}$.

7 Outflow speed and multiphase structure

In Section 6, we investigated the rate at which gas is expelled from galaxies and we saw that separating the cold component from the total outflow rate is important for an accurate comparison with the observations (Fig. 11). Here, we analyse the outflow speed of each component by studying its distribution of radial velocity $v_{r}$, so that we can separate inflows ($v_{r}<-v_{\rm vir}$), outflows ($v_{r}>v_{\rm vir}$), and hot gas in hydrostatic equilibrium or cold galactic gas ($|v_{r}|\ll v_{\rm vir}$).

Fig. 13 shows the distribution of radial velocity $v_{r}$ for the cold gas (cold ISM plus cold IGM), the hot ISM, and the hot CGM of the four galaxies listed in Table LABEL:table_vr. Galaxies g1.23e10 and g1.90e10 typify the bimodality of star formation in ultradwarves, which can be very inefficient (as in g1.23e10) or very intense (as in g1.90e10). Galaxies g3.21e11 and g1.12e12 are representative of massive dwarves/low-mass spirals and massive spirals, respectively, which constitute a more homogeneous population than ultradwarves, even though significant differences from one galaxy to another are observed.

The first result that emerges from Fig. 13 is that the presence of a hot ISM is linked to recent star formation. The passive ultradwarf g1.23e10 displays neither outflows nor a hot ISM. In contrast, a hot ISM is present in the other three galaxies, all of which display star formation and outflows.

Secondly, when dwarf galaxies do exhibit outflows (as in g1.90e10), they are normally dominated by the cold component. At higher masses, the hot dense component carries most of the gas with $v_{r}>200{\rm\,km\,s}^{-1}$. This is the same result that we had already found under a different guise in Figs. 11 and 12. Notice, however, that the bulk of the hot dense gas is not outflowing (all galaxies in Table LABEL:table_vr have $\bar{v}_{r}^{\rm hot\,ISM}<0.2v_{\rm esc}$). This is why we call this phase hot ISM and not wind.

Finally, the distributions of outflow speeds in g3.21e11 and g1.12e12 are not very different, but the escape speeds are. The vertical dashed lines in Fig. 13 show the escape speeds from $r_{\rm vir}$ computed with Eq. 9 but the argument that follows would still apply if we considered escape speeds from $r_{\rm g}$ instead. In g3.21e11, the gas fraction with $v_{r}>v_{\rm esc}$ is small (relative to the total mass of gas)777The observation that the mass in the wind is always small compared to the mass of the galaxy does not preclude the possibility that most of the baryons may be ejected. There are many more healthy people than sick ones. Yet, most of us will fall ill and die at some point. but still significant. In g1.12e12, no outflow is detected at $z=0$. In massive galaxies, SNe waste a lot of energy accelerating particles that will not escape. Since the efficiency of SNe that is relevant for the calculation in Eq. (8) refers to the kinetic energy of particles with $v_{r}>v_{\rm esc}$, this is another reason why $\epsilon_{\rm SN}$ decreases with $v_{\rm vir}$ in addition to the appearance of confining atmospheres at high masses discussed in Section 6 (comparing Fig. 13c and Fig. 13d shows that the hot CGM is more important in g1.12e12 than it is in g3.21e11).

8 The baryon cycle

Having discussed accretion (Section 4), star formation (Section 5) and outflows (Sections 6 and 7), we are now ready to address the fate of the ejected gas and close the loop of the galaxy baryon cycle.

Fig. 14 anticipates our most basic result. In dwarves with $M_{\rm vir}\sim 2-4\times 10^{11}\,M_{\odot}$, out of $100\,M_{\odot}$ expelled from the galaxy, $\sim 30\%$ escapes from the halo and never comes back, $\sim 20\%$ remains inside the halo until $z=0$, and $\sim 50\%$ accretes back onto the galaxy, although some of them will be blown out again. In massive galaxies with $M_{\rm vir}>3\times 10^{11}\,M_{\odot}$, the reaccreted fraction rises to $\sim 80\%$ and less than $10\%$ of the gas is expelled from the halo. Thus, reaccretion is statistically the single most likely outcome at all halo masses but especially at the highest. Our recycled fractions (green curve in Fig. 14) are consistent with the recycled fractions at $z=0$ from Fig. 9 of Anglés-Alcázar et al. (2017)).

We have stressed the possibility that some of the reaccreted gas may be ejected again because all figures in this paragraph are for a single ejection event. If a particle is ejected several times, the probability that it is expelled from the halo once will be larger than the $10\%$ to $30\%$ figure above.

To go beyond this final outcome, and study the pathways through which gas is expelled and reaccreted, we begin with our favourite galaxy g1.12e12 and consider all particles within the virial radius at the penultimate output time $z=0.05$ (Fig. 15), starting from those in the hot CGM.

The hot CGM that cools in the time interval $\Delta t=0.216\,$Gyr between $z=0.05$ and $z=0$ is about $1\%$ of the mass in the hot CGM (Fig. 15a). From this, one infers a cooling time of more than $20\,$Gyr, much longer than the one expected from the density, temperature and metallicity of the hot CGM (the cooling function by Sutherland & Dopita, 1993 gives $\sim 0.7{\rm\,Gyr}$). The difference arises because the cooling time that we infer from Fig. 15a is not $(3/2)nkT/{\cal C}$, where ${\cal C}$ is the power radiated per unit volume, but rather $(3/2)nkT/{\cal C-H}$, where ${\cal H}$ is the thermal power that SNe transfer to hot CGM per unit volume.

To test this explanation, we have computed the total power $L_{\rm X}$ radiated by CGM and we have compared it to the power $P_{\rm SN}=E_{\rm SN}\Psi_{\rm SN}{\rm SFR}$ from SNe, where $E_{\rm SN}=10^{51}\,$erg, and $\Psi_{\rm SN}=1/(120\,M_{\odot})$ (we consider the power from SNe only and neglect the contribution from pre-SN feedback because we assume that the latter is dissipated already at the level of the ISM). For galaxy g1.12e12 at $z=0$, we find $L_{\rm X}=1.2\times 10^{42}{\rm\,erg\,s}^{-1}$ and $P_{\rm SN}=1.6\times 10^{42}{\rm\,erg\,s}^{-1}$. Both quantities are likely to be overestimated: the NIHAO simulations use the cooling function by Shen et al. (2010), which considers the presence of an ultraviolet background and gives lower cooling rates than the one by Sutherland & Dopita (1993); the thermal power that SNe inject into the CGM will be lower than $P_{\rm SN}$ because some of the energy from SNe, too, may be dissipated at the level of the ISM and not reach the CGM. Furthermore, $P_{\rm SN}/L_{\rm X}$ varies stochastically from one snapshot to the next. However, within these uncertainties, our calculations support the conclusion that in massive spirals, such as g1.12e12, the cooling rate and the heating rate comparable, and thus balance each other off on average. This explains the very low cooling rate that we infer from Fig. 15a. We remark that this maintenance feedback is equivalent to the one by radio galaxies in clusters (e.g., Bîrzan et al., 2004, Rafferty et al., 2006, Cattaneo et al., 2009, and references therein), only spiral galaxies have more star formation and lower-mass black holes than elliptical galaxies, hence it is SNe that heat the gas.

Physically, SNe heat the hot CGM via two mechanisms: the hot ISM vents out of the galaxy, mixes with the hot CGM, and transfers its thermal energy to it; SN-driven outflows from the galaxy shock-heat the hot CGM. The cold CGM is located mainly in the plane of the disc and galactic wind mixes with the cold CGM only after radiating their thermal energy (in the descending part of the fountain). Hence, SN heating couples much more effectively to the hot CGM phase than to the cold CGM phase.

The prevalent motion in the cold CGM is accretion onto the galaxy (Fig. 15b). Particles that leave the cold ISM are, in their vast majority, particles that have been heated by SNe and transferred to the hot ISM (Fig. 15c). The hot ISM is a highly unstable phase that needs constant replenishment. By $z=0$, more than half of the particles in the hot ISM at $z=0.05$ have either cooled ($43\%$) or they have moved to the hot CGM ($13\%$; Fig. 15d). Hence, cooling is thrice more probable than becoming part of the hot CGM (twice more probable if we consider only cooling in situ, that is, within the ISM).

These results are based on instantaneous mass exchange rates for one galaxy in our sample, which is not a random object but the second most massive galaxy of the NIHAO series. To extend and generalise these findings to the entire NIHAO series, we have computed the mean lifetimes in the hot ISM, the cold CGM, the hot CGM and the IGM ($r>r_{\rm vir}$) for particles that entered these phases after been ejected from the cold ISM (Fig. 16). By construction, these lifetimes cannot be longer than the age of the Universe and they cannot be shorter than the time interval $\Delta t=0.216\,$Gyr between two consecutive outputs (dashed horizontal line).

Fig. 16 shows that the hot ISM is so short-lived that its lifetime is not adequately resolved by our analysis. To quantify the impact that this may have on our results, we have also considered one galaxy (g1.37e11) that was simulated with much closer output timesteps. Using $\Delta t=13\,$Myr allows us to resolve the time spent in the hot ISM, which is larger than $30\,$Myr by construction, since we model SN feedback as in Stinson et al. (2006).

In Fig. 16, the mean lifetimes for g1.37e11 are shown with open squares rather than open circles. The lifetime in the hot ISM is so short ($37\,$Myr) that the yellow square is outside the plotting range: as soon as gas is allowed to cool, it does it in just $7\,$Myr. The blue open square and the red open square lie a factor of two below the blue curve and the red curve, respectively, but their orders of magnitude are reproduced correctly. The black open square is consistent with the black curve.

Our findings for the blue open square and the red open square should not surprise us. Consider a particle that it is in a same phase at the consecutive output times $t_{1}$ and $t_{2}>t_{1}$; we conclude that this particle has been in that phase for a time $\geq t_{2}-t_{1}$. If a snapshot at an intermediate time $t_{3}$ becomes available, the additional information at $t_{3}$ may support this conclusion but also disprove it (if we find that at $t_{3}$ the particle is in a different phase). Hence, a finer time resolution can only shorten the mean time spent uninterruptedly in any phase, and this is what Fig. 16 shows. However, the integrated total time spent in any phase robust to this effect.

With the caveat above in mind, Fig. 16 shows that gas expelled from haloes remains at $r>r_{\rm vir}$ for a very long time (typically, until the simulation ends). In contrast, cold gas within haloes is usually reaccreted on a timescale of order $\sim 1\,$Gyr, which is shorter at high $M_{\rm vir}$. This behaviour has a simple physical explanation. Cold outflows are approximately ballistic. Hence, the reaccretion timescale $t_{\rm reaccr}$ is simply the dynamical time $r_{\rm vir}/v_{\rm vir}$ (in our cosmology, $r_{\rm vir}/v_{\rm vir}=1.9\,$Gyr at $z=0$) times a factor that depends on $v_{0}/v_{\rm vir}$ and concentration only, if we assume that the DM dominates the gravitational potential ($v_{0}$ is the initial speed at which gas is ejected)888In an NFW halo:

$t_{\rm reaccr}=2\frac{r_{\rm vir}}{v_{\rm vir}}\int_{0}^{x_{\rm max}}\frac{{\rm d}x}{\sqrt{2\frac{c}{f\!(c)}\frac{{\rm ln}(1+x)}{x}}},$

where $x_{\rm max}$ is found by solving the energy-conservation equation:

$\frac{1}{2}\left(\frac{v_{0}}{v_{\rm vir}}\right)^{2}=\frac{c}{f\!(c)}\left[1-\frac{{\rm ln}(1+x_{\rm max})}{x_{\rm max}}\right]$

and $f(c)={\rm ln}(1+c)-c/(1+c)$. . Differences due to concentration are small and the dynamical time is the same for all haloes at a given redshift. Thus, if $v_{0}/v_{\rm vir}$ were the same for all haloes, $t_{\rm reaccr}$ would be the same, too. If, on the contrary, $v_{0}$ were constant, then $t_{\rm reaccr}$ would be a strongly decreasing function of $v_{\rm vir}$. Its form would be more complicated than a simple power-law of $v_{\rm vir}$. However, if we assume an NFW gravitational potential with $c\sim 11-14$, a universal outflow speed implies $t_{\rm reaccr}\propto v_{\rm vir}^{-3}$ for $0.4<v_{\rm vir}/v_{0}<0.8$ and $t_{\rm reaccr}\propto v_{\rm vir}^{-2}$ for v_{\rm vir}/v_{0}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.8 ($v_{\rm vir}/v_{0}=0.8$ corresponds to $t_{\rm reaccr}\sim 0.4\,$Gyr at $z=0$). Fig. 13 shows an intermediate situation: $v_{0}$ increases with $v_{\rm vir}$ but not as rapidly as $v_{\rm vir}$. Hence, the reaccretion time decreases with $M_{\rm vir}$ but not as rapidly as assumed by Henriques et al. (2013); Henriques et al. (2015), who used $t_{\rm reaccr}\propto M_{\rm vir}^{-1}$ independent of redshift.

If we fit the blue circles in Fig. 16 with a power-law function of $M_{\rm vir}$, we find:

[TABLE]

(blue dashed line in Fig. 16), although it is possible that Eq. (11) may overestimate the reaccretion time for cold outflows by up to a factor of two (see above for a discussion of how time resolution can affect our results). Eq. (11) gives $t_{\rm reaccr}\sim 0.7{\rm\,Gyr}$ for $M_{\rm vir}=10^{12}\,M_{\odot}$, which becomes $t_{\rm reaccr}\sim 0.4{\rm\,Gyr}$ if we renormalise the blue dashed line so that it passes through the blue square. These reaccretion times are intermediate between those in Oppenheimer et al. (2010) and Christensen et al. (2016; $t_{\rm reaccr}\sim 1\,$Gyr), and those in Anglés-Alcázar et al. (2017; $t_{\rm reaccr}\sim 0.1-0.35\,$Gyr)999Anglés-Alcázar et al. (2017) used a time resolution of $0.1\,$Gyr, higher than the one used to compute the blue circles in Fig. 16 ($0.216\,$Gyr) but lower than the one used to compute the open square ($0.013\,$Gyr).. For comparison, the semiempirical model by Mitra et al. (2015) reproduces the observations with $t_{\rm reaccr}\propto 0.52(M_{\rm vir}/10^{12}\,M_{\odot})^{-0.45}(1+z)^{-0.32}\,$Gyr.

Discussing the lifetime in the hot CGM requires more subtlety. Hot outflows can be reaccreted only if they cool down and it is well known (since Rees & Ostriker, 1977, Silk, 1977 and Binney, 1977) that the cooling time is longer at higher masses. However, if we interpret the red curve in Fig. 16 as a measure of the radiative cooling time, then Fig. 16 gives results that conflict with our estimate for g1.12e12 (Fig. 15a). The reason is that gas can leave the hot CGM not only because it cools but also and more likely because it escapes from the halo (hot outflows). Hence, the typical $1-2\,$Gyr lifetime depicted by the red curve in Fig. 16 measures the time to leave the virial radius (a particle with an average speed of $100\,\rm km\,s^{-1}$ travels $100\,\rm kpc$ in about $1\,\rm Gyr$) more than the cooling time.

This conclusion may appear to conflict with our previous finding that most of the hot CGM has radial velocities lower than the escape speed (Fig. 13). However, the escape speed in Fig. 13 is the escape speed from the center of the halo, while most of the hot CGM is at several galactic radii.

Massive spirals such as g1.12e12 are the only ones that have proper CGM composed of hydrostatic rather than outflowing gas and are not representative of the average NIHAO galaxy in this regard. For these massive spirals, we have recomputed the mean lifetime in the hot CGM when we retain only particles that cool and return to the ISM101010We stress that this lifetime gives a cooling time and not a reaccretion time because after cooling gas must still fall back onto the galaxy. Hence, the reaccretion time will be longer.. We expected that this calculation would give us a longer time than the one we found by considering all particles. On the contrary, when we applied this analysis to g1.12e12 the survival time in the hot CGM decreased from $2\,$Gyr to $0.7\,$Gyr (in even greater disagreement with the cooling time that we infer from Fig. 15a). Our interpretation is that Fig. 15a gives the correct cooling time for the hot CGM of g1.12e12, but this cooling time is an average over two types a particles: a few particles that cool within the radiative cooling time $(3/2)nKT/{\cal C}\sim 0.7\,$Gyr expected if there were no heating and many more that never do.

The final step of our analysis is to study the different paths through which baryons can be ejected and reaccreted and to compute the masses processed through each of them. For this purpose, we introduce the concept of cycle. A cycle begins when a particle leaves the cold ISM. It finishes when the particle returns to the cold ISM (closed cycle) or when the simulation reaches $z=0$ (open cycle), whichever occurs first.

The frequency of a type of cycle (e.g., the cycle cold ISM $\rightarrow$ hot ISM $\rightarrow$ hot CGM $\rightarrow$ cold CGM $\rightarrow$ cold ISM) is quantified by the mass processed through the cycle divided the total mass $M_{\rm gal,acc}$ processed by the galaxy throughout the lifetime of the Universe. If all particles that accrete onto the galaxy went through this cycle once, then the cycle’s frequency would be unity. This is is impossible because some particles must form stars for feedback to occur. However, a cycle’s frequency could still be equal to one if half of the accreted particles formed stars and the other half went through the cycle twice. The cycle must be close for that to happen.

The hot ISM that vents out of galaxies may end up in either the hot CGM or the cold CGM, and the hot CGM may become cold CGM, but gas cannot move from the cold CGM to the hot ISM or the hot CGM (at least not without being reaccreted onto the galaxy)111111In reality, some particles do (Fig. 15a) but these transitions are not significant.. Neither is it physically possible that gas escapes from the galaxy to the ISM without passing through either the cold or the hot CGM, although it can happen in our analysis if the time resolution is inadequate. Hence, cycles can come in many different permutations, but there is a hierarchy of phases that allows to reduce them to four basic types (Fig. 17): i) the cold ISM $\rightarrow$ hot ISM $\rightarrow$ cold ISM cycle, ii) the cold ISM $\rightarrow$ hot ISM $\rightarrow$ cold CGM $\rightarrow$ cold ISM cycle (and its variant without the hot ISM), iii) the cold ISM $\rightarrow$ hot ISM $\rightarrow$ hot CGM $\rightarrow$ cold CGM $\rightarrow$ cold ISM cycle (and its variants without the hot ISM and/or the cold CGM), and iv) any cycle that passes through the IGM. This classification has been formulated for closed cycles but we can easily generalise it to include open cycles by dropping the requirement that particles should return to the cold ISM.

Fig. 18 shows the average number of type i, ii, iii and iv cycles that an average particle makes as a function of $M_{\rm vir}$. The solid curves are for closed cycles only. The dashed curve include open cycles. The yellow solid curve and the yellow dashed curve run close to each other. So do the blue solid curve and the blue dashed curve. The implication is that type i and type ii cycles are predominantly closed. They correspond to internal heating of the ISM by SNe and the galactic fountain, respectively. In contrast, the red solid curve and the red dashed curve are separated by about an order of magnitude. So are the black solid curve and the black dashed curve. Therefore, type iii and type iv cycles are predominantly open (gas that mixes with the hot CGM or leaves the halo is unlikely to be reaccreted).

Type ii cycles cold ISM $\rightarrow$ (hot ISM) $\rightarrow$ cold IGM $\rightarrow$ cold ISM are the most common at all but the highest masses (Fig. 18, blue solid curve). An average gas particle makes one or two of them after entering the galaxy (Anglés-Alcázar et al., 2017, too, find that an average gas particle is ejected and reaccreted twice). The intermediate hot-ISM phase is missing in most cases but especially at low masses (Fig. 19). Hence, type ii cycles trace cold material entrained by the less massive hot component (also see Fig. 3) rather than hot gas that has expanded and cooled. As the galaxy analysed with $\Delta t=0.013\,$Gyr (black square in Fig. 19) follows the same mean relation as the other NIHAO galaxies, we are confident that Fig. 19 depicts a real physical result rather than our inability to detect short spells in the hot ISM.

Type iv cycles are the second most common variety. Their median frequency is about half per particle (Fig. 18, black dashed curve). This result does not conflict with Fig. 14, which shows that only $10\%$ to $30\%$ of the ejected particles escape to $r>r_{\rm vir}$, because Fig. 14 refers to one ejection event, whereas Fig. 18 refers to the lifetime of the Universe. If we roll a dice once, the probability of an ace is one sixth. If we roll it several times, the probability to obtain an ace at least once will be larger.

Type i cycles are rare at low masses but become predominant at $M_{\rm vir}>6\times 10^{11}\,M_{\odot}$ (Fig. 18, yellow solid curve). Type iii cycles are rare at all masses but particularly at the lowest. Even in massive spirals, they affect at most one particle out of ten (Fig. 18, red dashed curve).

As we have seen that the time resolution of our analysis is a potentially important issue, we have compared the curves in Fig. 18 with our findings for g1.37e11, for which we have a snapshot every $13\,$Myr. Better time resolution allows us to detect a much greater number of type i cycles (compare the yellow squares and the yellow curves), but these cycles are so short that, even collectively, their contribution is not significant anyway. If it were, the hot ISM would contain a significant fraction of the ISM’s total mass. On the contrary, in the NIHAO simulations, the hot phase make up \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10\% of the mass ot the ISM on average, although the value can rise to \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}20\% for massive spirals such as the Milky Way. The frequencies of the cycles that matter most (those of type ii and iv) appear to be robust. So are the frequencies of type iii cycles when we lump closed and open cycles together. The only difference is the frequency of type iii closed cycles. Type iii cycles are more likely to be closed in the simulation analysed with $\Delta t=0.013\,$Gyr than they are on average in the simulations analysed with $\Delta t=0.216\,$Gyr (the red filled square is much closer to the red open square than the red solid line is to the red dashed line). However, type iii closed cycles are so rare that this difference concerns $2\%$ of the baryons at most.

Fig. 18 is very important because it summarises and quantifies the two main ways through which the expulsion of matter reduces the masses of galaxies: fountain and ejective feedback. In a typical galaxy, baryons are ejected and reaccreted once or twice and, in the process, $20\%$ to $80\%$ (thus, half on average) are blown out of the virial radius and never come back.

Fountain feedback acts by reducing the time gas spends in the cold ISM (the only phase where it may form stars). To quantify this effect, we have considered all gas particles that accrete onto a galaxy from when they were first accreted to their last time in the cold ISM, which can be when they were lasted ejected, when they formed stars, or when the simulation ended. Fig. 20 splits this time interval into two: the time spent in the cold ISM (black circles) and the time spent in closed cycles outside the cold ISM (red circles). In dwarf galaxies with $10^{10}\,M_{\odot}<M_{\rm vir}<10^{11}\,M_{\odot}$, an average gas particle spends $\sim 2.5\,$Gyr in the fountain and less than $1\,$Gyr in the cold ISM. Hence, it is available for star formation for only a quarter of the time. At $M_{\rm vir}>3\times 10^{11}\,M_{\odot}$, the time available for star formation becomes \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}50\% but the time spent in the fountain remains significant. The galaxy that we have analysed at intervals of $\Delta t=0.013\,$Gyr spends $\sim 10\%$ less time in the fountain than an average galaxy of the same mass (based on our analysis with $\Delta t=0.216$) but the difference is so small that it will not change anything in our conclusions.

9 Discussion and conclusion

The last decades have seen considerable progress in our understanding of the cosmological context in which galaxies form and evolve. These advancements are due to precise measurements of the cosmological parameters (Bennett et al., 2013; Planck Collaboration et al., 2014), galaxy surveys of unprecedented size and depth (SDSS, 2dF, DEEP2, VVDS, UltraVISTA, CANDELS), and hydrodynamic computer simulations (Hopkins et al., 2014; Schaye et al., 2015; Vogelsberger et al., 2014) but, even with today’s most powerful supercomputers, the range of scales that a simulation can explore is still limited. Hence, feedback – the way processes on the scale of stars and the ISM affect the evolution of galaxies at large – has emerged as the main outstanding problem in galaxy formation.

All simulations need subgrid models to capture the physics they cannot resolve, but some models are more physical than others. The FIRE collaboration (Hopkins et al., 2014) has run simulations that include type I and type II SNe, radiation pressure, stellar winds, photoionisation and photoelectric heating. We cannot claim that our simulations follow the interaction of stars with the ISM to a comparable level, but neither is that our objective. Our goal is not to study the physics of stellar feedback but rather to explore their effects on the formation and evolution of galaxies.

The NIHAO simulations (Wang et al., 2015) are perfectly suited for that purpose because they contain ninety galaxies that span the entire mass range from $M_{\rm vir}=4\times 10^{9}\,M_{\odot}$ to $M_{\rm vir}=4\times 10^{12}\,M_{\odot}$, and they have been shown to reproduce several properties of galaxies, which include stellar masses and SFRs (Wang et al., 2015), disc sizes and the Tully-Fisher relation (Dutton et al., 2017), galaxy velocity function (Macciò et al., 2016) and the metallicity in the CGM (Gutcke et al., 2016)

In this article, depletion times for Hi and H2 (Boselli et al., 2014) and mass-loading factors from CO and CH+ data (Cicone et al., 2014; Falgarone et al., 2017) have been added to the list. The NIHAO simulations’ ability to reproduce several observational constraints simultaneously gives us confidence that, even if the physics of stellar feedback are not resolved, their effects are described correctly. This confidence is important. Without it, there would be no reason to believe that what we learn from our simulations applies to the Universe.

Many of our results confirm findings that were already know from previous studies (Munshi et al., 2013; Hopkins et al., 2014; Shen et al., 2014; Mitchell et al., 2018; Anglés-Alcázar et al., 2017), but our work delineates a new coherent picture of stellar feedback as a multiscale and multiphase phenomenon. Stellar feedback operates on all scales from the ISM to several virial radii. Outflows have different components that we must separate for a meaningful comparison with observational data. We summarise our key results proceeding from the smallest scales to the largest ones.

On ISM scales, feedback reduces the efficiency of star formation via two mechanisms: it prevents the condensation of atomic gas into molecular clouds and it moves gas from the cold ISM to a warm/hot phase with $T>1.5\times 10^{4}\,$K. The first process, which we call anti-condensation feedback, is responsible for the long Hi’s depletion times in dwarf galaxies (Boselli et al., 2014 measures values longer than the Hubble time for galaxies with M_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{9}\,M_{\odot}). The specificity of anti-condensation feedback is that it does not expel or heat the gas. It simply maintains it to a temperature $T>1000\,$K at which it cannot reach the density required for star formation. In contrast, the warm/hot ISM with $T>1.5\times 10^{4}\,$K is not a significant mass reservoir in any galaxies, but even less so in dwarves, where gas vents out as soon as it is heated. Its significance is rather that it is the phase from which winds originate.

On galactic scales, galaxies with $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$ can expel more than $95\%$ of the accreted baryons but most of this gas ($50\%$ to $70\%$) cools very rapidly and falls back in a fountainlike motion on a timescale of the order of $1\,$Gyr, which is shorter at higher masses (Eq. 11). Fountain feedback reduces the masses of galaxies by a factor of two to four by reducing the time during which gas is available to form stars (in dwarves, gas spends thrice more time in the fountain than it does in the galaxy).

On the scale of the CGM, SN heating offsets the radiative losses and explains the hot CGM’s very long cooling time. This maintenance feedback is analogous to the one by supermassive black holes in X-ray groups and clusters (e.g., Cattaneo et al., 2009 and references therein).

Expulsion from the halo has probability of $10\%$ to $30\%$ per ejection event in massive spirals and dwarves, respectively. Hence, it is a less likely outcome than reaccretion through the fountain. However, after several ejection and reaccretion episodes $20\%$ to $70\%$ of the baryons that accrete onto galaxies escape from the virial radius (ejective feedback). These baryons seldom come back.

On megaparsec scales, galactic winds can divert the infalling gas while it is still several virial radii away and prevent its accretion onto galaxies. We call this mechanism pre-emptive feedback. The ability of SNe to reduce the gas inflow rate through a sphere of radius $r\sim 0.5-1r_{\rm vir}$ at low masses had already been remarked in cosmological simulations at z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}3 by Mitchell et al. (2018). However, we have been the first to show that outflows can disrupt inflows and affect the kinematics of gas on scales as large as $6r_{\rm vir}$. Our work shows that SNe and not photoionisation are responsible for the low baryonic fractions of haloes with M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3\times 10^{10}\,M_{\odot}.

The low $M_{\star}/M_{\rm vir}$ ratios of galaxies are the result of not one but all the above mechanisms put together. Pre-emptive feedback reduces the baryons that can accrete onto a halo with $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$ to half of the universal baryon fraction. Due to ejective feedback, only $30\%$ of these baryons are retained. Only a quarter of this $30\%$ are in the galaxy because the rest are in the fountain, and, due to anti-condensation feedback, the baryons in the galaxy form stars on a very long timescale (Fig. 10). If an average gas particle spends only $\sim 1\,$Gyr inside the galaxy and the star-formation timescale is $\sim 10\,$Gyr, then only a tenth of the baryons in the galaxy will be converted into stars. Only by multiplying all these factors do we find $M_{\star}/M_{\rm vir}\sim 0.004$, which is what we expect from abundance matching for a galaxy with $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$.

At even lower masses, star formation becomes sporadic. Ultradwarves are characterised by occasional bursts of star formation interspersed by long period of quiescence, during which there are not enough SNe for outflows of any consequence. In these galaxies, anti-condensation feedback and pre-emptive feedback are the main mechanisms.

Massive spirals accrete nearly $100\%$ of the universal baryon fraction and expel only a small fraction of their baryons. This is due to their deeper gravitational potential wells but also to the appearance of massive hot CGMa that confines winds within the central regions of DM haloes. The fountain is still operating but the mass it processes is at least four times lower, so that the average gas particle spends more time in the galaxy than it does in the fountain; and star formation is much more efficient in massive objects, which have a much higher molecular fraction ($t_{\rm sf}\sim 1\,$Gyr). All these differences combined explain why anticondensation, fountain, ejective and pre-emptive feedback are much less effective in massive spirals than they are in dwarves. Maintenance is the only process that becomes more important at high masses for the simple reason that low-mass haloes do not contain much hot gas, hence, it makes no difference whether it cools or not. In massive spirals, we also observe the appearance of a hot ISM, which is practically absent in dwarf galaxies.

Our picture of SN feedback is very different from the frequently portrayed one, in which the ejected gas mixes with the hot CGM, so that its reaccretion is determined by the hot CGM’s radiative cooling time. We find that only a small fraction of the baryons that accrete onto galaxies are ejected and reaccreted through this route121212The red solid curve in Fig. 18 suggests that this fraction is $<1\%$, although the real value could be larger if the black solid curve included a contribution from hot atmospheres that have spilled just a little outside $r_{\rm vir}$..

A physical picture of galactic winds requires considering their multiphase structure. We propose a model where winds have a hot and a cold component. Blastwaves from SNe form the hot ISM, which vents out of galaxies if it has enough pressure to overcome the depth of the gravitational potential well and the pressure of the surrounding CGM; this is the origin of the hot phase. The cold phase is mainly composed of cold gas that has been entrained by the hot wind. The cold phase is completely dominant at low masses, since there is no need of powerful winds to lift gas from the shallow potential wells of dwarf galaxies. At $M_{\star}\sim 10^{11}\,M_{\odot}$, however, the hot component is often dominant. It is only when we separate the cold phase from the hot phase that we are able to reproduce the outflow rates from IRAM (Cicone et al., 2014) and ALMA (Falgarone et al., 2017) over a mass range that spans two orders of magnitude in $M_{\star}$.

Separating cold outflows and hot outflows is important not only for a correct comparison with observations but also because of their different destinies. The hot component (gas that moves to $r>r_{\rm g}$ with $T>1.5\times 10^{4}\,$K) usually escapes from the halo and never comes back. The cold component falls back onto the disc and forms the fountain. Semianalytic models of galaxy formation could be improved by considering the presence of these two components and by using the NIHAO simulations to calibrate their importance:

[TABLE]

for $v_{\rm vir}>67{\rm\,km\,s}^{-1}$ and

[TABLE]

for $v_{\rm vir}<67{\rm\,km\,s}^{-1}$.

Eqs. (12) to (13) in conjunction with Eq. (11) for the cold component’s reaccretion timescale are the equations that we recommend as our reference model. We anticipate, however, that a semianalytic model based on these equations will be successful only if its cooling rate is as low as the one that we measure in the NIHAO simulations, that is, only if it also includes some form of maintenance feedback.

Acknowledgments

We thank G. A. Mamon for assistance in computing the stellar masses of M82, NGC 3256, NGC 3628, NGC 253 and NGC 2146. A. C. would also like to thank C. B. Brook for assistance with pynbody. We thank N. Bouché for his useful comments. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de) and the High Performance Computing resources at New York University Abu Dhabi.