Simulating the interstellar medium and stellar feedback on a moving mesh: Implementation and isolated galaxies

TL;DR

This paper presents the SMUGGLE model, a detailed stellar feedback implementation for the moving-mesh code AREPO, which accurately simulates the multiphase interstellar medium and gaseous outflows in Milky Way-like galaxies.

Contribution

The paper introduces the SMUGGLE model, a comprehensive stellar feedback scheme that self-consistently captures multiphase gas, outflows, and star formation regulation in high-resolution galaxy simulations.

Findings

Stellar feedback regulates star formation to observed levels.

Gaseous outflows with mass loading factors around unity are produced.

Multiphase interstellar medium with cold, warm, and hot gas phases is reproduced.

Abstract

We introduce the Stars and MUltiphase Gas in GaLaxiEs -- SMUGGLE model, an explicit and comprehensive stellar feedback model for the moving-mesh code arepo. This novel sub-resolution model resolves the multiphase gas structure of the interstellar medium and self-consistently generates gaseous outflows. The model implements crucial aspects of stellar feedback including photoionization, radiation pressure, energy and momentum injection from stellar winds and from supernovae. We explore this model in high-resolution isolated simulations of Milky Way-like disc galaxies. Stellar feedback regulates star formation to the observed level and naturally captures the establishment of a Kennicutt-Schmidt relation. This result is achieved independent of the numerical mass and spatial resolution of the simulations. Gaseous outflows are generated with average mass loading factors of the order of unity. Strong outflow activity is correlated with peaks in the star formation history of the galaxy with evidence that most of the ejected gas eventually rains down onto the disc in a galactic fountain flow that sustains late-time star formation. Finally, the interstellar gas in the galaxy shows a distinct multiphase distribution with a coexistence of cold, warm and hot phases.

Figures

Tables

Table 1. Table 1: Structural parameters of the Milky Way-type galaxy considered in this work. Columns indicate (from left to right): mass of the dark matter halo ( M halo subscript 𝑀 halo M_{\rm halo} ), circular velocity of the halo at the radius at which its average density is equal to 200 times critical density for closure ( v 200 subscript 𝑣 200 v_{200} ), halo concentration ( c 𝑐 c ), bulge mass ( M b subscript 𝑀 b M_{\rm b} ), bulge ( Hernquist 1990 sphere) scale length ( a 𝑎 a ), stellar disc mass ( M d subscript 𝑀 d M_{\rm d} ), stellar disc scale length ( r d subscript 𝑟 d r_{\rm d} ), stellar disc scale height ( h ℎ h ), gas disc mass ( M g subscript 𝑀 g M_{\rm g} ), gas disc scale length ( r g subscript 𝑟 g r_{\rm g} ), and gas fraction within R ⊙ = 8.5 ​ kpc subscript 𝑅 direct-product 8.5 kpc R_{\odot}=8.5\,{\rm kpc} (computed as the ratio between gaseous and stellar disc masses; f gas subscript 𝑓 gas f_{\rm gas} ).

$M_{\mathrm{halo}}$	$v_{200}$	$c$	$M_b$	$a$	$M_d$	$r_d$	$h$	$M_g$	$r_g$	$f_{\mathrm{gas}}$
($M_{\odot }$)	($\mathrm{km}{s}^{-1}$)		($M_{\odot }$)	(kpc)	($M_{\odot }$)	(kpc)	(pc)	($M_{\odot }$)	(kpc)	()
$1.53\times {10}^{12}$	169	12	$1.5\times {10}^{10}$	1.0	$4.73\times {10}^{10}$	3.0	300	$9\times {10}^9$	6.0	0.10

Table 2. Table 2: Parameters defining the resolution level at which the initial conditions are sampled. Columns indicate (from left to right): resolution level name, gravitational softening length for star particles ( ϵ ⋆ subscript italic-ϵ ⋆ \epsilon_{\star} ), minimum gravitational softening length for gas cells ( ϵ g subscript italic-ϵ g \epsilon_{\rm g} ; softening is adaptive and scales logarithmically with the cell size), target mass of gas cell ( m g subscript 𝑚 g m_{\rm g} ), mass of a stellar bulge particle ( m b subscript 𝑚 b m_{\rm b} ), and mass of a stellar disc particle ( m d subscript 𝑚 d m_{\rm d} ).

Resolution	${ϵ}_{\star }$	${ϵ}_g$	$m_g$	$m_b$	$m_d$
level	(pc)	(pc)	($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)
low	50.0	21.4	$9.0\times {10}^4$	$1.5\times {10}^5$	$1.2\times {10}^5$
intermediate	21.4	10.0	$1.1\times {10}^4$	$2.0\times {10}^4$	$1.5\times {10}^4$
high	7.2	3.6	$1.4\times {10}^3$	$2.3\times {10}^3$	$1.9\times {10}^3$

Table 3. Table 3: Parameters adopted in this works for physical processes implemented in our ISM and stellar feedback model with their fiducial value and description.

Parameter	Fiducial value	Units	Description
Star formation
$ϵ$	$0.01$	–	Star formation efficiency
${\rho }_{\mathrm{th}}$	$100$	${\mathrm{cm}}^{-3}$	Star formation density threshold
Supernova feedback
$f_{\mathrm{SN}}$	$1$	–	SN energy relative to fiducial value
$E_{\mathrm{SN}}$	${10}^{51}$	$\mathrm{erg}$	Fiducial energy per SN
$N_{\mathrm{ngb}}$	$64$	–	Effective neighbour number
$\Delta N_{\mathrm{ngb}}$	$1$	–	Neighbour number tolerance
Radiative feedback
${\kappa }_{\mathrm{IR}}$	$10\left(Z/Z_{\odot }\right)$	${\mathrm{cm}}^2{g}^{-1}$	Gas infrared opacity
$t_{\star ,\max }$	$5$	$\mathrm{Myr}$	Maximum stellar age for ionizing radiation
${\gamma }_{\star }$	${10}^3$	$L_{\odot }{M}_{\odot }^{-1}$	Stellar mass-to-light ratio
$T_{\mathrm{phot}}$	$1.7\times {10}^4$	$K$	Temperature of photoionized gas
$⟨h\nu ⟩$	$17$	$\mathrm{eV}$	Average ionizing photon energy

Equations

$$\begin{array}{rl}{\Lambda }_{\mathrm{mol}}=& 2.896\times 10^{-26}\{{\left(\frac{T}{125.215\mathrm{K}}\right)}^{-4.9202}+\\ & {{\left(\frac{T}{1349.86\mathrm{K}}\right)}^{-1.7288}+{\left(\frac{T}{6450.06\mathrm{K}}\right)}^{-0.3075}\}}^{-1}\times \\ & \left(0.001+\frac{0.10n_{\mathrm{H}}}{1+n_{\mathrm{H}}}+\frac{0.09n_{\mathrm{H}}}{1+0.1n_{\mathrm{H}}}+\frac{(Z/Z_{\text{sun}}{)}^2}{1+n_{\mathrm{H}}}\right)\times \\ & \left(\frac{1+(Z/Z_{\text{sun}})}{1+0.00143n_{\mathrm{H}}}\right)\times \exp \left(-{\left[\frac{T}{158000\mathrm{K}}\right]}^2\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^3.\end{array}$$

$$f_{\mathrm{ssh}}=(1-f){\left[1+{\left(\frac{n_{\mathrm{H}}}{n_0}\right)}^{\beta }\right]}^{{\alpha }_1}+f{\left[1+\left(\frac{n_{\mathrm{H}}}{n_0}\right)\right]}^{{\alpha }_2},$$

$${\Lambda }_{\mathrm{CR}}=-10^{-16}(0.98+1.65){\tilde{n}}_e{e}_{\mathrm{CR}}{n}_{\mathrm{H}}^{-1}\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^3,$$

$$e_{\mathrm{CR}}=\{\begin{array}{ll}9\times 10^{-12}& \text{ergcm−3}\text{for nH>0.01cm−3}\\ & \\ 9\times 10^{-12}\left(\frac{n_{\mathrm{H}}}{0.01{\mathrm{cm}}^{-3}}\right)& \text{ergcm−3}\text{for nH≤0.01cm−3},\end{array}$$

$$\begin{gathered} {\Lambda }_{\mathrm{phot}}=-1.3\times 10^{-24}{\tilde{e}}_{\nu }^{\mathrm{pe}}{n}_{\mathrm{H}}^{-1}\left(\frac{Z}{Z_{\text{sun}}}\right)\times \\ \left(\frac{0.049}{1+(x_{\mathrm{pe}}/1925{)}^{0.73}}+\frac{0.037(T/10^4\mathrm{K}{)}^{0.7}}{1+(x_{\mathrm{pe}}/5000)}\right)\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^3, \end{gathered}$$

$$\begin{array}{rl}{x}_{\mathrm{pe}}& \equiv \frac{{\tilde{e}}_{\nu }^{\mathrm{pe}}{T}^{0.5}}{{\Phi }_{\mathrm{PAH}}{\tilde{n}}_e{n}_{\mathrm{H}}},\end{array}$$

$$\dot{M}_{\star}=\left\{\begin{array}[]{ll}\,\,\,\,\,\,0&\rho<\rho_{\rm th}\\ \\ \displaystyle\epsilon\frac{M_{\rm gas}}{t_{\rm dyn}}&\rho\geq\rho_{\rm th},\end{array}\right.$$

$$t_{\mathrm{dyn}}=\sqrt{\frac{3\pi }{32G{\rho }_{\mathrm{gas}}}},$$

$${\alpha }_i=\frac{||\nabla \otimes {\mathbf{v}}_i|{|}^2+(c_{s,i}/\Delta x_i{)}^2}{8\pi G{\rho }_i},$$

$$||\nabla \otimes {\mathbf{v}}_i|{|}^2\equiv \sum _{i,j}{\left(\frac{\partial v_i}{\partial x_j}\right)}^2.$$

$$M_{\star ,i}=M_i\left[1-\exp \left(-\frac{{\dot{M}}_{\star }\Delta t}{M_i}\right)\right],$$

$$f_{{\mathrm{H}}_2}=1-\frac{3s}{4+s},$$

$$s=\frac{\log [1+\chi (0.6+0.01\chi )]}{0.6\tau },$$

$$\chi =2.3\frac{1+3.1(Z/Z_{\odot }{)}^{0.365}}{3},$$

$$\tau =0.067\left(\frac{Z}{Z_{\odot }}\right)\left(\frac{\Sigma }{{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2}}\right),$$

$$\Sigma =\rho \left(\frac{\rho }{||\nabla \rho ||}\right).$$

$$E_{\mathrm{SN}}=f_{\mathrm{SN}}{E}_{51},$$

$$v_{\mathrm{SN}}=\sqrt{\frac{2E_{\mathrm{SN}}}{M_{\mathrm{SN}}}},$$

$$p_{\mathrm{SN}}=M_{\mathrm{SN}}{v}_{\mathrm{SN}}=\sqrt{2E_{\mathrm{SN}}{M}_{\mathrm{SN}}}.$$

$$E_{\mathrm{SN},\mathrm{tot}}=f_{\mathrm{SN}}{E}_{51}(N_{\mathrm{SNII}}+N_{\mathrm{SNIa}}),$$

$$\begin{array}{rl}& p_{\mathrm{SN},\mathrm{tot}}=p_{\mathrm{SNII},\mathrm{tot}}+p_{\mathrm{SNIa},\mathrm{tot}}=\\ & \sqrt{2N_{\mathrm{SNII}}{E}_{\mathrm{SN}}{M}_{\mathrm{SNII},\mathrm{tot}}}+\sqrt{2N_{\mathrm{SNIa}}{E}_{\mathrm{SN}}{M}_{\mathrm{SNIa},\mathrm{tot}}}\end{array}$$

$$N_{\mathrm{SNII}}=M_{\star }{\int }_{M(t+\Delta t)}^{M(t)}\Phi (m)\mathrm{d}m,$$

$$M_{\mathrm{SNII},\mathrm{tot}}=M_{\star }{\int }_{M(t+\Delta t)}^{M(t)}M{f}_{\mathrm{rec}}(m,Z)\Phi (m)\mathrm{d}m.$$

$$N_{\mathrm{SNIa}}={\int }_t^{t+\Delta t}\mathrm{DTD}(t^{\prime })\mathrm{d}{t}^{\prime }.$$

$$\mathrm{DTD}(t)=\Theta (t-t_8)N_0{\left(\frac{t}{{\tau }_8}\right)}^{-s}\frac{s-1}{{\tau }_8},$$

$$M_{\mathrm{SNIa},\mathrm{tot}}=M_{\mathrm{SNIa}}{N}_{\mathrm{SNIa}}.$$

$$\Delta t_{\star }=\min \left(\Delta t_{\mathrm{grav}},\Delta t_{\mathrm{evol}}\right),$$

$$\Delta t_{\mathrm{evol}}=\min \left(\Delta t_{\mathrm{SNII}},\frac{t_{\mathrm{age}}}{300}\right)\mathrm{yr},$$

$$\Delta t_{\mathrm{SNII}}=\frac{{\tau }_8}{N_{\mathrm{SNII}}}.$$

$$r_{\mathrm{cool}}=28.4E_{51}^{2/7}{⟨n⟩}^{-3/7}f(Z)\mathrm{pc},$$

Full text

††pagerange: Simulating the interstellar medium and stellar feedback on a moving mesh: Implementation and isolated galaxies–A††pubyear: 2019

1 Introduction

The $\Lambda$ Cold Dark Matter ($\Lambda$CDM) cosmological model has emerged as the most successful theoretical framework to explain the formation and evolution of cosmic structures (Planck Collaboration XIII 2016). This model assumes that galaxies in the present-day Universe assembled from a nearly homogeneous initial state during almost 14 billion years of evolution. Invisible, non-baryonic, dark matter constitutes the backbone for this structure formation process. Dark matter haloes assemble through gravitational collapse allowing baryons to cool and collapse effectively at their centres where the density is the highest, and the potential is deepest. The potential wells associated with dark matter haloes host the formation of the first generation of stars, the formation of the first nascent galaxies, and eventually host the formation of more massive galaxies, as we are accustomed to seeing in the local Universe (White & Rees, 1978).

Understanding and modelling the formation of galaxies within their cosmological context is a daunting task owing to the large dynamical range in space and time as well as the wide range of physical mechanisms involved. For instance, the physics of a single supernova (SN) explosion occurs on sub-parsec scales, but megaparsec scales are required to capture the gravitational collapse and assembly of a halo such as that of the Milky Way. Even larger scales are needed to properly capture the galactic environment including large-scale tidal torques and the past merger and interaction history of galaxies. In other words, accurately simulating galaxy formation simultaneously demands modelling large volumes, while still resolving small scales in order to capture the relevant galaxy formation physics.

All currently existing cosmological simulations, both full volume and zoom-ins, rely on sub-resolution (or sub-grid) models. Sub-grid models are a set of effective numerical prescriptions that are needed to capture physical processes that occur on scales at or below the resolution of the simulation (i.e. those processes that are unresolved). For example, cosmological galaxy formation simulations are not able to resolve the formation of individual stars, nor stellar evolution. Therefore, models are required that describe both the formation process of stars and their subsequent stellar evolution. Nearly all physical processes involved in galaxy formation are implemented to some degree as sub-grid models including radiative cooling, star formation, stellar evolution, SN feedback, Active Galactic Nuclei (AGN) feedback, etc. Ideal sub-grid models are independent of numerical resolution above a minimum threshold resolution such that the model numerically converges with increasing resolution. However, a well converged model will also yield no new information once the convergent resolution requirement is achieved. It is therefore important to separately evaluate simulation limitations imposed by numerical resolution versus limitations imposed by inadequate modelling of physical processes. For galaxy formation, the most important physical processes beyond gravity and hydrodynamics are the formation of stars and their interactions with the surrounding interstellar medium (ISM).

It is observed that star formation is inefficient with less than $20-25$ per cent of a galaxy’s predicted total baryonic budget turning into stars, and with fractions as small as $1$ per cent estimated for low mass dwarf galaxies (Conroy & Wechsler, 2009; Moster et al., 2013; Behroozi et al., 2019). The remaining baryons may be kept in the form of gas associated with the disc in the ISM, dispersed into the more extended circumgalactic medium, or remain outside of their dark matter haloes (Bregman, 2007; Putman et al., 2012, and references therein). Most galaxies must have an active mechanism to prevent gas from cooling and turning into stars, with the exception of very low mass objects where reionization may prevent gas from accreting into the gravitational potential well of the dark haloes (Benson et al., 2002; Simpson et al., 2013). Several decades of research in this area have identified stellar feedback as one of the leading physical processes to explain this inefficiency in $L_{*}$ galaxies and below (e.g. Stinson et al., 2006; Sales et al., 2010; Ostriker & Shetty, 2011; Hopkins et al., 2014a; Kimm & Cen, 2014; Marinacci et al., 2014; Vogelsberger et al., 2014; Murante et al., 2015; Wang et al., 2015; Schaye et al., 2015; Martizzi et al., 2016; Li et al., 2017; Kim & Ostriker, 2018). Resolving the multiphase structure of the ISM and a detailed modelling of stellar feedback are therefore key targets in modern studies of galaxy formation.

In fact, the ISM structure in galaxies is complex, with hot, warm and cold phases of gas coexisting and interacting (e.g. McKee & Ostriker, 1977). Following this complex structure becomes computationally more expensive when the numerical resolution increases. In particular, the modelling of the cold dense regions is challenging since the time-steps demanded for hydrodynamical calculations can become rather short. A common solution to overcome this obstacle is to impose an effective equation of state for the dense gas instead of directly resolving the individual gas phases and processes within the ISM. Effective equation of state models have been implemented in several different ways. A widely-used approach is to treat the ISM as two-fluid gas assumed to be composed of cold clouds embedded into a hot and diffuse medium (Springel & Hernquist, 2003; Agertz et al., 2011). An imposed polytropic relation between density and temperature $T\propto\rho^{\gamma}$ (see also Agertz et al., 2011; Dalla Vecchia & Schaye, 2012) acts as a pressurization of the medium to prevent unresolved numerical fragmentation. Models of this nature are usually accompanied by a restricted minimum gas temperature comparable to warm ionized gas, $T\sim 10^{4}\,{\rm K}$, below which the cooling and further collapse of the cold gas phase is not considered.

Such an approach for the ISM, although numerically stable and with desirable convergence properties (e.g. Springel, 2000; Marinacci et al., 2014) may have unwanted consequences for the resulting structure of the simulated galaxies beyond simply not resolving the ISM complexity. For instance, Benítez-Llambay et al. (2018) demonstrated that the scale height of dark matter dominated discs in similar models is set purely by the combination of sound speed in the mid-plane and the circular velocity. In this case the minimum scale height of the disc becomes independent of numerical resolution, or in other words, while one might increase the numerical accuracy of the simulation, the vertical structure of the discs will nonetheless remain unresolved, resulting in discs that are thicker and kinematically hotter than observed in low mass galaxies. The overall structure of galaxies in such models is overly-smooth on scales of a few hundred parsecs – the scales of molecular clouds – as a result of the missing complexity in the gas phases, severely limiting the predictive power of these simulations on small scales. In this case the appearance of the ISM will not change even if the numerical resolution is increased dramatically, as it is intimately linked to the effective equation of state chosen to treat the ISM.

Efforts to build models of resolved multiphase ISM gas structure must go hand-in-hand with the development of a stellar feedback model that can regulate and prevent the runaway collapse of the cold ISM gas. Several energy and momentum injection channels have to be considered including contributions from SN explosions at the end of the lifetime of massive stars, but also the impact of their early ionizing radiation on the surrounding media (Stinson et al., 2013; Sales et al., 2014; Gupta et al., 2016; Emerick et al., 2018; Haid et al., 2018). Detailed modelling suggests that the feedback energy budget is dominated by radiation while the direct momentum injection budget is comparable between SN and radiation (Agertz et al., 2013). Radiation impacts the density and temperature of gas through photoionization and radiation pressure. Moreover, dust in the ISM can boost the impact of radiation pressure via photon trapping where a single photon will undergo multiple scatterings. All these processes should be included in models that aim to resolve the multiphase ISM.

What remains less well understood in this picture is the efficiency with which the available energy and momentum couples to the surrounding gas. This coupling efficiency determines how important each feedback channel is in regulating star formation and launching galactic outflows. Although analytic arguments and some numerical simulations suggest that radiation pressure might be responsible for the launching of outflows of several hundred ${\rm km\,s^{-1}}$ observed in galaxies (Veilleux et al., 2005; Heckman et al., 2000; Pettini et al., 2001), other authors find the effect of radiation pressure to be more modest (Krumholz & Thompson, 2012, 2013) and the overall impact of radiation to be determined mostly by photoionization effects (Sales et al., 2014; Rosdahl et al., 2015).

Idealized radiative transfer experiments in Sales et al. (2014) indicate that although radiation pressure has the ability to eventually push the gas to galactic-outflow speeds, the timescales required to propel the gas to these speeds are much longer than those of photoionization. As a result, stars will heat their surrounding gas via photoionization, driving gas expansion and lowering the surrounding gas density. This process further limits the number of photon absorption events and therefore limits the radiation pressure momentum coupling. The results from these idealized experiments are in agreement with other works in the literature that follow star formation and feedback on cloud scales (Rosdahl et al., 2015; Dale, 2017; Walch et al., 2012; Walch & Naab, 2015; Kim et al., 2018) as well as on idealized disc patches (Kannan et al., 2018). The relevance of radiation pressure is therefore questionable. There is general consensus that radiation pressure can dissipate local cold clouds, but uncertainties remain about its contribution to launching galactic scale outflows (see Krumholz, 2018). Furthermore, since the impact of radiation pressure might be considerably larger for very dusty and high density star forming regions, such as those characteristic of Ultra Luminous Infrared Galaxies (ULIRGs, Murray et al., 2010; Ochsendorf et al., 2014), the actual environment is also important.

Other sources of stellar feedback remain still under-explored and their coupling to SN and radiation feedback is not well understood. Examples include the impact of cosmic rays (e.g. Recchia et al., 2016; Pais et al., 2018; Jacob et al., 2018), stellar jets (Frank et al., 2014), high-energy photons from X-ray binaries (Kannan et al., 2016b), runaway stars (Kimm & Cen, 2014), among others. As such, caution should be exercised when interpreting the results of the models or the meaning of parameters currently used in numerical simulations since they are, by construction, only approximations to the real physics driving the structure and evolution of galaxies.

Similarly, even when the physics is well understood, there are still ad-hoc choices made in the numerical implementation of these processes. For instance, one might consider the total energy input per SN event in $10^{51}\,{\rm erg}$ as a relatively well constrained quantity. However, the number of neighboring resolution elements used to distribute that energy or the corresponding momentum can significantly change the overall impact of that SN event (e.g. Sales et al., 2010; Dalla Vecchia & Schaye, 2012). It is therefore important to benchmark new ISM models against both, numerical resolution as well as the robustness of internal numerical assumptions.

Current state-of-the-art zoom-in simulations have reached a resolution where it becomes desirable to remove the effective ISM approach and instead model the multiphase gas structure and related stellar feedback in more detail. Even full-volume simulations have reached mass resolutions of the order of $\sim 10^{5}\,\rm M_{\odot}$ and therefore are able to resolve ${L^{\star}}$ galaxies with millions of gas and stellar resolution elements. For example, the IllustrisTNG TNG50 simulation achieves a $8.5\times 10^{4}\,{\rm M_{\odot}}$ mass resolution and a minimum gas softening of $74\,{\rm pc}$ throughout a simulation box of $\simeq 50\,{\rm Mpc}$ on a side (Nelson et al., 2019; Pillepich et al., 2019). Similarly, the galaxies in the zoom-in simulations of the Auriga project (Grand et al., 2017) contain several million of resolution elements to describe the baryonic component of Milky Way (MW)-like galaxies, and tens of thousands for the surrounding low mass dwarfs. These resolutions are therefore, in principle, sufficient to resolve more details within the ISM phase, and hence allow for physically improved models of stellar feedback and galactic outflows.

Nevertheless, for the moment these simulations still rely on an effective equation of state model combined with a hydrodynamically decoupled wind-model, where winds are launched using phenomenological prescriptions. These limitations motivate our work here, where we aim to construct a new ISM model that captures the multiphase gas structure alongside a more explicit local stellar feedback model. Recent successful attempts in this direction have been presented (e.g. Hopkins et al., 2011; Agertz et al., 2011; Agertz et al., 2013; Hopkins et al., 2014a, 2018c), demonstrating the possibility to achieve a more detailed treatment of the gas and stars in galaxies, within the cosmological context and with high numerical resolution, on scales of MW galaxies and below (Wetzel et al., 2016; Oñorbe et al., 2015; Fitts et al., 2017). Improved and refined ISM models should ideally be numerically well posed (this, for example, also means that they should be implementable both in particle and mesh codes) and be feasible to execute at the resolution typically achieved by current state-of-the-art simulations. The goal of these novel ISM and stellar feedback models is to move the scales of numerical closure to smaller scales within the ISM. In contrast, effective models do not resolve structure within the ISM even if the numerical resolution is increased. The numerical closure scale is therefore essentially set by the ISM. Our goal is to go below this scale, to provide a model with a smaller scale for numerical closure that is better matched to the resolution of existing and upcoming simulations.

Here we therefore present the SMUGGLE model, a novel ISM and stellar feedback model for the moving mesh code Arepo. The resolution requirements of this new model are such that it can yield converged results at the resolution of current state-of-the-art cosmological simulations. Our paper is organized as follows. Section 2 gives a description of the implementation of the most important physical processes occurring in the ISM that are considered in our model. Section 3 describes the initial conditions of the simulations carried out in this work, and presents their main features such as star formation rate and histories (and their sensitivity to resolution) and the ISM multiphase structure. Section 4 shows the properties of the gaseous outflows generated by the model and the salient features of the gas in the simulations, while Section 5 studies the relative contribution of the different channels to the global momentum output of feedback. Finally, Section 6 compares the numerical implementation of the physical processes considered in SMUGGLE to ISM and feedback models existing in the literature, while Section 7 gives a brief summary of our results.

2 Methods

In this Section we describe the main ingredients of our ISM and stellar feedback model. Figure 1 provides a high-level schematic overview of the physical processes that we consider in our model. These main processes are: $i)$ gravity and hydrodynamics – modelled in Arepo with a standard oct-tree algorithm (Barnes & Hut, 1986) for gravity and a finite volume solver over an unstructured Voronoi mesh that is allowed to move freely with the flow for hydrodynamics (Springel, 2010; Pakmor et al., 2016)– $ii)$ gas heating and cooling mechanisms (see Sec. 2.1), which are needed to describe the emerging multiphase structure of the ISM, $iii)$ a stochastic implementation for the formation of stellar particles (the sites from which feedback, metal and mass enrichment are originating) from the gas phase (see Sec. 2.2), and $iv)$ stellar feedback processes from three main channels, SNe (Sec. 2.3), radiation (Sec. 2.4) and stellar winds (Sec. 2.5), which are important for regulating star formation and shaping the properties of the ISM.

Figure 2 shows further details of the modelling of the feedback channels and mass and metal return illustrating their salient features – in particular how feedback energy and momenta, as well as the mass returned and the metals synthesized from stars are coupled to the gas – and referring to the key equations that we used for their implementation. In the subsections below we give a more detailed description of the physical processes that we consider in this model and their numerical implementation.

2.1 Cooling and heating

We first describe the gas cooling mechanisms implemented in our model. The thermal state of gas is modified by a cooling and heating network that models a primordial mix of hydrogen and helium undergoing two-body processes, such as collisional excitation, collisional ionization, recombination, dielectric recombination and free–free emission (e.g. Katz et al., 1996), Compton cooling off CMB photons (Ikeuchi & Ostriker, 1986), and photoionization from a spatial-uniform UV background (Faucher-Giguère et al., 2009). This primordial network is complemented by high-temperature (T\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}10^{4}\,{\rm K}) metal line cooling based on a self-consistently updated gas metallicity field.

Net metal cooling rates are tabulated as a function of temperature, gas density and redshift based on cloudy calculations (Ferland et al., 1998) in the presence of UV background radiation while assuming solar metallicity. The net metal cooling rates are added to the primordial network, scaled by the gas total metallicity relative to solar. The details of these procedures are described in Vogelsberger et al. (2013). In the remaining part of this Section we only describe the implementation of the additional processes that we consider specifically for the SMUGGLE model.

An important aspect of our model is the creation of low-temperature (T\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10^{4}\,{\rm K}) gas. Additional cooling processes must be considered to reach these low temperatures. We allow gas to cool via low-temperature metal line, fine-structure and molecular cooling processes, which enables the gas to reach $\sim 10\,{\rm K}$ under a chemical and ionization equilibrium assumption. Cooling gas to low temperatures is particularly important for reaching the high densities necessary for the formation of a molecular gas phase and, ultimately, for star formation in our model (see Section 2.2). We account for these cooling processes by implementing the fit to the cloudy cooling tables presented in Hopkins et al. (2018c)

[TABLE]

here $T$ is the gas temperature, $n_{\rm H}$ (in ${\rm cm}^{-3}$) is the hydrogen number density, and $Z$ and $Z_{\sun}$ are the gas and the solar metallicities. To model self-shielding of the gas at high densities (n\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}10^{-3}\,{\rm cm}^{-3}), we multiply equation (1) by $1-f_{\rm ss}$, where $f_{\rm ss}$ is the self-shielding factor computed using the parameterization presented in Rahmati et al. (2013). More specifically, the factor varies as a function of redshift and gas density as

[TABLE]

where $n_{\rm H}$ is the hydrogen number density (in ${\rm cm}^{-3}$) of each cell and the parameters ($n_{0}$, $\alpha_{1}$, $\alpha_{2}$, $\beta$, $f$) encode the redshift dependence (see Rahmati et al., 2013, their Table A1). The self-shielding factor is also used to suppress the ionization and heating rates entering the primordial network and the normalization of the UV spectrum entering the determination of the contribution to the cooling of metal lines. Above $z=6$ no self-shielding correction is considered. More details about self-shielding can be found in Vogelsberger et al. (2013). Since the test simulations presented in Sec. 3 follow the evolution of an isolated, Milky Way-type galaxy (i.e. they are not cosmological), equation (2) adopts the redshift zero parameterization.

We also take into account cosmic ray and photo-electric heating. Both processes are thought to be important for the thermal balance and stability of the cold ($T\sim 50\,{\rm K}$) and warm ($T\sim 8000\,{\rm K}$) phases of the ISM (see e.g. Field et al., 1969; Wolfire et al., 1995). Cosmic ray heating is implemented using the form (see Guo & Oh, 2008)

[TABLE]

where $\tilde{n}_{e}$ is the electron number density in units of the hydrogen atom number density $n_{\rm H}$ (in ${\rm cm^{-3}}$) and $e_{\rm CR}$ is the cosmic ray energy density that is parametrized by

[TABLE]

in such a way that it becomes progressively less important for low density gas (see also Hopkins et al., 2018c). In equation (3) the first term represents the heating through hadronic losses while the second term is related to Coulomb interactions between cosmic rays and the gas. By putting together the last two equations the heating rate in dense gas is $\Lambda_{\rm CR}\simeq-2.37\times 10^{-27}n_{\rm H}^{-1}\,{\rm erg\,s^{-1}\,cm^{3}}$, a typical value used in ISM studies (see Tielens, 2010, equation 3.31).

For the photoelectric heating – i.e. the emission of electrons from dust grains, in particular polycyclic aromatic hydrocarbons (PAHs), caused by the photoelectric effect due to the interstellar radiation field – we employ the rate from Wolfire et al. (2003, equation 19)

[TABLE]

[TABLE]

where $\tilde{e}_{\nu}^{\rm pe}$ is the photon energy density, normalized to the Milky Way units, $\tilde{e}_{\nu}^{\rm pe}\equiv e_{\nu}^{\rm pe}/(3.9\times 10^{-14}\,{\rm erg\,cm^{-3}})$, which we set to one in our calculations, and $\Phi_{\rm PAH}$ is a factor that incorporates the uncertainties in the interaction rates between atoms and dust grains in the molecular regime. Observations of carbon density ratios in diffuse clouds suggest that $\Phi_{\rm PAH}=0.5$ (Jenkins & Tripp, 2001). Therefore, we fix this parameter to that value, which is also the fiducial value adopted in Wolfire et al. (2003). The factor $(Z/Z_{\odot})$ links the dust abundance in the gas to the metal abundance, which is a common choice in this type of modelling.

2.2 Star formation implementation

Cold dense gas will eventually be converted to stars. To convert gas cells into star particles we follow a probabilistic approach (e.g., Springel & Hernquist, 2003). Star particles, representing a coeval simple stellar population following a Chabrier (2001) initial mass function, are created stochastically according to the probability derived from the star formation rate ($\dot{M}_{\star}$) of each gas cell computed as follows

[TABLE]

where $\epsilon$ is an efficiency factor that we set to $0.01$, in line with observational determinations (see e.g. Krumholz & Tan, 2007), and $t_{\rm dyn}$ is the gravitational dynamical time of any given gas cell defined as

[TABLE]

with $\rho_{\rm gas}$ and $M_{\rm gas}$ indicating the gas density and mass, respectively. We note that other choices for the star formation efficiency parameter $\epsilon$ are possible. For instance, in Hopkins et al. (2018c) it is shown that the exact level of star formation is independent of the value of $\epsilon$, provided that star formation is indeed feedback regulated. We have not explicitly checked this in our model and instead we have chosen a value more in line with observational determinations (see also Smith et al., 2018). Also, Agertz et al. (2013) report that variations of a factor of two in the normalization of the Kennicutt-Schmidt relation Kennicutt (1998) are possible when passing from $\epsilon=0.01$ to $\epsilon=0.1$. We defer to future work for a detailed exploration of the sensitivity of the model to the values of the parameters (see Table 3 for a list of the key parameters of the model).

A first criterion for a gas cell to be eligible for star formation is to have a density above a specific density threshold $\rho_{\rm th}$. We set this parameter to $100\,{\rm cm^{-3}}$, which is in the range of average densities ($10^{2}-10^{3}\,{\rm cm}^{-3}$) of giant molecular clouds (Ferrière, 2001). Other choices of the density threshold above this value do not significantly impact the star formation history and the establishment of the Kennicutt-Schmidt relation in the simulations that we are going to present in Section 3. As for the star formation efficiency value, we postpone a more complete analysis of this aspect to future work.

We also add an additional criterion based on the virial parameter $\alpha$ to restrict star formation only to gravitationally bound regions, i.e. those that are prone to the onset of gravitational collapse (Semenov et al., 2017). Namely we compute for each cell $i$ (Hopkins et al., 2018c)

[TABLE]

and allow star formation only for cells in which $\alpha<1$, which indicates that the gas is not able to overcome gravitational collapse via gas motion and thermal support. In equation (9), $\boldsymbol{v}$ is the gas velocity, $c_{s}$ the sound speed, $\Delta x$ the cell size (defined as the radius of a sphere having the same volume as the cell), $\rho$ is the gas density and $G$ the gravitational constant, and $||\nabla\otimes\boldsymbol{v}_{i}||^{2}$ is defined as

[TABLE]

Once the star formation rate is determined, the stellar mass that would be formed in the current time-step $\Delta t$ by the cell $i$ is

[TABLE]

where $M_{i}$ is the mass of gas in the cell and $\dot{M}_{\star}$ the star formation rate determined in equation (7). Instead of forming stellar mass in a continuous way as prescribed by the previous equation, a stochastic approach, in which a gas cell is turned into a collisionless star particle with a probability consistent with the local SFR, is adopted. This probability can be easily found in equation (11) by noting that in a time $\Delta t$ a fraction $p\equiv 1-\exp(-\dot{M}_{\star}\Delta t/M_{i})$ of the cell mass will be converted to stars. Hence, to decide whether the conversion occurs, a uniformly distributed random variable $p^{\star}$ in the interval $[0,1]$ is extracted and compared with $p$. A cell is converted to a star particle if $p^{\star}<p$. The phase-space variables of the newly created star particle are inherited from those of the converted gas cell.

We note that with the procedure just described all the created star particles would have the mass of the parent gas cell. It is sometimes however desirable to change this behaviour. For instance, an equal mass for all stellar particles can be desirable. This can easily be accommodated in this scheme by multiplying $p$ with the ratio between the mass of the cell and the target mass for the star particle. Another typical case is to enforce that the stellar mass is within a predetermined factor from a target value, in analogy with the standard Lagrangian refinement criterion for gas cells adopted in Arepo (see Vogelsberger et al., 2012). As in these cases usually not all the gas cell mass is turned into a star, the conservative variables in the gas cell must be changed accordingly. We note that the scheme to generate stellar particles just described is also the standard approach in Arepo (see Vogelsberger et al., 2013).

Finally, our model can also base the computation of the star formation rate on molecular gas only, given the strong observational correlation between these two quantities (Bigiel et al., 2008; Leroy et al., 2013). This is achieved by multiplying $M_{\rm gas}$, appearing in equation (7), by the factor (see McKee & Krumholz, 2010; Krumholz & Gnedin, 2011):

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

However, for sufficiently dense gas, as it is the case for the simulations described below, $f_{\rm H_{2}}\sim 1$. $f_{\rm H_{2}}$ is enforced to be positive, meaning that if its value goes below zero then it is set to zero and no star formation occurs. This can represent a problem in the case of zero metallicity gas – a situation that occurs in cosmological simulations, in which the gas starts from a metal free condition – then no star will ever be formed. To circumvent this problem, a minimum metallicity floor of $Z=10^{-5}Z_{\odot}$ is imposed, thus allowing star formation to occur in gas with zero or very low metallicity values.

2.3 Feedback from supernovae

SN feedback is thought to play a crucial role in regulating star formation and the resulting ISM structure through the injection of momentum in the gas (e.g. Agertz et al., 2011; Agertz et al., 2013; Aumer et al., 2013; Guedes et al., 2011; Hopkins et al., 2011, 2014a; Kimm & Cen, 2014; Stinson et al., 2006, 2013). SN momentum injection is responsible for driving turbulence (Martizzi et al., 2016) and giving rise to galactic-scale outflows (Li et al., 2017). However, implementing SN feedback effectively is not trivial, especially in simulations with limited resolution (see e.g. Scannapieco et al., 2012). The primary difficulty associated with modelling SN feedback is to properly capture the early energy-conserving expansion phase of a SN remnant, the so-called Sedov-Taylor phase. During this phase, the momentum imparted to the ISM gas is generated by an overpressurized central gas bubble that expands into, sweeps up, and accelerates ambient material. The radial momentum of an expanding SN shell can reach a factor of $\sim 10$ above the initially injected momentum (see e.g. Martizzi et al., 2015). It is this boosted momentum injection that is largely responsible for the regulation of star formation (Ostriker & Shetty, 2011) or the driving of galactic winds (Kim & Ostriker, 2018). Therefore, an effective SN implementation that appropriately captures the total momentum and energy injection is important, especially if the early Sedov-Taylor blast phase expansion is not explicitly resolved.

Many approaches have been adopted to increase the efficiency of SN feedback. Common approaches are temporarily turning off the cooling of gas (Stinson et al., 2006) and selectively (Murante et al., 2010) or stochastically (Dalla Vecchia & Schaye, 2012) heat the gas affected by feedback. Both approaches allow sufficient momentum to build up to enable the launching of outflows. Additionally, some schemes inject momentum while hydrodynamically decoupling gas elements for a short period of time to facilitate the launching of a galactic wind (Springel & Hernquist, 2003; Oppenheimer et al., 2010; Vogelsberger et al., 2013; Davé et al., 2017; Valentini et al., 2017; Pillepich et al., 2018). In this Section, we describe our implementation of SN feedback that aims for a local treatment of energy and momentum injection into the ISM, regulating the buildup of cool/cold gas and generating gas outflows self-consistently.

2.3.1 SN energy and momentum budget at injection

We assume that the total energy injected in a single SN event is given by

[TABLE]

where $E_{51}=10^{51}\mathrm{ergs}$ and $f_{\mathrm{SN}}$ is a model parameter that encodes the SN feedback efficiency ($f_{\mathrm{SN}}\equiv 1$ in all our test runs). We assume that the blast wave velocity at explosion is

[TABLE]

such that the momentum carried by the blast wave at explosion is given by

[TABLE]

and $M_{\mathrm{SN}}$ is the ejecta mass per SN.

The stellar particles in our simulations do not represent single stars, but rather a stellar population. As such, for any given time step, $\Delta t$, multiple SN events may occur. To inject the proper amount of energy and momentum per timestep we must first compute the number of SN events and the total associated ejecta masses. We thus define for each stellar particle and for each time step $\Delta t$

[TABLE]

and

[TABLE]

where $N_{\mathrm{SNII}}$ and $N_{\mathrm{SNIa}}$ are the number of type II and type Ia SN in the time step and $M_{\mathrm{SNII,tot}}$ and $M_{\mathrm{SNIa,tot}}$ are the total ejecta mass associated with these events, respectively. The total number of SNII events and their associated ejecta mass are found by integrating over the initial mass function of the stellar particle as (Vogelsberger et al., 2013)

[TABLE]

[TABLE]

In equations (22) and (23), $M_{\star}$ is the mass of the star particle at birth, $\Phi(M)$ is the Chabrier (2001) initial mass function, $M(t)$ is the mass of a star that leaves the main sequence at an age $t$ and $f_{\mathrm{rec}}(M,Z)$ the amount of mass given back to the ISM by stellar evolution (which depends on stellar mass and metallicity, see e.g. Portinari et al., 1998). We set the minimum main sequence mass for a type II SN explosion to $8\,{\rm M_{\odot}}$ and adopt an IMF upper limit of $100\,{\rm M_{\odot}}$.

The calculation of SNIa events proceeds in a slightly different way. We parametrize the temporal distribution of SNIa events using a delay time distribution (DTD) and derive the number of type Ia SN events as (see Vogelsberger et al., 2013)

[TABLE]

The form of the DTD is poorly constrained. In this work we define it as follows

[TABLE]

where $\tau_{8}=40\,{\rm Myr}$ approximates the main sequence life time of an $8\,{\rm M_{\odot}}$ star, $N_{0}=2.6\times 10^{-3}$ SN ${\rm M_{\odot}^{-1}}$, $s=1.12$ and $\Theta$ is the Heaviside function that parameterizes the delay between the birth of the stellar population and the first SNIa event (see Maoz et al., 2012). This functional form agrees with theoretical models that link the SNIa rates to the orbital energy and angular momentum loss rate due to gravitational wave emission (Greggio, 2005) and is also consistent with the previous implementation of type Ia mass and metal return to the ISM adopted in cosmological simulation of galaxy formation with Arepo (Vogelsberger et al., 2013). Finally, each SNIa releases the same amount of ejecta, $M_{\mathrm{SNIa}}\simeq 1.37\,{\rm M_{\odot}}{\rm SN^{-1}}$ (Thielemann et al., 2003), and therefore the total mass return follows

[TABLE]

Finally, we would like to stress that since for stellar ages less than $t_{8}$ there are no SNIa events, the two types of SNe are occurring as two distinct, temporally-separated channels.

Since SN explosions are discrete events, our model aims to mimic their discrete nature using the following procedure (see also Hopkins et al., 2018c). We first impose a time-step constraint for each stellar particle based on its age (i.e. evolutionary stage) as

[TABLE]

where

[TABLE]

and $t_{\rm age}$ is the age of the star and

[TABLE]

In the previous equation $N_{\rm SNII}$ is the expected number of SN events for the stellar particle over $\tau_{8}$. For our choice of the IMF (Chabrier, 2001) and of the lower mass of type II SN progenitor ($8\,{\rm M}_{\odot}$) about $10^{-2}\,{\rm SN\,{\rm M}_{\odot}^{-1}}$ are expected.

Both $\tau_{8}$ and $N_{\rm SNII}$ in equation (29) are computed self-consistently by our stellar evolution model for each star particle, as the main sequence life time also depends on the metallicity of the star (see Vogelsberger et al., 2013). Imposing the limit derived in equation (27) ensures that the expectation value for the number of SN events per timestep is of the order of unity111There is of course a dependence on the mass of the stellar particle and, for the least resolved of our simulations, more than one SNII event may occur per time step. These are allowed (and dealt with self-consistently by our Poisson sampling) in order not to impose a too restrictive time step condition. For higher resolution simulations multiple SN explosions do not typically take place.. To determine whether a SN event takes place, we integrate equations (22) and (24) to calculate the expected number of SNe ($\lambda$) over the time step, $\Delta t_{\star}$. This value is taken as the expectation value of a Poisson distribution, which is sampled to obtain the actual number of discrete SN events per time step (usually either zero or one since $\lambda\ll 1$ at our fiducial resolution and because of the time step limit). More details about the Poisson sampling of SN events are given in Appendix A.

2.3.2 Accounting for PdV work in the Sedov-Taylor phase

The initial energy conserving Sedov-Taylor phase lasts until the SN approaches the cooling radius where the post-shocked gas reaches temperature of about $10^{6}\,{\rm K}$ at which point radiative losses become important. This cooling radius is given by (see, e.g. Cioffi et al., 1988; Hopkins et al., 2014b, 2018b)

[TABLE]

where $n$ is the average gas density within $r_{\rm cool}$ in ${\rm cm^{-3}}$, and $f(Z)$ is a function of gas metallicity defined in equation (33) below. If we do not fully resolve the cooling radius (i.e., if SN mass/energy is being returned on scales larger than $r_{\mathrm{cool}}$), as it is usually the case in our simulations, then we need to explicitly account for momentum that is generated by the $P{\rm d}V$ work by the hot post-shocked gas during the adiabatic Sedov-Taylor expansion phase of the SN blast.

Several schemes have been devised to account for this process (see e.g. Kimm & Cen, 2014; Gatto et al., 2015; Kimm et al., 2015; Kim & Ostriker, 2017; Rosdahl et al., 2017; Hopkins et al., 2018b; Smith et al., 2018). In our model, we do so by boosting the momentum imparted to each gas cell $i$ influenced by SN feedback as

[TABLE]

where $\tilde{w}_{i}$ is a weight function partitioning the energy and momentum injection among gas cells (see equations 35 and 36), $m_{i}$ is the mass of the gas cell, $\Delta m_{i}=\tilde{w}_{i}(M_{\mathrm{SNII,tot}}+M_{\mathrm{SNIa,tot}})$ and $p_{t}$ is the so-called terminal momentum, i.e. the final value of the momentum of the SN blast at $\approx$ the cooling radius, when the evolution of the blast wave transitions from the Sedov-Taylor phase to a momentum-conserving phase. The terminal momentum per SN is defined as

[TABLE]

and this expression is determined from high-resolution simulations of individual SN blast waves (see Cioffi et al., 1988, equation 4.7). In equation (32) $\left<n_{\rm H}\right>$ is the local gas hydrogen number density around the star (determined in an SPH-like fashion, see equations 34 and 43) and

[TABLE]

with $\left<Z\right>$ being the average SPH-weighted gas metallicity around the star particle computed in the same aperture as the local gas density. $Z_{\odot}=0.0127$ is the value we adopt for solar metallicity (Asplund et al., 2009). We note that to get the value for the terminal momentum, equation (32) is multiplied by $N_{\rm SNII}+N_{\rm SNIa}$, the number of SN events occurring at any given time step.

2.3.3 Supernova energy and momentum coupling

We spread the SN energy and momentum injection over a number of nearest gas particles using weight functions. To do so we identify for each star particle a predefined effective number of neighbours:

[TABLE]

Here, $h$ is a search radius (coupling radius), $W$ is the standard cubic spline SPH kernel (Monaghan & Lattanzio, 1985) and $\mathbf{r}_{i}$ and $\mathbf{r}_{s}$ are the position vectors of the $i$-th gas neighbour and of the star particle, respectively. Equation (34) defines an implicit relation for the search radius $h$, which is solved iteratively until the predetermined number of neighbours $N_{\rm ngb}$ is found. In our tests ‘ $N_{\rm ngb}=64$ and variations of $\pm 1$ (effective) neighbour are allowed.

Additionally, to avoid artificial numerical effects described in Sec. 2.3.4 we also impose a feedback limiter, namely, a maximum radius for the coupling of energy and momentum from the stars, $R_{\rm SB}$. The introduction of this second spatial scale requires a special differentiated handling of the weights for, on one hand, the redistribution of mass and metallicity from the stars (which must by definition be conserved) and, on the other hand, for the distribution of energy and momentum (which might not be occurring beyond $R_{\rm SB}$, see Sec. 2.3.4). In practice, we define the effective scale for coupling as $h_{\rm coupling}=\min(h,R_{\rm SB})$, where $R_{\rm SB}$ is introduced in equation (37) as the feedback limiter. If $h$ is larger than the feedback limiter radius, two sets of weights are computed for the two different values of the mass ($h$) and of the feedback ($h_{\rm coupling}$) coupling radii. Conversely, if $h$ is smaller than $R_{\rm SB}$, the two sets of weights have the same values. We now describe the procedure that we use to compute such weights.

Once $h_{\rm coupling}$ has been determined, weights are defined in such a way that each gas cell within $h_{\rm coupling}$ receives SN energy and momenta222Metals and mass returned from SNe are also coupled to the ISM, but with a set of weights computed within the unlimited coupling radius $h$. proportionally to the fraction of the $4\pi$ solid angle that it covers as seen from the star position. Namely, we define a weight $w_{i}$ such that (see also Hopkins et al., 2018c; Smith et al., 2018)

[TABLE]

where $|\mathbf{r}_{i}-\mathbf{r}_{s}|$ is the distance between the gas cell and the star particle and $A_{i}$ is the area of the gas cell defined as $A_{i}\equiv\pi\Delta x_{i}^{2}$ (recall that $\Delta x_{i}$ is the cell size). To ensure that the correct amounts of energy and momentum are imparted each cell receives a fraction

[TABLE]

of the amount of $p_{\mathrm{SN,tot}}$ and $E_{\mathrm{SN,tot}}$ determined in equations (21) and (20). Momenta are directed radially away from the star position. Feedback momentum and energy are injected in the reference frame in which the star particle is at rest, and the final momentum and total energy of the gas cell are later transformed to the appropriate values for the frame of reference adopted in the simulation. We highlight that gas particles receiving feedback from supernovae are never hydrodynamically decoupled after mass, energy and momentum are injected.

2.3.4 Maximum size for the SN coupling radius

The pressure associated with an expanding SN blast wave decreases with time or, equivalently, the radius reached by the shock wave during its expansion. The impact of a SN blast wave will have on the surrounding ISM gas becomes negligible (and actually the SN remnant ceases to exist as an individual entity and blends into the ISM) when the pressure associated with the SN ejecta becomes comparable to the ambient ISM pressure or alternatively and the velocity of the shock wave propagating into the ISM becomes equal to the turbulent velocity of the ambient gas.

At the resolution achieved in our simulations, however, we are unable to model single stars. Rather, each stellar particle formed represents a coeval single stellar population, or stellar cluster, in which multiple SN events occur as a result of the presence of associations of massive (namely OB) stars. This clustered SNe will inflate large bubbles in the galaxy’s ISM known as superbubbles (Weaver et al., 1977) with a size of a few hundreds of parsec across, thus influencing the ambient ISM on scales larger than those of an individual SN remnant.

The size of a superbubble can be estimated as (see Mac Low & McCray, 1988; Weaver et al., 1977)

[TABLE]

where $n_{\rm H}$ is the hydrogen number density (we have assumed a hydrogen mass fraction of $0.76$ for the estimation), $\Delta t$ is the life time of the OB association ($\sim 40\,{\rm Myr}$ for a lower SN progenitor mass of $8\,{\rm M}_{\odot}$) and $L_{\rm SN}$ is the mechanical SN energy, which can be defined as

[TABLE]

where $N_{\rm SN}$ is the total number of SNe in the single stellar population and $E_{\rm SN}$ is the energy of an individual SN event.

Beyond the distance estimated in equation (37) the energy/momentum injected by SNe will be unable to have a strong impact on the ISM properties, although the superbubble might be able to break out and vent material outside the star-forming disc when its size exceeds a few disc scale heights (see Mac Low & McCray, 1988, and the discussion of Figs 9 and 10 below). In general, typical values of $R_{\rm SB}$ are several hundreds of parsecs, depending on gas ambient density and numerical resolution. For instance, assuming $n_{\rm H}=1\;{\rm cm}^{-3}$, the corresponding superbubble radius is in the range $\simeq 0.3-1\;{\rm kpc}$ for the numerical resolution of the runs presented later in Sec. 3 (high to low resolution, respectively). We note, however, that the numerical convergence of our model improves when $R_{\rm SB}$ is kept constant instead of adjusting this value on the fly. Therefore, a flag in our feedback implementation allows the model to either follow equation (37) and adopt the limiter to density and resolution or, alternatively, to keep the value constant. Motivated by the better convergence behaviour, we choose the latter for the current work, adopting $R_{\rm SB}=0.86\;{\rm kpc}$ as our default value. We have tested that varying this in the range $0.3-3\;{\rm kpc}$ does not have a significant impact on our results.

2.4 Radiative feedback from young massive stars

Radiation, especially from young and massive stars, has an important impact on both, the thermal and the dynamical state of the ISM gas. It can alter its ionization state and therefore the gas temperature by photoionizing it. Moreover, each absorbed photon imparts additional momentum to the gas via radiation pressure. The imparted momentum can be further enhanced by multiple photon scattering and can then become relevant at high densities, such those found in giant molecular clouds, the birth clouds of stars.

This process is particularly important because of its timing. Radiative feedback can be responsible for the dispersal of such clouds before any SN goes off (Murray et al., 2010; Lopez et al., 2011; Walch et al., 2012), thereby rendering SN feedback more effective due to the reduced gas densities in which SNe subsequently explode. This process sometimes goes under the name of early stellar feedback (see e.g. Stinson et al., 2013). Some studies have also proposed radiative feedback as the key mechanism behind the launching of galactic-scale outflows (Murray et al., 2005), but this aspect is still strongly debated (Krumholz & Thompson, 2012, 2013). We describe our implementation of the radiative feedback, and specifically of the photoionization of gas and treatment of radiation pressure, in the subsections below.

2.4.1 Photoionization

Young massive stars are a copious source of ionising radiation. This radiation can in turn impact the ionisation state of the surrounding gas, leading to the emergence of Hii regions. In our model we capture the formation of such regions as follows.

We define the ionising photon rate emitted by a stellar particle as

[TABLE]

that is the luminosity of the star divided by the average photon energy emitted above $13.6\,{\rm eV}$ and we estimate the star particle luminosity given its mass $M_{\star}$ by assuming a mass-to-light ratio $\gamma_{\star}$. We choose $\langle h\nu\rangle=17\,{\rm eV}$ and $\gamma_{\star}=10^{3}\,{\rm L_{\odot}}/{\rm M_{\odot}}$ as our fiducial values, corresponding to the peak emission of a black-body Planck spectrum with temperature $T\sim 40,000\;\rm K$ (Rybicki & Lightman, 1986), consistent with massive OB stars.

The mass of gas that can be photoionised by a young stellar particle (i.e. the mass within the so-called Strömgren radius) is usually (much) smaller than the mass of gas contained in within $h_{\rm coupling}$, the scale at which the model couples the feedback energy. We therefore carry out the photoionization probabilistically. Cells are assigned a probability of being photoionized $p=n_{\star}/(\alpha_{\rm rec}n_{\rm H}^{2}V)$, where $\alpha_{\rm rec}\simeq 2.6\times 10^{-13}\,{\rm cm^{3}s^{-1}}$ is the hydrogen recombination rate, $n_{\rm H}=X\rho/m_{\rm p}$ is the average hydrogen number density of the cell – with $X$ being the hydrogen mass fraction and $m_{\rm p}$ the proton mass – $V$ is the cell volume and $n_{\star}=\tilde{w}_{i}N_{\star}$ is the rate of ionizing photons injected into the cell. In other words, the probability is given by comparing the number of recombinations expected in the cell to the total photon number emitted by the source scaled by the same solid angle weighting scheme adopted to model the supernova feedback.

A random number $p^{\prime}$ is then selected from a uniform distribution in the range [0,1] and cells where $p^{\prime}<p$ are tagged for photoionization. For those, we impose a temperature floor $T_{\rm phot}=1.7\times 10^{4}\,{\rm K}$ and disable their radiative cooling for a duration $t_{\rm off}$ equal to the star particle time step. The temperature floor is consistent with a gaseous medium made of hydrogen and photoionized by $\langle h\nu\rangle=17\,{\rm eV}$ photons. In addition to the probability criterion outlined in the previous Section, a gas cell has to fulfill two further requirements to be declared eligible for photoionization: (i) its $u_{\rm therm}$ must be less than $1.2\times u_{\rm phot}$ (i.e. the thermal energy per unit mass corresponding to $T_{\rm phot}$), (ii) must be able to cool (i.e $t_{\rm off}=0$). If any of these two conditions apply, the gas particle is considered to be already (photo)ionised – either because they are at a high enough temperature or are kept ionised by another stellar particle – and no action is taken. Finally, we point out that,while $t_{\rm off}$ is initialized to the star time step ($\Delta t_{\star}$), its value is updated according to the duration of the cell (individual) time step, so that cooling is disabled precisely for a duration equal to $\Delta t_{\star}$. The heating of eligible gas cells due to photoionization will generate some over-pressurization of the gas compared to their colder (and mostly neutral) surroundings. This will lead to the expansion of the photoionized regions, which in the spirit of the model should be capturing the expansion of the ${\rm H_{\sc II}}$ regions around young stellar systems, thus imparting momentum to the gas. However, please note that this momentum generation originates from a different physical mechanism than a direct transfer of momentum from the radiation field to the gas (i.e. radiation pressure), which is accounted for in SMUGGLE as described below.

2.4.2 Radiation pressure

Radiation coming from young stars has an additional impact on the dynamical state of the gas because of the pressure it can exert on it. This represents a source of momentum, which is important in dense, optically-thick regions, allowing a preprocessing of the gas environment in which SNe will subsequently explode (e.g. Agertz et al., 2013; Hopkins et al., 2014a; Stinson et al., 2013).

In our model, we account for this momentum source from young stars by considering an injection of momentum around each star particle of the form

[TABLE]

where $\Delta p$ is the total momentum injected by the star over the time step $\Delta t$, $L_{\star}$ is the star luminosity (see equation 39), $c$ the speed of light and $\tau_{\rm IR}=\kappa_{\rm IR}\Sigma_{\rm gas}$ is the optical depth of the gas to infrared radiation, which takes into account the multiple scattering of infrared photons in dense gas. $\kappa_{\rm IR}$ and $\Sigma_{\rm gas}$ are the opacity in the infrared band and the gas column density, respectively. We set $\kappa_{\rm IR}=10(Z/Z_{\odot})\,{\rm cm^{2}g^{-1}}$ (Hopkins et al., 2018c), although its precise value is uncertain and smaller $\kappa_{\rm IR}$ (by a factor of $2$) are also reasonable (Agertz et al., 2013). This would impact the total amount of momentum imparted by radiation. However, in our simulations radiation does not appear to be the dominant feedback channel (see Section 5). The total radiation momentum determined in equation (40) is coupled to the gas contained within $h_{\rm coupling}$ in the same way as it is done for SN feedback (see equations 35-36) and directed radially away from the star position. The same maximum value for the coupling radius adopted for SN feedback is adopted in both photoionization and radiation pressure calculations.

The last quantity to be determined to completely specify the model is an estimate of the column density of the gas $\Sigma_{\rm gas}$ within $h_{\rm coupling}$. We adopt a standard Sobolev approximation and compute for each star particle $j$ a Sobolev length as

[TABLE]

in which $\rho_{j}$ and $\nabla\rho_{j}$ are the gas density and the associated gradient determined by adopting a standard SPH approach (see equation 43 below and also Hopkins et al. 2018c). The column density is computed as

[TABLE]

and $\langle\rho_{s}\rangle_{j}$ is given by

[TABLE]

where $W$ is the cubic spline SPH kernel (see also equation 34), $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ are the position vectors of the $i$-th gas neighbour and of the star particle, respectively, $m_{i}$ the mass of the gas cell and $h$ appears in equation (34).

2.5 Feedback from OB and AGB stellar winds

Stars also contribute to feedback through stellar winds. For our purposes two classes of stars are important: massive, short-lived OB stars and asymptotic giant branch (AGB) stars. The former class helps in pre-processing the gas around young star particles before the onset of SNe (Matzner, 2002; Krumholz & Matzner, 2009) – and in particular in dispersing the dense gas clouds from which stars are born – a role similar to radiative feedback. In this sense they represent an additional early stellar feedback channel.

Moreover, AGB winds are an attractive source of feedback on their own, as they continue to act at later times than radiation and are associated to older stellar populations. Their mass return can further contribute to fuel later episodes of star formation in galaxies with their own subsequent injection of associated energy and momentum. It is worth noting that a detailed accounting of the mechanical power carried by winds of massive stars shows that their momentum injection rate per unit stellar mass formed is comparable to that of the SNe (Agertz et al., 2013), with the difference that most of this power is released before the first SN explosion (t\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}5\,{\rm Myr}) and rapidly declines thereafter. Therefore, stellar winds represent a feedback channel that should be explicitly modelled in galaxy formation simulations, given its potential in enhancing the effect of subsequent SN explosions because of the reduction in the ambient density.

We implement these two contributions in our model as follows. In a first step we compute the mass loss due to these two processes. For OB stars we parametrize the cumulative mass loss per unit stellar mass as

[TABLE]

where (see Hopkins et al., 2018c)

[TABLE]

In the previous expressions $\tilde{Z}$ is the metallicity of the stellar particle in solar units (and capped to 1.5), time is expressed in Myr, and for $t>100\,{\rm Myr}$ $M_{\rm closs}=h(100)$. Over a time step $\Delta t$ the mass lost by a star particle with an initial mass equal to $M_{\star}$ is $M_{\rm loss}=M_{\star}[m_{\rm closs}(t+\Delta t)-m_{\rm closs}(t)]$. Moreover, OB stellar winds leave each stellar particle with the same chemical composition as the particle itself. We consider OB stars wind injection only for stellar life times smaller than the one associated with the minimum mass of type II SNe progenitors ($8\,{\rm M_{\odot}}$). Finally, AGB stars are treated by the stellar evolution model already present in Arepo (Vogelsberger et al., 2013), so we consider those values for their contribution to the metal and total mass return.

To compute the energy and momentum injection from stellar winds we set (Hopkins et al., 2018c)

[TABLE]

where $M_{\rm loss}$ is the mass loss determined above and

[TABLE]

The total momentum injection rate is thus determined as

[TABLE]

The mass, metallicity, momentum and energy are injected as in the SN case in the rest frame of the star and the final values of momentum and total energy of the gas cell are then transformed back into the reference frame of the simulations. The same maximum values to the coupling radii also apply. Finally, we would like to note that differently from SNe, the feedback injection from stellar winds is a continuous process and is implemented accordingly. However, OB wind injection and the associated feedback is only performed if the returned mass over a given time step is larger than $10^{-4}$ times the mass $M_{\star}$ of the stellar particle at birth. If this is not the case, the mass loss is accumulated until the first time step at which the mass return threshold for this channel is reached. Our implementations of cooling, star formation, and stellar feedback are inspired by those in the FIRE-2 simulations and largely follow, as highlighted by the references throughout this Section, the methods laid out in (Hopkins et al., 2018b), but with a number of modifications as summarized in detail in Section 6.

3 Isolated Galaxies

3.1 Initial conditions and model setup

We test our new feedback implementation for isolated late-type galaxies. Our main goal is to explore the ability of the new model to regulate star formation, and to investigate its impact on the properties of the ISM of the galaxy.

We construct the initial conditions following the technique described in Springel et al. (2005, see also ). We set up an equilibrium compound galaxy model representative of the Milky Way consisting of a dark matter halo, a bulge and a stellar and gaseous discs. Both the halo and the bulge are modelled with a Hernquist (1990) profile. The stellar and gaseous discs are exponential in the radial direction. While the stellar disc follows a ${\rm sech}^{2}$ distribution in the vertical direction, the vertical profile of the gaseous disc is computed self-consistently to ensure hydrostatic equilibrium at the beginning of the simulation. The initial gas temperature is set to $10^{4}\,{\rm K}$.

For the structural parameters of these components we selected the same values that were used to construct the Milky Way galaxy model presented in Hopkins et al. (2012). This model has a total mass of $1.6\times 10^{12}\,{\rm M_{\odot}}$ of which $\simeq 1.5\times 10^{10}{\rm M_{\odot}}$ are contained in the bulge, $\simeq 4.73\times 10^{10}\,{\rm M_{\odot}}$ and $\simeq 9\times 10^{9}\,{\rm M_{\odot}}$ in the stellar and gaseous discs, respectively, resulting in a disc gas fraction of $f_{\rm gas}\simeq 10$ per cent within $R=8.5\,{\rm kpc}$. The gas in the disc has a metallicity equal to the solar value of $0.0127$ (see Asplund et al., 2009). A summary of the structural parameters is listed in Table 1.

Three different resolution levels (low, intermediate and high) are created, differing in the number of element used to sample each galaxy component. At the lowest resolution level we put $N_{\rm bulge}=10^{5}$ particles in the bulge, $N_{\rm disc}=4\times 10^{5}$ particles in the stellar disc and $N_{\rm gas}=10^{5}$ gas cells in the gaseous disc. Each time we increase the resolution level, we do so by increasing the particle number by a factor of 8. The parameters defining each resolution level are listed in Table 2. For efficiency reasons, the simulations model the dark matter halo as a static gravitational field. Although modelling the dark matter in a static way is an approximation, in the sense that the dark matter halo will not react to changes occurring in the system, the results are only marginally affected by this choice. For instance, the final stellar mass formed in our simulations with a static versus live halo agree within $\approx 10$ per cent.

Finally, we enclose the system in a cubic volume with a side of $857~{}{\rm kpc}$ in which we inserted, with the procedure described in Section 9.4 of Springel (2010), a background grid consisting of low-density gas cells to avoid vacuum boundary conditions. The coarsest level of resolution of the background grid corresponds to that of a grid composed by $16^{3}$ cells. Cells are allowed to be refined (de-refined) whenever they are above (below) a factor of two (half) of a predetermined target mass equal to the average gas mass before the introduction of the grid. The de-refinement procedure is not applied to the background grid cells, which are identified as having a volume larger than 10 per cent of the average gas cell volume in the gaseous disc. The resulting initial conditions are then evolved for $\sim 1~{}{\rm Gyr}$ with the schemes described in Section 2. In Table 3 we list the fiducial values adopted for the numerical parameters associated with the different physical processes modeled in this work.

3.2 Star formation rates and the Kennicutt-Schimdt relation

An essential aspect of any stellar feedback model is its ability to regulate star formation, which otherwise will be a runaway process occurring on much shorter time-scales (or at much larger rates) than observed in “normal” star-forming galaxies. In this Section, we will investigate the ability of our new model to achieve star formation regulation.

Figure 3 (left-hand panel) presents the star formation history (SFH) for our Milky Way run. The solid lines in the figure give the evolution of star formation in the three different resolution runs: low (blue), intermediate (red) and high (green). For comparison, we also include two runs with the Springel & Hernquist (2003) model one with the basic model implementation (i.e. without galactic winds, dark orange) and one including decoupled galactic winds (purple, at intermediate resolution) as well as a run without any feedback source (black, at intermediate resolution). The Springel & Hernquist (2003) run including winds has a weak outflow strength as compared to the parametrizations used in the Illustris (Vogelsberger et al., 2014) and Auriga (Grand et al., 2017) projects. This is to avoid that very strong stellar feedback prematurely terminates star formation given that in our isolated setup the galaxy lacks a fresh supply of gas due to cosmological gas accretion. All the lines show the instantaneous star formation rate, computed as the sum of the star formation rates of each individual gas cell at a given time $t$.

Comparisons of the colored lines with the no-feedback run in Fig. 3 indicate that the star formation rate in the simulated Milky Way reaches a self-regulated state due to the presence of feedback, avoiding a large peak in the beginning of the simulation and the subsequent decline shown in the no-feedback case (black curve). This decline is due to the rapid gas consumption during the early stages of the no-feedback run. Instead, our new feedback model shows a remarkably constant and stable SFR around $\sim 3\,{\rm M_{\odot}}\,{\rm yr^{-1}}$, in agreement with observations of this galaxy type. There is a small residual dependence on resolution at early times (see for example $t\leq 0.2$ Gyr), but passed the initial stages and settling of the disc the three resolution levels agree quite well, which is also clear from the total integrated stellar mass formed, shown on the right panel.

Most importantly, the predicted level of star formation broadly agrees with the results obtained with the Springel & Hernquist (2003) model, which is the most commonly adopted model in arepo thus far and which relies on an effective equation of state. Achieving a similar averaged SFR together with self-regulation but without having to resort to the imposition of an equation of state to treat the gas is one of the main accomplishments of our new model. The inclusion of decoupled galactic winds in the Springel & Hernquist (2003) run contributes to further lowering the SFR of that model by about 35 per cent. Finally, we note that the inclusion of feedback suppresses star formation by factors as large as $\sim 10$, especially at early times.

Figure 3 (right-hand panel) is the cumulative version of the previous plot and shows the total amount of stellar mass formed by the different simulations as a function of time. In line with the results discussed for the instantaneous star formation rate, by the end of the simulations the stellar mass formed in our new feedback model agrees within $\sim 13$ per cent for the different resolution levels, within $\simeq 18$ per cent for the total mass produced by the Springel & Hernquist (2003) model without winds and shows an overall $\sim 3-4$ suppression factor compared to the no-feedback run. The addition of kinematically decoupled winds to the Springel & Hernquist (2003) model pushes the difference in cumulative mass with SMUGGLE  to within a factor of two.

A related assessment of the performance of our model to regulate star formation can be made by analyzing the Kennicutt-Schmidt relation (Kennicutt, 1998) connecting the projected density of gas, $\Sigma_{\rm gas}$ to the star formation rate density (SFRD). This is shown in Fig. 4, with the left panel highlighting different spatial locations for our simulated galaxy at $t=0.7$ Gyr and the right panel showing the disc-averaged values (up to the radius enclosing 90 per cent of the star formation rate) at different $7$ Myr time intervals. Blue, red and green colors indicate the different resolution levels as before, together with orange/purple and black used for the Springel & Hernquist (2003) without/with galactic winds and no-feedback runs, respectively. The simulated relations have been obtained by considering the gas in a region centred on the disc and defined by the conditions $R<50\,{\rm kpc}$, $|z|<5\,{\rm kpc}$. For the left panel, selected gas cells have been binned into concentric annuli with width $0.1\,{\rm kpc}$ over which the gas and star formation rate surface densities have been estimated.

Our simulations including feedback compare well with the original observational sample by Kennicutt (1998) for normal star-forming galaxies (empty squares), whereas the non-feedback run strongly over-predicts star formation. To guide the eye, we include with dashed and dash-dotted lines the best fit of the relation obtained by Kennicutt (1998, dashed) and Bigiel et al. (2008, dash-dotted). The dotted lines show constant depletion times of $10$, $1$ and $0.1\,{\rm Gyr}$ from bottom to top. As a result of feedback, star formation is self regulated and we find good agreement between our new model predictions and the observations, regardless of the numerical resolution.

In other words, the resulting star formation in our model is rather inefficient with gas being converted into stars with a typical time scale of $1-2\,{\rm Gyr}$, rather than on a free fall time. Interestingly, in the left panel, there is also evidence of a break in the simulated relations at \,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10\,{\rm M_{\odot}}{\rm pc^{-2}} (Bigiel et al., 2008). We note that the Springel & Hernquist (2003) model successfully reproduces these trends, albeit completely free of or with minimal scatter and show a less pronounced downturn of the relation below $\sim 10\,{\rm M_{\odot}}{\rm pc^{-2}}$. These trends hold independently of the inclusion of galactic winds in the model, which has only a minimal impact on the recovered slope and normalization of the simulated Kennicutt-Schmidt relation. Instead, our new feedback and star formation prescriptions show a level of scatter that compares well with observational estimates. This increased scatter is not surprising given the fact that the Springel & Hernquist (2003) model was specifically formulated to explicitly reproduce the median observed relation.

The large impact of feedback on the depletion time-scale $\tau$ can also be seen in more detail in Fig. 5, following the same color coding as before. $\tau$ is defined as the ratio between the gas mass and the instantaneous SFR in the galaxy, and acts as a proxy for the typical time over which a galaxy will consume its gas reservoir due to star formation processes. To compute both quantities a spatial cut defined by the conditions $R\leq 40\,{\rm kpc}$ and $|z|\leq 1.5\,{\rm kpc}$ has been adopted. Depletion time scales in the new feedback runs are up to $10$ times larger compared to the no-feedback case and reasonably converged with resolution, especially once the runs have settled into self-regulation (t\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}0.3-0.4 Gyr). The resulting timescales after feedback is included agree well with observational estimates (Sancisi et al., 2008; Bigiel et al., 2011; Kennicutt & Evans, 2012; Leroy et al., 2013; Putman, 2017).

Interestingly, results from the SMUGGLE model are in line with those from the Springel & Hernquist (2003) model without winds, where $\tau\simeq 3\,{\rm Gyr}$ and is relatively constant over time (dark orange line). The inclusion of galactic winds in the Springel & Hernquist (2003) model (purple) further pushes up $\tau$ above $\simeq 4\,{\rm Gyr}$. This is caused by a combination of the ISM pressurization due to the adoption of an effective equation of state and of the gas ejection due to the galactic wind presence. We hasten to add that the agreement between the depletion time scale in SMUGGLE and in the no wind implementation of the Springel & Hernquist (2003) model is by no means straightforward, as in our model star formation suppression and self-regulation is achieved self-consistently by following the different gas phases, whereas in the Springel & Hernquist (2003) approach this is entirely due to the pressurization of gas resulting from the implemented effective equation of state for all gas elements with densities above the star formation density threshold (which represent an equilibrium two-phase ISM model where the hot gas phase is heated by supernova energy input). For the typical parameters of our simulations the star formation density threshold of the Springel & Hernquist (2003) model is $\sim 0.1\,{\rm cm}^{-3}$. As we will see in Sec. 2.5, the generation of winds in our model also follows self-consistently from the gas treatment whereas they need to be explicitly added in the Springel & Hernquist (2003) model.

3.3 ISM structure

As discussed above, an important feature of our ISM model is the natural emergence of a multiphase ISM during the evolution of the simulation. This subsection is devoted to the presentation of the properties and the structure of this multiphase ISM. Fig. 6 shows the phase diagram of all the gas cells in the high resolution Milky Way simulation at $t=0.7\,{\rm Gyr}$. The figure presents the phase diagrams in terms of a temperature versus density plot and is constructed as a two-dimensional histogram. Darker shades correspond to increasingly larger gas mass values.

To guide the eye, we highlight the multiphase structure of the gas in Fig. 6 through dashed lines and labels to indicate the different regions. A large fraction of the ISM is composed of warm gas at a fairly constant temperature of $\sim 10^{4}\,{\rm K}$, covering a density range from about $10^{-2}\,{\rm cm}^{-3}$ to $\sim 10^{2}\,{\rm cm}^{-3}$ (d), and a colder phase, starting off at $10^{-2}\,{\rm cm}^{-3}$ and reaching low (\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}100\,{\rm K}) temperatures and high densities (e and f). Indeed, only gas which is sufficiently dense ($n>10^{-2}\,{\rm cm}^{-3}$) can effectively self-shield from the UV background and effectively cool and become dense enough to form stars (e). All the material present in the region of the phase diagram in between these two loci is thermally unstable, and populates these region of the diagram because of stellar feedback.

Two other features of the diagram are interesting as well. The first is the continuation of the warm gas phase locus to higher density (\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}1\,{\rm cm}^{-3}) and with approximately the same $\simeq 10^{4}\,{\rm K}$ temperature (c). This results from the photoionization of the gas described in Section 2.4.1. The second one is the presence of gas at low density (\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10^{-2}\,{\rm cm}^{-3}) and relatively high ($\sim 10^{5}-10^{6}\,{\rm K}$) temperatures, coincident with the hot gas phase (a). This is largely the result of the ejection of hot gas from within the star forming disc due to the momentum and energy delivered by SN feedback (b).

A more quantitative determination of the contribution of each gas phase to the ISM is presented in Fig. 7, which shows the gas density probability density function (PDF) for our high resolution Milky Way galaxy at $t=0.7\,{\rm Gyr}$. The left and right panels correspond to the mass- and volume-weighted versions of this figure, respectively. Gas belonging to the disc is selected in the region $R<40\,{\rm kpc}$ and $|z|<1.5\,{\rm kpc}$. The PDF is plotted for all gas cells in the disc region (solid black lines) and for the cold ($T<2\times 10^{3}\,{\rm K}$, dashed blue lines), warm ($2\times 10^{3}\,K<T<4\times 10^{5}\,{\rm K}$, dash-dot yellow lines) and hot ($T>4\times 10^{5}\,{\rm K}$, dotted red lines) ISM phases, to better highlight the different contribution to the PDFs and their relative importance.

The multiphase nature of the ISM can be identified in both representations in Fig. 7 through the broad distribution of possible gas densities. A typical $n\sim 1\,{\rm cm^{-3}}$ marks the transition between the warm and cold dominated gas phases. At high densities ($n\geq 10^{2}\;\rm cm^{-3}$), cold gas dominates in both distributions. Across all densities, a significant fraction ($\simeq 32$ per cent) of the total mass is in this phase, although it occupies very little volume ($\simeq 0.5$ per cent of the total). For densities \,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10^{4}\,{\rm cm}^{-3} the hot phase, heated by SN feedback starts to appear, being more dominant in the volume-weighted distribution than in mass-weighted one, especially at very low densities (\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10^{-3}\,{\rm cm}^{-3}). However, the global contribution of the hot phase is always very small both in a mass- ($\simeq 0.1$ per cent) and in a volume-weighted ($\simeq 1.6$ per cent) sense.

4 Gas dynamics and galactic-scale fountain flows

An essential aspect of the SMUGGLE model is the natural generation of gaseous outflows, which are locally correlated with star formation and are also fully coupled to the hydrodynamic equations. Figure 8 (left-hand panel) displays the measured gas outflow (blue solid line) and inflow (green dashed line) rates onto the disc as a function of time for the high resolution Milky Way simulation. To compute the rates only gas at a distance $R\leq 40\,{\rm kpc}$ and contained in two slabs of thickness $300\,{\rm pc}$ and located at a distance $\pm 2\,{\rm kpc}$ from the midplane of the galaxy has been considered. For each snapshot, after removing from the gas the average vertical gas velocity $\left<v_{z}\right>$, each cell contained in the slab contributes to the total rate as $mv_{z}/\Delta z$, where $v_{z}$ is the outflow/inflow velocity with respect to the disc plane.

Our implementation of stellar feedback is able to launch outflows at a rate of a few ${\rm M_{\odot}}\,{\rm yr}^{-1}$ throughout the entire simulation. The outflow rate is larger at early times, reaching a maximum of $\sim 10\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ at $t\simeq 0.35-0.4\,{\rm Gyr}$, and subsequently declining. Particularly noticeable is the bursty appearance of the mass outflow rate that cycles through its maxima and minima on a time scale of $\sim 100\,{\rm Myr}$. Comparing these outflows rates to the instantaneous SFR of the galaxy333The SFR per bin is determined within the same galactocentric distance cut and at heights $|z|<2\,{\rm kpc}$. we compute the mass loading factor parameter $\beta$ (right-hand panel), which results in $\beta\sim 1$ approximately constant with time. The solid red and dashed black lines show the raw mass loading factor obtained from the snapshots and a smoothed version estimated over a time-scale of $50\,{\rm Myr}$, respectively.

The left panel of Fig. 8 also includes the evolution of gas inflows with time (green dashed curve). Peaks in the values of the inflow rates are usually correlated, albeit offset, with respect to peaks in the outflow rate or $\beta$ hinting at a galactic fountain origin of the accreting gas. Therefore, ejected gas in our runs is not entirely lost for the galaxy but instead cycles back in a galactic fountain flow, a cycle that has been shown to be important for the metallicity and angular momentum build up of the disc and hot corona (e.g. Marinacci et al., 2010; Übler et al., 2014; Christensen et al., 2016; Anglés-Alcázar et al., 2017).

As a result of capturing this complexity in the gas phases, the structure of the gas in the galaxies simulated with SMUGGLE is also visually very different from previous models in arepo. We show images of the evolution of the gas column density (Fig. 9) and density-weighted gas temperature (Fig. 10) of the fiducial Milky Way galaxy run with our new feedback model at high resolution (first two rows) and of the same galaxy run with the Springel & Hernquist (2003) model (third and fourth row). Each figure shows face-on and edge-on projections taken at different times ($t=0.14,0.35,1\,{\rm Gyr}$, from left to right). Each panel is $40\,{\rm kpc}$ across and in projection depth, with a total number of $512^{2}$ pixels for a spatial resolution of $\approx 80\,{\rm pc}$.

Stars tend first to form along the spiral pattern in the gaseous disc, where the density is highest. Particularly noticeable is the effect of stellar feedback in these regions, where cavities of low-density gas are carved by the combined action of radiation feedback and SN explosions. As time progresses, the gas distribution becomes increasingly structured, with dense gas regions in which star formation is active that are surrounded by cavities of lower density gas. Low-density bubbles are propelled from the star-forming disc into the halo region by the previous generation of stars. The whole process reaches self-regulation at later times. The feedback-driven gaseous outflows are more easily seen in the edge on projections, in which the amount of gas above and below the disc – and also the disc thickness – tends to increase as a function of time.

In contrast, none of these features are present in the simulations run with the Springel & Hernquist (2003) model. In this case, the use of an effective equation of state is able to regulate star formation at the price of losing information on the ISM structure. In particular, the appearance of the ISM is smoother and no SN cavities are visible. Furthermore, since the momentum injection from stellar feedback is not directly modelled (we recall that in the most basic implementation of the Springel & Hernquist 2003 model supernova feedback acts only as a heating source for the hot phase of the ISM), gaseous outflows are not generated, as can be clearly appreciated in the edge-on projections. This is the reason why in cosmological applications the Springel & Hernquist (2003) model has to be complemented with a prescription for hydrodynamically decoupled galactic winds (see also their Section 4).

Similar trends are also visible in the temperature plots (Fig. 10), for which it is worth noticing that the densest gas regions feature the lowest temperatures, whereas the gas cavities carved by feedback are filled with shock-heated hot ($\sim 10^{6}\,{\rm K}$) gas, as naturally expected in this multiphase gas media. This hot gas is eventually vented outside the disc in a superbubble blowout. These features are smoothed over in the case of the Springel & Hernquist (2003) simulation, which also develops high temperature gas region in the disc midplane as a result of the imposed effective equation of state.

The development of a gas circulation cycle over galactic scales can be further inspected in the kinematical face-on and edge-on maps shown in Fig. 11, constructed following the same specifications from the previous two figures. The highest vertical velocities \,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}150\,{\rm km~{}s^{-1}} are reached in the low-density gas cavities generated by stellar (and in particular SN) feedback, which is best appreciated in the face-on projection (top row) of Fig. 11. This is consistent with the picture that the gas originally contained in such cavities has been displaced by the momentum injected by star particles and accelerated to velocities compatible with its ejection above or below the disc mid-plane. The fact that a large-scale outflow is generated in the SMUGGLE simulation is clearly visible in the edge-on projections, in which most of the gas that is a few ${\rm kpc}$ above or below the disc plane has outflow velocities in excess of $150\,{\rm km~{}s^{-1}}$.

However, the kinematics of the gas is more complex than a simple ejection occurring across the whole star-forming disc. The outflow is strongest near the galactic centre, where star formation is most effective, but at late times it can be seen that gas is also accreting onto the disc. As argued before in the discussion of Fig. 8 it implies that outflows do not completely escape the galaxy, but rather create a galaxy-wide gas circulation compatible with a galactic fountain flow (Bregman, 1980). The accreted (recycled) gas thus becomes the fuel that sustains star formation at later times.

These results, combined with the discussion of Fig. 8, suggest that the produced galactic-scale winds are still rather weak overall ($\beta\sim 1$) and do not seem to be able to produce baryonic loss from the halo, at least not at the halo mass scale ($M_{\rm tot}\sim 10^{12}\,M_{\odot}$) and redshift ($z\sim 0$) probed by these calculations. However, cosmological simulations indicate that strong outflows that lead to baryonic loss are needed to produce sensible galaxies in the $\Lambda$CDM cosmology (e.g. Schaye et al., 2015; Muratov et al., 2015; Pillepich et al., 2018; Tollet et al., 2019). Testing the ability of this model of doing so below $M_{\rm tot}\sim 10^{12}\,M_{\odot}$ (i.e. the dwarf galaxy mass scale) and in cosmological applications is therefore needed to assess its potential in producing realistic galaxies in the full cosmological context. We plan to address this important issue in future work.

5 Contributions of the different feedback channels

The spatial distribution of newly formed stars in our model sets the location and conditions under which feedback is injected. Thus, we start our analysis by showing in Fig. 12 face-on (top) and edge-on (bottom) stellar light maps for the fiducial Milky Way. Light is assigned to stellar populations assuming a Chabrier (2001) IMF using the starburst99 stellar population synthesis models. Simple line-of-sight attenuation is included assuming a constant dust-to-metal ratio with a Milky Way attenuation curve. Younger stellar populations (identified by their blue colors) are not uniformly distributed, but rather clustered in regions of recent star formation activity. The gas distribution (identified by the darkened features) is also distinctly non-uniform. The edge-on panels show clear dust lanes, with a large fraction of the dense, star-forming gas being constrained to a disc plane. Newly born stars are mostly born within a very thin disc plane owing to the low (\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}10^{3}\,{\rm K}) temperatures reached by the dense ISM gas from which they are born.

In SMUGGLE, consistent with other models in the literature, the clustering of newly formed stars is the result of tying the active star formation regions to cold and relatively high-density (n\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}100\,{\rm cm}^{-3}) gas. Indeed, in the simulations performed with the Springel & Hernquist (2003) model, that features a lower ($n\sim 0.1\,{\rm cm}^{-3}$) star formation density threshold, the distribution of younger stellar particles is more uniform, albeit dependent on gas density. The presence of galactic outflows in the Springel & Hernquist (2003) model does not significantly alter this result.

The clustering of star formation might have an impact on the injection of supernova feedback and the response of the gas to it. Some studies (e.g. Buck et al., 2019) also suggests that the strength of the clustering signal of young stellar particles can be used to derive constraints on the most appropriate value of the star formation density threshold to use in cosmological simulations. The detailed study of the clustering properties of star formation, their dependence on model parameters (particularly on the star formation density threshold) and the implications for stellar feedback is outside the scope of the present paper and will be addressed in future work.

As stars are born and stellar evolution takes over they inject energy and momentum back to the surrounding media. We keep track of the cumulative momentum injected into the gas by all feedback channels (radiative, OB/AGB winds and SNe). This enables some interesting comparisons among the different channels and allows us to gauge their relative contributions to the star formation self regulation and outflows. Additionally, as in the case of radiative and SN feedback, one can evaluate their “boost” factors, i.e. the ratio between momentum actually injected into the gas compared to the one expected.

In Fig. 13 we display the ratios between cumulative momenta injected into the gas by different feedback channels for our high resolution simulation. The different panels show the ratios between OB/AGB winds and radiative feedback (left-hand panel), OB/AGB winds and SNe (central panel), and radiative feedback and SNe (right-hand panel). The plots are computed as cumulative distribution of the ratio of the total momenta injected by each stellar particle in their life-time in each feedback channel (blue solid line). A variable number of constant-width bins has been used to compute the distribution of all the ratios. In particular, the bin width for ratios involving OB/AGB winds is 0.1, whereas it is fixed to 1 for the SN and radiative feedback case. For plots including SN momentum, only stars older than $5\,{\rm Myr}$ are considered. This is to exclude the nominal delay of $\sim 5\,{\rm Myr}$ between the creation of the star particle and the first SN event in order to not bias the analysis. We note, however, that because of our discrete sampling of the SN events there is a residual dependence of the tail of the distribution on the adopted age cut (see also the discussion below).

Several interesting trends are noticeable. First, the left panel of Fig. 13 indicates that stellar winds (OB/AGB) and radiation feedback provide approximately the same amount of momentum to the gas within a $1\sigma$ dispersion444This is computed in analogy with a Gaussian distribution by considering half the range between $16-$th and $84-$th percentile of the distribution. of $\simeq 1.66$. Second, as indicated by the middle panel, SN momentum clearly dominates over that from OB/AGB winds. Indeed, the momentum imparted by stellar winds is on average about 10 per cent or less than the one delivered by SNe and the distribution around the mean is quite narrow ($1\sigma$ dispersion is $\simeq 0.09$). We note, however, that a small minority of stars have a tail extending to momentum ratios up to p_{\rm AGB}/p_{\rm rad}\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}15. We have checked that this is caused by young stars (age \,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle<}}{{\sim}};}40\,{\rm Myr}), for which there is a delay in the onset of the first SN explosions and the energy and momentum injected by stellar winds dominates for a short period of time from the creation of the star particle and before the first SN event.

We now turn our attention to the ratio between the momentum injected by radiation compared to SN, with the distribution shown in the rightmost panel of Fig. 13. The behaviour is similar to the stellar winds analyzed before, with radiation being largely subdominant to the SN channel. Radiation can deliver on average a factor of $\sim 0.71$ times the momentum of SNe with a $1\sigma$ dispersion of $\simeq 0.64$. However, the estimate for the mean is biased by the long tail of stellar particles extending up to p_{\rm SN}/p_{\rm rad}\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}100. By using the median of the distribution, the momentum imparted by the radiation is merely $4$ per cent of the one injected by SNe. Once again the extended tail at large $p_{\rm rad}/p_{\rm SN}$ values is due to very young stars for which the momentum budget is dominated by radiation because there has not been enough time for SN to set in. Moreover, even by using a more conservative cut on stellar ages ($40\,{\rm Myr}$) leaves a tail of particles for which values of $p_{\rm rad}/p_{\rm SN}\sim 15$ can be reached. This signals that the momentum imparted by the radiation can locally play an important role, especially in high gas density regions. However, if one looks at the overall budget, this analysis suggests that SNe dominate the momentum feedback input in our simulations.

The results discussed above appear to be quite robust to numerical resolution at which the simulations are performed. Although the median values of the imparted momentum ratios between different feedback channels change with resolution, the conclusions on their relative contribution to the overall momentum output hold also at lower resolution levels. In particular, the median $p_{\rm AGB}/p_{\rm rad}$ ratio increases from 0.36 to 1.53 passing from low to high resolution, the median $p_{\rm AGB}/p_{\rm SN}$ ratio is almost constant across levels (decreasing only from 0.1 to 0.07 from low to high resolution), and the median $p_{\rm rad}/p_{\rm SN}$ ration decreases from 0.23 at low resolution to 0.04 at high resolution. From these findings it is evident that, for increased resolution, the radiation pressure channel becomes on average less important for the overall momentum output, whereas the ratio between the momentum imparted to the gas by OB/AGB winds and supernovae remains approximately constant. The radiation pressure is the most sensitive feedback channel to resolution changes, because its effects are more strongly influenced by the detailed structure of the ISM given the explicit dependence of the imparted momentum on the gas column density (see eq. [40]).

The coupling and propagation of feedback in the ISM can cause the effective momentum achieved in the gas to exceed the one originally injected, a phenomenon referred to as “boost” factor. Figure 14 shows the cumulative distribution of these boost factors for SN (top) and radiation (bottom) momentum in the high resolution Milky Way simulation. As in Figure 13, this is computed cumulatively with respect to the total momenta injected by each stellar particle during their entire life-time. Logarithmically-spaced bins are used to compute the histograms for both quantities.

In the SN case, the boost factor $\beta$ indicates the ratio between the momentum imparted to the gas after the Sedov-Taylor phase and the value injected at the explosion. We find on average a $\sim 7.1$ momentum enhancement for SNe. The dispersion around the mean value is $\simeq 4.6$, implying that larger values of the boost factor are possible, but that the probability of very large values is rather small – only approximately less than $1$ per cent of stars have \beta\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}15 – and that, therefore, $\beta$ has only a very weak dependence on the environmental conditions of the gas, i.e. density and metallicity (although metallicity is constant and equal to $Z_{\odot}$ in our simulations), surrounding the star particle. We would like to note that these values for the momentum boost factor are in line with previous high-resolution simulation studies (Cioffi et al., 1988; Kim & Ostriker, 2015; Martizzi et al., 2015; Thornton et al., 1998; Walch & Naab, 2015).

On the other hand, for the radiative feedback case, the boost factor accounts for multiple infrared scattering, which increases the amount of momentum transferred to the gas via radiation pressure by $(1+\tau_{\rm IR})$. This quantity is proportional to the gas optical depth to infrared radiation, which is in turn depending quite strongly on the physical conditions (and in particular the column density) of the gas. Therefore, contrary to the SN case, the momentum boost is more affected and gets significantly larger in denser gas. On average, the momentum boost in our simulations is $\simeq 16.37$, with a median of $\simeq 3.86$, which are reasonable values based on previous work (Agertz et al., 2013). Compared to the SN case there is a larger scatter ($\simeq 16.14$) because of the (stronger) environmental dependence. For the environments analyzed here, namely a Milky-Way like disc at $z=0$, values for the $(1+\tau_{\rm IR})$ boost factor can reach $>100$ in very rare cases, but no more than $\sim 5$ per cent of the stars have (1+\tau_{\rm IR})\,\lower 3.01385pt\hbox{;\stackrel{{\scriptstyle\textstyle>}}{{\sim}};}100.

Based on this analysis, we conclude that all feedback channels included here play a role in regulating star formation, with boost factors of about $7.1$ and $16.37$ that seem to be common for SN and radiation pressure, respectively. However, the momentum input and generation of galactic-scale fountain flows for the scale of Milky Way galaxies analyzed here is largely dominated by the SNe channel only. We acknowledge, however, that according to the gas physical conditions, which may vary in galaxies of different mass and redshifts, the relative contribution to the global feedback energy and momentum input among these channels may vary, a topic that we plan to investigate in future work.

6 Comparison to previous work

The SMUGGLE ISM and feedback model is not the first attempt towards a more explicit description of the ISM physics in simulations of galaxy formation. Several of the processes included here (for example the parameterization for molecular cooling, heating process in the ISM and stellar winds) are in common with other models previously introduced in the literature. However, our adopted numerical implementation is different in several aspects. We will first focus on a comparison with the cases of the FIRE-2 (Hopkins et al., 2018c) and the Agertz et al. (2013) ISM treatments, being among the closest match to our newly developed model. The main differences can be summarized as follows:

• •

hydrodynamic solver: Arepo adopts a moving-mesh finite volume scheme versus mesh-less finite mass employed by gizmo (FIRE-2) and the AMR technique in ramses (Agertz et al., 2013). This makes a fundamental difference in the way feedback energy and momentum are distributed, especially for updating the fluid conserved variables in our model;

• •

feedback energy and momentum: one of the differing aspects of our model is the explicit justification of a maximum coupling radius for the SN feedback energy and momentum (see Sec. 2.3.4), which comes through superbubble evolution results. Instead, in (Hopkins et al., 2018c) the maximum coupling radius is treated as a fixed parameter and set to $\simeq 1-2\,{\rm kpc}$. The distribution of feedback among the neighboring cells can also differ because of a different choice in the computations of weights or in the number of neighbours. For instance, in the computation of the weights $\tilde{w_{i}}$ (see equation 36), we do not take into account a vector/tensor correction term (see e.g. Hopkins et al., 2018b; Smith et al., 2018) that would ensure linear momentum conservation to machine precision. This choice was made for reasons of computational efficiency/simplicity and it would be interesting to compare the outcome of the two approaches. Instead, in Agertz et al. (2013) energy and momentum are spread over the 26 neighbours surrounding each stellar particle (basically the closest neighbours of the cell containing the stellar particle, which, however, is excluded from the feedback loop). In our model we choose an SPH-like kernel with 64 cells, but the weights may be adjusted based on the maximum coupling radius, as discussed in Sec. 2. We also note that in order to have an efficient feedback loop regulating star formation, Agertz et al. (2013) may also adopt a non-thermal energy reservoir approach (based on the turbulent feedback model by Teyssier et al., 2013) that in practice is comparable to a delayed cooling set-up. We do not include such energy source here;

• •

radiation feedback: our model adopts a stochastic photoionization scheme based on the calculation of the total number of ionizing photons that are then distributed to the neighbouring cells surrounding a star particle based on a solid angle weighting scheme (see equations 35 and 36). The probability of ionization is then computed as the ratio between the number of ionizing photons and that of the recombinations occurring within the cell. This is different from the FIRE-2 implementation in which cells are sorted by distance and become photoionized in such order until the total available numbers of ionizing photons is accounted for. Moreover, we do not consider the long-range radiation component in our model and we also use a prescribed mass-to-light ratio and a fixed time over which radiative feedback occurs. This is a different choice with respect to FIRE-2 in which these processes are included as the output of starburst99 calculations. Note that we follow photoionization and radiation pressure effectively as a single-band calculation and do not track separately the 5-band approximation (photo-ionizing, photo-electric, UV, optical, IR) as done in FIRE-2. On the other hand, the model in Agertz et al. (2013) does not explicitly take into account the impact of photoionization, but includes radiation pressure. Their boost due to multiple scattering of the infrared photons is tracked only for the lifetime of the birth clouds of the stars ($3$ Myr, instead of the $5$ Myrs considered here based on the age of the OB stars) and the infrared opacity in Agertz et al. (2013) is half of our value ($\kappa_{\rm IR}=5(Z/Z_{\odot}){\rm cm}^{2}\,{\rm g}^{-1}$)555A multiplicative pre-factor of 2 is used in Agertz et al. (2013) for both direct and infrared momentum injection to account for grid smearing effects; see their equation (5)..

• •

Photo-electric heating: SMUGGLE does not include an explicit treatment of this processes, whereas FIRE-2 does. Comparison to radiative transfer calculations have shown that the inclusion of photo-electric heating effects does not significantly impact the dynamics and global properties of galaxies (Hopkins et al., 2018a), however, it may play a relevant role in the properties of the ISM (e.g. Forbes et al., 2016). We defer to future work a detail study of this in SMUGGLE.

• •

AGB winds: while we adopt the FIRE-2 parameterization for the wind speed, the mass loss due to the AGB phase is directly taken from our stellar evolutionary model (see below);

• •

stellar evolution: besides differing on the adopted metal production yields, the IMF adopted in SMUGGLE and in Agertz et al. (2013) follows a Chabrier (2001) distribution666However, the metal enrichment outcome of the stellar evolution can differ, as only oxygen and iron are taken as proxy for gas metal enrichment in (Agertz et al., 2013)., whereas FIRE-2 uses a Kroupa (2001) IMF instead. There are also subtle differences in the mass ranges for type II supernova progenitors. For instance, in our model we take $8\,{\rm M}_{\odot}$ as lower limit, whereas $6\,{\rm M}_{\odot}$ is used in FIRE-2. This has a direct influence on the number of SNII events per unit stellar mass formed, which is nearly twice as high in FIRE-2 compared to our model. Furthermore, this also translates in an average mass per SN event of $\sim 16\,{\rm M}_{\odot}$ in SMUGGLE  and $\sim 10\,{\rm M}_{\odot}$ for FIRE-2. Finally, the sampling of type II supernova mass return is also treated differently. In particular, in FIRE-2 each supernova returns the same ejecta mass, which is equal to the global IMF average of the mass return of all supernova events. In our model, each supernova has a different ejecta mass depending on the age of the stellar particle. Indeed, the ejecta mass is computed as an IMF average over the local star particle time step only. Basically, this translates to a progressively lower ejecta mass as a function of the star particle age;

• •

star formation criteria: while for our model and in FIRE-2 this is a combination of local gas density and virial parameter conditions (see Sec. 2.2 and equation 9), Agertz et al. (2013) use a density threshold criterion with metallicity dependence to discriminate between star-forming and non-starforming cells. This follows the form $n_{*}\simeq 25(Z/Z_{\odot})^{-1}\,{\rm cm}^{-3}$ and represents the transition between the neutral and molecular dominated phase of the interstellar gas according to Gnedin & Kravtsov (2011). Morevoer, only gas with $T<10^{4}\,{\rm K}$ is eligible for star formation in their model.

We hasten to add that other efforts have been made to attempt to produce a resolved ISM medium in Arepo (see Smith et al., 2018). However, there are several major differences with respect to our new model. Most importantly, radiative feedback and OB/AGB winds are not taken into account in the Smith et al. (2018) treatment. Other differences include the scheme adopted to distribute SN feedback energy and momentum. Smith et al. (2018) identify the gas cell containing the star particle (which they denote as host cell) and they assign $5$ per cent of the energy budget to it. The remaining part is then given to the nearest neighbours of the host cell by employing a modified version of the solid angle weighting scheme that we present here. Smith et al. (2018) also adopts a Kroupa (2001) IMF, in contrast with the Chabrier (2001) in SMUGGLE. Furthermore, each supernova event (sampled from a Poisson distribution) releases $10\,{\rm M}_{\odot}$ of which $2\,{\rm M}_{\odot}$ are metals or, in other words, a constant metallicity of $Z=0.02$ is assumed for the supernova ejecta. No other metal enrichment channel is considered. Lastly, the star formation criterion demands only a gas density threshold, whereas our model imposes an extra check on the virialization state of the gas, and additionally, an explicit dependence on the molecular fraction (see equation 12).

Although the impact of these specific differences in the outcome of each model is unclear, it is a useful reminder of the many unconstrained numerical choices that are necessarily made in each attempt to follow the physics of the ISM at sub-kpc resolution level. A detailed evaluation of our results and their dependence on variations of some of these choices will be presented in future work.

Finally, we note that important physical processes in the ISM are not taken into account in our model or are included with rather crude approximations. These processes include, for example, cosmic rays transport and physics (which are important as a source of pressure support and as an additional feedback channel to launch galactic-scale outflows), magnetic fields (which might provide additional pressure support to the gas and determine the transport properties of relativistic particles and thermal energy), dust production, destruction and evolution (important for both the molecular chemistry of the ISM and the reprocessing of radiation fields), thermal conduction (which operates at the transition layers between the supernova-heated hot gas and the cold phase). Also, while we include some form of feedback due to radiation, its transport and interaction with the gas are treated in a rather simplistic way. A more detailed implementation is available in Arepo for the majority of these physical processes (Pakmor & Springel, 2013; Kannan et al., 2016a; Pfrommer et al., 2017; Kannan et al., 2019; McKinnon et al., 2018); their inclusion in our model is a line of research that we intend to pursue. For example, efforts are currently under way to combine the SMUGGLE model with a realistic treatment of radiation transport and dust physics (Kannan et al., in prep.).

7 Summary and conclusions

Modelling the formation of galaxies and their evolution is a challenging task due to the large dynamic range and the variety of physical processes that concur in shaping galaxies. Numerical simulations are the ideal tool to tackle this non-linear and complex problem. Cosmological simulations have recently achieved important successes in modelling of Milky Way-like galaxies (Hopkins et al., 2014a; Grand et al., 2017; Wang et al., 2015) as well as in reproducing realistic populations of galaxies within the cosmological context (Vogelsberger et al., 2014; Schaye et al., 2015; Springel et al., 2018). However, these simulations have to rely on effective physical prescriptions, with different degrees of sophistication, to include important physical processes below their resolution scale – the so-called sub-grid physics. For the galaxy formation problem, sub-grid models are particularly important for the treatment of the interstellar medium (ISM) and stellar feedback processes. In this paper we haved presented the SMUGGLE model, a novel numerical implementation for the treatment of the ISM gas and stellar feedback in the moving mesh code Arepo.

This new model does not employ an effective equation of state for the ISM, but instead aims at resolving the different gas phases in the ISM and self-consistently implements a local stellar feedback scheme that generates gaseous outflows. The energy and momentum feedback inputs, and the associated generation of outflows, are tied to the location and characteristics of newly created stars rather than by imposing a prescribed outflow mass loading set based on the more global properties of the simulated galaxy. Specifically, our model includes the following physical processes: low-temperature molecular cooling; cosmic ray and photo-electric heating; star formation based on a high density threshold, virial parameter and molecular gas criteria; mechanical feedback from supernovae taking into account the momentum boost due to $P{\rm d}V$ work during the (unresolved) Sedov-Taylor expansion phase; radiative feedback from young, massive stars both in the form of photoionization and radiation pressure accounting for multiple infrared scattering; and the energy and momentum injection from OB and AGB stellar winds.

We have explored this novel implementation with simulations of an isolated non-cosmological Milky Way-like disc. Our model yields a star formation rate of $\sim 3\,{\rm M_{\odot}\,yr^{-1}}$, consistent with observational data and reproduces the observed Kennicutt-Schmidt relation. The regulation of star formation increases the gas depletion time-scale from a few hundred ${\rm Myr}$ in the no feedback case to a few ${\rm Gyr}$ once feedback processes are included, which is in line with observations. These results are largely independent of resolution. As a consequence of feedback, gaseous outflows across the star-forming disc are generated with mass loading factors of order unity. The gas outflow rates are quite bursty on a time-scale of $\simeq 100\,{\rm Myr}$, which is reflected in the behaviour of the mass loading factor. Contrary to the outflow rate, however, the mass loading factor remains approximately constant. The majority of the ejected gas cycles back onto the disc, as the similarity between gas outflow and inflow rates demonstrates.

Most importantly, the model is also able to produce a distinct multiphase ISM in the galactic disc with three gas phases coexisting. There is still a residual dependence in the amount of the hot gas present in the simulation and the numerical resolution. In particular, higher resolution corresponds to larger amounts of hot gas. This is likely due to the fact that the cooling of the gas is overestimated with coarser numerical resolution. The creation of the multiphase medium in the simulated galaxies is due to the interplay between gas cooling, star formation and the different feedback channels implemented in the model. Although all feedback processes considered here are important towards regulation of star formation, in the Milky Way-type galaxy simulated in this work, the dominant channel is represented by supernova feedback. Indeed, compared to radiative feedback and stellar winds, supernovae inject into the ISM about a factor of $10$ times more momentum. We note, however, that these ratios might change in more dense environments or at high redshifts.

Our new model is intended to pave the way towards the next generation of cosmological large volumes simulations together with associated zoom-ins of individual objects. While the future of galaxy formation simulations will require the direct treatment of additional physics such as radiative transfer, effects of cosmic rays, etc., the computational power needed to achieve such goals on large cosmological scales are beyond reach in the near future. Our model aims to bridge the gap between detailed physical modelling of the ISM requiring pc-resolution scale (e.g. Rosdahl et al., 2015; Simpson et al., 2016; Raskutti et al., 2016; Kannan et al., 2018; Krumholz, 2018; Emerick et al., 2018; Haid et al., 2018) and the current cosmological volume simulations. The application of our model to such volume simulations will increase the predictive power of these simulations, in particular for galactic structure, outflows and circum-galactic gas predictions. Such theoretical forecasts over a large number of galaxies will be directly contrasted with observations of galaxy population offering further constraints to galaxy formation models within the $\Lambda$CDM scenario.

Acknowledgments

The authors wish to thank Rahul Kannan, Eve Ostriker, Andrey Kravtsov, Lars Hernquist and Philip Hopkins for stimulating discussions about this work. FM acknowledges support through the Program “Rita Levi Montalcini” of the Italian MIUR. LVS is thankful for financial support from the Hellman Foundation as well as NSF and NASA grants, AST-1817233, HST-AR-14552 and HST-AR-14553. MV acknowledges support through an MIT RSC award, a Kavli Research Investment Fund, NASA ATP grant NNX17AG29G, and NSF grants AST-1814053 and AST-1814259. The simulations were performed on the joint Harvard-MIT computing cluster, supported by the Faculty of Arts and Sciences of Harvard University and the MIT Kavli Institute for Astrophysics and Space Research, and the HiPerGator computing cluster at the University of Florida. Preparatory work was made possible by the award of the XSEDE grants AST-140082, AST-140078, AST-150007 and AST-160006. All the figures in this work were created with the matplotlib graphics environment (Hunter, 2007).

Appendix A Discrete supernova rates

Supernovae are discrete and rare777In our stellar evolutionary model, each stellar particle produces $\simeq 3\times 10^{-4}\,{\rm SN\,M_{\odot}^{-1}\,Myr^{-1}}$ of type II supernovae for a main sequence lifetime of an $8\,{\rm M_{\odot}}$ star of $40\,{\rm Myr}$. Type Ia supernovae are even rarer with a maximum rate of $\simeq 7.8\times 10^{-8}\,{\rm SN\,M_{\odot}^{-1}\,Myr^{-1}}$ at a stellar age of $40\,{\rm Myr}$ (see equation 25). events, as such, they are sampled in our model as originating from a Poisson distribution

[TABLE]

where $p$ is the probability of having $n$ supernova events given an expected value $\lambda$ over the time step $\Delta t$.

Figure 15 presents the cumulative number of SNII (top) and SNIa (bottom) as a function of stellar age to illustrate the accuracy of our discrete sampling. Each plot shows a two-dimensional histogram (blue shades) counting the number of stellar particles of a given age yielding a cumulative number of SN per unit stellar mass formed (see figure caption for details on binning). The red solid line is the average number of cumulative SN events per age bin with black dashed lines showing the standard deviation around this mean value. The grey dashed lines are the number of SN events per unit mass expected from the stellar evolution model adopted in this work, with the associated standard deviation represented by the grey dotted lines.

It is immediately obvious that the two supernova types have two different distributions as a function of the stellar age. In particular, type II supernovae display a very steep rise in the first $40\,{\rm Myr}$ after which they plateau. This is a direct consequence of the fact that we set the minimum mass for a star to explode as a type II supernova to $8\,{\rm M_{\odot}}$, that in turn implies a stellar lifetime of $\sim 40\,{\rm Myr}$ (with a very slight dependence on metallicity). In other words after about that age, all supernova II events have exploded and no further events can occur. Type Ia supernova instead have a different trend: events start only after $40\,{\rm Myr}$ and then follow the trend described in equation (24). Please note that our sampling method is able to recover both the average number of supernova events as a function of age and their scatter around the mean value. Given the lower number of type Ia events, the signatures of discrete supernova explosions are clearly visible in the bottom panel as the horizontal stripes in the histogram.

It is worth noting that the level of scatter around the mean value, i.e. its amplitude, is expected to increase with resolution. This is because progressively less massive stellar particles (with increasing resolution) tend, on average, to produce less supernova events per particle, thus increasing the scatter. This can be easily shown as follows. Let us define

[TABLE]

as the expected number of supernova events at $1\sigma$ above/below the average – the last equality holds because for a Poisson distribution with expected value $\lambda$, $\sigma^{2}=\lambda$. The amplitude of the scatter will be

[TABLE]

or per unit mass

[TABLE]

in which $n=N/m_{\star}$ and $m_{\star}$ is the average mass of a stellar particle. The number of supernova events (of both types) per unit mass $n$ is independent of resolution and determined only by the stellar evolutionary model adopted. However, for increasing resolution $m_{\star}$ decreases and the amplitude of the scatter around the expected number of supernova events per unit mass increases accordingly.