Snap, Crackle, Pop: sub-grid supernova feedback in AMR simulations of disc galaxies
We compare 5 sub-grid models for supernova (SN) feedback in adaptive mesh refinement (AMR) simulations of isolated dwarf and L-star disc galaxies with pc resolution. The models are thermal dump, stochastic thermal, “mechanical” (injecting energy or momentum depending on the resolution), kinetic, and delayed cooling feedback. We focus on the ability of each model to suppress star formation and generate outflows. Our highest-resolution runs marginally resolve the adiabatic phase of the feedback events, which correspond to 40 SN explosions, and the first three models yield nearly identical results, possibly indicating that kinetic and delayed cooling feedback converge to wrong results. At lower resolution all models differ, with thermal dump feedback becoming inefficient. Thermal dump, stochastic, and mechanical feedback generate multiphase outflows with mass loading factors , which is much lower than observed. For the case of stochastic feedback we compare to published SPH simulations, and find much lower outflow rates. Kinetic feedback yields fast, hot outflows with , but only if the wind is in effect hydrodynamically decoupled from the disc by using a large bubble radius. Delayed cooling generates cold, dense and slow winds with , but large amounts of gas occupy regions of temperature-density space with short cooling times. We conclude that either our resolution is too low to warrant physically motivated models for SN feedback, that feedback mechanisms other than SNe are important, or that other aspects of galaxy evolution, such as star formation, require better treatment.
keywords:galaxies: formation, galaxies: evolution, methods: numerical
In our CDM Universe, most of the mass is made up of dark matter. On large scales baryons trace the dark matter and its gravitational potential. Baryonic gas falls into galaxies at the centres of dark matter haloes, where it cools radiatively and collapses to form stars. By naive gravitational arguments, star formation should be a fast affair, consuming the gas over local free-fall times. However, from observations we know that it is a slow and inefficient process, taking free-fall times, depending on the scale under consideration (e.g. Zuckerman & Evans, 1974; Krumholz & Tan, 2007; Evans et al., 2009). Also, while observers have a notoriously hard time confirming the existence of gas flowing into galaxies (Crighton et al., 2013), they instead routinely detect oppositely directed outflows at velocities of hundreds of km/s (see review by Veilleux et al., 2005).
To understand the non-linear problem of galaxy formation and evolution, theorists use cosmological simulations of dark matter, describing the flow and collapse of baryonic star-forming gas either with directly coupled hydrodynamics or semi-analytic models. Strong feedback in galaxies is a vital ingredient in any model of galaxy evolution, be it hydrodynamical or semi-analytic, that comes even close to reproducing basic observables, such as the star formation history of the Universe, the stellar mass function of galaxies, the Kennicutt-Schmidt relation, rotational velocities, and outflows (e.g. Vogelsberger et al., 2013; Dubois et al., 2014; Hopkins et al., 2014; Schaye et al., 2015; Wang et al., 2015; Somerville & Davé, 2015).
In order to capture the inefficient formation of stars, the first generation of galaxy evolution models included core collapse (or type II) supernova (SN) feedback, where massive stars () end their short lives with explosions which inject mass, metals, and energy into the inter-stellar medium (ISM). In early hydrodynamical simulations, the time-integrated type II SN energy of a stellar population, ergs per SN event, was dumped thermally into the gas neighbouring the stellar population (Katz, 1992). However, such thermal dump feedback had little impact on star formation, resulting in an over-abundance of massive and compact galaxies. This so-called over-cooling problem is partly numerical in nature, and a result of low resolution both in time and space. As discussed by Dalla Vecchia & Schaye (2012), the energy is injected into too large a gas mass, typically resulting in much lower temperatures than those at work in sub-pc scale SN remnants. The relatively high cooling rates at the typical initial temperatures attained in the remnant, of K, allow a large fraction of the injected energy to be radiated away before the gas reacts hydrodynamically, resulting in suppressed SN blasts and hence weak feedback. Gas cooling is, however, also a real and physical phenomenon, and while it is over-estimated in under-resolved simulations, a large fraction of the energy in SN remnants may in fact be radiated away instead of being converted into large-scale bulk motion (Thornton et al., 1998).
A number of sub-resolution SN feedback models have been developed over the last two decades for cosmological simulations, with the primary motivation of reproducing large-scale observables, such as the galaxy mass function, by means of efficient feedback. The four main classes of these empirically motivated SN feedback models are i) kinetic feedback (Navarro & White, 1993), where a fraction of the SN energy is injected directly as momentum, often in combination with temporarily disabling hydrodynamical forces (Springel & Hernquist, 2003), ii) delayed cooling (e.g. Gerritsen, 1997; Stinson et al., 2006), where radiative cooling is turned off for some time in the SN remnant, iii) stochastic feedback (Dalla Vecchia & Schaye, 2012), where the SN energy is re-distributed in time and space into fewer but more energetic explosions, and iv) multiphase resolution elements that side-step unnatural ‘average’ gas states at the resolution limit (Springel & Hernquist, 2003; Keller et al., 2014).
In principle, a physically oriented approach to implementing SN feedback with sub-grid models is desirable. The goal is then to inject the SN blast as it would emerge on the smallest resolved scale, by making use of analytic models and/or high-resolution simulations that capture the adiabatic phase, radiative cooling, the momentum driven phase, and the interactions between different SN remnants. However, these base descriptions usually include simplified assumptions about the medium surrounding the SN remnant, and fail to capture the complex inhomogeneities that exist on unresolved scales and can have a large impact on cooling rates. In addition, even if the SN energy is injected more or less correctly at resolved scales, it will generally fail to evolve realistically thereafter because the multi-phase ISM of simulated galaxies is still at best marginally resolved. Hence there remains a large uncertainty in how efficiently the SN blast couples to the ISM. This translates into considerable freedom, which requires SN feedback models to be calibrated to reproduce a set of observations (see discussion in Schaye et al., 2015).
The most recent generation of cosmological simulations has been relatively successful in reproducing a variety of observations, in large part thanks to the development of subgrid models for efficient feedback and the ability to calibrate their parameters, as well as the inclusion of efficient active-galactic nucleus (AGN) feedback in high-mass galaxies. However, higher-resolution simulation works (e.g. Hopkins et al., 2012; Agertz et al., 2013) suggest that SNe alone may not provide the strong feedback needed to produce the inefficient star formation we observe in the Universe.
Attention has thus been turning towards complementary forms of stellar feedback, which provide additional support to the action of SNe. Possible additional feedback mechanisms include stellar winds (e.g. Dwarkadas, 2007; Rogers & Pittard, 2013; Fierlinger et al., 2016), radiation pressure (e.g. Haehnelt 1995; Thompson et al. 2005; Murray et al. 2010, but see Rosdahl et al. 2015), and cosmic rays (e.g. Booth et al., 2013; Hanasz et al., 2013; Salem & Bryan, 2014; Girichidis et al., 2016).
None the less, SN explosions remain a powerful source of energy and momentum in the ISM and a vital ingredient in galaxy evolution. For the foreseeable future a sub-resolution description of them will remain necessary in cosmological simulations and even in most feasible studies of isolated galaxies. The true efficiency of SN feedback is still not well known, and hence we do not know to what degree we need to improve our SN feedback sub-resolution models versus appealing to the aforementioned complementary physics.
Rather than introducing a new or improved sub-resolution SN feedback model, the goal of this paper is to study existing models, using controlled and relatively inexpensive numerical experiments of isolated galaxy discs modelled with gravity and hydrodynamics in the Eulerian (i.e. grid-based) code Ramses (Teyssier, 2002). We use those simulations to assess each model’s effectiveness in suppressing star formation and generating galactic winds, the main observational constraints we have on feedback in galaxies.
We study five subgrid prescriptions for core-collapse SN feedback in isolated galaxy discs. We explore the ‘maximum’ and ‘minimum’ effects we can get from SN feedback using these models, and consider how they vary with galaxy mass, resolution, and feedback parameters where applicable. The simplest of those models is the ‘classic’ thermal dump, where the SN energy is simply injected into the local volume containing the stellar population. Three additional models we consider have been implemented and used previously in Ramses. These are, in chronological order, kinetic feedback, described in Dubois & Teyssier (2008) and used in the Horizon-AGN cosmological simulations (Dubois et al., 2014), delayed cooling, described in Teyssier et al. (2013), and mechanical feedback, described in Kimm & Cen (2014) and Kimm et al. (2015). In addition, for this work we have implemented stochastic feedback in Ramses, adapted from a previous implementation in the smoothed particle hydrodynamics (SPH) code Gadget, described in Dalla Vecchia & Schaye (2012, henceforth DS12).
The layout of this paper is as follows. First, we describe the setup of our isolated galaxy disc simulations in §2. We then describe the SN feedback models in §3. In §4 we compare results for each of these models using their fiducial parameters in galaxy discs of two different masses, focusing on the suppression of star formation and the generation of outflows. In §5 we compare how these results converge with numerical resolution, both in terms of physical scale, i.e. minimum gas cell size, and also in terms of stellar particle mass. In Sections 6 - 8 we take a closer look at the stochastic, delayed cooling, and kinetic feedback models respectively, and study how varying the free parameters in each model affects star formation, outflows and gas morphology. The reader can skip those sections or pick out those of interest, without straying from the thread of the paper. We discuss our results and implications in §9, and, finally, we conclude in §10.
Before we introduce the SN feedback models compared in this paper, we begin by describing the default setup of the simulations common to all runs.
We run controlled experiments of two rotating isolated disc galaxies, consisting of gas and stars, embedded in dark matter (DM) haloes. The main difference between the two galaxies is an order of magnitude difference in mass, both baryonic and DM. We use the AMR code Ramses (Teyssier, 2002), which simulates the interaction of dark matter, stellar populations and baryonic gas, via gravity, hydrodynamics and radiative cooling. The equations of hydrodynamics are computed using the HLLC Riemann solver (Toro et al., 1994) and the MinMod slope limiter to construct variables at cell interfaces from the cell-centred values. We assume an adiabatic index of to relate the pressure and internal energy, appropriate for an ideal monatomic gas. The trajectories of collisionless DM and stellar particles are computed using a particle-mesh solver with cloud-in-cell interpolation (Guillet & Teyssier 2011; the resolution of the gravitational force is the same as that of the hydrodynamical solver).
2.1 Initial conditions
The main parameters for the simulated galaxies and their host DM haloes are presented in Table 1. We focus most of our analysis on the lower-mass galaxy, which we name g9. It has a baryonic mass of , with an initial gas fraction of , and it is hosted by a DM halo of mass . We also compare a less detailed set of results for feedback models in a more massive galaxy, g10, similar (though somewhat lower) in mass to our Milky-Way (MW), with , , and . Each simulation is run for Myr, which is orbital times (at the scale radii) for both galaxy masses, and enough for star formation and outflows to settle to quasi-static states.
The initial conditions are generated with the MakeDisk code by Volker Springel (see Springel et al., 2005; Kim et al., 2014), which has been adapted to generate Ramses-readable format by Romain Teyssier and Damien Chapon. The DM halo follows an NFW density profile (Navarro et al., 1997) with concentration parameter and spin parameter (Macciò et al., 2008). The dark matter in each halo is modelled by one million collisionless particles, hence the g9 and g10 galaxies have DM mass resolution of and , respectively. The initial disc consists of stars and gas, both set up with density profiles which are exponential in radius and Gaussian in height from the mid-plane (scale radii of kpc for g9 and kpc for g10, and scale heights one tenth of the scale radius in both cases). The galaxies contain stellar bulges with masses and scale radii both one tenth that of the disc. The initial stellar particle number is , a million of which are in the disc and the remainder in the bulge. The mass of the initial stellar particles is and for the g9 and g10 galaxies, respectively, close to the masses of stellar particles formed during the simulation runs, which are shown in Table 1. While contributing to the dynamical evolution and gravitational potential of the rotating galaxy disc, the initial stellar particles do not explode as SNe. This initial lack of feedback results in over-efficient early star formation and a subsequent strong feedback episode which typically then suppresses the star formation to a semi-equilibrium state within a few tens of Myr (see e.g. star formation rate plots in §4.2). To overcome this shortcoming, future improvements should include sensible age assignments to the initial stellar particles, which could be used to perform SN feedback right from the start of the simulations.
The temperature of the gas discs is initialised to a uniform K and the ISM metallicity is set to and for the g9 and g10 galaxies, respectively. The circum-galactic medium (CGM) initially consists of a homogeneous hot and diffuse gas, with , K, and zero metallicity. The cutoffs for the gas discs are chosen to minimize the density contrast between the disc edges and the CGM. The square box widths for the g9 and g10 galaxies are and kpc, respectively, and we use outflow (i.e. zero gradient) boundary conditions on all sides.
The same initial conditions and similar simulation settings were used in Rosdahl et al. (2015), where we studied stellar radiation feedback in combination with thermal dump SNe. The main differences from the setup of the previous work, apart from not including stellar radiation, is that here we include a homogeneous UV background, we form stellar particles that are about a factor three more massive, and the previous work included a bug, now fixed111Thanks to Sylvia Ploeckinger for finding and fixing the issue., in metal cooling, where the contribution of hydrogen and helium was double-counted at Solar metallicity222We have checked and verified that the metal cooling bug has a negligible effect on the results of both this and our previous work.. The most significant of these changes is the larger stellar particle mass, which boosts the efficiency of thermal dump SN feedback in suppressing star formation and, to a lesser extent, in generating outflows.
2.2 Adaptive refinement
Each refinement level uses half the cell width of the next coarser level, starting at the box width at the first level. Our simulations start at level , corresponding to a coarse resolution of cells, and adaptively refine up to a maximum level 14, corresponding to an effective resolution of cells. This corresponds to an optimal physical resolution of pc and pc in the less and more massive galaxies, respectively. Refinement is done on mass: a cell is refined if it is below the maximum refinement level, if its total mass (DM+stars+gas) exceeds (see mass values in Table 1), or if its width exceeds a quarter of the local Jeans length
2.3 Gas thermochemistry
The gas temperature and the non-equilibrium ionization states of hydrogen and helium are evolved with the method presented in Rosdahl et al. (2013), which includes collisional ionization/excitation, recombination, bremsstrahlung, di-electronic recombination, and Compton electron scattering off cosmic microwave background photons. We include hydrogen and helium photo-ionization and heating of diffuse gas from a redshift zero Faucher-Giguère et al. (2009) UV background, and enforce an exponential damping of the UV radiation above the self-shielding density of .
Above K, the contribution to cooling from metals is added using Cloudy (Ferland et al., 1998, version 6.02) generated tables, assuming photoionization equilibrium with a redshift zero Haardt & Madau (1996) UV background. Below K, we use fine structure cooling rates from Rosen & Bregman (1995), allowing the gas to cool radiatively to K.
2.4 Star formation
We use a standard star formation (SF) model which follows a Schmidt law. In each cell where the hydrogen number density is above the star formation threshold, , gas is converted into stars at a rate , where is the gas (mass) density and is the star formation efficiency per free-fall time, , where is the gravitational constant. Stellar populations are represented by collisionless stellar particles that are created stochastically using a Poissonian distribution (for details see Rasera & Teyssier, 2006), which returns the stellar particle mass as an integer multiple of (see Table 1). We use in this work (e.g. Krumholz & Tan, 2007). In future work we will consider how varying the details of star formation affects the efficiency of SN feedback, but that is beyond the scope of the present paper. The stellar particle masses are given in Table 1, and are equal to the SF density threshold, , times the volume of a maximally refined gas cell333We do not allow more than of the cell gas to be removed when forming stars. Thus, stellar particles actually do not form below a density of ..
2.5 Artificial Jeans pressure
To prevent numerical fragmentation of gas below the Jeans scale (Truelove et al., 1997), an artificial ‘Jeans pressure’ is maintained in each gas cell in addition to the thermal pressure. In terms of an effective temperature, the floor can be written as , where we have set K (and is the aforementioned star formation threshold), to ensure that the Jeans length is resolved by a constant minimum number of cell widths at any density – 7 and 3.5 cell widths in the smaller and larger galaxy simulations, respectively (see Eq. 3 in Rosdahl et al., 2015). The pressure floor is non-thermal, in the sense that the gas temperature which is evolved in the thermochemistry is the difference between the total temperature and the floor – therefore we can have .
3 SN feedback
Supernova feedback is performed with single and instantaneous injections of the cumulative SN energy per stellar population particle. Each stellar particle has an energy and mass injection budget of
respectively, where is the fraction of stellar mass that is recycled into SN ejecta444Note that we will neglect the mass that ends up in stellar remnants of SNe., is the average stellar mass of a type II SN progenitor, and, as a reminder, is the mass of the stellar particle. We assume a Chabrier (2003) initial mass function (IMF) and set and , giving at least ergs per particle in the g9 galaxy and ergs in the g10 galaxy. We neglect the metal yield associated with stellar populations, i.e. the stellar particles inject no metals into the gas, and the metallicity of the gas disc stays at roughly the initial value of solar (which is negligibly diluted due to mixing with the pristine CGM). The time delay for the SN event is 5 Myr after the birth of the stellar particle.
The model for SN energy and mass injection, and how it affects the galaxy properties and its environment, is the topic of this paper. We explore five different SN models, which we now describe.
3.1 Thermal dump feedback
This is the most simple feedback model, and one which is well known to suffer from catastrophic radiative losses at low resolution (e.g. Katz, 1992). The (thermal) energy and mass of the exploding stellar particle are dumped into the cell hosting it, and the corresponding mass is removed from the particle. Unless the Sedov-Taylor phase is well resolved in both space and time, the thermal energy radiates away before it can adiabatically develop into a shock wave. The primary aim of each of the SN models that follow is to overcome this ‘overcooling’ problem.
Note that in SPH simulations, the energy in thermal dump feedback is typically distributed over neighbouring gas particles, whereas in our implementation all the energy is injected into a single cell. Consequently, in SPH simulations with similar resolution, the amount of gas that is heated is typically larger. This can lead to lower temperatures and larger radiative losses in SPH, but in the case of strong density gradients around the feedback event, it can also enhance feedback efficiency if SPH particles with a low density receive part of the SN energy.
3.2 Stochastic thermal feedback
While the other SN models described in this paper existed previously in Ramses and have been described and studied individually in previous publications, we have for this work adapted to Ramses the stochastic SN feedback model presented in Dalla Vecchia & Schaye (2012, a.k.a. DS12), which has so far only been used in SPH. The idea is to heat the gas in the cell hosting the stellar particle to a temperature high enough that the cooling time is long compared with the sound crossing time across the cell. The energy then has the chance to do significant work on the gas before being radiated away, and overcooling is reduced.
As argued in DS12, a single SN energy injection should heat the gas enough that the ratio between the cooling time and sound crossing time across a resolution element is . Given a local gas density , a physical resolution , and assuming cooling is dominated by Bremsstrahlung (true for K), an expression from DS12 (their eq. 15) can be used to derive an approximate required temperature increase, , to enforce this minimum ratio and thus avoid catastrophic cooling, resulting in the condition that
Some specified time delay after the birth of the stellar population particle ( Myr in this work), it injects its total available energy, , into the hosting gas cell. Since may be smaller than what is needed for the required temperature increase , the feedback event is done stochastically, with a probability
where is the gas mass of the host cell (including the SN ejecta) and
is the required specific energy, with the Boltzmann constant, the proton mass, and the mean particle mass in units of 555We use , assuming the gas to be ionized.. When a stellar particle is due to inject SN energy, is calculated via Eq. (4). If , the available energy is sufficient to meet the cooling time constraint and it is simply injected into the host cell. On the other hand, if , a random number between 0 and 1 is drawn: only if is the energy injected into the cell, otherwise no energy injection takes place. With this approach, the energy averaged over the whole simulation box and sufficiently long time-scales is close to the available SN energy budget, as we have confirmed in our runs – it is just distributed unevenly in space and time in order to overcome the cooling catastrophe. Note that the mass (and metal) yield from the stellar particles is not subject to the stochastic process, but is always injected into the host cell, regardless of whether or not the energy injection takes place666This results in slight cooling of gas in those cells where stellar particles inject mass but not energy..
We note that our stochastic feedback implementation in AMR differs significantly from the original SPH implementation from DS12 in two regards. First, whereas the probability for a feedback event varies with the local density in AMR, the event probability is a constant over a simulation run in SPH, since each candidate for energy injection is a gas particle and all gas particles typically have identical mass. Second, thermal dump explosions in SPH are normally injected into a number of (typically ) neighbouring gas particles, and the objective of the stochastic feedback model is then to reduce the number of neighbours receiving the feedback energy in each event (in effect making it more similar to our implementation of thermal dump feedback). However, in AMR, the energy is only released into the gas cell hosting the exploding stellar particle, so our stochastic feedback model redistributes feedback events such that some stellar particles explode with boosted energy, and some not at all.
Naively, our stochastic feedback implementation may be presumed to re-distribute SN energy towards lower gas densities, as the probability for each SN event scales inversely with the density via . This is not the case: there is no (average) re-distribution over densities, since the injected energy scales inversely with the probability, and hence the average energy injected per ‘candidate’ feedback event, at any density, is
We study the effects of varying the parameter in Sec. §6. For the comparison of feedback models, we use the fiducial value of K, because: i) it is roughly the minimum value given by Eq. (3) using our star formation density threshold, ii) in our comparison of values in Sec. §6 we find this gives a similar star formation and Kennicutt-Schmidt relation as higher values (though note that higher values of do produce stronger outflows), and iii) it is the fiducial value used in DS12 and in the EAGLE simulations (Schaye et al., 2015), allowing us to qualitatively compare the efficiency of our AMR version of the stochastic feedback model to that of previous SPH works. Assuming a resolution of pc and that stellar particles typically explode close to the star-formation density threshold of , Eq. (3) gives K, so our chosen fiducial value is well above the estimated requirement from DS12777Indeed, with , , and pc, the sound crossing time is yr (eq. 9 in DS12) and the cooling time is yr (eq. 13 in DS12).. Using the same values, Eq. (4) shows that the probability for feedback events is below unity at densities , i.e. slightly above our adopted star formation density threshold, for both the g9 and g10 galaxies we simulate (the lower resolution in g10 is exactly counter-balanced by the larger stellar particle mass).
3.3 Delayed cooling thermal feedback
Another widely used method for overcoming the numerical overcooling problem in galaxy formation simulations is to turn off radiative cooling in SN heated gas for a certain amount of time. This method, usually referred to as delayed cooling, has been used in SPH simulations (e.g. Gerritsen & de Blok, 1999; Stinson et al., 2006; Governato et al., 2010), giving both a strong suppression in star formation and enhancement in outflows.
Teyssier et al. (2013) introduced an AMR version of this method, which we use in this paper. Here, a specific energy tracer, , is stored on the grid in the form of a passive scalar, and typically associated with an unresolved turbulent energy. It is advected with the gas and decays on a time-scale as
For every feedback event, the SN energy of the stellar particle, , is injected as thermal energy into the host cell, as in thermal dump feedback, but at the same time it is added to the non-thermal energy density in the same cell. As long as the local “turbulent velocity” is above a certain limit in a given cell, , radiative cooling is disabled in that location, mimicking the non-thermal nature of turbulent energy. When the local turbulent velocity has fallen below , via decay, diffusion, and mixing, radiative cooling is enabled again.
The main free parameter in the model is , which determines how quickly the turbulent energy disappears. is also an adjustable parameter, but it has more or less the same effect as , so we keep it fixed at (corresponding to about of the injected specific energy of a SN, or about th of the velocity in its unloaded remnant). The value of can be motivated by an underlying physical mechanism, e.g. the crossing time over a few cell widths, after which the resolved hydrodynamics should take over the unresolved advection of energy. The appendix of Dubois et al. (2015) derives an expression for the choice of an appropriate , given the local SN rate, density, and resolution (their eq. A8), for which our g9 simulation settings (, , , pc) give Myr. However, in this paper we follow the literature (Teyssier et al., 2013; Roškar et al., 2014; Mollitor et al., 2015; Rieder & Teyssier, 2016), and use a much larger fiducial value of Myr for the g9 galaxy, and Myr in low-resolution versions of g9 and in the g10 galaxy. Assuming decay dominates over diffusion and mixing, and assuming the SN remnants travel at () km/s, our fiducial Myr corresponds to a delay length scale of () kpc. We explore variations of in §7 (including values close to that derived by Dubois et al. 2015).
The disadvantage of delayed cooling is that while overcooling is in part a numerical problem, radiative cooling is a real and physical process, without which stars would not form at all. By neglecting radiative cooling altogether, even if for a relatively short time, delayed cooling is likely to result in over-efficient type II SN feedback (but, at the same time, it perhaps compensates for the neglect of other feedback processes which may be important in galaxy evolution). In addition, delayed cooling can result in the gas populating parts of the temperature-density diagram where the cooling time is short, which may yield unrealistic predictions for absorption and emission diagnostics.
3.4 Kinetic feedback
We use the kinetic feedback model presented in Dubois & Teyssier (2008). Here, the trick to overcoming numerical overcooling is to skip the unresolved Sedov-Taylor phase and directly inject the expected collective result of that phase for a stellar population, which is an expanding momentum-conserving shock wave (or snowplow). Note, however, that the injected kinetic energy may subsequently be converted into thermal energy if shocks develop.
SN mass and momentum is injected into gas within a bubble radius of the exploding stellar particle. The free parameters for the method are , the fraction of which is released in kinetic form, , the radius of the bubble, and , the sub-resolution mass loading factor of the Sedov-Taylor phase, describing how much mass, relative to the stellar mass, is redistributed from the cell at the bubble centre to the bubble outskirts.
The redistributed mass consists of two components: one is the SN ejecta, , removed from the stellar particle, the other is the swept up mass, , removed from the central host cell (no more than of the central cell mass is removed, hence for individual feedback events at relatively low densities it may happen that is smaller than ). The total wind mass is thus , which is redistributed uniformly (i.e. uniform density) to all cells inside the bubble.
The kinetic energy, , is likewise distributed to the bubble cells, but with an injected velocity (directed radially away from the stellar particle) that increases linearly with distance from the centre, such as to approximate the ideal Sedov-Taylor solution:
where is the mass added to the cell, is the position of the centre of the cell relative to the stellar particle, is a bubble normalisation constant888The normalisation constant is the volume-weighted average distance from the centre, for each volume element in the bubble. In the ideal case of infinitely small cells, the factor is . required to ensure that the total redistributed energy is equal to , and
is the unnormalised wind velocity, where we used Eq. (1) for the latter equality. Note that this is the velocity of the added mass, i.e. each cell gains momentum
so if the mass already in the cell is substantial compared to the added mass, the resulting velocity change can be small. The injection is performed in the mass-weighted frame of the SN particle (with ) and host cell (with ). The remaining thermal energy, , is then distributed uniformly between the bubble cells.
In this work, we use fiducial parameters , , and pc, a size comparable to galactic super-bubbles (note that it is also comparable to the initial scale height of the stellar and gas disc in our simulations, which is pc and pc for the g9 and g10 galaxies, respectively). These values give a velocity for the gas ejected from the central cell (from Eq. 9) of km/s. Our choice of implies that there have been neither radiative losses nor momentum cancellation from the set of unresolved individual SNe inside the bubble. We explore the effects of a smaller bubble and higher mass loading in §8.
3.5 Mechanical feedback
This model was introduced to the Ramses code by Kimm & Cen (2014, see also ), and an analogue SPH scheme was earlier described independently in Hopkins et al. (2014). Here, momentum is deposited into the neighbour cells of a SN hosting cell, with the magnitude adaptively depending on whether the adiabatic phase of the SN remnant is captured by this small bubble of cells and the mass within it, or whether the momentum-conserving (snowplow) phase is expected to have already developed on this scale. In the first case, the momentum is given by energy conservation, while in the latter case, the final momentum, which depends via the cooling efficiency on the density and metallicity, is given by theoretical works (Blondin et al., 1998; Thornton et al., 1998).
In a single SN injection event, SN momentum (and any excess energy) is distributed over the nearest neighbours (sharing at least two vertices) of the SN host cell. The number of such cells can vary, depending on the cell refinement, but given the extreme limit where all the neighbours are at a finer level (i.e. half the cell width) of the SN host cell, the maximum number of neighbours is . When a neighbouring cell is at the same level as the host, or one level coarser, it is given an integer weight corresponding to how many of the finer level cell units it contains ( if sharing a plane with the host, if sharing a line). The SN host cell has a weight of , so the total number of cell units receiving direct SN energy injection is .
The goal is to inject into each neighbour cell a momentum , corresponding to that generated during the energy conserving phase if that is resolved, but to let converge towards that of the momentum-conserving snowplow phase in the limit that the energy conserving phase is unresolved. In each SN neighbour cell, this limit (energy vs momentum conserving) depends on the local mass-loading, i.e. the ratio of the local wind mass, to the SN ejecta given to that cell,
with , , the mass initially contained in the SN host cell, and the initial neighbouring gas mass, with and the gas neighbour cell gas density and host cell width, respectively.
The momentum injected into each neighbour cell, radially from the source, is
Here, the upper expression represents the resolved energy conserving phase, and comes from assuming the (final) cell gas mass of receives a kinetic energy (we ignore for the moment). The lower expression represents the asymptotic momentum reached in the snowplow phase, with is the total SN energy (i.e. ) in units of erg, the local hydrogen number density in units of , and . The Solar metallicity form of the expression was derived from analytic arguments, and confirmed with numerical experiments, in Blondin et al. (1998), and the dependency was added in the numerical work of Thornton et al. (1998).
The phase transition ratio, , is found by equating the snowplow expression in Eq. (12) with , where is the fraction of the SN energy assumed to be in kinetic form at the transition. This gives
where we used Eq. (1), , and normalised to typical values for the g9 galaxy in the latter equality. The function
ensures a smooth transition between the two expressions in Eq. (12).
If the momentum injection results in removal of total energy in a cell, due to cancellation of velocities, the surplus energy (initial minus final) is added to the cell in thermal form. As it has no preferred direction, the SN host cell receives only thermal energy. In Kimm & Cen (2014) and Kimm et al. (2015), due to wrong book-keeping of the surplus energy, the thermal energy injection during the adiabatic phase was overestimated (by a factor ) in regions where the swept up mass is large compared to the SN ejecta (by a similar factor), but the correct momentum and energy was used during the snowplow phase and the adiabatic phase with little mass loading (). This bug has since been corrected.
If we assume that all the cells receiving the SN energy have the same refinement level (i.e. the same cell width) and density and that the density is at least as high as the threshold for star formation, then the initial mass of the neighbour dominates and we can get a rough estimate for the local mass loading,
Here, the last equality, which comes from comparing with Eq. (13) and normalising to the g9 simulation parameters, shows that mechanical feedback events are marginally resolved, with the momentum injection being done using the upper expression in Eq. (12) for , but switching to the final snowplow momentum, i.e. the lower expression in Eq. (12), for higher gas densities. For the g10 galaxy, where the resolution is lower ( pc), the stellar mass higher (), and the metallicity higher (), the SN blasts are slightly worse resolved, with (at ) for the same assumptions. Here the effects of lower resolution and higher metallicity, towards worse-resolved SN blasts, are counter-weighted by the higher stellar particle mass.
4 SN feedback model comparison
We begin by comparing all SN feedback models using the fiducial settings. Later in this paper we will study each feedback model in more detail and show how the results vary with the values of the free parameters. We focus on star formation, outflows, and galaxy morphologies. Unless stated otherwise, our analysis will be restricted to the lower-mass g9 galaxy.
4.1 Galaxy morphologies
In Fig. 1, we show the total hydrogen column density face-on and edge-on at the end of the Myr run for each feedback model. The maps illustrate how the gas morphology is affected by the SN feedback models. Without feedback (top left panel), the galaxy becomes clumpy, containing dense star-forming clouds which accrete gas, thus creating large ‘holes’. The gas outside the thin edge-on disc is diffuse and featureless.
Compared to the no feedback case, thermal dump feedback (top right panel) significantly changes the gas morphology, reducing the gas clumpiness and thickening the disc. In fact, comparing to other panels in Fig. 1, it has here a very similar morphological effect as the stochastic and mechanical feedback models (middle left and bottom right panels, respectively). We will come back to this similarity in later sections.
Delayed cooling (middle right panel) and kinetic feedback (bottom left), on the other hand, produce quite different morphologies from other models in Fig. 1. Delayed cooling diffuses the gas more, with less obvious spiral structure, and the disc becomes thicker, indicating increasing feedback efficiency. In stark contrast, kinetic feedback results in a very thin disc plane, and thin, well defined spiral filaments. Judging qualitatively from these images of column density, delayed cooling appears most efficient in terms of smoothing out the gas, thickening the disc, and creating outflows, while kinetic feedback visually appears weakest, with relatively dense and filamentary gas in the disc and low column densities out of the disc plane. However, as we will see in what follows, kinetic feedback actually has the strongest and fastest (but relatively diffuse) outflows. We note that these distinct features of kinetic feedback are sensitive to the radius of momentum and mass injection, i.e. the parameter. As we will argue in §4.5, with our fiducial bubble size of pc, the momentum injection is essentially hydrodynamically decoupled from the galactic disc, and as we show in §8, a considerably smaller bubble leads to kinetic feedback behaving similarly to thermal dump, stochastic, and mechanical feedback.
4.2 Star formation
The feedback efficiencies can be quantified and compared via the star formation rates (SFRs), which we show for the g9 galaxy in Fig. 2. The SFRs are calculated by binning the stellar mass formed over time intervals of Myr. They vary by almost two orders of magnitude, depending on the feedback model utilised, and one order of magnitude at the end of the simulation runtime, by which time the rate of evolution has settled down after the initial collapse of the disc (due to radiative cooling and lack of initial feedback) and burst of star formation around the Myr mark.
Focusing on the star formation around Myr, we find that the feedback models separate roughly into the same three groups as in our assessment of the morphologies. Thermal dump, stochastic, and mechanical feedback all perform almost identically in terms of star formation, indicating that thermal dump is not strongly affected by overcooling (see §4.6). The SFR is suppressed by about a factor of compared to the no feedback case (labelled NoFB in the plot). This may seem an inefficient suppression, compared to the inferred average efficiency of star formation observed in the Universe, but it should be kept in mind that the star formation model already has a built in sub-resolution efficiency of only (see §2.4). We comment further on the choice and effect of in the discussion (§9.1).
At Myr, kinetic feedback has a SFR fairly close to those three aforementioned models. The difference is that the star formation rate has not stabilised, but is declining steadily. As we will see, this is due to the strong outflow removing gas from the star-forming ISM.
Delayed cooling is by far the most effective at suppressing star formation. The feedback from the initial peak in the SFR almost blows apart the gas disc, but once it has settled again the SFR stabilises around , though it remains somewhat bursty. The final SFR at Myr is almost an order of magnitude lower than for the other feedback models.
4.2.1 Star formation in the more massive galaxy
In Fig. 3 we show the SFRs in the ten times more massive (and lower resolution) g10 galaxy simulations.
Due to the combination of the deeper gravitational potential, stronger (metal) cooling, lower resolution, and the SN events happening at higher gas densities (typically by dex) we find more differences between the feedback models in their ability to suppress star formation than for the g9 galaxy. Thermal dump feedback is weak, with the star formation stabilising at the same rate as for the case of no feedback. With stochastic and mechanical feedback the star formation is suppressed by about a factor of two compared to thermal dump, with mechanical feedback being somewhat stronger. Kinetic feedback shows the same qualitative behaviour as in the lower-mass galaxy, with an initially high SFR that declines steadily due to gas outflows. Again, delayed cooling gives SFRs that are much lower than for the other models.
4.2.2 The Kennicutt-Schmidt relation
In the local Universe, SFR surface densities, , are observed on large scales to follow the Universal Kennicutt-Schmidt (KS) relation, , where is the gas surface density (Kennicutt, 1998). We plot in Fig. 4 the relation between SFR and gas surface densities in our simulations at Myr for the different feedback models, and compare it with the empirical relation shown as a solid line (normalised for a Chabrier IMF, see Dalla Vecchia & Schaye 2012). In this plot we include results from both the low-mass g9 and high-mass g10 galaxies, in order to show a wide range of surface densities, and to demonstrate how the feedback efficiency changes for each model with galaxy mass, metallicity, and physical resolution. Results for the g9 galaxy are shown with smaller opaque symbols, while the g10 galaxy is represented by larger and more transparent symbols. The gas and SFR surface densities are averaged along annulli around the galaxy centre, with equally spaced azimuthal bins of pc, and we only include gas within a height of kpc from the disc plane ( kpc in the case of the g10 disc).
All feedback models, and even the case of no feedback, produce slopes in the KS relation in rough accordance with observations at gas surface densities substantially above the ‘knee’ at , though the slopes tend to be slightly steeper than observed. The similarity to the observed slope is in large part a result of the built-in star formation model, . However, even though all simulations have the same sub-resolution local star formation efficiency of , the normalisation varies by about an order of magnitude, with delayed cooling being most efficient at suppressing the star formation for any given gas surface density, owing to the large scale height of the disc. At high gas surface densities (), all methods predict too high , except delayed cooling which predicts too low values.
For the lower-mass g9 galaxy (smaller solid symbols), thermal dump, stochastic, mechanical feedback, and delayed cooling are all similar in the KS plot, though delayed cooling does not produce as high gas surface densities as the other models. Kinetic feedback has significantly higher SFR surface densities for given gas surface densities (but relatively low maximum gas surface densities), owing to the very thin disc produced by the almost decoupled injection of momentum.
For the more massive g10 galaxy, which was simulated with lower resolution, the picture is quite different (large transparent symbols). With thermal dump feedback, the SFR surface densities shift significantly upwards and the relation is quite similar to the no feedback case. Stochastic feedback, and to a lesser extent, mechanical feedback, also shift upwards, away from the observed relation. For kinetic feedback, the relation is however almost unchanged in the more massive galaxy (except for low gas surface densities, where it is higher), but consistently remains about a factor two above the observed relation. With delayed cooling, the gas surface densities become much higher than in the lower mass galaxy, but the SFR surface densities are significantly lower than observed.
For delayed cooling, we can calibrate the available free parameter to improve the comparison to observations. Halving the delayed cooling time-scale in the g10 galaxy, to the same value as used for the g9 galaxy, results in a KS relation which is very close to the observed one. For the other models, we cannot calibrate the feedback parameters to close in on the observed relation, and other measures are required, such as increasing the feedback energy per unit stellar mass. Another option is to reduce the star formation efficiency parameter, , in which case a fair match to observations can be produced, but at the cost of making the feedback insignificant compared to the no feedback case in terms of morphology, total SFR, and outflows (the feedback essentially all becomes captured inside ).
Galactic outflows are a vital factor in delaying the conversion of gas into stars. Feedback processes in the ISM are thought to eject large quantities of gas from the galaxy, some of the gas escaping the gravitational pull of the galactic halo altogether. Most of the gas, however, is expected to be ejected at velocities below the halo escape velocity and to be recycled into the disc. Galactic outflows are routinely detected in observations (e.g. Steidel et al., 2010; Heckman et al., 2015), and while the outflow speed of cold material can be fairly accurately determined, other properties of the outflows are not well constrained, including the mass outflow rate, the fraction of gas escaping the halo, the density, and thermal state of the gas.
Outflows are often characterised in terms of the mass loading factor, which is the ratio of the outflow rate and the rate of star formation in the galaxy. Its definition is somewhat ambiguous, as it depends on the geometry and distance from the galaxy at which the outflows are measured, which is hard to determine in observations. Observational works have inferred outflow mass loading factors well exceeding unity (see e.g. Bland-Hawthorn et al., 2007; Schroetter et al., 2015), and many theoretical models require mass loading factors of in sub-L galaxies to reproduce observable quantities in the Universe (e.g. Puchwein & Springel, 2013; Vogelsberger et al., 2013; Barai et al., 2015; Mitra et al., 2015).
It is therefore important to consider outflow properties when evaluating SN feedback models. Models that produce weak or no outflows, with mass loading factors well below unity, could be at odds with current mainstream theories of galaxy evolution (although it is not known whether SN feedback is directly responsible for outflows – e.g. cosmic rays could play a major role; Booth et al. 2013; Hanasz et al. 2013; Salem & Bryan 2014; Girichidis et al. 2016).
In Fig. 5 we compare the time-evolution of outflows from the g9 galaxy with the different SN feedback schemes999We show outflow plots for the g9 galaxy only, but we comment on outflows in the more massive g10 galaxy (which have similar properties) at the end of this subsection.. We measure the gross gas outflow (i.e. ignoring inflow) across planes parallel to the galaxy disc, at a distance of kpc in the left panels and further out at kpc in the right panels. The top row of panels shows the mass outflow rate across those planes (), the middle row shows the mass loading factor (), and the bottom row shows the mass-weighted average of the outflow velocity perpendicular to the outflow plane ().
In terms of outflows kpc above the disc plane (left panels of Fig. 5), kinetic feedback is strongest, with and slightly above unity. Delayed cooling produces a substantially lower outflow rate, but since the SFR is also much lower, the mass loading factor is higher. The other feedback models give much lower outflow rates, and have mass loading factors . At a larger distance from the disc of kpc (right panels of Fig. 5), the situation is quite similar. All models except for kinetic feedback have declining outflow rates and mass loading factors, owing to the strong initial starburst, which can be seen to result in an outflow rate peaking around Myr.
In the bottom panels of Fig. 5 we compare the average outflow velocities to the DM halo escape velocity101010The escape velocity estimate ignores the contribution of baryons. Hence, it is an underestimate, that is likely non-negligible close to the disc, but insignificant at kpc.,
where (e.g. Mo & Mao, 2004), which has been marked with horizontal grey solid lines. Close to the disc, the average velocity for kinetic feedback is marginally higher than escape, but slowly declining due to the declining SFR. For the other feedback models, the outflow velocity is well below escape. Ten times further out, the mean outflow velocities are considerably higher for all feedback models. This is to some degree due to the initial starburst, which ejects high-velocity outflows early in the simulation runs, and to some degree due to gas at lower velocities not having reached kpc and thus not contributing to the velocity average. In any case, these (average) velocities are still below the escape velocity, again with the exception of kinetic feedback. This implies that for all models except kinetic feedback, most of the outflowing gas will not escape to infinity, but instead fall back on the galaxy where it will eventually produce stars. Kinetic feedback does give the gas high enough velocity so that in principle it can escape the halo entirely, while in practice this may be complicated by CGM and IGM gas which stands in the way and needs to be swept out.
The main message of Fig. 5, however, is not the escape velocity, but the low mass loading factors for thermal dump, stochastic, and mechanical feedback, far below the order unity inferred from observations of local galaxies. As before, we see a strong similarity between the results produced by thermal dump, stochastic, and mechanical feedback.
In Fig. 6, we study how the outflow properties kpc from the disc scale with the local gas surface density. The panels show, from top to bottom, gross outflow rate per unit area (), mass loading factor , and mass-weighted gross outflow velocity. We split the face of the disc into a kpc wide grid of pc squares (that is, squares), and extract the outflow properties in each square. In the plots, the outflow properties are binned by the gas surface density, and the shaded regions show the logarithmic standard deviation in each bin.
For each model, the general trend is that higher gas surface densities correspond to higher outflow rates and velocities, but lower mass loading factors. The outflow velocities are noticeably higher than those (at Myr) in Fig. 5. The reason is that Fig. 5 shows the mass-weighted average of all outflowing gas, while Fig. 6 is restricted to a kpc wide square plane directly above and below the disc, hence capturing the more collimated part of the outflows. Kinetic feedback clearly stands out as having the highest outflow velocities, peaking close to km/s (at the peak surface densities), with little scatter. With kinetic feedback, almost all of the outflowing gas directly above or below the disc is moving faster than the escape velocity (indicated with a horizontal solid line). The other feedback models produce outflow rates and velocities that are alike (within roughly a factor of two), and much lower than for kinetic feedback, with the exception of delayed cooling, which has a massive, slow outflow, and the highest mass loading factor.
In Fig. 7 we show images of slices along the -plane (at ) at Myr, with each set of two panels showing the hydrogen density (left) and temperature (right) for a given feedback model. Delayed cooling and kinetic feedback clearly stand out here. The former yields dense and cold ( K) outflows. The outflows for the latter are diffuse, hot (), and the most extended (which is expected since Fig. 6 shows they are by far the fastest). The remaining three feedback models produce qualitatively similar multiphase ( K) outflows. The clear distinction between the most effective feedback model, i.e. delayed cooling, and the other, less effective models, in the outflow properties, could be used in future work as an observational probe into how accurately those models represent actual feedback in galaxies.
4.3.1 Outflows from the more massive galaxy
For the more massive g10 galaxy, the outflow differences, which we do not plot, are qualitatively similar to those in the above analysis. Kinetic feedback gives the highest mass loading factor, which is again of order unity both at and kpc. All the other models give similar mass loading factors kpc above the disc as for g9, but in contrast to g9 the mass loading drops by orders of magnitude at kpc (the biggest drop occurring for delayed cooling), due to the stronger gravitational pull. The outflow velocities are slightly higher for all models, but they are still much (marginally for kinetic feedback) lower than the ( km/s) escape velocity.
4.4 Gas properties
In Fig. 8 we compare gas properties for runs with the different SN feedback models, using phase diagrams of gas temperature versus density at Myr. The colour scheme represents the mass of gas in each temperature-density bin. The mass-weighted mean density in each diagram is represented by a solid vertical line, while the mass-weighted mean temperature in each density bin is shown by solid red curves. Star-forming gas is enclosed by a dotted box, while gas with temperatures below the artificial Jeans temperature (which has been subtracted from the ‘thermal’ temperature plotted here) is indicated by the shaded diagonal region in the bottom right corner of each diagram.
We continue to see the same qualitative picture as before: delayed cooling and kinetic feedback stand out, while the remaining three models look similar. Delayed cooling yields by far the lowest mean density, , which is well below the star formation density threshold of , and almost two orders of magnitude below the mean density without feedback. The other models all yield mean densities near .
Delayed cooling produces the highest mean temperatures at intermediate densities of with a lot of gas at temperatures K, which is probably unphysical given the short radiative cooling times in this regime. Curiously, in more diffuse gas, delayed cooling produces by far the lowest temperatures, K, while in the same density regime the gas is typically at K with other feedback models (even in the absence of feedback). From identical phase diagrams excluding the galaxy disc, we have confirmed that this diffuse gas is primarily ‘CGM’ gas outside the disc. In the case of delayed cooling, these outflows, i.e. gas outside the disc, span a wide range of densities, , in stark contrast to the other feedback models, where the CGM gas reaches maximum densities of a few times . With kinetic feedback, the CGM has a clear bi-modality not seen for other models, with some of the gas following an adiabat starting around , K, and extending towards much lower densities, and the remainder at K and spanning densities of (the latter is flowing ’diagonally’ from the disc, i.e. at a steep angle from the axis of disc rotation).
All feedback models produce hot and relatively diffuse gas, populating the region K, in Fig. 8. One might expect this to belong to the outflowing CGM. However, if we exclude the disc (out to kpc in height and kpc in radius), all of this hot diffuse gas disappears from the phase diagrams, indicating that it in fact belongs to the ISM. In the case of thermal dump, stochastic, and mechanical feedback, the outflowing CGM is indeed warm to hot, but dilute, with , while kinetic feedback produces the aforementioned bi-modality in the outflowing CGM, and delayed cooling produces circum-galactic outflows that are predominantly at temperatures between and K.
4.5 Impact of feedback on the local environment
A major factor in any feedback model is how efficiently it clears away those dense clouds where stars can form. When dense regions are quickly cleared by early SN explosions in a stellar cluster, this can also boost the efficiency of subsequent SN explosions which then take place at lower densities where cooling is less efficient and the momentum obtained in the SN remnant increases (Eq. 12; see also Kim et al., 2014; Martizzi et al., 2015). SNe exploding in the diffuse ISM have been suggested to prevent the formation of star-forming clouds (e.g. Iffrig & Hennebelle, 2015), to maintain the hot volume-filling ISM, and to generate fast outflows (e.g. Ceverino & Klypin, 2009).
In Fig. 9 we show the gas densities at which stellar particles are born (dashed curves) and the densities at which the SN events take place (solid curves) for each of the feedback models in the g9 galaxy. Focusing first on the star formation densities, they are almost indistinguishable for all the models, and differ only for the case of no feedback, for which the stars typically form at significantly higher densities. This shows, as we have already seen from the morphological comparison in Fig. 1, that all the feedback models are efficient at preventing and/or destroying star-forming clouds in this g9 galaxy, and almost all the stellar particles are formed within one dex of the star formation density threshold of . For the no feedback case, the clouds can collapse to higher densities, impeded only by the density-dependent pressure floor.
For the densities at which the SN events take place, there are larger differences between the feedback models. Thermal dump, stochastic, and mechanical feedback are similar, with a non-negligible per cent of the SN energy injected below ( per cent below ), and SN events consistently taking place at lower densities than star formation.
Delayed cooling stands out, in having more SNe than those aforementioned models at densities , but fewer SNe for lower densities. This comes from the efficiency of delayed cooling in diffusing and thickening the ISM disc, resulting in a mass-weighted density distribution of the ISM (not shown) which peaks at , about a dex lower than for thermal dump, stochastic, and mechanical feedback, but a volume-weighted distribution (also not shown) which peaks at , about a dex higher than for those other models. Delayed cooling hence smooths out not just the density peaks in the ISM, but also the density throughs, such that stars form at lower average densities, but SNe explode at higher minimum densities than for the other models.
Standing out much more distinctly, kinetic feedback SNe explode in gas with almost exactly the same densities (and even higher) as the stars are formed, with almost no SNe exploding at lower densities. This signature of kinetic feedback suggests a decoupling between the injected momentum and the immediate environment surrounding the exploding stellar particle. Instead of quickly dispersing the sites of star formation, the gas is gradually transported away from those sites, out of the ISM, without coupling to the immediately surrounding gas. This explains the thin bubble-free gas disc (Fig. 1) and the gradual decline of the SFR (Fig. 2), which is due to slow gas depletion. These distinct features are however strongly dependent on the size of the bubble, , into which the SN momentum and mass is injected. The fiducial size pc results in this decoupling between the SNe and the surrounding gas. In §8 we experiment with a smaller bubble size ( pc) and find a very different behaviour for kinetic feedback, with results resembling those for thermal dump, stochastic, and mechanical SNe, i.e. much lower outflow rates, a flatter star formation rate with time, and a thicker disc.
4.6 Similarity of three models in the low-mass galaxy
For the low-mass g9 galaxy (but not for the more massive g10 galaxy), we have seen that the results for thermal dump, stochastic, and mechanical feedback are near identical in terms of morphologies, star formation, and gas properties.
In Eq. (4) we showed that the probability for stochastic feedback events is above unity at the density threshold for star formation () in the g9 galaxy, and in Eq. (15) we saw that, also at , mechanical feedback remains in the resolved adiabatic phase. In addition, we found in §4.5 that SN events do typically take place at densities close to . Hence there appears to be no significant numerical overcooling issue in the g9 runs, and it is then no surprise to find similar results for thermal dump, stochastic, and mechanical feedback. It can be argued that the adiabatic phase of thermal dump feedback is resolved. Note that this may imply that delayed cooling and kinetic feedback, for the fiducial parameters we have chosen, converge to wrong results for the effects of SN feedback.
For the more massive g10 galaxy, this is not the case: these aforementioned feedback models give quite different results (Fig. 3), and thermal dump does not do much to suppress star formation compared to the no feedback case. Purely from Equations 4 (stochastic probability) and 15 (mechanical feedback phase), one might expect the adiabatic phase to be resolved here as well, since the changes in stellar mass and resolution, compared to the g9 galaxy, cancel out. However, in part due to stronger metal cooling, but more importantly due to the larger gravitational potential of the g10 galaxy, stars form and explode at significantly higher densities than in the g9 galaxy. Hence the probability for stochastic feedback events becomes lower than unity (on average in the stochastic feedback run), most mechanical feedback events become pure snowplow momentum injections, and thermal dump becomes a victim of numerical overcooling.
Kim et al. (2014) derived a density limit at which the momentum created by single thermal dump type II SN explosions is resolution-converged with grid-hydrodynamics. They found that convergence is maintained with a cell width , where . Taking we would need a resolution of pc for converged thermal dump feedback. That is a considerably higher resolution than ours ( pc), so naively we would expect overcooling to be significant for the g9 thermal dump simulation. In light of the above finding, that thermal dump feedback appears more or less resolved in the g9 galaxy, the lack of resolution according to Kim et al. (2014) is mitigated by the fact that each stellar particle in our g9 simulations releases the equivalent of 40 type II SN explosions instantaneously, instead of one, as assumed in Kim et al. (2014).
4.7 SN model comparison summary
For the low-mass g9 galaxy, the results for thermal dump, stochastic, and mechanical feedback are near identical in terms of morphologies, star formation, and gas properties. This is an indication that the adiabatic phase of SN explosions is resolved. For the more massive g10 galaxy, thermal dump is significantly weaker than stochastic and mechanical feedback in suppressing star formation (though not so much in generating outflows).
Delayed cooling is by far most efficient at suppressing star formation and yields results closest to the observed Kennicutt-Schmidt relation (at least for our assumed star formation efficiency.
Thanks to a large fiducial ‘bubble radius’ of pc for momentum and mass injection, kinetic feedback has the highest outflow rates and a mass loading factor of order unity. Delayed cooling follows, with weaker outflow rates but a slightly higher (but declining) mass loading factor. The other feedback models produce much lower outflow rates and mass loading factors than those two more efficient models. In the more massive ( MW) g10 galaxy, the mass loading factor is similar to that of g9 close to the galaxy plane, but drops by orders of magnitude ten times further out at kpc, for all models except kinetic feedback, which maintains a mass loading factor of unity.
The feedback models producing the lowest outflow rates and mass loading factors produce hot and dilute outflows, while delayed cooling yields distinctively cold and dense outflows. For kinetic feedback the outflows have a clear bi-modal phase structure, with relatively cold and dense outflows close to the disc, and hot and diffuse outflows further out following a temperature-density relation suggesting adiabatic cooling.
Given the large fiducial bubble radius, which effectively decouples feedback from the ISM, kinetic feedback produces by far the fastest outflow, some of it above the escape velocity. All other models produce outflow velocities about an order of magnitude lower, well below the escape velocity.
5 Resolution convergence
Resolution convergence is an important factor in assessing SN feedback models. Ideally, the effects of feedback should remain constant if the resolution is increased, or at least if it is varied within reasonable limits, i.e. within a factor of a few111111Convergence with a dramatic change in resolution, on the other hand, is usually not a desired goal. With much lower resolution, a sub-grid model becomes meaningless as the structures with which the model is supposed to interact become completely unresolved and a ‘lower level’ of sub-grid physics must take over the current ones. With much higher resolution, some of the real physics become resolved, and the sub-grid model becomes irrelevant (though ideally it should then converge towards a ‘first-principles’ methodology).. In practice such constancy, while desirable, is hard to obtain without making significant sacrifices, such as disabling physical processes like hydrodynamical interactions or radiative cooling in the ISM. A second best choice is a small and predictable change with varying resolution, so the feedback parameters can be easily calibrated for different setups in order to achieve “weak convergence” (Schaye et al., 2015).
In this section we aim to understand how and to what extent measurable galaxy properties change with resolution for the different feedback models. For this purpose, we use a lower-resolution version of the g9 galaxy, which we call g9lr. The setup is identical to g9, except that the minimum cell width is two times larger, i.e. pc, and the particle mass (in the initial conditions as well as for new stellar particles) is eight times higher (i.e. for stellar particles, for DM particles). For simplicity, and because the Jeans length is already resolved by cell widths in the higher-resolution runs (and hence by cell widths in the lower-resolution ones), we leave the pressure floor and star formation threshold unchanged.
In the left column of panels in Fig. 10, we plot, for each feedback model, the ratios of SFRs (upper panels) and mass loading factors kpc from the disc (lower panels) for g9lr over g9 runs.
For delayed cooling, the resolution has a significant effect on the relative SFR, although it should be kept in mind that the SFR for delayed cooling is quite small in the first place (and hence the absolute change is low compared to the SFRs of other models). For other models, the SFRs change insignificantly with resolution, though there is a systematic tendency of slightly lower SFR with lower resolution.
The lower resolution has a larger effect on the outflow rates, shown in the bottom left panel of Fig. 10. Decreasing the resolution systematically increases the outflow rates for mechanical feedback (by up to a factor four), thermal dump, stochastic feedback, and delayed cooling, (by roughly a factor two), while for kinetic feedback the outflow rate is almost unchanged. The outflow rate increase is likely connected to the winds being launched on larger scales, due to the larger cell width and mass (note that the specific energy, i.e. the ratio between the SN energy and receiving gas mass is similar to the higher resolution run, since the particle mass is times larger). The obvious exception is kinetic feedback, where the energy is distributed within a bubble radius which we have kept constant, and indeed the outflow rate remains unchanged.
We also consider the effect of ‘SN’ resolution, where we lower the grid resolution (and that of the initial conditions particles), just as in the g9lr runs, but keep the mass of newly formed stellar particles fixed compared to the g9 runs. These runs, which we call g9lr_m, give another measure of resolution convergence for the feedback models, as the specific energy per feedback event (i.e. SN energy over the local gas mass) is reduced by a factor of eight compared to g9 runs, while it was kept constant in Fig. 10. We compare those to the g9 runs in the right column of panels in Fig. 10, showing the ratios of SFRs and mass loading factors for the feedback models. For the SFR, the largest difference occurs for thermal dump feedback, which becomes much less efficient at suppressing star formation rates due to the lower particle masses. The adiabatic phase becomes severely unresolved, and there is no mechanism built into the model to compensate. Other feedback models maintain similar average star formation with the lower specific energies. The mass loading factors increase somewhat if we decrease the spatial resolution and the specific SN energies (bottom right panel of Fig. 10), but at Myr the increase is smaller than for fixed specific energies (bottom left panel). Delayed cooling is an exception, showing a decrease in mass loading at some time intervals, but an increase in others, which is likely just caused by the relatively large variations in the SFRs.
With thermal dump feedback, resolution non-convergence is a well known problem. With lower resolution, a larger gas mass is heated in a single feedback event, resulting in lower initial temperatures given a fixed SN energy. This would normally lead to higher SFRs, but in the left panels of Fig. 10 the effect is counter-balanced by the use of more massive stellar particles, slightly increasing the feedback efficiency due to the higher SN energy per feedback event.
For the case of stochastic feedback, the fairly good convergence of the SFR with both spatial resolution and stellar particle mass is not a big surprise, since the stochasticity is built in to ensure that the heated gas receives a fixed specific energy, regardless of cell size, density, and particle mass. The outflow rates are relatively poorly converged for stochastic feedback and low stellar particle mass, which likely comes from the aforementioned tendency for the outflow rate to increase with lower resolution, and hence larger launching scales, even if the specific SN energy is constant. Mechanical feedback was shown by Kimm & Cen (2014) to converge well with resolution in terms of the final momentum reached, and indeed the whole point of the model is to maintain the same momentum injection regardless of whether the momentum buildup is captured numerically or not. While we confirm that the convergence is good for the SFR, the outflow rates are not very well converged, neither in terms of spatial resolution nor stellar particle mass. For delayed cooling, resolution convergence is not an obvious property, but if cooling is turned off long enough, the energy (i.e. cooling) losses should become insignificant and hence independent of the resolution. The fact that the SFR (and to a lesser extent the mass loading factor) is poorly converged for delayed cooling hints that cooling losses are still significant, but we remind that the absolute change is actually lower than for the other feedback models. As long as kinetic feedback is sufficiently decoupled from the ISM surrounding the SN explosion due to the use of a large SN bubble radius, it is not surprising to see good resolution convergence, since the main effect of the feedback is then simply to slowly deplete the disc of gas mass.
In summary, except for thermal dump feedback, the feedback models converge fairly well in terms of SFRs. However, none except for kinetic feedback (with a bubble radius larger than the disc height) converge well in terms of the outflow mass loading factor, with the mass loading factor decreasing with higher resolution.
6 Varying the temperature jump for stochastic feedback
We now examine the effect of varying the stochastic heating, parameter, . We do not go lower than the fiducial value of K, because already here the feedback is not very stochastic: in the g9 galaxy the average probability for SN candidates (which, as a reminder, is the ratio of the available SN energy of the stellar particle to the energy required to heat the host gas cell by ) is ( in g10). Lower would lead to order unity probabilities for SN explosions, converging towards thermal dump feedback.
We thus consider three values for in addition to the fiducial one, each time increasing the injected specific energy by half a dex, i.e. , and K.
Fig. 11 shows maps of the hydrogen column density for and K. For the case of K, there is not a significant difference from the fiducial run (the middle left panel in Fig. 1), though the gas disc becomes slightly thicker and clumpier. For K, the difference is much clearer. Here the face-on disc is more diffuse overall, but at the same time it contains more dense clumps and filaments, especially close to the centre and at the disc outskirts. This is due to the increased stochasticity of SN explosions: the low probability for SN events, on average , allows dense star-forming clumps to live longer before they are hit by the first SN explosion121212We note that Dalla Vecchia & Schaye (2012) advise against such high values of that probabilities for feedback events become , which is clearly the case here.. This effect is amplified at the outskirts, where there is relatively little star formation and thus a low rate of SN explosions per unit volume. Looking at the disc edge-on, there is much more structure in the CGM for our maximum value of .
In Fig. 12 we zoom out and look at the large-scale outflows in the ‘maximum’ case of K. The outflows are dramatically different from the fiducial setting for stochastic feedback shown in the top right panel of Fig. 7: they are much denser, clumpier, (mostly) hotter, and more extended. Indeed, looking at the dashed lines in Fig. 13, we see that the outflow rate at kpc increases almost linearly with , and is about times higher than the SFR (solid) at the end of the run for K. The average outflow velocity (not shown) increases by a factor at kpc for the highest considered (compared to the lowest), but is unaffected at kpc.
Varying has a much weaker effect on the SFR than on the outflow, as shown in Fig. 13. The ‘first’ two increases in have almost no effect on the SFR, while the highest value produces an initially higher SFR which then declines gradually, much like in the case of kinetic feedback, and ends up significantly lower than for lower values. As with kinetic feedback, the decline in SFR is likely due to the strong outflow depleting the galaxy of fuel for star formation. Also, due to the lower heating probability, star-forming clumps can continue to form stars longer without having strong SN events disrupting them.
We recall that the SN energy is not directly redistributed to lower gas densities with stochastic feedback (see Eq. 6). On the contrary, we find that increasing indirectly results in the SN energy being deposited at higher gas densities due to the increased clumpiness of the gas (not shown; K results in energy being deposited at dex higher densities, compared to the fiducial case). Hence the stochasticity increases feedback efficiency purely by increasing the injected energy of a given SN event at a given density (while the total SN energy over the galaxy is unchanged), and hence also the local cooling time.
In summary, increasing the stochasticity of feedback by increasing strongly increases the outflows. The SFRs are insensitive to these stochasticity variations except for the highest value of considered. With very large stochasticity, the disc also becomes increasingly clumpy.
7 Varying the time-scale in delayed cooling feedback
Of the feedback models we have compared in §4, with their fiducial setup parameters, delayed cooling feedback suppresses star formation most strongly and yields the highest mass loading factors for the outflows (see Figures 2, 3, 4, and 5).
Since delayed cooling feedback in its fiducial setup is strong, we examine here what happens if we reduce the value of its free parameter, which is the cooling delay time-scale, . The fiducial value is Myr, so here we compare with two runs with significantly lower delay time-scales of and Myr. We find that the results are highly sensitive to variations in . In Fig. 14 we show the g9 SFRs and gross outflow rates across planes kpc from the disc, for those variations. As expected, a shorter delay time-scale strongly decreases the feedback efficiency, producing higher SFRs and lower outflow rates. However, even for the shortest delay time-scale considered here, the SN feedback is still more efficient in terms of suppressing star formation than any of the other feedback models (and second only to kinetic feedback in terms of outflow rates). Note that the mass loading factor is particularly sensitive to the value of the free parameter. At Myr it decreases from for Myr to for Myr (both at and kpc from the disc). The outflow velocities, which we do not show here, are unaffected at kpc, and are fairly insensitive to at kpc (few tens of percent increase in outflow velocity from smallest to highest ).
We analysed the KS relation for the same variations in (not shown). Decreasing produces a KS relation which looks increasingly like that of mechanical feedback (for the g9 galaxy) in Fig. 4, with a similar slope and maximum gas surface density of for Myr. Morphologically, the galaxy with the shortest dissipation time-scale also looks quite similar (though a bit more diffuse) as the runs with mechanical, stochastic, and thermal dump feedback in Fig. 1.
While the outflow rates decrease for shorter delay times-cales, the morphological and phase structure of the outflows remain qualitatively similar.
8 Variations in kinetic SN feedback
In Fig. 15, we show variations in the SFRs and gross mass outflow rates at kpc for the g9 galaxy, for variations in the free kinetic feedback parameters, which are the bubble radius (fiducially pc) and the sub-resolution mass loading parameter (fiducially set to ).
Decreasing the bubble radius, , to about two times the minimum cell width (of pc), has a significant effect on the SFR and dramatically suppresses the gas outflow rate, which both (and also the outflow speed) become similar to those produced by thermal dump, stochastic, and mechanical feedback. Morphologically (not shown), the runs with the small bubble radius also resemble the runs with these other feedback models.
The smaller bubbles are less efficient at driving large-scale outflows, because a larger fraction of the energy is now deposited into dense ISM gas. For the same reason, the SFR at early times decreases for smaller bubbles. At late times the SFR is higher because there remains substantially more gas in the galaxy as a consequence of the weaker large-scale winds.
Increasing the sub-resolution mass loading, , by an order of magnitude has little effect on star formation (Fig. 15), outflows (rates and speeds), and morphology (not shown). This is because most SN events occur quite close to the density threshold for star formation (see Fig. 9), meaning that usually there is not much more mass available in the host cell than roughly matches the stellar particle mass. Hence, an initial mass loading of more than is not possible, since it would mean removing more mass from the host cell than is available for re-distribution to the SN bubble. To investigate the effect of lower stellar particle masses, or star formation happening well above the density threshold, we performed runs where we used a three times lower stellar particle mass, . Here increasing reduces the outflow rates significantly (as it does to a smaller extent in Fig. 15 for the smaller bubble radius). Moreover, the feedback becomes more coupled to the disc, which becomes much thicker. This can be understood from the fact that with a higher mass loading, i.e. a larger ejected mass, the velocity of the ejected gas must decrease to conserve momentum, and hence it is less likely to escape from the galaxy.
We also looked at the KS relation for these variations in kinetic feedback (not shown). Here, again, varying has a negligible effect, while the smaller bubble radius gives a relation very similar to e.g. mechanical feedback.
The effects of the SN feedback models we have studied are sensitive to the resolution and/or mass of the simulated galaxies. For our lower-mass g9 galaxy, thermal dump, stochastic, and mechanical feedback give very similar results in terms of SFRs, the Kennicutt-Schmidt relation, and outflow mass loading factors, and we have argued that this is an indication that the adiabatic phase of thermal SN feedback is ‘resolved’ in this galaxy, with the caveat that we inject the equivalent of individual SN explosions instantaneously. For our (ten times) more massive g10 galaxy, however, this no longer applies, and thermal dump feedback is significantly weaker than any of the other models.
If we compare our results to observations of SFRs in the local universe, delayed cooling comes closest to reality. Observed SFRs for local star-forming galaxies of similar stellar masses to g9 are typically in the range (Figure 11 in Chang et al., 2015). In Fig. 2, we see that delayed cooling is the only feedback model giving SFRs matching those observations, while the other models give values well above the upper limits. For stellar masses similar to our g10 galaxy (), Chang et al. (2015) give SFRs in the one-sigma range of , which actually falls slightly below the delayed cooling SFR in Fig. 3. While such a comparison to observations gives some indication of what is required to produce a realistic suppression of star formation, we stress that comparing the results from our isolated and somewhat idealised setup to observed star formation rates has many caveats. Perhaps most significantly, the gas fractions in our galaxies are high compared to those of local galaxies (see e.g. the compilations in Bahé et al., 2016; Sales et al., 2016).
Such caveats matter less if we consider the Kennicutt-Schmidt relation. Here, all models except delayed cooling produce a slope that is too steep compared to observations and/or SFR surface densities that are significantly too high for a given gas surface density. Delayed cooling does relatively well in terms of slope but for our fiducial parameter choice it undershoots the SFR surface density. However, unlike for the other models, this can be calibrated towards a reasonable fit to the observed KS relation by tuning the delayed cooling time-scale.
In addition, delayed cooling and ‘decoupled’ kinetic feedback are the only models able to produce mass loading factors exceeding unity, while for other models with their fiducial parameters the mass loading factors are at best an order of magnitude below unity.
While one can argue that these factors make delayed cooling for these observables empirically most successful, the model produces unrealistic thermal signatures in the ISM and CGM, where large amounts of gas occupy a region in temperature-density space where the cooling time is very short. Moreover, the convergence between thermal dump, stochastic, and mechanical feedback suggest the adiabatic phase is resolved and hence that the results from delayed cooling and kinetic feedback may be unphysical.
Mechanical feedback could be argued to be most physically motivated, being based on analytic calculations and high (sub-pc) resolution simulations of the final momentum attained in various environments in terms of gas density and metallicity. Yet it does not come near observations in suppressing star formation and produces weak outflows in our simulations. Perhaps this discrepancy can be traced to the idealistic assumptions made when deriving the mechanical momentum (Eq. 12). While the momentum is realistic and converged for a homogeneous medium, such a medium is not a good description of a multi-phase and porous feedback-regulated star-forming ISM, and may lead to an under-prediction of the generated momentum (e.g. Kimm et al., 2015). We also neglect the preprocessing of the local environment by stellar radiation, which lowers the surrounding gas densities and has been shown to increase the momentum injection from stellar populations, typically by a factor (e.g. Geen et al., 2015).
It was recently reported by Gentry et al. (2016) that idealised experiments of SN feedback in the literature have under-predicted the final momentum by up to an order of magnitude, due to i) the neglect of multiple successive SN explosions, ii) a lack of resolution, and iii) a preference for Eulerian hydrodynamical solvers, which are argued to suffer from over-mixing and hence over-cooling. If these results are confirmed, we may have much more momentum to play with.
Another culprit may be the simplified setup of our simulations. For control and to reduce the computational cost, we used isolated galaxy simulations, assuming an initial state of the galaxy and its dark matter halo rather than letting it evolve naturally in a cosmological context. Ignoring environmental factors such as mergers and gas accretion, may change the behaviour of SN feedback.
Finally, there are other feedback processes at play in galaxies that we have neglected, such as AGN (thought to be important at high, larger than MW masses; e.g. Crain et al. 2015), radiation pressure (e.g. Rosdahl et al., 2015), and cosmic rays (e.g. Booth et al., 2013; Hanasz et al., 2013; Salem & Bryan, 2014; Girichidis et al., 2016; Simpson et al., 2016). The efficiency and interplay of each of those processes is not well constrained, but they likely provide an additional suppression of the SFR on top of SN feedback. If they turn out to play an important role in galaxy evolution, some empirically calibrated sub-resolution models for SN feedback may be re-interpreted to also represent other feedback processes at play in galaxies.
9.1 Dependence on factors other than SN feedback
Galaxy evolution involves a complex interplay of many physical processes, and the SN feedback efficiencies that we have reported may be sensitive to factors other than the SN feedback models.
One of the most important choices in our simulations aside from the implementation of feedback is the local star formation efficiency, which we have set to . While this is a normal choice in simulations of galaxy evolution, Agertz & Kravtsov (2015) argued that such a low value artificially suppresses the effects of feedback, and that higher local SF efficiencies of are needed to allow feedback to self-regulate star formation. We have experimented with in the g9lr galaxy (while keeping the other parameters fixed). While a full analysis is beyond the scope of the present work, it is worth mentioning the qualitative effects. For all SN feedback models considered, it results in a high early SFR peak, followed by a dip and then convergence towards similar but somewhat lower SFRs compared to those obtained for the fiducial . However, the SFR surface density moves in the wrong direction, i.e. away from the observed KS relation, both in terms of slope and normalisation.
We also investigated the effect of reducing , to find what calibration is required to reproduce the observed KS relation. A value of , i.e. ten times lower than our fiducial value, produces a good match to the observed relation, but at the cost of making the self-regulating effect of SN feedback negligible. In other words, the parameter takes over as the dominant feedback mechanism when set to a very low value.
While our non-comprehensive probes of the effect of a varying numerical star formation efficiency are discouraging, we have thus far not studied variations in in combination with other factors. For example, in combination with a high local star formation efficiency, early feedback may play a significant role by suppressing runaway star formation for the Myr from the birth of the first stars until the onset of SN explosions in a given star-forming region (e.g. Hopkins et al., 2012; Stinson et al., 2013; Agertz & Kravtsov, 2015). Also the locality of star formation may be important for the efficiency of feedback, i.e. it may matter whether the distribution of star formation, and hence SN explosions, is scattered smoothly in space and time, or happens in short localised bursts (e.g. Federrath & Klessen, 2013; Hopkins et al., 2013; Walch et al., 2015).
We will investigate the role of in more detail in future work, where we combine higher local SF efficiencies with more stringent SF criteria and early feedback in the form of stellar radiation.
There are many more variations which may or may not affect the feedback efficiency, and again we have made limited explorations, which we summarise below.
We varied the density threshold for star formation, , by a factor of ten in either direction, for stochastic feedback. Higher does give lower initial SFRs, but after Myr the values are the same to within a few percent for two orders of magnitude variations in . Outflow rates remain nearly unchanged.
In the more massive galaxy we ran identical simulations with Solar metallicity, instead of the default Solar metallicity. The reduced metallicity has a marginal effect on the outflows and on the SFRs in the case of delayed cooling and kinetic feedback. However, for thermal dump, stochastic, and mechanical feedback, it reduces the SFRs by a few tens of percent and increases the outflow mass loading factor by about a factor of two.
Switching from non-equilibrium hydrogen and helium thermochemistry to an equilibrium assumption in the g9 galaxy significantly increases outflow rates for thermal dump (factor ), stochastic feedback (factor ), and mechanical feedback (factor ), while the SFRs are only marginally reduced (and there is almost no effect on either SFRs or outflows with kinetic feedback or delayed cooling). In a forthcoming paper, we will explore the physical effects of equilibrium versus non-equilibrium thermochemistry in the context of SN feedback.
With thermal dump and stochastic feedback we scaled the Jeans pressure floor by a factor of three in each direction (i.e. three times higher and lower pressure at a given gas density). The disc morphology becomes slightly but noticeably smoother with a higher Jeans pressure, but the effect on star formation and outflows is negligible. A very large increase of the pressure floor suppresses the effect of SN feedback with any model, since it overtakes the role of SN feedback in suppressing the collapse of gas to star-forming densities.
With mechanical feedback, we explored the effect of runaway OB stars in the g9 galaxy, giving a kick to each newborn stellar particle with a random direction and random speed of . This has negligible effect on the SFRs, but the outflow rates are enhanced by almost a factor of ten, due to SN explosions taking place at significantly lower densities on average.
Also with mechanical feedback, we simulated individual erg SN explosions, stochastically sampled for each stellar particle over the Myr lifetimes assumed for massive stars in the (Chabrier) IMF (see details in Kimm et al., 2015). The effect is similar to the above test with runaway OB stars: the SFRs are marginally higher, while the outflow rates are increased significantly, though somewhat less than for the runaway OB stars. The reason for the increased outflow rates is that for a given particle, early SN explosions can clear away the dense environment, leading to late SN explosions taking place at reduced densities. We also combined runaway OB stars and individual SN explosions in the same run (again for the g9 galaxy). This produces SFRs marginally higher than individual SNe only and outflows marginally higher than runaway OB stars only, i.e. it does not give an extra boost to the outflows.
Changing the time delay for stochastic SN feedback (zero to Myr) has similar effects as runaway OB stars: The SFRs are only marginally affected, while an increased delay increases the outflow mass loading factor (it is halved for zero delay and doubles for a Myr delay).
9.2 Comparison with stochastic heating in GADGET
The stochastic feedback model included in this work is based on the scheme introduced in DS12 and used in the Gadget code. Similarly to this work, DS12 explore the effects of their stochastic feedback model using rotating isolated galaxy discs of two masses, the main difference being that their lower mass disc is roughly ten times less massive than our g9 disc (while their higher mass disc is comparable in mass to our g10). They compare different values for stochastic heating and find K to be efficient in suppressing star formation (in the massive disc) and generating strong outflows (in both discs). In DS12 this value, which we also choose as our fiducial value for stochastic heating, gives outflow mass loading factors (at of the virial radii) of and for the low- and high-mass galaxies, respectively. This is times larger than our mass loading factors, pointing to a significant difference between our simulations and DS12 in the efficiency of stochastic feedback in suppressing star formation and generating outflows.
In Appendix A we discuss the differences between our simulations and those of DS12 and attempt to more closely reproduce their setup. We conclude that there is a major difference between the two versions (AMR and SPH) of stochastic feedback, with the SPH version being much more efficient at generating outflows, for a given SFR. Pinpointing the reason(s) for this difference, however, and whether it is due to subtle differences in setup or resolution, or to more fundamental differences between AMR and SPH, remains a challenge for follow-up work.
We used simulations of isolated galaxy discs with the Ramses code to assess sub-resolution models for SN feedback in AMR simulations, in particular their efficiency in suppressing star formation and generating outflows. We focused our analysis on a dwarf galaxy, ten times less massive than the MW, using a spatial resolution of pc and a stellar (DM) mass resolution of (), but also included a more limited analysis of a MW mass galaxy (using a resolution of pc, stellar particles and DM particles).
We studied five SN feedback models: i) thermal dump of SN energy into the host cell of the star particle, ii) stochastic thermal feedback, where the SN energy is re-distributed into fewer but more energetic explosions, iii) kinetic feedback, where momentum is deposited directly into a bubble around the star particle, iv) delayed cooling, where cooling is suppressed temporarily in the expanding SN remnant, and v) mechanical feedback, which injects energy or momentum depending on the resolution. Three of those models can be calibrated with adjustable parameters, which are the minimum local heating temperature for stochastic feedback, the bubble size and local mass loading for kinetic feedback, and the cooling suppression time for delayed cooling). The mechanical feedback model has no free parameters (once the SN energy has been decided) and the injected momentum is based on analytic derivations and high-resolution simulations of cooling losses in expanding SN blasts. We compared the results produced using these models with their fiducial settings, and for those models with adjustable parameters we studied the effects of parameter variations. Our main results are as follows.
For our low-mass, high-resolution galaxy, thermal dump, stochastic, and mechanical feedback produce nearly identical results (Figs. 2 and 5). We showed that at our current resolution and star formation densities, stochastic feedback is actually not that stochastic, and mechanical feedback is still mostly in the adiabatic phase. Hence those feedback models are converged in that setup, and thermal dump feedback adequately resolves the energy injection (by multiple SNe in a single event). For our more massive galaxy, stochastic and mechanical feedback become significantly stronger than thermal dump feedback, but are still weak compared to delayed cooling and kinetic feedback (Fig. 3).
Strong outflows are not easily generated in our AMR simulations. Mass loading factors of unity or above require extreme measures, such as turning off cooling for a prolonged time, or kinetic feedback that is in effect hydrodynamically decoupled due to the bubble radius exceeding the disc height (Fig. 5). The outflows produced by delayed cooling and kinetic feedback are very distinct, the former being cold, dense, and slow, while the latter are hot, diffuse, fast, and featureless (Figs. 6, 7, and 8). The other models produce slow and remarkably similar outflows at intermediate densities and temperatures.
Save for thermal dump feedback, all models do well in terms of resolution convergence when considering SFRs, while, with the exception of kinetic feedback, they produce significantly higher outflow rates at lower resolution (Fig. 10).
Although a direct comparison is difficult, stochastic feedback appears to produce much weaker outflows than in the similar disc runs with the original SPH version of the model of DS12. This discrepancy is perhaps a result of subtle setup differences between our discs and those of DS12, but we cannot rule out a more fundamental AMR versus SPH difference. Stochastic feedback does become efficient at generating massive outflows in our AMR discs if we use very high values for the stochastic heating temperature (up to K), but this comes at the cost of strong stochasticity due to low SN probabilities (Figs. 13 and 11).
The major handle on the generation of outflows appears to be how well the SN feedback model circumvents gas cooling, directly or indirectly. Delayed cooling is the only model which succeeds at generating outflows with mass loading factors exceeding unity (Fig. 5), at reproducing the observed main sequence SFRs (Figs. 2 and 3), and the Kennicutt-Schmidt relation (with appropriate calibration; Fig. 4). The other models fail to produce SN feedback strong enough to reproduce these observations. This is discouraging, as delayed cooling retains too much energy for too long, which in reality is partly lost to radiative cooling, while the other feedback models are arguably more physically motivated. Moreover, for the low-mass galaxy we argued that thermal dump, stochastic, and mechanical feedback converge because we resolve the adiabatic phase of the feedback events. This implies that in this case delayed cooling and kinetic feedback yield incorrect answers. In particular, delayed cooling results in gas occupying regions of temperature-density space where the cooling time is very short, which compromises predictions for observational diagnostics.
Possible reasons for the disconnect between observations and our results are: i) a lack of additional feedback physics, such as radiation feedback or cosmic rays, ii) an incomplete setup, i.e. an insufficiently realistic description of galaxies, iii) other aspects of the subgrid physics, such as star formation, are unrealistic (e.g. Agertz & Kravtsov, 2015; Semenov et al., 2016), iv) overcooling on galactic scales is still an issue at our resolution, even if different feedback models converge to the same results, and a significantly higher resolution is required.
The current analysis will serve as a foundation for future studies of feedback in galaxies, where we will use a sub-set of these models to study the interplay of SN feedback with different sub-grid methods for star formation and with feedback in the form of stellar radiation.
We thank Jérémy Blaizot, Léo Michel Dansac, Julien Devriendt, Sylvia Ploeckinger, and Maxime Trebitsch for useful discussions, and the anonymous referee for constructive comments. This work was funded by the European Research Council under the European Unionâs Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement 278594-GasAroundGalaxies. TK was supported by the ERC Advanced Grant 320596 “The Emergence of Structure during the Epoch of Reionization”. The simulations were in part performed using the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. We also acknowledge PRACE for awarding us access to the ARCHER resource (http://www.archer.ac.uk) based in the UK at the University of Edinburgh (PRACE-3IP project FP7 RI-312763).
- Agertz & Kravtsov (2015) Agertz O., Kravtsov A. V., 2015, AJ, 804, 18
- Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, AJ, 770, 25
- Bahé et al. (2016) Bahé Y. M. et al., 2016, MNRAS, 456, 1115
- Barai et al. (2015) Barai P., Monaco P., Murante G., Ragagnin A., Viel M., 2015, MNRAS, 447, 266
- Bland-Hawthorn et al. (2007) Bland-Hawthorn J., Veilleux S., Cecil G., 2007, Astrophys. Space Sci., 311, 87
- Blondin et al. (1998) Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, AJ, 500, 342
- Booth et al. (2013) Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 777, L16
- Ceverino & Klypin (2009) Ceverino D., Klypin A., 2009, AJ, 695, 292
- Chabrier (2003) Chabrier G., 2003, Publ. Astron. Soc. Pacific, 115, 763
- Chang et al. (2015) Chang Y.-Y., van der Wel A., da Cunha E., Rix H.-W., 2015, ApJS, 219, 8
- Crain et al. (2015) Crain R. A. et al., 2015, MNRAS, 450, 1937
- Crighton et al. (2013) Crighton N. H. M., Hennawi J. F., Prochaska J. X., 2013, ApJ, 776, L18
- Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
- Dubois et al. (2014) Dubois Y. et al., 2014, MNRAS, 444, 1453
- Dubois & Teyssier (2008) Dubois Y., Teyssier R., 2008, A&A, 477, 79
- Dubois et al. (2015) Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502
- Dwarkadas (2007) Dwarkadas V. V., 2007, AJ, 667, 226
- Evans et al. (2009) Evans N. J. et al., 2009, ApJS, 181, 321
- Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, AJ, 703, 1416
- Federrath & Klessen (2013) Federrath C., Klessen R. S., 2013, AJ, 763, 51
- Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, Publ. Astron. Soc. Pacific, 110, 761
- Fierlinger et al. (2016) Fierlinger K. M., Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2016, MNRAS, 456, 710
- Geen et al. (2015) Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2015, MNRAS, 454, 4484
- Gentry et al. (2016) Gentry E. S., Krumholz M. R., Dekel A., Madau P., 2016, eprint arXiv:1606.01242
- Gerritsen (1997) Gerritsen J. P. E., 1997, Ph.D. thesis
- Gerritsen & de Blok (1999) Gerritsen J. P. E., de Blok W. J. G., 1999, A&A, 342, 655
- Girichidis et al. (2016) Girichidis P. et al., 2016, ApJ, 816, L19
- Governato et al. (2010) Governato F. et al., 2010, Nature, 463, 203
- Guillet & Teyssier (2011) Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756
- Haardt & Madau (1996) Haardt F., Madau P., 1996, AJ, 461, 20
- Haehnelt (1995) Haehnelt M. G., 1995, MNRAS, 273, 249
- Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, ApJ, 777, L38
- Heckman et al. (2015) Heckman T. M., Alexandroff R. M., Borthakur S., Overzier R., Leitherer C., 2015, AJ, 809, 147
- Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581
- Hopkins et al. (2013) Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647
- Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488
- Iffrig & Hennebelle (2015) Iffrig O., Hennebelle P., 2015, A&A, 576, A95
- Katz (1992) Katz N., 1992, AJ, 391, 502
- Keller et al. (2014) Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013
- Kennicutt (1998) Kennicutt R. C., 1998, AJ, 498, 541
- Kim et al. (2014) Kim J.-h. et al., 2014, ApJS, 210, 14
- Kimm & Cen (2014) Kimm T., Cen R., 2014, AJ, 788, 121
- Kimm et al. (2015) Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900
- Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, AJ, 654, 304
- Macciò et al. (2008) Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
- Martizzi et al. (2015) Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504
- Mitra et al. (2015) Mitra S., Davé R., Finlator K., 2015, MNRAS, 452, 1184
- Mo & Mao (2004) Mo H. J., Mao S., 2004, MNRAS, 353, 829
- Mollitor et al. (2015) Mollitor P., Nezri E., Teyssier R., 2015, MNRAS, 447, 1353
- Murray et al. (2010) Murray N., Quataert E., Thompson T. A., 2010, AJ, 709, 191
- Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, AJ, 490, 493
- Navarro & White (1993) Navarro J. F., White S. D. M., 1993, MNRAS, 265, 271
- Puchwein & Springel (2013) Puchwein E., Springel V., 2013, MNRAS, 428, 2966
- Rasera & Teyssier (2006) Rasera Y., Teyssier R., 2006, A&A, 445, 1
- Rieder & Teyssier (2016) Rieder M., Teyssier R., 2016, MNRAS, 457, 1722
- Rogers & Pittard (2013) Rogers H., Pittard J. M., 2013, MNRAS, 431, 1337
- Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188
- Rosdahl et al. (2015) Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 4553
- Rosen & Bregman (1995) Rosen A., Bregman J. N., 1995, AJ, 440, 634
- Roškar et al. (2014) Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837
- Salem & Bryan (2014) Salem M., Bryan G. L., 2014, MNRAS, 437, 3312
- Sales et al. (2016) Sales L. V. et al., 2016, MNRAS, 464, 2419
- Schaye et al. (2015) Schaye J. et al., 2015, MNRAS, 446, 521
- Schroetter et al. (2015) Schroetter I., Bouché N., Peroux C., Murphy M. T., Contini T., Finley H., 2015, AJ, 804, 83
- Semenov et al. (2016) Semenov V. A. et al., 2016, AJ, 826, 200
- Simpson et al. (2016) Simpson C. M., Pakmor R., Marinacci F., Pfrommer C., Springel V., Glover S. C. O., Clark P. C., Smith R. J., 2016, AJ, 827, L29
- Somerville & Davé (2015) Somerville R. S., Davé R., 2015, Annu. Rev. Astron. Astrophys., 53, 51
- Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
- Springel & Hernquist (2003) Springel V., Hernquist L., 2003, Astrophys. Supercomput. using Part. Simulations, 208, 273
- Steidel et al. (2010) Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N. A., Bogosavljević M., Rudie G. C., Rakic O., 2010, AJ, 717, 289
- Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
- Stinson et al. (2013) Stinson G. S., Brook C., Macciò A. V., Wadsley J., Quinn T. R., Couchman H. M. P., 2013, MNRAS, 428, 129
- Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
- Teyssier et al. (2013) Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068
- Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, AJ, 630, 167
- Thornton et al. (1998) Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998, AJ, 500, 95
- Toro et al. (1994) Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25
- Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman J. H., Howell L. H., Greenough J. A., 1997, AJ, 489, L179
- Veilleux et al. (2005) Veilleux S., Cecil G., Bland-Hawthorn J., 2005, Annu. Rev. Astron. Astrophys., 43, 769
- Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
- Walch et al. (2015) Walch S. et al., 2015, MNRAS, 454, 238
- Wang et al. (2015) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83
- Zuckerman & Evans (1974) Zuckerman B., Evans N. J., 1974, AJ, 192, L149
Appendix A Stochastic feedback comparison with DS12
DS12 introduced a model for stochastic SN feedback in the SPH code Gadget, on which we base our AMR version of stochastic feedback. In simulations of two isolated disc galaxies and using the same fiducial stochastic heating temperature difference as us of K, they find mass loading factors of and , at of the virial radius, for galaxies of baryonic masses and , respectively. For a lower-mass galaxy roughly ten times lower in mass than theirs and a higher-mass galaxy roughly similar to theirs, we find mass loading factors at comparable distances that are orders of magnitude lower than their values.
While this may be due to the different hydrodynamical solvers, i.e. SPH and AMR, a direct comparison is complicated by the fact that in a number of ways, our simulations are set up differently from those in DS12.
One difference is that in the disc simulations of DS12, the CGM initially has zero gas density, while we are forced in AMR hydrodynamics to use positive nonzero density everywhere. We can rule this out as a major issue though: we ran simulations where we changed the initial CGM density by factors of ten in each direction, and found our simulation results remain unchanged.
The remaining trivial setup differences are the different pressure floors, cooling functions, star formation criteria, and time delay between the birth of stellar particles and their SN events. To assess the significance of those differences, we have performed a g10 run, and also of a galaxy similar to the low-mass galaxy in DS12, with our fiducial stochastic feedback, but with a setup closer to that in DS12, i.e. with zero metallicity, equilibrium cooling, , an (almost) identical pressure floor,
and a time delay of Myr from the birth and SN events of stellar particles (a factor of longer than our fiducial time delay).
Although this brings us closer to the simulation settings of DS12, we emphasize that it is still not an ideal comparison. For one thing, the pressure floor is applied slightly differently. In DS12 it is an actual floor, while in our simulations, it is added to the thermal pressure which is evolved separately. Also, in DS12 the pressure floor is only applied above the density threshold for star formation, of , while we apply it everywhere (though in the absence of metals, gas does not cool below the floor below the density threshold for star formation). Third, the star formation law in DS12 is pressure-dependent, and hence very different from ours. Fourth, it is not obvious whether the resolution is comparable between the AMR and SPH runs131313While the number and mass of DM particles and gravitational ‘softening’ scales are similar between our simulations and those of DS are comparable, it becomes more tricky for both the gas mass and physical resolution around SN events. In our AMR runs, the mass resolution of the gas elements receiving the SN energy is highly variable and depends linearly on the gas density, while in SPH the gas (and stellar) mass resolution is fixed ( and for the low-mass and high-mass galaxies in DS12, respectively). At the same time, the physical resolution is fixed in AMR (assuming the highest refinement level), but variable in SPH (depending, again, on the density). Still, we find no reason to expect that any of these discrepancies should lead to order of magnitude differences in the mass loading factors.
In Fig. 16, we show the SFRs and outflow rates at of the virial radius for our massive g10 galaxy and stochastic feedback, both for our fiducial simulation settings (green curves) and for these alternative settings described above to mimic those of DS12 (orange curves). With arrows at the right side of the plot we indicate the SFR and outflow rate at Myr in the corresponding run in DS12 (their g12 galaxy), also with stochastic feedback and K. By comparing the green and orange solid curves we see that applying the DS12 simulation settings has the effect of reducing the SFR by a factor and reducing the outflow rate by a factor of few. Clearly, this does not help in bringing the mass loading factor closer to the one found in DS12. The SFR in our run mimicking the DS12 simulation settings is actually close to that in DS12, and the suppression is likely mostly due to the strong pressure floor. However, even with a similar SFR as in DS12, the outflow rate is more than two orders of magnitude lower.
Fig. 17 shows the same kind of comparison, but now for a galaxy with a baryonic mass of , similar to the low-mass galaxy in DS12. We use the exact setup parameters as for the lowest-mass g8 galaxy described in Rosdahl et al. (2015), except that as with the other galaxies in this paper, the mass of formed stellar particles is times the cell volume at maximum resolution, and we now include a UV background. We performed one run with the fiducial settings described in §2, and another run with the DS12 settings described above for the g10 galaxy.
Comparing to the two arrows from DS12, our SFRs are a nearly an order of magnitude higher and not highly sensitive to the setup, except it is more bursty with our fiducial setup. The outflow rate, however is quite sensitive to the setup, being nearly an order of magnitude lower than that of DS12 for our fiducial setup, but nearly identical to that of DS12 when we mimic their setup. Hence, the mass loading factor is more than an order of magnitude lower in our fiducial runs and a factor of lower in our runs mimicking the DS12 settings. As with the g10 galaxy, though not as significantly, SN feedback is less efficient at regulating our galaxy and generating outflows than in DS12, hinting towards subtle ‘non-linear’ differences in simulation settings that remain undetected, or more fundamental differences between the hydrodynamical solvers (or both). We stress, however, that we must treat these hints with caution. While differences between SPH and AMR in terms of feedback efficiency and the generation of winds may exist, we cannot rule out other causes for now, and a further exploration must be left for future work.