# Towards simulating star formation in turbulent high-z galaxies with mechanical supernova feedback

Taysun Kimm, Renyue Cen, Julien Devriendt, Yohan Dubois, Adrianne Slyz

► To cite this version:

Taysun Kimm, Renyue Cen, Julien Devriendt, Yohan Dubois, Adrianne Slyz. Towards simulating star formation in turbulent high-z galaxies with mechanical supernova feedback. Monthly Notices of the Royal Astronomical Society, 2015, 451, pp.2900-2921. 10.1093/mnras/stv1211. insu-03644881

HAL Id: insu-03644881

https://insu.hal.science/insu-03644881

Submitted on 20 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Towards simulating star formation in turbulent high-z galaxies with mechanical supernova feedback

Taysun Kimm,<sup>1★</sup> Renyue Cen,<sup>1</sup> Julien Devriendt,<sup>2,3</sup> Yohan Dubois<sup>4,5</sup>
and Adrianne Slyz<sup>2</sup>

Accepted 2015 May 27. Received 2015 May 26; in original form 2015 January 21

## ABSTRACT

To better understand the impact of supernova (SN) explosions on the evolution of galaxies, we perform a suite of high-resolution (12 pc), zoom-in cosmological simulations of a Milky Way-like galaxy at $z=3$ with adaptive mesh refinement. We find that SN explosions can efficiently regulate star formation, leading to the stellar mass and metallicity consistent with the observed mass–metallicity relation and stellar mass–halo mass relation at $z \sim 3$. This is achieved by making three important changes to the classical feedback scheme: (i) the different phases of SN blast waves are modelled directly by injecting radial momentum expected at each stage, (ii) the realistic time delay of SNe is required to disperse very dense gas before a runaway collapse sets in, and (iii) a non-uniform density distribution of the interstellar medium (ISM) is taken into account below the computational grid scale for the cell in which an SN explodes. The simulated galaxy with the SN feedback model shows strong outflows, which carry approximately 10 times larger mass than star formation rate, as well as smoothly rising circular velocity. Although the metallicity of the outflow depends sensitively on the feedback model used, we find that the accretion rate and metallicity of the cold flow around the virial radius is impervious to SN feedback. Our results suggest that understanding the structure of the turbulent ISM may be crucial to assess the role of SN and other feedback processes in galaxy formation theory.

**Key words:** galaxies: formation – galaxies: high-redshift – galaxies: ISM.

## 1 INTRODUCTION

Outflows are common in actively star-forming galaxies. Multiwave-
length images of a local starburst galaxy, M82, exhibit a clear bipolar
outflow, with X-ray emitting gas at the centre (e.g. McCarthy, van
Breugel & Heckman 1987; Strickland & Heckman 2009). Lyman-
break galaxies (LBGs) at high redshifts ($z \sim 2-3$) ubiquitously
reveal significant blueshifted absorption of low-ionization metal
lines, such as C\,{\sc ii} $\lambda$1334 (Shapley et al. 2003; Steidel et al. 2010).
These outflows often carry a large amount of gas that is comparab-
le to or even several times larger than the mass of newly formed
stars (Martin 1999, 2005; Steidel et al. 2010; Newman et al. 2012;
Chisholm et al. 2014), indicating that stellar feedback is the main
process shaping galaxy properties. However, little is yet known

about which types of stellar feedback mechanisms are most im-
portant in regulating star formation in the baryon-driven picture of
galaxy formation

Studies show that a variety of galactic phenomena can naturally
be explained by the presence of feedback from stars. Shetty &
Ostriker (2012) and Kim, Ostriker & Kim (2013) demonstrate that
energetic supernova (SN) explosions can drive turbulence in the in-
terstellar medium (ISM), resulting in a low star formation efficiency
of ~1 per cent per dynamical time-scale of a galaxy, consistent with
galactic scale observations (Kennicutt 1998; Bigiel et al. 2008;
Evans et al. 2009). The evolution of the mass–metallicity relation
(MZR; Erb et al. 2006; Maiolino et al. 2008; Zahid et al. 2013;
Steidel et al. 2014) also seems to require that metals are not instan-
taneously recycled to form the next generation stars, but instead
are first blown away from the galaxy (Dalcanton 2007; Finlator &
Davé 2008). Finally, semi-analytic calculations of galaxy forma-
tion suggest that, in order to match the galaxy luminosity functions

<footnote n="1"><sup>1</sup>Princeton University Observatory, Peyton Hall, 4 Ivy Lane, Princeton, NJ 08544-1001, USA</footnote>

<footnote n="2"><sup>2</sup>Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="3"><sup>3</sup>Observatoire de Lyon, UMR 5574, 9 avenue Charles André, F-69561 Saint Genis Laval, France</footnote>

<footnote n="4"><sup>4</sup>Sorbonne Universités, UPMC Univ. Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75005 Paris, France</footnote>

<footnote n="5"><sup>5</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote>*E-mail: kimm@astro.princeton.edu</footnote>
below $L_*$ (e.g. Bell et al. 2003) or the stellar mass to halo mass relation (e.g. Guo et al. 2010; Moster et al. 2010), the regulation of star formation by SN needs to be very effective so that only a small fraction (≈20 per cent) of baryons turns into stars (e.g. Cole et al. 2000; Khochfar et al. 2007; Somerville et al. 2008, references therein).

Attempts have been made to understand the formation of galaxies with SN feedback in a fully cosmological context using numerical methods (Cen & Ostriker 1992; Katz 1992; Navarro & White 1993). However, Navarro & Steinmetz (2000) find that simulated galaxies with kinetic feedback (KFB) rotate too fast, inconsistent with the local Tully–Fisher relation (Tully & Fisher 1977). The angular momentum catastrophe is attributed in part to the fact that dynamical friction of infalling gas is overestimated as standard smoothed particle hydrodynamics (SPH) cannot capture hydrodynamic instabilities at contact discontinuities accurately and gas comes in cold and clumpy (Navarro & Steinmetz 2000; Sijacki et al. 2012). However, simulations with grid-based codes and weak SN feedback also find that galaxies are too compact (Joung, Cen & Bryan 2009; Kimm et al. 2011a; Hummels & Bryan 2012), as gas cools excessively and forms stars at the galaxy centre. How to overcome the overcooling problem (Katz 1992) and robustly include the effects of SNe are important tasks for all numerical galaxy formation studies.

In order to ensure that energy from SNe is not radiated away artificially, several authors introduced the cooling suppression model in which gas near young stars is assumed to be adiabatic for several to tens of million years (Governato et al. 2010; Guedes et al. 2011). With this feedback model, Guedes et al. (2011) find that a realistic disc galaxy with rising rotation curves can be produced within a Λ cold dark matter paradigm if stars are permitted to form only in dense environments ($n_H = 5 \text{ cm}^{-3}$), so that star formation and outflows are strongly clustered. They claim that the collective winds are efficient at removing low-angular momentum gas that is accreted at high redshifts (see also Brook et al. 2011). Oh et al. (2011) analyse the simulation of a dwarf galaxy by Governato et al. (2010), and find that the simulated inner slope of mass density as well as the shape of rotation curves are consistent with those of THINGS galaxies (Walter et al. 2008), demonstrating that strong feedback is crucial to reproduce the observations. Yet the time-scale imposed for the adiabatic phase is controversial, given that it relies on the density of the ISM, which is largely underresolved in the simulations.

Recently, other forms of stellar feedback have been suggested as important mechanisms to control the growth of galaxies. Walch et al. (2012) and Dale et al. (2014) show that photoionization is capable of dispersing a small cloud of mass $10^4 \text{ M}_\odot$ in ~3 Myr by generating an overpressurized H II region (see also Vázquez-Semadeni et al. 2010; Sales et al. 2014). Dale et al. (2014) point out that winds from massive stars can create a cavity of ~10 pc at the centre of $10^4–10^5 \text{ M}_\odot$ clouds, but they alone cannot unbind the gas from the system (see Geen et al. 2015, for a uniform case). One also expects that radiation pressure from ionizing photons has a similar effect as stellar winds (for metal-rich gas), given that their momentum injection rates are comparable (Leitherer et al. 1999). However, Sales et al. (2014) assert that the inclusion of radiation pressure by ionizing photons on top of photoionization has a subdominant effect, because gas near a massive star is transparent to the Lyman-limit photons once ionized and cannot be accelerated continuously (cf. Wise et al. 2012b). The presence of dust can change the conclusion though, as it can also absorb the ionizing radiation. More importantly, when gas is optically thick to infrared (IR) photons, the latter can in principle be multiply scattered supplying 10–100 times the

original momentum of the ionizing radiation (Murray, Quataert & Thompson 2010). Hopkins and collaborators conduct systematic studies to gauge the relative importance of photoheating, radiation pressure, and SN explosions in isolated discs (Hopkins, Quataert & Murray 2011, 2012b,a), and find that radiation pressure with a momentum boost of 10–100 is required to drive violent galactic outflows in actively star-forming galaxies (cf. Rosdahl et al. 2015). A similar conclusion is reached by Agertz et al. (2013) and Aumer et al. (2013).

It should be noted, however, that several inaccurate assumptions are made concerning SN explosions in previous cosmological simulations. First, even though the initial explosion energy of $10^{51}$ erg is used, as guided by observations (e.g. Nomoto et al. 1993), the momentum transfer from SN to the ISM is likely to be underestimated due to limited resolution. Because a large amount of energy is deposited in a very concentrated form, SN blast waves must undergo an adiabatic phase during which radial momentum should increase by up to an order of magnitude (≈3 × $10^5$ km s⁻¹M⊙), compared to the initial ejecta momentum (Taylor 1950; Sedov 1959; Chevalier 1974; Cioffi, McKee & Bertschinger 1988; Blondin et al. 1998). Without modelling this phase properly, hydrodynamic simulations underestimate the momentum transfer from SN. Thornton et al. (1998) performed a parameter study of the explosion in uniform media using one-dimensional hydrodynamic simulations, and showed that the final radial momentum is a weak function of the background density ($\propto \rho^{-0.12}$) and metallicity ($\propto Z^{-0.14}$; see also Geen et al. 2015; Kim & Ostriker 2015). The final momentum transferred to an inhomogeneous medium turns out to be surprisingly similar (Iffrig & Hennebelle 2015; Kim & Ostriker 2015; Martizzi, Faucher-Giguère & Quataert 2015), as it has a weak dependence on density. Although this has been appreciated in studies of the generation of turbulence in the ISM (e.g. Shetty & Ostriker 2008), only recently a few cosmological simulations begun to take into account the momentum at different phases of SN expansions (Kimm & Cen 2014; Hopkins et al. 2014). Secondly, each star particle is often assumed to release SN energy in a single event at a fixed time delay of ~10 Myr (Dubois & Teyssier 2008). This may overestimate the impact of SN, as clustered explosions can work together to create galactic fountains. But at the same time, it is likely that the regulation of star formation is mistimed when the delay for SN explosions is fixed to a single time such as 10 Myr. Aumer et al. (2013) reported that in the absence of feedback mechanisms active before SNe, to regulate star formation in their simulations, they needed to use an early time delay of 3 Myr. This time delay appeared to leave galaxy star formation histories (SFHs) unchanged when radiation pressure was included in the simulations. Lastly, even though up-to-date cosmological simulations can now employ a spatial resolution as high as 10 pc, this is not enough to resolve the turbulent structure of the ISM where the volume-filling density is smaller than the mean (e.g. Federrath & Klessen 2012). Iffrig & Hennebelle (2015) demonstrate that the expansion of an SN out to some fixed radius depends on the local environments of the explosion. This implies that determining the mass swept up from the host cell of an SN may not be trivial. In this respect, questions remains regarding how momentum should be injected on resolved scales.

The aim of this paper is to understand the impact of SN by realistically modelling the momentum injection from SN. Can the momentum input from SN alone possibly explain the observed relations, such as the MZR or stellar mass to halo mass relation? Why do the classical feedback schemes fail to reproduce these? In order to address these questions, we perform a suite of zoom-in, cosmological hydrodynamic simulations with adaptive mesh refinement
<figure><img src="image_2.png" /><figcaption>Figure 1. Composite images of the Nut galaxy at $z = 7, 5$, and 3 from the NutMFBmp run. The pink, blue, and green colours represent the projected distributions of density, temperature, and metallicity, respectively. The image measures 2.715 Mpc (comoving) on a side. The white circle denotes the virial radius of the host halo of the Nut galaxy (41 kpc at $z = 3$, physical), which is fed by roughly three large-scale filaments. It can be seen that the metal-enriched SN-driven winds extend out to several virial radii of the halo at $z = 3$.</figcaption></figure>

(AMR). We describe our simulation set up and a variety of feedback models tested in this work in Section 2. Our results on the evolution of stellar mass, metallicity, and kinematic properties are presented in Section 3. We discuss the effects of each prescription and why some feedback models are more efficient at suppressing star formation in Section 4. Finally, we summarize our findings in Section 5.

## 2 SIMULATIONS

We use the tree-based Eulerian AMR hydrodynamics code, RAMSES (Teyssier 2002), to study the impact of SN feedback on the evolution of a Milky Way (MW)-type galaxy. The initial conditions for the simulations are identical to those used for the NUT suite presented in Powell, Slyz & Devriendt (2011) and Geen, Slyz & Devriendt (2013), which were generated using the MPGRAFIC software (Prunet et al. 2008; Bertschinger 2001) with Wilkinson Microwave Anisotropy Probe 5 cosmological parameters ($\Omega_m = 0.258, \Omega_\Lambda = 0.742, \Omega_b = 0.045, \sigma_8 = 0.8, H_0 = 72 \text{ km s}^{-1} \text{ Mpc}^{-1}$; Dunkley et al. 2009). The computational domain represents a periodic box of $9 \text{ h}^{-1} \text{ Mpc}$ on a side, covered with $128^3$ cells. Within this volume, three nested grids are placed in a spherical region of radius 2 Mpc (comoving) encompassing a dark matter halo of $M_{vir} = 5.5 \times 10^{11} \text{ M}_\odot$ at $z = 0$. The corresponding mass resolution of a dark matter particle is $5.5 \times 10^4 \text{ M}_\odot$. Only the spherical region is allowed to be further refined. Fig. 1 shows the composite images of density, temperature, and metallicity around the Nut halo at three different redshifts ($z = 3, 5$, and 7) from a run with SN feedback.

The Poisson equation is solved using a multigrid method (Guillet & Teyssier 2011) for grids with level 7 $\le l v \le 9$, while a conjugate gradient method is used for more refined grids ($l v \ge 10$; Teyssier 2002). The Euler equations are integrated in supercomoving coordinates (Martel & Shapiro 1998) with an HLLC Riemann solver (Toro, Spruce & Speares 1994). We adopt the typical courant factor of 0.7.

Metal-dependent radiative cooling is modelled with cooling curves from Sutherland & Dopita (1993), down to $\sim 10^4$ K. Gas cooler than $10^4$ K can lose energy further via metal fine structure transitions, following Rosen & Bregman (1995). A uniform ultraviolet background is turned on at $z = 8.5$ using Haardt & Madau

(1996), and modelled as a heating term in the energy equation. However, gas denser than $n_H = 0.01 \text{ cm}^{-3}$ is assumed to be self-shielded from the radiation (Faucher-Giguère et al. 2010; Rosdahl & Blaizot 2012). The simulations start with gas of metallicity $Z = 2 \times 10^{-5}$ at $z = 499$, as the primordial SN explosions in minihaloes are not resolved in our simulations (e.g. Whalen et al. 2008; Wise et al. 2012b).

In order to better resolve the ISM in galaxies, we employ a maximum physical spatial resolution of $12 \text{ pc}^1$ ($=12.5 \text{ Mpc}/2^{20}$) by triggering 10 more levels of refinement when baryonic plus dark matter mass in each cell becomes greater than $4.4 \times 10^5 \text{ M}_\odot$. In addition, we also enforce that gas denser than $n_H = 1, 42.5, 340 \text{ cm}^{-3}$ is always resolved at least to resolutions of 48, 24, and 12 pc cell, respectively. At the final redshift ($z = 3$), the main halo is resolved with $3.0-3.7 \times 10^6$ computational cells, depending on the feedback model. We find that 59, 14, and 5 per cent of the gas mass inside $0.1 R_{vir}$ ($\approx 3.6 \text{ kpc} \sim 2 R_e$) is resolved on 48, 24, and 12 pc resolution elements in the case of the NutMFBmp run. The main galaxy in the NutTFB run uses more high-resolution cells (40, 32, and 24 per cent of the gas mass is resolved on at least 48, 24, and 12 pc cells, respectively), as the ISM becomes denser than that of MFBmp due to the overcooling.

Dark matter haloes are identified with the HaloMaker (Tweed et al. 2009). The virial radius of each halo is defined as the radius within which the mean density is equal to the critical density of the universe ($\rho_{crit} = 3H(z)^2/8\pi G$) times the virial overdensity ($\Delta_{crit} \equiv 18\pi^2 + 82x - 39x^2$; Bryan & Norman 1998), where $x \equiv \Omega_m/(\Omega_m + a^3\Omega_\Lambda) - 1$ and $a$ is the scale factor.

### 2.1 Star formation

We model star formation as a stochastic process (e.g. Katz 1992; Rasera & Teyssier 2006), based on a Schmidt law (Schmidt 1959),

$$ \dot{\rho}_{\text{star}} = \epsilon_{\text{ff}} \rho_{\text{gas}} / t_{\text{ff}}, \quad (1) $$

<footnote n="1">¹ We perform a convergence test by comparing the MFBmp run and its twin version with a higher resolution (6 pc) down to $z = 5.5$, while keeping other parameters fixed, except for the mass of a star particle (203 $M_\odot$). We find that the stellar mass of the main galaxy is converged within 10 per cent at $z = 5.5$.</footnote>
where $\rho_{\text{gas}}$ is the gas density, $t_{\text{ff}} = \sqrt{3\pi/32G\rho_{\text{gas}}}$ is the free-fall time of the gas, and $\epsilon_{\text{ff}}$ is the efficiency of star formation per free-fall time. Specifically, we examine the probability of forming a star particle in a cell denser than a certain threshold density ($n_{\text{SF}}$) in each fine time step. We adopt $n_{\text{SF}} = 42.5 \text{ cm}^{-3}$, motivated by the Larson-Penston density ($\rho_{\text{LP}} \simeq 8.86 c_s^2 / \pi G \delta x^2$) above which isothermal gas undergoes a runaway collapse (Larson 1969; Penston 1969; Gong & Ostriker 2011, 2013), where $c_s$ is the sound speed and $\delta x$ is the scale of the gas cloud. For a typical temperature of the cold ISM ($\sim 30$ K) and the finest size of our computational grid (12 pc), this yields $n_H = 42.5 \text{ cm}^{-3}$. Although the threshold density is slightly lower than the typical mean density of a giant molecular cloud (GMC, $n_H \sim 100 \text{ cm}^{-3}$), we note that the majority ($\gtrsim 70$ per cent) of stars form in such environments, as star formation is generally a slow process (see below Kennicutt 1998; Krumholz & Tan 2007; Evans et al. 2009).

The number of stars formed ($N_p$) is determined by drawing a random number from a Poisson distribution with a mean of

$$ \lambda \equiv \epsilon_{\text{ff}}(\rho \Delta x^3 / m_{*,\text{min}}) (\Delta t_{\text{sim}} / t_{\text{ff}}), \quad (2) $$

where $\Delta x$ is the size of a computational cell, and $m_{*,\text{min}} = \alpha N_p n_{\text{SF}} m_H \Delta x_{\text{min}}^3 / X_H$ is the minimum mass of a star particle. Here, $m_H$ is the hydrogen mass, $X_H (=0.76)$ is the mass fraction of hydrogen for the primordial composition, and $\alpha$ is a parameter that controls the minimum mass, which we use $\alpha = 0.264$ to achieve $m_{*,\text{min}} = 610 \text{ M}_\odot$. We assume that 2 per cent ($\epsilon_{\text{ff}} = 0.02$) of the dense gas ($n_H \gtrsim 42.5 \text{ cm}^{-3}$) turns into stars per free-fall time (Krumholz & Tan 2007).

## 2.2 Feedback models

We examine six different feedback models to investigate the impact of core-collapse SN on the evolution of galaxies. We use the Chabrier initial mass function (IMF; Chabrier et al. 2005) throughout this study. We assume that 31 per cent of the total mass is lost via Type II SN explosion. For a progenitor mass of $19.1 \text{ M}_\odot$, this means that a star particle of $610 \text{ M}_\odot$ would host 10 SNe. The corresponding specific energy per $M_\odot$ at the time of explosion is $1.6 \times 10^{49}$ erg. The metallicity of the ejecta is taken to be 0.05, which is equivalent to a yield of 0.016 for a simple stellar population with $1 \text{ M}_\odot$.

Feedback models may be distinguished by how the energy or momentum is put into the ISM. Note that we do not turn off the radiative cooling (Mori et al. 1997; Stinson et al. 2006) or hydrodynamically decouple winds from the ISM (Oppenheimer & Davé 2006; Scannapieco et al. 2006). We describe the details of each model below.

### (i) Reference run (NutCO)

The reference run (NutCO) does not include any energy or momentum feedback, but gas and newly processed metals from SNe are added to their host cell, once a star particle becomes older than 10 Myr (so called 'metal feedback'). Since there is no direct energy input in this case, star formation is essentially controlled by the cosmic inflow and uniform background UV heating. Note that we use a higher minimum mass for star particles ($m_{*,\text{min}} = 2310 \text{ M}_\odot$) than other runs to reduce the computational cost.

### (ii) Thermal feedback (NutTFB)

The total energy of $10^{51}$ erg per SN is simply added to the host cell of the star particle, along with ejected gas and metals. Each star particle is assumed to generate a single SN event at 10 Myr after the birth.

This can increase the temperature of a cell up to $\sim 3 \times 10^8$ K in a low-density medium (i.e. $M_{\text{gas}} \ll M_{\text{ej}}$). In principle, the overpressurized gas can drive the expanding motion depositing $\sim 30$ per cent of the total into the kinetic energy (Chevalier 1974). However, as is well known in the literature, the atomic and metal cooling processes can rapidly radiate the internal energy away before the SN-driven blast wave sweeps up the surrounding medium in dense environments or in simulations where the cooling radius is underresolved (Katz 1992; Navarro & White 1993; Abadi et al. 2003; Slyz et al. 2005; Hummels & Bryan 2012; Kimm & Cen 2014). To circumvent this issue, some studies assume that gas is adiabatic for several tens of million years (Mori et al. 1997; Thacker & Couchman 2001; Stinson et al. 2006; Governato et al. 2007; Teyssier et al. 2013). Our thermal feedback (TFB) model does not include the prescription to properly address excessive radiative losses.

### (iii) Kinetic feedback with a small mass loading (NutKFB)

We also test the KFB scheme of Dubois & Teyssier (2008), in which the initial explosion energy of $E_{\text{SN}} = 10^{51}$ erg is modelled by increasing the kinetic energy$^2$ of the surrounding gas. An important parameter in this model is a mass-loading factor ($\eta_w \equiv M_{\text{load}}/M_{\text{ej}}$), which accounts for how much gas is entrained with the SN ejecta. Notice that a mass-loading factor ($\eta$) is often defined as the mass ratio between the outflow rate ($\dot{M}_{\text{out}}$) and the star formation rate (SFR, $\dot{M}_{\text{star}}$), which is different from what we use in this study. Empirical determinations of the mass-loading factor suggest that it varies depending on the mass of a galaxy (Martin 2005; Rupke, Veilleux & Sanders 2005; Weiner et al. 2009; Martin et al. 2012; Newman et al. 2012; Arribas et al. 2014; Chisholm et al. 2014), but is roughly of the order of $\eta_w \sim 10$. Motivated by this, we adopt $\eta_{w, \text{max}} = 10$ (i.e. $\eta \sim 3$). Notice, however, that $\eta_w$ can be as small as zero if there is little gas in the host cell of SN. The velocity of the neighbouring gas is then determined by $v = \sqrt{2 E_{\text{SN}}/(1 + \eta_w) \dot{M}_{\text{ej}}}$, hence the input momentum can be as 3.3 times as high as the initial momentum if there is enough gas to entrain from the host cell. Again, SNe are modelled as a single event at 10 Myr per individual star particle, but this is evaluated at every coarse time step ($\sim 0.16$ Myr) to reduce the message passing interface communications between computing cores.

### (iv) Mechanical feedback (NutMFB)

Kimm & Cen (2014) introduced a physically based feedback model, which can properly explain the transfer of momentum at all stages of a Sedov-Taylor blast wave (from the adiabatic to snowplough phase; see also Hopkins et al. 2014). The main difference from the KFB model is that the input momentum is calculated according to the stage of the Sedov-Taylor blast wave. Although the momentum available during the free expansion phase is only $\approx 4.5 \times 10^4$ km s$^{-1}$ M$_\odot$, it increases as more gas is swept up by the shock ($p \propto \sqrt{m}$). When the energy lost at the shell becomes significant ($\Delta E/E \gtrsim 0.3$) and the pressure inside the bubble becomes comparable to that of the pre-shock region, the expansion of the shell becomes momentum-conserving. At this point, the shell mass gain is much higher than the decrease in velocity resulting in the net increase of the total momentum by a factor of $\sim 7-10$ at the end of the adiabatic phase. This factor can be smaller if the number of SNe is larger or the background density of gas is higher or the metallicity is higher, as the radiative cooling becomes more significant.

<footnote>$^2$ When the momentum added is cancelled out, the kinetic energy is used to increase the temperature to conserve the total energy. This is also the case for the MFB runs.</footnote>
Specifically, the input momentum during the snowplough phase is calculated as (Blondin et al. 1998; Thornton et al. 1998, see also Geen et al. 2015; Kim & Ostriker 2015)

$$p_{\text{SN,snow}} \approx 3 \times 10^5 \text{ km s}^{-1} \text{M}_\odot E_{51}^{16/17} n_H^{-2/17} Z'^{-0.14}, \quad (3)$$

where $E_{51}$ is the number of SNe, $n_H$ is the hydrogen number density, and $Z'$ is the metallicity in solar units, with a minimum of 0.01 ($Z' \equiv \max[Z/Z_\odot$, 0.01]). Whether the expansion is in the adiabatic or momentum-conserving phase is determined by the mass swept up from the host and neighbouring cell. This is done by calculating the mass ratio ($\chi$) for each 48 neighbours (=N$_{nbor}$) per explosion as

$$\chi \equiv dM_{\text{swept}}/dM_{\text{ej}}, \quad (4)$$

where

$$dM_{\text{ej}} = (1 - \beta_{\text{sn}})M_{\text{ej}}/N_{\text{nbor}}, \quad (5)$$

$$dM_{\text{swept}} = \rho_{\text{nbor}} \left(\frac{\Delta x}{2}\right)^3 + \frac{(1 - \beta_{\text{sn}})\rho_{\text{host}} \Delta x^3}{N_{\text{nbor}}} + dM_{\text{ej}}, \quad (6)$$

$\Delta x$ and $\rho_{\text{host}}$ is the size and gas density of the host cell, respectively, and $N_{\text{nbor}} = 48$ is the number of neighbouring cells (see fig. 15 in Kimm & Cen 2014). Here, $\beta_{\text{sn}}$ determines the mass fraction of ($M_{\text{ej}} + \rho_{\text{host}}\Delta x^3$) to be left in the host cell. We take $\beta_{\text{sn}} = 4/52$ to distribute the mass as evenly as possible to the host and neighbouring cells when the cells are on the same level. If $\chi$ is greater than the transition mass ratio,

$$\chi_{\text{tr}} \equiv 69.58 E_{51}^{-2/17} n_H^{-4/17} Z'^{-0.28}, \quad (7)$$

we input $p_{\text{SN, snow}}$ (equation 3), whereas momentum appropriate to the adiabatic phase is injected otherwise,

$$p_{\text{SN}} = \begin{cases} p_{\text{SN,ad}} = \sqrt{2\chi} M_{\text{ej}} f_e E_{\text{SN}} & (\chi < \chi_{\text{tr}}) \\ p_{\text{SN,snow}} & (\chi \ge \chi_{\text{tr}}) \end{cases}, \quad (8)$$

where $f_e = 1 - \frac{\chi^{-1}}{3(\chi_{\text{tr}}-1)}$ is used to smoothly connect the two regimes. Note that the local properties of gas at which an SN explodes can be significantly different from those of gas at which the star particle is formed, as we allow for the ISM to dynamically evolve for 10 Myr before the SN explosion.

(v) *Mechanical feedback with multiple explosions* (NutMFBm) This model is based on the same mechanical feedback (MFB) scheme as NutMFB, but differs in that SNe from a star particle are distributed in time, as in the real Universe. In cosmological simulations with finite resolution, an individual star particle represents a coeval star cluster, which is usually assumed to have a single explosion at some fixed age. Our models (NutCO, TFB, KFB, MFB) also assume that SNe explode at 10 Myr, meaning that tens of SNe explode simultaneously. In reality, some SN explosions commence at as early as 3 Myr, whereas an 8 M$_\odot$ star evolves off the main sequence only at ~40 Myr. The multiple SNe model is designed to mimic the discrete nature of the multiple explosions by randomly sampling the lifetime of massive stars for each particle. To do so, we use the inverse sampling method for a polynomial fit to the integrated rate of SN Type II (Fig. 2, left-hand panel),

$$n(<t) = \sqrt{a(Z) + b(Z)} t_6 + c(Z), \quad (9)$$

where $t_6 \equiv t/10^6$ yr, $a(Z) = -0.268 + 0.139Z' - 0.575Z'^2$, $b(Z) = 0.042 + 0.022Z' + 0.0790Z'^2$, $c(Z) = -0.086 - 0.170Z' + 0.187Z'^2$, and $Z' = \max(Z, 0.008)/0.02$. We use the population synthesis code STARBURST99 (Leitherer et al. 1999) to compute the SN rates as a function of time and metallicity. Fig. 2

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left: integrated SN rates based on the population synthesis code STARBURST99 (solid) and the corresponding polynomial fits (dashed lines, equation 9). Different colours indicate different metallicities, as shown in the legend. Right: an example of the random sampling of discrete SN explosions for a 1000 M$_\odot$ star particle.</figcaption></figure>

(right-hand panel) shows an example of the random sampling and corresponding mass evolution of a 1000 M$_\odot$ star particle. Note that, unlike NutKFB and MFB, we evaluate the explosion at every fine time step (several hundreds to thousand years) for a better accuracy.

(vi) *Mechanical feedback with multiple explosions in a porous ISM* (NutMFBmp)

Particular attention is paid to how momentum is distributed to the neighbouring cells in this model. In NutMFB and NutMFBm, all of the gas in the host cell of SN is re-distributed uniformly to its surroundings. This may be a reasonable assumption for an ISM with a smooth structure or for an ISM with a low Mach number where the volume-filling density is not very different from the average density (Vazquez-Semadeni 1994; Federrath & Klessen 2012). However, high-z galaxies are normally turbulent (e.g. Förster Schreiber et al. 2009), and our resolution elements may not fully resolve the porous structure of the turbulent ISM. If an SN explodes in a porous medium, energy will preferentially propagate through low-density channels (Iffrig & Hennebelle 2015; Kim & Ostriker 2015; Martizzi et al. 2015), unlike SN bubble expansion in a uniform medium. In this regard, it is likely that too much gas is entrained from the host cell of the SN in simulations with limited resolution (including ours), leading to the underestimation of the initial gas outflow velocity in the neighbouring cells. For example, if an SN explodes in a cell with $n_H = 10 \text{ cm}^{-3}$ resolved to 24 pc and a neighbouring cell has a density of $1 \text{ cm}^{-3}$, the initial wind velocity would be about 30 km s$^{-1}$ for a uniform ISM. However, if the ISM were highly turbulent with a Mach number of 10, as is often observed in the Galactic star-forming clouds (e.g. Roman-Duval et al. 2010), the volume-filling density would be roughly 10 times smaller for a collapsing cloud (Federrath & Klessen 2012; Molina et al. 2012) and the wind velocity would therefore be roughly two times faster than the estimate with the uniform case. Note that the increase in the wind velocity is partially attributed to the fact that the input momentum would be 30 per cent larger in the porous case, as the background density is lower. In other words, underresolving the host cell of SN can artificially smooth out the winds, reducing the initial velocity of the neighbouring gas, especially along low-density channels.
<table><caption>**Table 1.** Summary of NUT42 simulation parameters and physical ingredients. From left to right, columns are as follows: simulation name, type of SN feedback, minimum grid size (physical units), mass of a dark matter particle, minimum mass of a star particle, threshold density for star formation, star formation efficiency per free-fall time, final redshift of each simulation, time delay of a core-collapse SN, and remarks. $\eta_{w, \max}$ is the maximum mass-loading factor with respect to the ejecta mass, used in Dubois & Teyssier (2008, DT08). $f_{w, \text{host}}$ indicates the fraction of the mass entrained from the host cell of an SN.</caption><thead><tr><th>Simulations</th><th>SN II feedback</th><th>Δx<sub>min</sub><br>(pc, physical)</th><th>m<sub>DM</sub><br>(M<sub>☉</sub>)</th><th>m<sub>star, min</sub><br>(M<sub>☉</sub>)</th><th>n<sub>th</sub><br>(H/cm<sup>3</sup>)</th><th>ε<sub>ff</sub></th><th>z<sub>end</sub></th><th>t<sub>delay</sub><br>(Myr)</th><th>Remarks</th></tr></thead><tbody><tr><td>NutCO</td><td>Metals only</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>2310</td><td>42</td><td>0.02</td><td>3</td><td>10</td><td></td></tr><tr><td>NutTFB</td><td>Thermal</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>610</td><td>42</td><td>0.02</td><td>3</td><td>10</td><td></td></tr><tr><td>NutKFB</td><td>Kinetic (DT08)</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>610</td><td>42</td><td>0.02</td><td>3</td><td>10</td><td rowspan="3">η<sub>w, max</sub> = 10</td></tr><tr><td>NutMFB</td><td>Mechanical</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>610</td><td>42</td><td>0.02</td><td>3</td><td>10</td></tr><tr><td>NutMFBm</td><td>Mechanical</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>610</td><td>42</td><td>0.02</td><td>3</td><td>Realistic</td></tr><tr><td>NutMFBmp</td><td>Mechanical</td><td>12</td><td>5.5 × 10<sup>4</sup></td><td>610</td><td>42</td><td>0.02</td><td>3</td><td>Realistic</td><td>Porous ISM (f<sub>w, host</sub> = 0.1)</td></tr></tbody></table>

In an attempt to examine the effect of the porous ISM and estimate its possible consequence on the impact of SNe, we allow for a faster wind by reducing the amount of gas mass entrained from the SN cell. We simply assume that only 10 per cent of the gas from the host cell of SN ($f_{w, \text{host}} = 0.1$) is directly involved in the expansion. This means that equation (6) reads

$$dM_{\text{swept}} = \rho_{\text{nbor}} \left(\frac{\Delta x}{2}\right)^3 + \frac{f_{w,\text{host}} \rho_{\text{host}} \Delta x^3}{N_{\text{nbor}}} + dM_{\text{ej}}. \quad (10)$$

We emphasize that this number is rather arbitrarily chosen, hence the comparison with observations should be taken with caution. Nevertheless, we note that it can clearly demonstrate the importance of detailed modelling of momentum injection and address the question as to whether the momentum budget from SN suffices to regulate star formation.

The parameters used in the six runs are summarized in Table 1.

## 3 RESULTS

In this section, we compare the baryonic properties of simulated galaxies to evaluate the effect of different SN feedback models. Our zoom-in region contains 99 dark matter haloes<fnref n="3" /> with $M_{\text{vir}} \ge 10^8$ M☉ at $z=3$, including the progenitor of an MW-like galaxy. Our analysis mostly focuses on the most massive galaxy, but we also present results for other smaller galaxies to see a general trend of the physical properties as a function of galaxy mass.

The centre of a galaxy is often taken as the densest region of a dark matter halo in large-scale simulations (e.g. Tweed et al. 2009). But it is not necessarily the case in high-resolution simulations where mergers are more frequent and a galaxy is often very clumpy. Starting from the stellar plus gas distributions within $0.5 R_{\text{vir}}$, we redefine the galaxy centre by iteratively calculating the centre of baryonic mass within its half-mass radius ($r_{\text{eff, m}}$) until converged. This ensures that the centre of a galaxy is normally taken as the centre of the most massive clump. The galaxy stellar mass is then computed by summing up the mass of the star particles within $0.2 R_{\text{vir}}$ from the centre.

The main galaxy in our simulation is hosted by a dark matter halo with $5.5 \times 10^{11}$ M☉ at $z=0$ (Kimm et al. 2011a). At $z=3$, its virial mass is $\approx 10^{11}$ M☉, which is the typical host halo mass of Lyman alpha emitters (LAE; Gawiser et al. 2007). As a comparison, LBGs are thought to reside in more massive haloes with $M_{\text{halo}} \gtrsim 10^{11.5}$

(Giavalisco & Dickinson 2001; Arnouts et al. 2002; Adelberger et al. 2005; Lee et al. 2009; Bielby et al. 2013).

Using the two Millennium simulations (Springel et al. 2005; Boylan-Kolchin et al. 2009), Fakhouri, Ma & Boylan-Kolchin (2010) measured that a halo of mass $10^{11}-10^{12}$ M☉ would typically have $\sim 2, 10$, and 20 mergers of the ratio 1:3:3, 1:33, and, 1:100 at $3 < z < 9$, respectively. The Nut halo undergoes six minor mergers (mass ratio of 4–20), and 15 very minor (mass ratio of 20–100) mergers at $3 < z < 10$, hence it has a relatively quiet merger history (Fig. 3). Note that the only major merger occurs at $z = 2.2$ (Kimm et al. 2011a). In Fig. 4, we show the projected distributions

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Halo merger tree of the Nut galaxy. The leftmost branch corresponds to the merger history of the main galaxy. The digits shown in the right-hand side of a merger event indicate a merger ratio ($\gamma = M_{\text{cen}}/M_{\text{sat}}$) below 20. The ratio for the very minor merger ($20 \le \gamma \le 100$) is omitted. Major ($\gamma \le 4$) and minor mergers are displayed as thick and thin solid lines, respectively. A bigger symbol indicates a more massive halo. Note that actual galaxy mergers take place later than the halo merger.</figcaption></figure>

<footnote n="3">³ Subhaloes are not included in this sample.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Projected distributions of the density (first and third columns), metallicity (second column), stellar density (fourth column), and SFR density (fifth column) of the Nut galaxy in runs with different feedback models, as indicated in the top-left corner of each panel. The dotted circles in the first two columns denote the virial radius of the dark matter halo. The last three columns show the central region of the halo. Note that all images are produced along the same direction (y-axis) for direct comparisons. The stellar component is markedly different from each other depending on the feedback model adopted. The spin direction of the galaxy is particularly sensitive to the choice of the feedback model. The strong outflow in the NutMFBmp run prevents the formation of a well-ordered gaseous disc at $z = 3$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of the galaxy stellar mass (left) and SFH (right) of the Nut galaxy. The abscissa indicates the age of the universe ($t_H$). Different colour-codings correspond to different feedback models, as indicated in the legend. The right-hand panel shows SFHs calculated using star particles located within 0.2 virial radius of a dark matter halo at $z=3$. Note that this includes stars formed in all progenitors. It can be seen that mechanical SN feedback is more effective at suppressing star formation than either thermal or KFB. Star formation is very episodic in the MFBmp run, as SNe efficiently destroy star-forming clouds.</figcaption></figure>

of gas density, metallicity, stellar density, and SFR density of the Nut halo and galaxy at $z=3$ in each run.

## 3.1 Star formation histories

We begin our investigation of the evolution of galaxy properties with the cooling run (NutCO) in which energy from SNe is neglected. Because there is no pressure support from stars and the UV background radiation is self-shielded above $0.01 \text{ H cc}^{-1}$ (Faucher-Giguère et al. 2010; Rosdahl & Blaizot 2012), star-forming gas clumps experience a runaway collapse converting a large fraction of the total baryons into stars at all redshifts (Fig. 5, the left-hand panel). The growth of the stellar component closely follows the assembly history of its host dark matter halo (shown as a dashed line), indicating that star formation is limited by gas accretion. The main progenitor of the Nut galaxy forms stars at a rate of $\sim 10 \text{ M}_\odot \text{ yr}^{-1}$ in the last $\sim 1.5 \text{ Gyr}$ (Fig. 5, right-hand panel), resulting in a final stellar mass of $M_{\text{star}} = 1.9 \times 10^{10} \text{ M}_\odot$ at $z=3$. As a consequence of the efficient star formation, the galaxy ends up being gas-poor at $z=3$ ($f_{\text{gas}} = M_{\text{gas}}/[M_{\text{gas}} + M_{\text{star}}] \approx 0.15$), inconsistent with the observations at $z \sim 2-3$ that galaxies in a similar mass range are typically gas rich ($f_{\text{gas}} \gtrsim 0.5$; e.g. Mannucci et al. 2009; Tacconi et al. 2010).

Fig. 6 shows the mass of the progenitor dark matter haloes in which each star particle is formed. The plot demonstrates that most of the stars populating the Nut galaxy $z=3$ in the NutCO run are formed in massive haloes. The amount of stars formed between $3 \le z \le 6$ (when the halo mass is $10.5 \le \log M_{\text{halo}}/\text{M}_\odot \le 11$) accounts for half of the total stellar mass, which reflects that more massive galaxies show a higher level of star formation activity (Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Salim et al. 2007; González et al. 2010; Whitaker et al. 2012; Salmon et al. 2015). Thus, feedback must regulate star formation at least in haloes within this mass range, if not larger, in order to be consistent with the observational findings that only a small fraction ($\lesssim 10\%$ per cent)

of baryons is turned into stars at $z \sim 3$ (Gawiser et al. 2007; Guo et al. 2010; Moster et al. 2010; Behroozi, Wechsler & Conroy 2013; Moster, Naab & White 2013).

The inclusion of thermal SN explosions (NutTFB) has a weak impact on the SFH of our intermediate-size galaxy. This is not unexpected, given that thermal energy is susceptible to severe

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The mass of dark matter haloes in which each star particle is born. The y-axis indicates the amount of stars formed in haloes of different masses. A larger amount of stars form in more massive haloes, consistent with the observed main star formation sequence. Star formation should be suppressed not only in small haloes but also in massive haloes to be consistent with the empirical results that only a small fraction of baryons turns into stars.</figcaption></figure>
artificial radiative losses. Nevertheless, the galaxy stellar mass from the NutTFB run is reduced by 31 per cent ($M_{\text{star}} = 1.3 \times 10^{10} M_{\odot}$), compared with that from the cooling run. Fig. 5 (the right-hand panel) shows that the suppression of star formation occurs mostly at an early epoch ($z \gtrsim 4$), suggesting that TFB is more efficient in dark matter haloes with a shallower potential.

There is a clear trend that star formation is more suppressed in smaller haloes with KFB (Figs 5 and 6). This is essentially because the ejecta momentum from SN explosions with a slight boost ($p \lesssim 3p_{\text{ej}}$) is significant enough to disperse small star-forming clouds present in dwarf-sized galaxies. However, as the halo becomes larger and more gas is accreted, the galaxy forms a denser and more massive gas cloud at its centre, leading to efficient star formation. The resulting stellar mass of the Nut galaxy at $z = 3$ in the NutKFB run becomes greater ($1.5 \times 10^{10} M_{\odot}$) than that in the NutTFB run.

On the other hand, the use of the appropriate momentum from the adiabatic to snowplow phase ensures that star formation is notably suppressed in the NutMFB run. The galaxy stellar mass in this run is 2.4 times smaller than that of the NutCO run at $z = 3$. It is worth noting that the mass difference is even more significant at high redshifts ($z \gtrsim 5$), indicating again that the regulation of star formation is more effective in small-mass haloes. The amount of stars formed in haloes with $M_{\text{halo}} \lesssim 10^{10.5} M_{\odot}$ is reduced by an order of magnitude (Fig. 6). However, once a giant gas clump forms ($M_{\text{core}} \sim 10^8 M_{\odot}$) as a result of mergers of massive gas clumps ($M_{\text{gas}} \sim 10^{5-6} M_{\odot}$) at the centre of the halo at $z \sim 4.8$, even momentum from SNe fail to disrupt the cloud (Fig. 7). We will discuss the cause of the inefficiency in more detail in Section 4.

In previous cosmological simulations, owing to finite resolution, assuming that each star particle of mass $10^4-10^6 M_{\odot}$ experiences a single explosion was unavoidable. It is now feasible to model discrete, multiple SNe, at least in zoom-in simulations, and to investigate their effect (NutMFBm). There are two important differences in MFBm compared to the aforementioned runs. First, since SNe start occurring as early as 3 Myr after the formation of a star particle in the NutMFMM run, the self-regulation process starts dispersing dense gas clouds earlier than in the runs with a single event (Nut-TFB, KFB, and MFB). On the other hand, as SN explosions are distributed over time, they become less disruptive in dense environments. This is precisely the motivation of using a ‘superbubble’ feedback in some studies (e.g. Sharma et al. 2014). Figs 5 and 6 indeed show that the stellar mass grows faster than the NutMFB run, indicating that the feedback is not as efficient as the single explosion case (NutMFB).

The last model (NutMFBmp) explores the possibility that SN explosions channel preferentially through a lower density medium which may not be fully resolved in our simulations. Note that the momentum input from SNe is not boosted artificially in this run, but calculated according to equation (8) with an updated average density, as a smaller amount of mass from the cell hosting the SN is now entrained ($f_{\text{host}} = 0.1$, equation 10). We find that the subtle change in modelling momentum deposition can significantly alter SFHs of the main galaxy from NutMFBmp (Fig. 5, red lines). As the outflow becomes faster along lower density channels and momentum injection is slightly increased ($\sim 30$ per cent<fnref n="4" />, if the density

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Formation of a massive stellar core due to the overcooling triggered by mergers of gas clumps in the NutMFB run at $z \approx 4.8$. The top and middle panels display the projected gas density at different redshifts. The numbers in the parentheses are the age of the universe at each redshift. The circles indicate the clumps that lead to a formation of a giant clump of $1.4 \times 10^8 M_{\odot}$ at $z = 4.71$ ($M_{\text{star}} = 8 \times 10^7 M_{\odot}$, $M_{\text{gas}} = 6 \times 10^7 M_{\odot}$) for which cooling is accelerated. The clumps with the same colours merge first into a clump before the white and yellow circles merge at $z = 4.77$. Also included in the bottom panels are the stellar density distributions at $4.71 \le z \le 4.80$. The escape velocity of the massive clump is $\sim 110 \text{ km s}^{-1}$. Note that winds with $v > 110 \text{ km s}^{-1}$ are difficult to generate inside a dense, gas-rich core via momentum input from SN. A possible solution to prevent the onset of the overcooling is by having early, efficient (SN) feedback.</figcaption></figure>

is reduced by a factor of 10), SN explosions ensure to remove gas from dense star-forming regions continuously and suppress gas accretion on to the galaxy centre (see Section 4). As a result, the galaxy stellar mass is reduced by more than an order of magnitude at $z = 3$ ($M_{\text{star}} = 7.4 \times 10^8 M_{\odot}$), compared to that of the cooling run. Notice that we are not claiming that the outflow along lower density channels would easily disrupt the massive gas clouds with $m_{\text{gas}} \gg 10^7 M_{\odot}$ which already suffer from the overcooling problem. Rather, we find that the feedback scheme is able to prevent the formation of overcooled structures in the first place through the continuous outflows and keep the SFR low.

We also note that star formation in the NutMFBmp run becomes episodic, as a small number of gas clumps govern the total SFR (Fig. 4). When these clumps are disrupted by the effective feedback, star formation drops rapidly. Such bursty star formation is also seen in Hopkins et al. (2014, fig. 9). At lower redshifts ($z < 1$) when the galaxy becomes more massive and its specific SFR drops, many star-forming clouds are likely to form in the galactic disc (see fig. 1 in Hopkins et al. 2014, for example) and SFHs become smoother.

It is worth mentioning that the role of the outflow along lower density channels in regulating star formation may vary under certain circumstances. For example, using a set of idealized simulations with an initially turbulent medium, Hennebelle & Iffrig (2014) show that SNe exploding in lower density environments enhance star formation (their model D) than the case where SNe always explode in a dense region (their model C1, fig. 13). They argue that this is

<footnote n="4">⁴ Note that only a small fraction of cells (10 per cent) would have 30 per cent increase in the SN momentum by assuming that the volume-filling density is 10 times smaller than the mean density of the host cell of SN (see Section 4). The total momentum input from SNe in the galaxy is increased only by $\sim 4$ per cent.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Stellar mass fraction of galaxies in simulations with different feedback models at $z = 3$. Note that galaxies that are not accreted on to the main halo are also included. Different symbols and colours correspond to different feedback schemes. Empirical determinations based on the abundance matching technique are shown as different lines, as indicated in the legend (Behroozi et al. 2013; Moster et al. 2013). The shaded regions denote $1\sigma$ and $2\sigma$ uncertainties for the Moster et al. (2013) and Behroozi et al. (2013) results, respectively. Independent measurements from the clustering analysis of $z \sim 3$ star-forming galaxies selected from the VIMOS ultra deep survey (Durkalec et al. 2015) and Ly$\alpha$ emitters at $z \sim 3.1$ (Gawiser et al. 2007) are shown as a blue pentagon and a cyan diamond symbol, respectively.</figcaption></figure>

likely because most of the SN energy is carried away with diffuse gas, compressing the neighbouring dense clouds (model D). This may sound contradictory to our findings, but we note that the different conclusions may be reached due to different initial conditions and time-scales involved. First, while dense clouds develop simultaneously over the entire computational domain in Hennebelle & Iffrig (2014), gas is accreted smoothly through the cosmic web over several gigayears in our simulations. Thus, the diffuse gas accelerated by SN explosions is less likely to encounter the collapsing clouds and trigger star formation in our case. Secondly, although the diffuse hot gas may not destroy nearby star-forming clouds, it can reduce the gas cooling and accretion on to the galaxy centre by interacting with halo gas, leading to a lower star formation activity on cosmological time-scales.

In Fig. 8, we examine the stellar mass fraction of the galaxies ($M_{\text{star}}/f_{\text{bar}}M_{\text{halo}}$) in our simulations. Note that galaxies outside the virial radius of the main halo are also included. For comparison, we include the empirical determinations of the stellar mass fraction based on the abundance matching technique (Guo et al. 2010; Moster et al. 2010; Behroozi et al. 2013; Moster et al. 2013). Independent measurements based on the clustering analysis of 3022 star-forming galaxies selected from the VIMOS Ultra Deep Survey (Durkalec et al. 2015, $z \sim 3$) and 162 Ly$\alpha$ emitters at $z \sim 3.1$ (Gawiser et al. 2007) are also included as a blue pentagon and a cyan diamond with error bars, respectively. We stress that the comparison between the empirical results and our simulations should be taken with caution, given that the simulated galaxy sample is quite small.

Nevertheless, several interesting features can be gleaned from this plot. First, a large fraction of baryons is converted into stars in massive haloes in the cooling run, while the fraction is small in haloes of a few times $10^8 M_{\odot}$, as the UV background heating prevents gas from collapsing. We confirm that a similar trend is found when the fraction of total baryonic mass fraction is plotted instead of the stellar mass fraction. Secondly, there is a general trend that the stellar mass fraction is lower in smaller haloes in the runs with feedback. Thirdly, the use of more realistic time delays does reduce stellar mass in the small haloes ($M_{\text{halo}} \lesssim 10^{9.5} M_{\odot}$). Fourthly, the stellar fraction in small haloes agrees well between MFBmp and MFBmp. This means that the fraction of SN host cell mass entrained is irrelevant, as long as the momentum injection from SNe is high enough to unbind the star-forming cloud. Finally, the MFBmp run produces an amount of stars consistent with the empirical sequence, indicating that momentum input from SNe can provide enough pressure support to regulate star formation. Whether this momentum can be efficiently transferred to the turbulent medium, as postulated in the MFBmp run, remains to be investigated, however. It can be seen that star formation in galaxies embedded in lower-mass haloes is more suppressed in general. Mechanical SN feedback is more efficient than other feedback schemes at regulating star formation. This plot demonstrates that momentum input from SNe alone may be able to match the empirical sequence (the NutMFBmp run), provided that the ISM is highly turbulent.

## 3.2 Metal enrichment

It is well established that the gas-phase metallicity in galaxies increases with galaxy luminosity or mass (Garnett 2002; Pilyugin, Vílchez & Contini 2004; Tremonti et al. 2004; Kewley & Ellison 2008). The MZR is also observed in galaxies at higher redshifts (Erb et al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Yuan, Kewley & Richard 2013; Zahid et al. 2013; Steidel et al. 2014; Wuys et al. 2014), but a systematic offset is found such that more distant galaxies are more metal-poor for a given mass. This suggests that metals are not instantly recycled to enrich the ISM, but possibly blown away through galactic winds (Dalcanton 2007).

Stellar metallicity is an important indicator of the overcooling problem. Wise et al. (2012a) show that the metallicity of stars is very sensitive to the strength of stellar feedback. We also confirm that the metal enrichment is directly affected by a galaxy's ability to regulate star formation (Fig. 9). In the cooling run, the stellar metallicity quickly approaches the solar value ($Z = 0.02$), as newly synthesized metals are recycled instantaneously within their birth clouds. Visual inspection of metallicity distributions within the dark matter halo shows that the intergalactic medium (IGM) is barely metal-enriched (Fig. 4, the top second panel). This demonstrates that galaxy interactions are not an efficient process to re-distribute metals in this intermediate mass haloes. The galaxies in runs with SN feedback start with a lower metallicity when the halo mass is small, but again most of them approach solar metallicity as soon as a massive stellar core which feedback cannot destroy forms in the galaxy (see Fig. 7). The only exception is the NutMFBmp run in which strong outflows prevent the formation of such a massive core by continuously dispersing gas clumps (Fig. 4).

Fig. 10 compares the MZR of the simulated galaxies from different feedback runs at $z = 3$. We compute the gas-phase metallicity by assuming the solar abundance $[(O/H)_{\odot} = 8.69$; Asplund et al. 2009. Our minimal run (NutCO) shows a steep MZR in small galaxies ($M_{\text{star}} \lesssim 10^7 M_{\odot}$), as star formation is mitigated by the strong UV background radiation. The MZR becomes nearly
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of stellar metallicity in the main galaxy. Different colour-codings indicate the runs with different feedback models. Also included as black or red dotted lines are the corresponding gas metallicities. Stellar metallicity is a good indicator of the overcooling problem, as demonstrated by Wise et al. (2012a). We find that the stellar metallicity increases rapidly as soon as a dense core forms at the galaxy centre (see also Fig. 7). Only the efficient feedback model (MFBmp) forms stellar populations with a significantly subsolar metallicity.</figcaption></figure>

flat, however, for more massive galaxies, again due to efficient recycling of metals, in contradiction to observations. The MZR has a slope close to the observed one only when SN feedback is included.

The normalization of the MZR in simulations is directly affected by the strength of SN feedback. While observations at high redshifts ($z \sim 3$) show that galaxies with $10^9 \lesssim M_{\text{star}} \lesssim 10^{9.5} M_{\odot}$ are metal-poor ($\sim 0.1-0.3 Z_{\text{sun}}$) (Maiolino et al. 2008; Mannucci et al. 2009), simulations with weak SN feedback (thermal and kinetic) predict metal-rich dwarf galaxy populations with $\sim 0.3-1.0 Z_{\text{sun}}$. Only the galaxies with star formation regulated by strong SN feedback (e.g. MFBmp) reasonably follow the observed sequence.

It is worth noting that the simulated main galaxies with different feedback, except the one from the MFBmp run, show similar gas-phase metallicities at $z = 3$. This is essentially because the gas metallicity saturates at around the solar value. Although the metallicity of the ejecta adopted in this study is $2.5 Z_{\text{sun}}$, the pristine gas is accreted continuously on to the galaxy at high redshifts, diluting the enriched gas. Furthermore, since stellar mass increases together with the metallicity (see also Agertz & Kravtsov 2015, fig. 6), the gas-phase metallicities of the main galaxies seem to agree with the observational data within the current measurement uncertainties. Given the comparison of stellar mass fraction, this means that the gas-phase metallicity of the massive galaxies ($M_{\text{star}} \gtrsim 10^{10} M_{\odot}$) may not be a good indicator for the overcooling problem, although dwarf galaxies with masses $10^8 \lesssim M_{\text{star}} \lesssim 10^9 M_{\odot}$ will still be very useful to put a strong constraint on galaxy formation models.

The origin of the MZR is often attributed to the fact that star formation efficiency is lower in smaller galaxies. Some authors argue that outflows may not be necessary to explain the MZR as long as the stellar fraction is a strong function of galaxy mass (Finlator & Davé 2008; Calura et al. 2009; Mannucci et al. 2009). However,

<figure><img src="image_11.png" /><figcaption>**Figure 10.** The MZR at $z = 3$ in the Nut simulations. Different symbols and colours correspond to the runs with different feedback models. For comparison, we include the observed MZRs drawn from nine galaxies at $z \sim 3.5$ (Maiolino et al. 2008) and eight LBGs at $z \sim 3.1$ (Mannucci et al. 2009). The solid line is the fit to the results from Maiolino et al. (2008), while the local MZR derived from the local SDSS galaxies is shown as a dashed line (Kewley & Ellison 2008). A fraction of simulated galaxies exhibits very enhanced metallicity compared to their mean MZR sequence, because SNe blow out star-forming gas temporarily and the gas metallicity is governed by the enriched outflow. Note that without outflows (NutCO) the slope of the MZR is shallower than the observations. The effectiveness of SN feedback determines the offset from the observed MZR at $z \sim 3$.</figcaption></figure>

our simulation (NutCO) suggests that the fraction is always high and not a strong function of galaxy (or halo) mass without feedback (NutCO, see Fig. 8). This is essentially because gas experiences a runaway collapse, shifting the density distribution of star-forming gas towards a very dense regime. Note that the typical density of star-forming clouds without feedback is set by the force resolution of numerical simulations, and we find that this is roughly $n_{\text{H}} \sim 10^{3.5} \text{ cm}^{-3}$ in our runs with 12 pc resolution (see also fig. 12 of Hopkins et al. 2012a for higher resolution runs where the typical density is larger $n_{\text{H}} \gtrsim 10^4 \text{ cm}^{-3}$). Thus, for a star formation efficiency of $\epsilon_{\text{eff}} = 0.02$, one would expect that more than 80 per cent of the gas cloud is converted into stars in 100 Myr. Since the Hubble time is much longer than this for $z = 3$ galaxies, metallicity evolves similarly as a closed box model in the absence of feedback, except that the ‘equilibrium’ metallicity is smaller than the metallicity of the stellar ejecta ($Z_{\text{ej}} = 0.05$) due to the infall of pristine gas. Finlator & Davé (2008) already pointed out that in the absence of outflows the capability of reproducing the MZR depends on the gas consumption time-scale (see their fig. 13).

### 3.3 Galaxy size, circular velocity, and morphology

An important prediction from tidal torque theory is that late accretion carries larger specific angular momentum ($j$) than the gas inflow at early times (Peebles 1969; Doroshkevich 1970; White 1984; Catelan & Theuns 1996). Because of this, the spin parameter of a dark matter halo ($\lambda' = j_{\text{dm}}/\sqrt{2}R_{\text{vir}}V_c$; Bullock et al. 2001) can be kept nearly constant across different redshifts (van den Bosch et al. 2002; Peirani, Mohayaee & de Freitas Pacheco 2004; Bett
<figure><img src="image_12.png" /></figure>

**Figure 11.** Evolution of the half-mass radius ($R_{\text{eff,m}}$) of the main galaxy at $3 \le z \le 15$ ($\Delta z \approx 0.5$). The galaxy at $z=3$ is shown as a bigger symbol than others. Different colour-codings and symbols are the same as in Fig. 10. Galaxies grow with time, as late accretion supplies gas with larger specific angular momentum. Strong feedback tends to produce a more extended stellar component. Note that the half-mass radii can also increase temporarily due to galaxy mergers. The late increase in radius at late times in the run with TFB is also the outcome of serendipitous mergers (see the text).

et al. 2007). Numerical simulations indeed confirm that $j$ of dark matter near the virial radius is a factor of 2–3 larger than the net $j$ of the halo (Kimm et al. 2011b; Pichon et al. 2011; Stewart et al. 2013; Danovich et al. 2015).

Baryons gain angular momentum by the same large-scale tidal torques as dark matter (e.g. Kimm et al. 2011a), and thus galactic discs are expected to grow by late gas inflows. Fig. 11 indeed shows that the half-mass radius of the stars ($R_{\text{eff,m}}$) in the main galaxy gets larger with time regardless of the feedback model. One can also see that the strength of SN feedback generally leads to a larger $R_{\text{eff,m}}$. The NutMFB and NutMFBmp runs produce galaxies with $R_{\text{eff,m}} = 0.83$ and 1.5 kpc, while the half-mass radii of the main galaxy in NutCO, NutKFB, and NutMFBm runs are 0.43, 0.40, 0.31 kpc, respectively.

An exception is the TFB run, which produces a rather extended stellar component with $R_{\text{eff,m}} = 0.93$ kpc. We find that the relatively large half-mass radius originates from the serendipitous mergers at $5 < z < 6$. In this run, the merging satellite galaxies are more massive than those in other runs with feedback (Fig. 8), and we find that the orbital angular momenta of them are coincidently aligned with the spin of the galactic disc. Consequently, the merger remnant becomes a more extended clumpy disc (Fig. 4, fourth column), even though there is still a significant central concentration of mass at its centre. However, it is very unlikely that such serendipitous mergers occur only in certain feedback models. Although we have not made dwarves a focus of this study, we confirm that dwarf galaxies with $10^7 \lesssim M_{\text{star}} \lesssim 10^9 M_{\odot}$ from the NutTFB run are in fact smaller ($R_{\text{eff,m}} \sim 100$ pc) than the runs with MFB ($R_{\text{eff,m}} \sim 500$ pc).

It is also interesting to note that the galaxy size is closely linked with the SFH and metal enrichment in stars. The most evident example can be seen from the NutMFB run at $4.5 \lesssim z \lesssim 5.5$. This corresponds to the epoch during which stellar mass and metallicity

<figure><img src="image_13.png" /></figure>

**Figure 12.** Circular velocities of the main galaxy from a suite of Nut simulations. The circular velocity is measured as $V_c(<r) = \sqrt{GM(<r)/r}$, where $M(<r)$ is the total mass (baryons+dark matter) inside radius $r$. Different lines correspond to runs with different feedback models, as indicated in the legend. As known in the literature, the runs with weak SN feedback show a peaked velocity curve towards the galaxy centre. The run with MFB in a porous ISM (MFBmp) suppresses the central peak in the velocity curve.

increase rapidly indicating the onset of severe runaway collapse and overcooling (Fig. 7). Accordingly, the galaxy gets significantly smaller during this period ($8.5 \le \log M_{\text{star}}/M_{\odot} \le 9.5$).

Another useful galaxy quantity is the circular velocity ($V_c(<r) = \sqrt{GM(<r)/r}$). Observed bright spiral galaxies are characterized by the smoothly varying baryonic rotation curves (e.g. Kassin, de Jong & Weiner 2006). However, Fig. 12 shows that the maximum circular velocity appears at very small radii ($r \sim 100$ pc) in the runs in which the gas-to-stellar conversion efficiency of the main galaxy is higher than the empirical estimates (CO, TFB, KFB, MFB, and MFBmp), confirming earlier results (Navarro & White 1993; Abadi et al. 2003; Governato et al. 2007; Scannapieco et al. 2009; Hummels & Bryan 2012; Agertz & Kravtsov 2015). The peak near the centre sets in earlier and the maximum circular velocity is larger in simulations with less effective feedback. The velocity curve in the NutMFB run is smoothly rising at $z=7$, but the peak at the centre begins to emerge at $z \lesssim 5$ again due to the runaway gas collapse (Fig. 7). Only the run with the mechanical SN feedback based on the assumption of a porous ISM (NutMFBmp) does not show the peaked circular velocity curve. Note, however, that the rotation curve of the galaxy turns out to be very different from that of the circular velocity, indicating that there is no stellar disc with a well-ordered motion (see Fig. 4 for the stellar density plot in the fourth column).
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Top: fraction of stars supported by rotation ($f_{\text{rot}}$) in the main galaxy as a function of the Hubble time. We simply assume that a star is rotationally supported if more than half of its total kinetic energy is accounted for by rotational motion. The different colours and symbols correspond to different feedback models, as indicated in the legend. Bottom: $f_{\text{rot}}$ of the stars younger than 10 Myr. The main galaxy in the runs having more stars than the empirical estimates (CO, TFB, KFB, MFB, and MFBm) is marginally rotation-supported. Interestingly, we find that the most effective feedback (MFBmp) leads to the formation of a spheroidal galaxy.</figcaption></figure>

In order to see the kinematic properties of the simulated galaxies, we quantify the fraction of stars supported by rotation ($f_{\text{rot}}$) in Fig. 13. We first determine the spin direction of the total stellar component within 0.2 virial radius and the rotation velocity ($v_{\text{rot}}$) of each star particle along the azimuthal direction. Then the mass fraction of stars whose kinetic energy ($E_{\text{kin},*} = mv^2/2$) is dominated by the rotational component ($E_{\text{rot},*}/E_{\text{kin},*} > 0.5$) is calculated within the inner radius ($r \le 0.2R_{\text{vir}}$), where $E_{\text{rot},*} = mv_{\text{rot}}^2/2$. $f_{\text{rot}}$ may be viewed as a proxy for the galaxy morphology in the sense that a purely rotating disc would have $f_{\text{rot}} = 1$, while a stellar system with only radial motions would have $f_{\text{rot}} = 0$. Fig. 13 (top panel) shows that the simulated galaxies with a high stellar-to-halo mass ratio are marginally supported by rotation. Approximately, 55–70 per cent of the stars have orbits that are more or less aligned ($\theta < 45^\circ$) with the spin axis of the galaxy, giving rise to geometrically thick discs (Fig. 4, fourth columns)<fnref n="5" />.

On the other hand, a large fraction of stars (~80 per cent) in the main galaxy from MFBmp is supported by random motions, resulting in a spheroidal morphology rather than a disc (Fig. 4). This

seems to be in conflict with the previous claim that strong feedback preferentially blows out gas with low angular momentum (e.g. Governato et al. 2007; Brook et al. 2011), but consistent with recent high-resolution simulations by Roškar et al. (2014) and Hopkins et al. (2014). Using the same grid-based code as ours (RAMSES), Roškar et al. (2014) performed idealized simulations varying the optical depth to IR photons to investigate the impact of radiation pressure on galaxy properties. They find that the model with strong stellar feedback, which matches the empirical stellar-to-halo mass relation, leads to the formation of a spheroidal galaxy. Hopkins et al. (2014) perform cosmological simulations to study the formation of disc galaxies with a sophisticated stellar feedback model, including radiation pressure and SN feedback. They show that the morphology of the MW progenitor at $z = 3.4$ is irregular due to strong outflows (see their fig. 1). Somewhat encouragingly, despite the absence of the well-ordered disc at $z \sim 3$, they find that smooth gas accretion at later times leads to the formation of a realistic, extended disc. In this regard, the spheroidal morphology we obtain from the MFBmp run may not be problematic, but we will need to run the simulation down to $z = 0$ to confirm this.

Although a significant fraction (~30 per cent) of stars has orbits misaligned with the overall spin of the galaxy, the stars younger than 10 Myr form out of a well-organized, gaseous disc ($f_{\text{rot}} \gtrsim 0.8$) in the CO, TFB, KFB, and MFBm runs. This can also be seen in the projected SFR density distributions (Fig. 4, fifth columns). One may wonder why $f_{\text{rot}}$ of the whole stellar population does not increase notably for $3 \le z \lesssim 6$, given that $f_{\text{rot}}$ of the young component is high (~0.9). This can be attributed to two factors. First, the young stars form in the core of the galaxy (~300 pc) where there exists a massive bulge of ~2–4 × $10^9 M_{\odot}$. Since the motion of the young stars is mainly governed by the gravity of the bulge, the interaction with the spherical structure is likely to isotropize their motions. More importantly, the kinematically cold young stars are susceptible to external disturbances, such as harassment or mergers. In order to check this, we follow the evolution of $f_{\text{rot}}$ of the young stars ($t \le 10$ Myr) formed at different redshifts ($z = 4.5–4.8$, $\Delta z = 0.05$) in the NutTFB galaxy during which the spin vector of the (young) stellar component is nearly unchanged (~15°), while their $f_{\text{rot}}$ drops to ~0.5. We find that the evolutionary patterns of $f_{\text{rot}}$ in the young population with different ages are well synchronized, following the merging histories of star clusters and/or galaxies, irrespective of their age. This means that the young stars are continuously transformed into a bulge by the dynamical disturbances.

We find that star-forming gas shows irregular morphologies when SN feedback is strong (see Fig. 4, MFBmp) or when the galaxy is small ($z \gtrsim 7$ in the MFB runs). This is essentially because the feedback keeps blowing out/away a substantial amount of gas before a kinematically cold gaseous disc forms. Let us suppose that there exists a well-defined, clumpy galactic disc in the x-y plane. If star formation is very active, it would drive strong outflows along the z direction destroying the star-forming clouds. Even though the gas is initially coplanar, the spin direction of the returning gas would no longer be aligned with the galactic disc, randomizing the direction of the spin direction of the subsequently collapsing star-forming clouds. This effect is likely to be more significant if bursty star formation takes place in a small number of gas clouds, as in the MFBmp run (Fig. 4). However, if a gas-rich disc becomes massive enough to behave like a buffer region in which the returning galactic fountain gas is first mixed with the disc gas or other fountain gas, star formation in the well-defined disc is likely to continue, increasing the fraction of kinematically settled discs (Kassin et al. 2012).

<footnote n="5">⁵ Note that we do not reorient the galaxy to the edge-on configuration in Fig. 4 to demonstrate that the spin direction of the simulated galaxies may be significantly altered by the feedback model.</footnote>
## 3.4 Galactic outflows and inflows

We now turn our attention to the properties of gas outside the main galaxy. The inflow or outflow rate is measured by computing the flux at radius $r$ as

$$ \dot{M}_{\text{inf}}(r) = \int d\Omega 4\pi r^2 \rho_{\text{gas}}(r; \Omega) v_{\text{rad}}(r; \Omega) \Theta(-v_{\text{rad}}), $$

$$ \dot{M}_{\text{out}}(r) = \int d\Omega 4\pi r^2 \rho_{\text{gas}}(r; \Omega) v_{\text{rad}}(r; \Omega) \Theta(v_{\text{rad}}), $$

where $\Omega$ is the solid angle, $v_{\text{rad}}$ is the radial velocity with positive meaning outflows, and $\Theta$ is the Heaviside step function. Since we are mainly interested in the effects of SN explosions by the central galaxy, we exclude the contribution from the satellites by ignoring the gas within their tidal radius.

We find that the cosmic gas inflow is impervious to SN explosions (Fig. 14, panel a). Although strong SN feedback (MFBmp) suppresses star formation by an order of magnitude, it cannot completely destroy the dense gas filaments around the virial radius (see Fig. 4, first columns), and both the inflow rate and the net infalling velocity are only mildly reduced (panel d). Other runs with weak feedback show a comparable inflow rate as in the NutCO run. This suggests that, if a significant fraction of baryons is missing from the host dark matter halo of an MW-like galaxy (McGaugh et al. 2010, cf. Werk et al. 2014), it is unlikely to be caused by the suppression of smooth gas accretion, but rather by stellar feedback blowing gas out of the halo.

The more effective a feedback model is at reducing star formation, the larger the outflow rate it produces (panel b). While the runs with TFB and KFB blow gas out of the virial sphere at a rate of $1-2 M_{\odot} \text{yr}^{-1}$, the MFB runs show more significant outflows of $1-5 M_{\odot} \text{yr}^{-1}$ (MFB and MFBm) or $5-10 M_{\odot} \text{yr}^{-1}$ (MFBmp). The outflow rates measured at smaller radii ($0.5 R_{\text{vir}}$) tend to be larger than those at the virial radii (panel c), indicating that the outflow is launched more or less ballistically. Although we do not include the figure in this paper, we find that a large fraction ($\sim 60\%$ per cent) of baryons escapes from the halo in the MFBmp run, resulting in $(M_{\text{star}} + M_{\text{gas}})/M_{\text{vir}} = 0.071$ at $z=3$ (see Table 2). Note that this is a factor of three lower than that in the cooling run ($(M_{\text{star}} + M_{\text{gas}})/M_{\text{vir}} = 0.238)^6$.

Due to the enhanced outflow and reduced star formation in the MFBmp run, a high mass-loading factor ($\eta = \dot{M}_{\text{out}}/\dot{M}_{\text{star}}$) of $5-10$ is observed for galaxy stellar masses of $10^8 \lesssim M_{\text{star}} \lesssim 10^9 M_{\odot}$ (Fig. 15). This is roughly two orders of magnitude larger than those of weaker feedback runs ($\eta \sim 0.1$, TFB or KFB). Here, we use the outflow rate at the virial radius, and thus the mass-loading factor could be even higher if it is measured at smaller radii. We note that such a high wind efficiency in dwarf-sized galaxies is required in large-scale cosmological simulations to match the stellar mass

functions (Oppenheimer et al. 2010). Hopkins et al. (2012b) also show using idealized simulations with parsec-scale resolutions that the mass loading is $\eta \sim 10$ for their Sbc galaxy where the maximum circular velocity of the dark matter halo ($86 \text{ km s}^{-1}$) is comparable to our main galaxy at $3 \lesssim z \lesssim 6$. The comparison is by no means sufficient to demonstrate the success of the models, however, and future works will have to directly compare to observable signatures of the outflows, such as low-ionization absorption lines (Martin 2005; Rupke et al. 2005).

We also find that the physical properties of the outflowing gas are generally similar in all feedback runs. Fig. 14 (panels e, f, and i) shows the flux-weighted velocity and hydrogen number density of the outflow in different runs. The galactic winds (panels e, f) are faster ($150-300 \text{ km s}^{-1}$) than the inflow (panel d), and their density at $0.5 R_{\text{vir}}$ is comparable to the mean baryonic density of a virialized halo ($\approx 178 \rho_c(z) \Omega_b / \Omega_m$), regardless of the feedback model used (panel i), where $\rho_c(z)$ is the critical density of the universe. This suggests that it is the opening angle of the outflow that determines the overall outflow rate. To substantiate this, we compare the distribution of the outflowing gas between NutTFB and NutMFBmp in Fig. 16. The plot shows the HEALPIX map (Górski et al. 2005) of hydrogen number density on a sphere of $0.5 R_{\text{vir}}$ at $z=3.9$ when the outflow velocity of the two runs is very similar. For clarity, both the infalling gas and the gas belonging to satellites are removed. It is evident from the figure that the outflow rate in the MFBmp run is larger, as it subtends a larger solid angle than the one in the TFB.

Note that the typical density shown in panel (i) (Fig. 14) is sensitive to the radius at which it is measured. For example, it can increase up to $0.1-1 \text{ cm}^{-3}$ in the MFBmp run at $0.3 R_{\text{vir}}$, as small gas clumps that are blown away from the galaxy can contribute to the outflows (Fig. 16, the bottom-right panel). These gas clumps with $n_H > 0.01 \text{ cm}^{-3}$ move outwards more slowly ($\sim 40 \text{ km s}^{-1}$) than the mean outflow velocity ($\sim 90 \text{ km s}^{-1}$) at $0.3 R_{\text{vir}}$, and thus we expect that many of these clumps do not escape from the halo (see also Cooper et al. 2008). On the other hand, the runs with weak feedback (NutTFB or NutKFB) cannot form many outflowing clumps, as SN explosions cannot transfer the correct amount of momentum to the ISM.

Fig. 14 (panel h) shows that the outflowing gas is metal-rich, but generally subsolar at $z \ge 3$ in the main halo. This means that SN ejecta entrains a considerable amount of gas, which is at least $1.5$ times the ejecta mass, even in the weak feedback runs (TFB and KFB), given that the assumed metallicity of the SN ejecta is $2.5$ solar. Thus, a small mass-loading factor of $\eta \sim 0.1$ seen in Fig. 15 arises simply because not all SN explosions drive large-scale winds. Because of the significant entrainment, the metallicities of the outflow are comparable to the gas-phase metallicity of the main galaxy (Fig. 10). For example, the main galaxy in the MFBmp run with a mass-loading factor of $\eta \sim 10$ reveals a significantly smaller metallicity of $\sim 0.1 Z_{\odot}$ than the ejecta metallicity. The cooling run displays somewhat complex behaviour, as the outflow in this case is driven by mergers and interactions with satellite galaxies.

As accretion is dominated by the cold filamentary gas in our main halo at $z \ge 3$ (Kimm et al. 2011a), the metallicity of the inflow is found to be an order of magnitude smaller ($0.005-0.03 Z_{\odot}$) than that of the outflow. The different metallicity between the inflow and outflow demonstrates that the outflow is not efficiently mixed with the cold flows (see Fig. 4, second columns). Indeed, the metallicities turn out to be different even for the MFBmp run where the outflow rate is comparable to the infall rate. Recently, Crighton, Hennawi & Prochaska (2013) analyse the absorption features of a $z = 2.44$

<footnote>$^6$ Faucher-Giguère, Kereš & Ma (2011) also found a similar feature in their SPH simulation without feedback, and attributed this to the delayed gas infall due to photoionization by the UV background. However, we find that the main halo in the NutCO run shows the high fraction even before the UV background is instantaneously turned on at $z=8.5$. We argue that the high baryon fraction is due to the fact that baryons and dark matter particles settle into a dark matter halo in a different way (Kimm et al. 2011a). Whereas gas agglomerates at the centre of a halo once accreted, collision-less dark matter particles pass through the central region first and populate the outskirts or regions beyond the virial radius of the halo depending on their angular momentum. Thus, it may not be surprising that the baryon fraction is higher than $f_b$ if one measures the fraction in the most baryon-rich region (i.e. halo centre).</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Properties of gas inflow and outflow around the main galaxy in the simulations. Different feedback models are shown as different colour-codings and line styles. The top panels (a–c) display the accretion ($dM_{\text{inf}}/dt$) or outflow rate ($dM_{\text{out}}/dt$) measured at the virial or half of the virial radius. The middle (d–f) and bottom panels (g–h) show the flux-weighted radial velocity ($V_r$) and metallicity ($Z_{\text{gas}}$) of the inflow and outflow, respectively. Note that we exclude the gas within the tidal radius of satellite galaxies, hence the properties represent the smooth component. Also included in the last panel (i) is the flux-weighted hydrogen number density of inflow (dashed lines) or outflow (solid lines) at the virial radius or $0.5 R_{\text{vir}}$, respectively (see the legend). For comparison, the mean baryonic density of a virialized halo ($\approx 178 \rho_c(z) \Omega_2 / \Omega_\text{m}$) is shown as a grey line, where $\rho_c(z)$ is the critical density of the universe. We find that the gas accretion rate is more or less similar regardless of the feedback model, indicating that the cold gas filaments are impervious to SN feedback (see Fig. 4, first columns). The outflow is stronger in the runs in which star formation is more suppressed. However, the flux-weighted velocity and number density of the outflowing component at $0.5 R_{\text{vir}}$ is largely indistinguishable among the different feedback models. The outflow seen in the run without feedback (NutCO) is not feedback-driven, but induced by mergers. Finally, we find that the inflowing gas is roughly 10 times more metal-poor than outflows.</figcaption></figure>

system with the background quasar at $z = 2.66$, and find a cool ($T < 20\,000\,\text{K}$) component with $0.007–0.015$ solar metallicity. Their neutral column density ($N_H = 10^{19.50 \pm 0.16} \text{ cm}^{-2}$) as well as the low metallicity are surprisingly similar with the properties of the inflow in our simulations, supporting their claim that it is a direct detection of the cold flow around an LBG.

# 4 DISCUSSION

## 4.1 Overcooling problem and a viable solution

In the previous section, we have shown that a large amount of gas is converted into stars at the centre of a galaxy if the feedback
<table><caption>**Table 2.** Simulated properties of the most massive galaxy embedded in a $1.2 \times 10^{11}$ M$_{\odot}$ halo at $z = 3$. From left to right, each column corresponds to the model name, the total stellar mass of a galaxy inside $0.2 R_{\text{vir}}$ of a dark matter halo ($M_{\text{star}}$), SFR averaged over the last 10 Myr ($\dot{M}_{\text{star}}$), the gas-to-stellar mass ratio ($M_{\text{gas}}/M_{\text{star}}$) within the virial radius, the stellar-to-halo mass ratio ($M_{\text{star}}/M_{\text{halo}}$) within the virial radius, the fraction of baryons inside a virial radius ($M_{\text{bar}}/M_{\text{halo}}$), the mean metallicity of stars ($Z_{\text{star}}$), the mean metallicity of gas ($Z_{\text{gas}}$) within $0.2 R_{\text{vir}}$, outflow rate measured at $0.5 R_{\text{vir}}$ ($\dot{M}_{\text{out}}$), and the half-mass radius of the main galaxy.</caption><thead><tr><th>Model</th><th><i>M</i><sub>star</sub><br>(M<sub>☉</sub>)</th><th><i>M</i>&#775;<sub>star</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th><i>M</i><sub>gas</sub><sup>vir</sup>/<i>M</i><sub>star</sub><sup>vir</sup></th><th><i>M</i><sub>star</sub><sup>vir</sup>/<i>M</i><sub>halo</sub></th><th><i>M</i><sub>bar</sub><sup>vir</sup>/<i>M</i><sub>halo</sub></th><th><i>Z</i><sub>star</sub><br>(Z<sub>☉</sub>)</th><th><i>Z</i><sub>gas</sub><br>(Z<sub>☉</sub>)</th><th><i>M</i>&#775;<sub>out</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th><i>r</i><sub>eff, m</sub><br>(kpc)</th></tr></thead><tbody><tr><td>NutCO</td><td>1.9 &times; 10<sup>10</sup></td><td>6.2</td><td>0.45</td><td>0.147</td><td>0.238</td><td>1.03</td><td>0.65</td><td>0.02</td><td>0.43</td></tr><tr><td>NutTFB</td><td>1.3 &times; 10<sup>10</sup></td><td>4.3</td><td>0.69</td><td>0.096</td><td>0.184</td><td>0.95</td><td>0.63</td><td>0.50</td><td>0.93</td></tr><tr><td>NutKFB</td><td>1.5 &times; 10<sup>10</sup></td><td>5.7</td><td>0.58</td><td>0.115</td><td>0.195</td><td>0.98</td><td>0.50</td><td>0.03</td><td>0.40</td></tr><tr><td>NutMFB</td><td>7.8 &times; 10<sup>9</sup></td><td>5.5</td><td>1.31</td><td>0.061</td><td>0.144</td><td>0.62</td><td>0.58</td><td>0.67</td><td>0.81</td></tr><tr><td>NutMFBm</td><td>9.8 &times; 10<sup>9</sup></td><td>2.4</td><td>1.14</td><td>0.062</td><td>0.148</td><td>0.82</td><td>0.30</td><td>0.38</td><td>0.31</td></tr><tr><td>NutMFBmp</td><td>7.4 &times; 10<sup>8</sup></td><td>0.6</td><td>10.2</td><td>0.006</td><td>0.071</td><td>0.06</td><td>0.07</td><td>1.20</td><td>1.5</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** The mass-loading efficiency of winds in the main galaxy at redshifts $3 \le z \le 12$. The loading factor is computed by dividing the outflow rate measured at the virial radius with the SFR averaged over 10 Myr. The strong SN feedback model (MFBmp) leads to a high $\dot{M}_{\text{out}}$/SFR of 5–10 when the galaxy stellar mass is greater than $10^8$ M$_{\odot}$. Such a high mass loading is normally required to reproduce the stellar mass functions in large-scale hydrodynamic simulations (e.g. Oppenheimer et al. 2010).</figcaption></figure>

cannot regulate star formation effectively. The overcooling problem is a well-known issue and reported in many earlier works based on thermal SN feedback (Katz 1992) or a KFB model (Dubois et al. 2014). In this work, we demonstrated that the use of the momentum available at the end of the adiabatic phase (NutMFB and NutMFBm) cannot eliminate the overcooling problem, even though it better suppresses star formation in small galaxies. More importantly, we find that SFHs of galaxies predicted from hydrodynamic simulations are very sensitive to the details of how momentum from SNe is distributed (see also Hennebelle & Iffrig 2014).

In order to understand why some feedback models are more effective at blowing gas away than others, we compute the probability distribution functions (PDFs) of the density at which stars primarily form ($n_{\text{H, form}}$, top panel) and the density at which SNe explode ($n_{\text{H, SN}}$, middle panel) in Fig. 17. The PDFs of $n_{\text{H, form}}$ are obtained by computing the actual probability ($\propto M_{\text{gas}}/t_{\text{ff}}$) of star formation

in the ISM of the main galaxy from $\sim 50$ snapshots at $3 \le z \le 10$. The sites of SN explosions are chosen as regions where young stars with the age $9 \le t \le 10$ Myr are located for the TFB, MFB, and CO runs, or by taking into account the realistic SN rates based on the metallicity and age of each star in the case of the MFBm and MF-Bmp runs at $3 \le z \le 10$ (Leitherer et al. 1999). Also included in the bottom panel is the mass-weighted PDF of the ISM, which we simply calculate for gas within a radius $r_{\text{ISM}} \equiv \min(0.2 R_{\text{vir}}, 2 \text{ kpc})$. The figure shows that in the absence of SN feedback, the majority of stars are born in a very dense medium ($n_{\text{H, form}} > 10^3 \text{ cm}^{-3}$). Because the stars do not undergo violent dynamical events (except for mergers), most of the metals are simply distributed around their birth cloud (Fig. 4).

The inclusion of TFB lowers the density peak of the PDF for star-forming clouds to $n_{\text{H, form}} \sim 10^3 \text{ cm}^{-3}$, and allows $\approx 10$ per cent of SNe to emerge at densities $n_{\text{H, SN}} \le 1 \text{ cm}^{-3}$. However, most of the SNe explode in the dense ISM ($n_{\text{H, SN}} \sim 10^3 \text{ cm}^{-3}$). Note that the shell formation radius at this density is only $\approx 2$ pc for metal-poor gas ($Z = 0.01 Z_{\odot}$; Thornton et al. 1998). Given that the radius should be resolved by at least three computational cells to capture the final momentum from SN (Kim & Ostriker 2015), it is not surprising that TFB fails to regulate star formation in the simulation with 12 pc resolution. Only SNe that explode in an ISM with $n_{\text{H}} \le 2 \text{ cm}^{-3}$ may approximate the Sedov–Taylor solution with 12 pc resolution, but this requires a very stringent refinement that a cell should not contain more than 110 M$_{\odot}$ of gas.

The KFB scheme by Dubois & Teyssier (2008) is originally designed to alleviate the overcooling problem by distributing a boosted momentum of $p_{\text{SN}} = \sqrt{(1 + \eta_w)2E_{\text{SN}}M_{\text{ej}}}$ to the neighbouring cells of an SN, where $\eta_w = \min(\eta_{\text{w, max}}, M_{\text{host}}/M_{\text{ej}})$ and $M_{\text{host}}$ is the gas mass of the host cell of an SN. However, if the adiabatic phase of SN expansion is not resolved, then radial momentum cannot build up and so the impact of SNe is essentially limited by the initial momentum injected to the ISM. Furthermore, the assumption used in the model that the mass loading only comes from the gas in the host cell of an SN reduced the estimation of the loading factor ($\eta_{\text{w}}$) in a low-density medium or in regions where there are multiple SNe, compared to the MFB run where neighbouring gas is also taken into account to calculate the swept-up mass. As a result, the sites of star formation and SN explosion in the KFB run turn out to be quite similar to those of the TFB run.

The most notable difference in the MFB runs (MFB, MFBm, and MFBmp) from others is that a significant fraction of SNe explodes in low densities of $n_{\text{H}} \le 1 \text{ cm}^{-3}$ (40, 26, and 78 per cent, respectively). The PDFs of $n_{\text{H, SN}}$ are bimodal in these runs, which reflect
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Differences in the opening angle of the outflow between the thermal SN feedback (top) and MFB run with the assumption of a highly porous ISM (bottom). The panels display the distribution of hydrogen number density on a sphere of 0.5 $R_{vir}$ (=13.3 kpc, left) at $z = 3.9$ when the flux-weighted outflowing velocity and density are very similar in the two runs. The right-hand panels correspond to the density maps at smaller radii (0.3 $R_{vir}$ = 8.0 kpc). The gas belonging to satellites galaxies or infalling gas is removed for clarity. The plot demonstrates that the outflow in the MFBmp run is stronger than the TFB run, because it subtends a larger solid angle than the former. At smaller radii, the NutMFBmp run drives winds composed of many gas clumps with $n_H \gtrsim 0.1 \text{ cm}^{-3}$.</figcaption></figure>

the fact that the low-density gas ($n_H \sim 10^{-2} \text{ cm}^{-2}$) becomes more volume-filling in the main galaxy. This low-density gas has a wide distribution of temperature ranging from $10^4$ to a few times $10^7$ K, indicating that the multiphase ISM is naturally generated with SN explosions (Cox & Smith 1974; McKee & Ostriker 1977). Kim et al. (2013) simulate SN-regulated star formation by injecting the final momentum and unlike us, they find that the ISM is composed of cold and warm ($T \lesssim 10^4$ K) gas. We attribute this difference to the fact that SNe in their simulations are assumed to emerge instantaneously as soon as a dense gas cloud of $n_H = 100 \text{ cm}^{-3}$ forms, whereas the time delay of 3–40 Myr allows SNe to explode not only in the dense environment, but also in a low-density medium in our simulations.

Fig. 17 demonstrates that a realistic time delay for SNe is crucial for preventing gas from collapsing to very high densities ($n_H > 10^4 \text{ cm}^{-3}$, the MFBm and MFBmp runs). This is essentially because more frequent explosions keep stirring dense gas. The bottom panel indicates that, in the case of MFBm, the mass fraction for gas with densities in the range $100 \lesssim n_H \lesssim 10^4 \text{ cm}^{-3}$ increases, lowering to $n_H \lesssim 300 \text{ cm}^{-3}$ the typical density at which an SN explodes. However, the gas mass contained in a 12 pc cell with a density of $100 \text{ cm}^{-3}$ ($\approx 5500 \text{ M}_\odot$) already exceeds the mass swept up by the end of the adiabatic phase ($\approx 1300 \text{ M}_\odot$, for 0.1 $Z_\odot$). Thus, the velocity gained from an individual SN would be only $\sim 2 \text{ km s}^{-1}$ if the final momentum of $2.4 \times 10^5 \text{ M}_\odot \text{ km s}^{-1}$ (equation 3) is equally distributed to 48 neighbouring cells of the same density. Given that the local escape velocity ($\sim 8 \text{ km s}^{-1}$) is greater than this, it is very unlikely that a single SN would disrupt the entire cloud. One way to unbind the gas is for collective winds to continually lift up the gas. This may be plausible because the gas experiencing a runaway collapse will always produce young stars. On the other hand, the ram pressure due to gas accretion on to a collapsing gas cloud is likely to counterbalance the expanding motion. The MFBm run suggests that such collective events do not occur frequently enough to unbind the gas. We note, however, that this may simply be due to the fact that the structure of the ISM in our simulation is intrinsically smooth at our resolution limits

(12–48 pc), and that fast outflows along low-density channels are not resolved.

Small-scale turbulence simulations show that the structure of the ISM is complex (e.g. Padoan, Jones & Nordlund 1997; Ostriker, Stone & Gammie 2001; Federrath & Klessen 2012; Molina et al. 2012), and that the volume-filling gas density is smaller than the mean density by $\exp[-0.5 \ln(1 + \beta b^2 M^2/(1 + \beta))]$, where $\beta$ is the ratio of thermal to magnetic pressure, $b$ is the turbulence forcing parameter ($0.3 \le b \le 1$), and $M$ is the Mach number. As aforementioned, it is likely that our simulations underestimate the extent to which SNe channel through a lower density medium on scales of 12–48 pc (Iffrig & Hennebelle 2015; Kim & Ostriker 2015), i.e. our resolution scale, hence, too much gas is likely to be entrained along with the SN ejecta from the host cell. Reducing the mass loading from the host cell, as in the MFBmp run, increases the outflow velocity in low-density neighbours, and the input momentum during the snowplough phase increases slightly ($\propto n_H^{2/17}$).

Fig. 18 illustrates an example of how much we expect the outflow velocity of each cell to be boosted when we assume that $f_{w,host} = 0.1$ (i.e. only 10 per cent of the gas from the host cell is entrained) instead of $f_{w,host} = 1$. We use the ISM structure of the main galaxy at $z = 3$ from the cooling run to compare the velocity for the case where $f_{w,host} = 0.1$ ($v_{f0.1}$) to that where $f_{w,host} = 1$ ($v$). We assume that young stars with $3 \le t \le 40$ Myr are SNe. The top panel shows that $v_{f0.1}$ can be a factor of $\sim 10$ larger than $v$. Note that such a significant increase is rare because it can only occur around a very sharp edge of a dense clump where there is a rapid spatial variation of the density ($\rho_{\text{nbor}}/\rho_{\text{host}} \sim 0.01$). We still find, however, that $\approx 10$ per cent of the gas cells that are directly affected by SN explosions would have a velocity more than double their value in the case in which the ISM structure is assumed to be smooth on the resolution scale. This occurs in gas cells where $\rho_{\text{nbor}}/\rho_{\text{host}} \sim 0.1$. Of the 10 per cent, the fraction of cells accelerated more than the typical escape velocity of the local GMCs (i.e. $10 \text{ km s}^{-1}$) increases from 25 to 65 per cent for $0.1 \lesssim n_H \lesssim 10 \text{ cm}^{-3}$. Thus, there is a higher probability that the gas is blown away through the lower density channels. On the other hand, in dense regions with $n_H \gtrsim 100 \text{ cm}^{-3}$, the velocity change
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Local environments at which star formation occurs (top) or SNe explode (middle). The PDFs of the density at which stars form ($n_{\text{H, form}}$) are obtained by computing $M_{\text{gas}}/t_{\text{ff}}$ of the ISM in the main galaxy at $3 \le z \le 10$. Similarly, we compute the PDFs of the density at which SNe explode ($n_{\text{H, SN}}$) using potential SN particles in the main galaxy (see the text). Note that although there is no energy input, the NutCO run also has SNe that release metals at 10 Myr after the birth of a star particle. Also included in the bottom panel is the mass-weighted PDF of the gas within a small radius [$r = \min(0.2 R_{\text{vir}}, 2 \text{ kpc})$]. The realistic time delay of SNe commencing at 3 Myr suppresses the formation of very dense gas ($n_{\text{H}} \gtrsim 10^4 \text{ cm}^{-3}$, MFBm and MFBmp), allowing SNe to explode at lower density ($n_{\text{H}} \lesssim 300 \text{ cm}^{-3}$) than the weak feedback runs (TFB or KFB). More SNe emerge in lower density environments when SN feedback becomes stronger.</figcaption></figure>

is not significant in general, as the local ISM structure is more or less smooth in the NutCO run and the mass contained in the neighbour is already significant, compared to that entrained from the host cell of SN. Even though only a small fraction of gas gains a substantial velocity, the comparison between MFBm and MFBmp indicates that this plays an important role in dispersing the star-forming gas (Fig. 17). We do not find any sign of overcooling in the MFBmp run. This leads us to conclude that the propagation of SN momentum through lower density channels may be crucial in modelling star formation in galaxies, providing a viable solution to the overcooling problem.

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Expected increase in the outflow velocity of gas when the ISM is assumed to be highly porous. The top panel shows the ratio between the outflow velocity estimated assuming that all the gas in the host cell of the SN is entrained (v) and the velocity in the case where the volume-filling density is 10 times lower than the mean density ($v_{f0.1}$). Approximately, 10 per cent of the cells that are directly affected by the SN explosions would have a twice faster outflow along low-density channels than the fiducial case. The bottom panel represents the velocities of the cells with $v_{f0.1} > 2v$. Black and orange points denote v and $v_{f0.1}$ in physical units, respectively. The knee present around $n_{\text{H}} \sim 1 \text{ cm}^{-3}$ corresponds to the transition between the adiabatic and snowplough phase.</figcaption></figure>

## 4.2 Comparison to previous work

To single out the impact from SNe, our simulations do not include any form of radiation pressure. However, other feedback processes come into play before the first SN ($t < 3$ Myr), and the role of SNe in this work may be somewhat exaggerated. Photoionization heating can reduce the density by an order of magnitude around young stars by overpressurizing the medium (Krumholz & Tan 2007; Sales et al. 2014). Although the determination of the IR optical depth is not straightforward (Krumholz & Thompson 2012; Davis et al. 2014; Rosdahl & Teyssier 2015), radiation pressure can also drive additional turbulence in optically thick regions. Including the additional feedback mechanisms may further reduce star formation than the model with SN feedback alone (e.g. Agertz et al. 2013) or the importance of SNe may simply be mitigated, as the early processes substitute their role. Bearing this in mind, we compare our results with previous studies in this section.

Using an AMR code, ENZO (Bryan et al. 2014), Hummels & Bryan (2012) performed 16 cosmological simulations of MW-like galaxies with a few hundreds of parsec resolution. They included SN feedback by heating the host cell of SN. Because the thermal energy is quickly radiated away, more than 50 per cent of the baryons ends up turning into stars, leading to a too peaked rotation curve (see also Joung et al. 2009). The resulting SFH is found to be smooth, in agreement with the results from the NutTFB run. Although there are several differences in numerical methods (e.g. hydrodynamic solver) between the two studies, this suggests that adopting a higher resolution does not solve the angular momentum
catastrophe (Navarro & Steinmetz 2000; Governato et al. 2004) so long as the spatial resolution is larger than the shell formation radius.

To compensate for the finite resolution (~50–300 parsecs) of cosmological simulations, several studies adopted the cooling suppression model (Stinson et al. 2006; Governato et al. 2007, 2010; Avila-Reese et al. 2011; Brooks et al. 2011; Guedes et al. 2011; Teyssier et al. 2013; Shen et al. 2014). As gas is prevented from collapsing on to the centre of dwarf to MW-sized galaxies by the overpressure from SNe, galaxies in these studies showed a smoothly rising rotation curve in general, as in observations. Because they only presented the results at low redshifts ($z \le 1$), it is difficult to make a direct comparison with our work, but SFHs in the MW-sized galaxies appear rather smooth at high redshifts ($z \gtrsim 2$; Guedes et al. 2011, cf. Shen et al. 2014 for more episodic SF in dwarf galaxies), compared to our MFBmp run or Hopkins et al. (2014). Hummels & Bryan (2012) also examined the model in which the radiative cooling is turned off for 50 Myr in a $10^{12}$ M$_\odot$ halo, and found that the central peak in the rotation curve is reduced, but still too many stars formed (~60 per cent), compared to the local $M_{\text{star}}-M_{\text{halo}}$ relation (e.g. Moster et al. 2010).

Recently, Hopkins et al. (2011, 2012b) investigated the effects of momentum transfer from the UV-IR photons on the ISM in four different types of isolated galaxies (Small Magellanic Cloud – SMC-like dwarf, MW-like spiral, gas-rich star-forming galaxies hosted by an ~$10^{11}$ M$_\odot$ halo, and star-bursting disc in a $10^{12}$ M$_\odot$ halo) with an entropy-conserving SPH code (Springel et al. 2005). They carried out a series of simulations, varying the feedback model, with parsec scale resolution, and concluded that the radiation pressure is mainly responsible for the generation of galactic winds in the actively star-forming, gas-rich galaxies, whereas SN explosions play a major role in relatively quiescent galaxies (SMC and MW-like galaxies). However, their SN feedback model is based on the initial energy of $10^{51}$ erg, which is likely to have underestimated its impact. Our Nut galaxy is similar to their gas-rich, star-forming galaxy (Sbc), but we find that the SN explosions are capable of driving a strong outflow if the ISM is assumed to be porous (MFBmp). As mentioned, the wind efficiency ($\dot{M}_{\text{out}}/\dot{M}_{\text{star}}$) from the MFBmp run (Fig. 15) is comparable or slightly larger than the star-forming galaxy (Sbc) in Hopkins et al. (2012b). In a subsequent paper, Hopkins et al. (2014) incorporate the final momentum from the SN blast wave in a dense medium, in a very similar spirit to our MFB scheme, and concluded that both SNe and radiation pressure are essential for reproducing the observed relations, such as $M_{\text{star}}-M_{\text{halo}}$ relation or Kennicutt-Schmidt law (Kennicutt 1998; Bigiel et al. 2008).

By including SN, stellar winds, and momentum transfer from radiation in the RAMSES code (Teyssier 2002), Agertz et al. (2013) systematically examined their relative importance in an isolated cloud of $10^6$ M$_\odot$ and a disc galaxy embedded in a halo of mass $10^{12}$ M$_\odot$. Because of their finite resolution, they used the local size and mass relation of cluster/clumps to estimate the surface density, which is distinct from Hopkins et al. (2011), who directly measured it from their simulations. Agertz et al. (2013) pointed out that the inclusion of early feedback is important for disrupting the dense gas in star-forming regions. Applying the radiation pressure model to cosmological simulations, Agertz & Kravtsov (2015) argue that reproducing a realistic disc galaxy requires two conditions: (i) radiation pressure should be included on top of the energy variable that dissipates away slowly on a 10 Myr time-scale, (ii) star formation should be modelled, such that a large fraction ($\epsilon = 10$ per cent) of gas is converted into stars on a local dynamical time-scale. They find

that both the $M_{\text{star}}-M_{\text{halo}}$ relation and the MZR are well accounted for if these conditions are met simultaneously. It is interesting to note that when a low star formation efficiency ($\epsilon = 1$ per cent) is used, their radiation pressure model could not slow down star formation in the MW-like galaxy and converted 50 per cent of the baryons into stars. The need for a high $\epsilon$ appears to suggest that star formation should be bursty for radiation pressure to control the growth of galaxies. Given that our MFBmp run achieved such bursty SFH with a relatively low star formation efficiency ($\epsilon = 2$ per cent), the requirement for a high $\epsilon$ may be model-dependent.

Several authors also studied the role of radiation pressure, and come to a different conclusion. Using the ART code (Kravtsov, Klypin & Khokhlov 1997), Ceverino et al. (2014) argue that radiation pressure from the ionizing photons alone provide enough energy to unbind star-forming gas clumps, making the role of multiply scattered IR photons to be less significant. However, the results are not straightforward to compare with other studies, because they model the radiation pressure as a non-thermal term that does not dissipate away until the density becomes lower than the star formation threshold ($n_H = 1 \text{ cm}^{-3}$). Agertz et al. (2013) tested a similar model by using a 'feedback energy variable', and showed that the resulting SFHs are quite sensitive to the dissipation time-scale as well as the fraction of the energy that is deposited into the non-thermal term. By contrast, Aumer et al. (2013) claim that a high IR optical depth of ~20 is requisite at $z > 2$ to form a realistic disc galaxy in SPH simulations. They could reasonably match not only the $M_{\text{star}}-M_{\text{halo}}$ relation but also the size and metallicity evolution in 16 dwarf to intermediate-size galaxies resimulated from the Aquarius project (Springel et al. 2008) and some simulations presented in Oser et al. (2010). The different conclusions on the role of radiation pressure by UV to IR photons call for further tests based on hydrodynamics calculations fully coupled with radiative transfer (Krumholz & Thompson 2012; Davis et al. 2014; Rosdahl & Teyssier 2015).

## 5 CONCLUSIONS

We have investigated the impact of SN feedback on the evolution of the progenitor of an MW-like galaxy, embedded in a $10^{11}$ M$_\odot$ halo at high redshifts ($z \ge 3$). For this purpose, we have carried out a suite of zoom-in cosmological hydrodynamic simulation with high spatial (12 pc) and mass ($m_{\text{dm}} = 5.5 \times 10^4$ M$_\odot$, $m_{\text{star}} = 610$ M$_\odot$) resolution using the AMR code, RAMSES, varying the SN feedback model. Our main results can be summarized as follows.

(i) Our simulations confirm that too many stars would form when the feedback from SN cannot drive strong outflows. The cooling run without SN feedback (NutCO) turns 85 per cent of baryons into stars inside $0.2 R_{\text{vir}}$ of the main halo ($M_{\text{halo}} \approx 10^{11}$ M$_\odot$), resulting in a much smaller gas-to-stellar mass ratio (Fig. 5) than observations at high redshifts. The run with thermal feedback (NutTFB) also shows the overcooling problem, as the cooling radius of the SN blast wave is still underresolved even in our zoom-in simulations. Inclusion of the SN explosion in a kinetic form with a small mass loading of $\eta_w \le 10$ (NutKFB) does not greatly help either, as shock-heated gas is again subject to artificial radiative losses and the subsequent generation of radial momentum during the adiabatic phase is not guaranteed. Nevertheless, the increase in the input momentum ($\sqrt{1+\eta_w}$) used in the KFB scheme helps suppress star formation in low-mass haloes ($M_{\text{halo}} \lesssim 10^{10}$ M$_\odot$) by an order of magnitude, compared to the cooling run (Fig. 8).

(ii) The MFB models (NutMFB, MFBm, and MFBmp) that transfers the correct amount of radial momentum from the adiabatic to
snowplough phase is able to substantially suppress star formation in small haloes ($M_{\text{halo}} \lesssim 10^{10} M_{\odot}$). The galaxy stellar mass of the haloes is reduced by roughly two orders of magnitude from that found in the cooling run (Fig. 8).

(iii) The prediction of SFHs in the main galaxy depends on the details of how momentum from SN is distributed. The model in which SNe emerge relatively late (10 Myr, MFB) could not prevent a runaway collapse of a massive cloud at the galaxy centre during frequent mergers of clumps (at $z \lesssim 5$, Fig. 7). This led to a rapid increase in galaxy stellar mass from $z \sim 5$. The final stellar mass at $z=3$ was overpredicted, compared with the typical observational estimates from the clustering analysis of LAEs or the abundance matching technique.

(iv) Inclusion of a realistic time delay for SNe, commencing at 3 Myr, efficiently suppresses the build-up of very dense gas ($n_{\text{H}} \gtrsim 10^4 \text{cm}^{-3}$) in the main galaxy (Fig. 17, MFBm). But because gas outflow velocities become smaller, distributing SN explosions in time makes a galaxy more susceptible to the overcooling problem than the case where at least 10 SNe explode simultaneously (MFB). As a result, star formation is not strongly regulated from high redshifts ($z \gtrsim 5$), and follows that of the weak feedback models (TFB or KFB, Fig. 5).

(v) The observed relations, such as the MZR and stellar mass to halo mass ratio, are reproduced by the SN feedback model with the realistic time delay, provided that the ISM is highly porous (MFBmp). The latter assumption is motivated by the observation that galaxies at high redshifts are more turbulent than their local counterpart (Förster Schreiber et al. 2009). If the ISM is highly turbulent, the volume-filling density of the medium would be systemically smaller than the mean. This means that taking the uniform gas density at the computational grid scale (as in MFBm) is likely to underestimate the wind velocities along low-density channels by entraining too much gas from the host cell of SN. We find that reducing the volume-filling density to 10 per cent (equivalent to the supersonically turbulent ISM with the Mach number of $\sim 10$) suffices to prevent the overcooling issue, yielding a mass loading that is 10 times larger than the SFR ($\eta \sim 10$, Fig. 15), as postulated in the momentum-driven winds (Oppenheimer et al. 2010). Although the degree of turbulence is arbitrarily chosen, this demonstrates that the momentum input from SNe alone may be able to regulate star formation, as observed. However, future works based on resolved SN feedback will be required to draw a firm conclusion on the relative role of different stellar processes on the evolution of galaxies though.

(vi) The MFBmp model shows stronger outflow than the weak feedback models (e.g. TFB, Fig. 14 b and c). The difference originates from the opening angle of the outflow on the sky (Fig. 16), not the density or velocity, which are found to be insensitive to the feedback prescriptions (Fig. 14 e, f, and i). On the other hand, metallicity of the outflow is lower in the model with stronger feedback and generally follows that of the galactic gas. The fact that outflow metallicity ($0.07 \lesssim Z/Z_{\odot} \lesssim 0.7$) is much lower than the metallicity of the ejecta ($2.5 Z_{\odot}$) indicates that the outflow is efficiently mixed with the ISM before it is blown away. The strong outflow in this MFBmp model leads to the formation of a spheroid with $M_{\text{star}} = 7.8 \times 10^8 M_{\odot}$ rather than a well-ordered stellar disc at $z=3$ (Fig. 13).

(vii) The cold filamentary inflow around the Nut halo ($M_{\text{halo}} \approx 10^{11} M_{\odot}$) is impervious to SN feedback. The accretion rate at $R_{\text{vir}}$ is found to be more or less similar, irrespective of the strength of the feedback (Fig. 14 a). The metallicity of the inflow ($Z = 0.005-0.02 Z_{\odot}$) is an order of magnitude smaller

than that of the outflow, and not very sensitive to the feedback model.

## ACKNOWLEDGEMENTS

We are grateful to Sam Geen, Chang-Goo Kim, and Eve Ostriker for fruitful discussions, and Romain Teysnier for making his code RAMSES publicly available. We also thank the anonymous referee for suggestions, which improved the clarity of the paper. Computing resources were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, the DiRAC facility jointly funded by STFC and the Large Facilities Capital Fund of BIS, and the Horizon-UK programme through DiRAC-2 facilities. The research is supported in part by NSF grant AST-1108700 and NASA grant NNX12AF91G and in part by grant Spin(e) ANR-13-BS05-0005 of the french ANR. JD and AS's research is supported by funding from Adrian Beecroft, the Oxford Martin School and the STFC.

## REFERENCES

* Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 597, 21

* Adelberger K. L., Steidel C. C., Pettini M., Shapley A. E., Reddy N. A., Erb D. K., 2005, ApJ, 619, 697

* Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18

* Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

* Arnouts S. et al., 2002, MNRAS, 329, 355

* Arribas S., Colina L., Bellocchi E., Maiolino R., Villar-Martín M., 2014, A&A, 568, A14

* Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

* Aumer M., White S. D. M., Naab T., Scannapieco C., 2013, MNRAS, 434, 3142

* Avila-Reese V., Colín P., González-Samaniego A., Valenzuela O., Firmani C., Velázquez H., Ceverino D., 2011, ApJ, 736, 134

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Bell E. F., McIntosh D. H., Katz N., Weinberg M. D., 2003, ApJS, 149, 289

* Bertschinger E., 2001, ApJS, 137, 1

* Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215

* Bielby R. et al., 2013, MNRAS, 430, 425

* Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

* Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ, 500, 342

* Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150

* Brook C. B. et al., 2011, MNRAS, 415, 1051

* Brooks A. M. et al., 2011, ApJ, 728, 51

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Bryan G. L. et al., 2014, ApJS, 211, 19

* Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240

* Calura F., Pipino A., Chiappini C., Matteucci F., Maiolino R., 2009, A&A, 504, 373

* Catelan P., Theuns T., 1996, MNRAS, 282, 436

* Cen R., Ostriker J. P., 1992, ApJ, 399, L113

* Ceverino D., Klypin A., Klimek E. S., Trujillo-Gomez S., Churchill C. W., Primack J., Dekel A., 2014, MNRAS, 442, 1545

* Chabrier G., Baraffe I., Allard F., Hauschildt P. H., 2005, preprint (astro-ph/0509798)

* Chevalier R. A., 1974, ApJ, 188, 501

* Chisholm J., Tremonti C. A., Leitherer C., Chen Y., Wofford A., Lundgren B., 2014, preprint (arXiv:1412.2139)

* Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252

* Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
Cooper J. L., Bicknell G. V., Sutherland R. S., Bland-Hawthorn J., 2008, ApJ, 674, 157

Cox D. P., Smith B. W., 1974, ApJ, 189, L105

Crighton N. H. M., Hennawi J. F., Prochaska J. X., 2013, ApJ, 776, L18

Daddi E. et al., 2007, ApJ, 670, 156

Dalcanton J. J., 2007, ApJ, 658, 941

Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694

Danovich M., Dekel A., Hahn O., Ceverino D., Primack J., 2015, MNRAS, 449, 2087

Davis S. W., Jiang Y.-F., Stone J. M., Murray N., 2014, ApJ, 796, 107

Doroshkevich A. G., 1970, Astrophysics, 6, 320

Dubois Y., Teyssier R., 2008, A&A, 477, 79

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dunkley J. et al., 2009, ApJ, 701, 1804

Durkalec A. et al., 2015, A&A, 576, L7

Elbaz D. et al., 2007, A&A, 468, 33

Erb D. K., Shapley A. E., Pettini M., Steidel C. C., Reddy N. A., Adelberger K. L., 2006, ApJ, 644, 813

Evans N. J., II et al., 2009, ApJS, 181, 321

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Faucher-Giguère C.-A., Kereš D., Dijkstra M., Hernquist L., Zaldarriaga M., 2010, ApJ, 725, 633

Faucher-Giguère C.-A., Kereš D., Ma C.-P., 2011, MNRAS, 417, 2982

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Finlator K., Davé R., 2008, MNRAS, 385, 2181

Förster Schreiber N. M. et al., 2009, ApJ, 706, 1364

Garnett D. R., 2002, ApJ, 581, 1019

Gawiser E. et al., 2007, ApJ, 671, 278

Geen S., Slyz A., Devriendt J., 2013, MNRAS, 429, 633

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Giavalisco M., Dickinson M., 2001, ApJ, 550, 177

Gong H., Ostriker E. C., 2011, ApJ, 729, 120

Gong H., Ostriker E. C., 2013, ApJS, 204, 8

González V., Labbé I., Bouwens R. J., Illingworth G., Franx M., Kriek M., Brammer G. B., 2010, ApJ, 713, 115

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Governato F. et al., 2004, ApJ, 607, 688

Governato F., Willman B., Mayer L., Brooks A., Stinson G., Valenzuela O., Wadsley J., Quinn T., 2007, MNRAS, 374, 1479

Governato F. et al., 2010, Nature, 463, 203

Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Guo Q., White S., Li C., Boylan-Kolchin M., 2010, MNRAS, 404, 1111

Haardt F., Madau P., 1996, ApJ, 461, 20

Hennebelle P., Iffrig O., 2014, A&A, 570, A81

Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950

Hopkins P. F., Quataert E., Murray N., 2012a, MNRAS, 421, 3488

Hopkins P. F., Quataert E., Murray N., 2012b, MNRAS, 421, 3522

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hummels C. B., Bryan G. L., 2012, ApJ, 749, 140

Iffrig O., Hennebelle P., 2015, A&A, 576, A95

Joung M. R., Cen R., Bryan G. L., 2009, ApJ, 692, L1

Kassin S. A., de Jong R. S., Weiner B. J., 2006, ApJ, 643, 804

Kassin S. A. et al., 2012, ApJ, 758, 106

Katz N., 1992, ApJ, 391, 502

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

Khochfar S., Silk J., Windhorst R. A., Ryan R. E., Jr, 2007, ApJ, 668, L115

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

Kim C.-G., Ostriker E. C., Kim W.-T., 2013, ApJ, 776, 1

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Devriendt J., Slyz A., Pichon C., Kassin S. A., Dubois Y., 2011a, preprint (arXiv:1106.0538)

Kimm T., Slyz A., Devriendt J., Pichon C., 2011b, MNRAS, 413, L51

Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73

Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304

Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155

Larson R. B., 1969, MNRAS, 145, 271

Lee K.-S., Giavalisco M., Conroy C., Wechsler R. H., Ferguson H. C., Somerville R. S., Dickinson M. E., Urry C. M., 2009, ApJ, 695, 368

Leitherer C. et al., 1999, ApJS, 123, 3

McCarthy P. J., van Breugel W., Heckman T., 1987, AJ, 93, 264

McGaugh S. S., Schombert J. M., de Blok W. J. G., Zagursky M. J., 2010, ApJ, 708, L14

McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148

Maiolino R. et al., 2008, A&A, 488, 463

Mannucci F. et al., 2009, MNRAS, 398, 1915

Martel H., Shapiro P. R., 1998, MNRAS, 297, 467

Martin C. L., 1999, ApJ, 513, 156

Martin C. L., 2005, ApJ, 621, 227

Martin C. L., Shapley A. E., Coil A. L., Kornei K. A., Bundy K., Weiner B. J., Noeske K. G., Schiminovich D., 2012, ApJ, 760, 127

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680

Mori M., Yoshii Y., Tsujimoto T., Nomoto K., 1997, ApJ, 478, L21

Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191

Navarro J. F., Steinmetz M., 2000, ApJ, 538, 477

Navarro J. F., White S. D. M., 1993, MNRAS, 265, 271

Newman S. F. et al., 2012, ApJ, 761, 43

Noeske K. G. et al., 2007, ApJ, 660, L43

Nomoto K., Suzuki T., Shigeyama T., Kumagai S., Yamaoka H., Saio H., 1993, Nature, 364, 507

Oh S.-H., Brook C., Governato F., Brinks E., Mayer L., de Blok W. J. G., Brooks A., Walter F., 2011, AJ, 142, 24

Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265

Oppenheimer B. D., Davé R., Kereš D., Fardal M., Katz N., Kollmeier J. A., Weinberg D. H., 2010, MNRAS, 406, 2325

Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

Ostriker E. C., Stone J. M., Gammie C. F., 2001, ApJ, 546, 980

Padoan P., Jones B. J. T., Nordlund Å. P., 1997, ApJ, 474, 730

Peebles P. J. E., 1969, ApJ, 155, 393

Peirani S., Mohayaee R., de Freitas Pacheco J. A., 2004, MNRAS, 348, 921

Penston M. V., 1969, MNRAS, 144, 425

Pichon C., Pogosyan D., Kimm T., Slyz A., Devriendt J., Dubois Y., 2011, MNRAS, 418, 2493

Pilyugin L. S., Vílchez J. M., Contini T., 2004, A&A, 425, 849

Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671

Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R., Gottlieber S., 2008, ApJS, 178, 179

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492

Rosdahl J., Blaizot J., 2012, MNRAS, 423, 344

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, preprint (arXiv:1501.04632)

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837

Rupke D. S., Veilleux S., Sanders D. B., 2005, ApJS, 160, 115

Sales L. V., Marinacci F., Springel V., Petkova M., 2014, MNRAS, 439, 2990

Salim S. et al., 2007, ApJS, 173, 267

Salmon B. et al., 2015, ApJ, 799, 183

Scannapieco C., Tissera P. B., White S. D. M., Springel V., 2006, MNRAS, 371, 1125
Scannapieco C., White S. D. M., Springel V., Tissera P. B., 2009, MNRAS, 396, 696

Schmidt M., 1959, ApJ, 129, 243

Sedov L. I., 1959, Similarity and Dimensional Methods in Mechanics. Academic Press, New York

Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., 2003, ApJ, 588, 65

Sharma P., Roy A., Nath B. B., Shchekinov Y., 2014, MNRAS, 443, 3463

Shen S., Madau P., Conroy C., Governato F., Mayer L., 2014, ApJ, 792, 99

Shetty R., Ostriker E. C., 2008, ApJ, 684, 978

Shetty R., Ostriker E. C., 2012, ApJ, 754, 2

Sijacki D., Vogelsberger M., Kereš D., Springel V., Hernquist L., 2012, MNRAS, 424, 2999

Slyz A. D., Devriendt J. E. G., Bryan G., Silk J., 2005, MNRAS, 356, 737

Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, MNRAS, 391, 481

Springel V. et al., 2005, Nature, 435, 629

Springel V. et al., 2008, MNRAS, 391, 1685

Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289

Steidel C. C. et al., 2014, ApJ, 795, 165

Stewart K. R., Brooks A. M., Bullock J. S., Maller A. H., Diemand J., Wadsley J., Moustakas L. A., 2013, ApJ, 769, 74

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Strickland D. K., Heckman T. M., 2009, ApJ, 697, 2030

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Tacconi L. J. et al., 2010, Nature, 463, 781

Taylor G., 1950, Proc. R. Soc. A, 201, 159

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Thacker R. J., Couchman H. M. P., 2001, ApJ, 555, L17

Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998, ApJ, 500, 95

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Tremonti C. A. et al., 2004, ApJ, 613, 898

Tully R. B., Fisher J. R., 1977, A&A, 54, 661

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

van den Bosch F. C., Abel T., Croft R. A. C., Hernquist L., White S. D. M., 2002, ApJ, 576, 21

Vazquez-Semadeni E., 1994, ApJ, 423, 681

Vázquez-Semadeni E., Colín P., Gómez G. C., Ballesteros-Paredes J., Watson A. W., 2010, ApJ, 715, 1302

Walch S. K., Whitworth A. P., Bisbas T., Wünsch R., Hubber D., 2012, MNRAS, 427, 625

Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt R. C., Jr, Thornley M. D., Leroy A., 2008, AJ, 136, 2563

Weiner B. J. et al., 2009, ApJ, 692, 187

Werk J. K. et al., 2014, ApJ, 792, 8

Whalen D., van Veelen B., O'Shea B. W., Norman M. L., 2008, ApJ, 682, 49

Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, ApJ, 754, L29

White S. D. M., 1984, ApJ, 286, 38

Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012a, MNRAS, 427, 311

Wise J. H., Turk M. J., Norman M. L., Abel T., 2012b, ApJ, 745, 50

Wuyts E. et al., 2014, ApJ, 789, L40

Yuan T.-T., Kewley L. J., Richard J., 2013, ApJ, 763, 9

Zahid H. J., Geller M. J., Kewley L. J., Hwang H. S., Fabricant D. G., Kurtz M. J., 2013, ApJ, 771, L19

This paper has been typeset from a *TEX/LATEX* file prepared by the author.