LES of isolated disc galaxies

Large-eddy simulations of isolated disc galaxies with thermal and turbulent feedback

H. Braun, W. Schmidt, J. C. Niemeyer, and A. S. Almgren
Institut für Astrophysik, Universität Göttingen, Friedrich-Hund Platz 1, D-37077 Göttingen, Germany
Center for Computational Sciences and Engineering, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
E-mail: hbraun@astro.physik.uni-goettingen.deE-mail: schmidt@astro.physik.uni-goettingen.de
Accepted 2014 June 3. Received 2014 June 3; in original form 2013 December 18
Abstract

We present a subgrid-scale model for the Multi-phase Interstellar medium, Star formation, and Turbulence (MIST) and explore its behaviour in high-resolution large-eddy simulations of isolated disc galaxies. MIST follows the evolution of a clumpy cold and a diffuse warm component of the gas within a volume element which exchange mass and energy via various cooling, heating and mixing processes. The star formation rate is dynamically computed from the state of the gas in the cold phase. An important feature of MIST is the treatment of unresolved turbulence in the two phases and its interaction with star formation and feedback by supernovae. This makes MIST a particularly suitable model for the interstellar medium in galaxy simulations. We carried out a suite of simulations varying fundamental parameters of our feedback implementation. Several observational properties of galactic star formation are reproduced in our simulations, such as an average star formation efficiency 1 per cent, a typical velocity dispersion around in star-forming regions, and an almost linear relationship between the column densities of star formation and dense molecular gas.

keywords:
methods: numerical - galaxies: ISM - stars: formation - turbulence
pagerange: Large-eddy simulations of isolated disc galaxies with thermal and turbulent feedbackLarge-eddy simulations of isolated disc galaxies with thermal and turbulent feedbackpubyear: ????

1 Introduction

Stars are a product of a complex sequence of competing and interacting processes on a vast range of spatial and temporal scales that concentrate initially dilute gas into compact cores. It is not yet fully understood how the interplay of all of the processes involved, such as gravitational collapse, cooling, turbulence, magnetism, and stellar feedback, leads to the observed properties of the interstellar medium and stars in galaxies. A recent review on the properties of star formation was presented by Kennicutt & Evans (2012). Of particular interest is the mechanism regulating the observed low efficiency of star formation. Measures of the star formation efficiency are gas depletion - or consumption - time-scales , which relate star formation to the available gas supply. Wong & Blitz (2002); Evans (2008); Bigiel et al. (2008); Blanc et al. (2009) infer  Gyr in local disc galaxies from CO- -, and UV-measurements with resolutions down to 200 pc. Comparable measurements of gas-rich galaxies by Daddi et al. (2010); Tacconi et al. (2013); Saintonge et al. (2013) and others indicate a significantly shorter time  Gyr, corresponding to a relative gas consumption rate per free fall time . According to the KS relation (Schmidt, 1959; Kennicutt, 1998, and others), the star formation rate is well correlated with the local gas supply. The power-law slope of measured KS relations depends, however, on the tracers used, the resolution achieved, and other observational limitations (e.g. Onodera et al., 2010; Lada et al., 2010; Leroy et al., 2013). Recent observations show a good linear correlation between star formation rate and dense/molecular gas (Gao & Solomon, 2004; Lada et al., 2010; Bigiel et al., 2011). Evans et al. (2009); Murray (2011) showed that the local efficiency in individual actively star-forming molecular clouds is much greater than the efficiency on galactic scales, and their depletion time-scale is considerably shorter ( Myr). This implies that molecular clouds convert a sizable fraction 0.1-0.4 of their mass into stars during their lifetime (a few 10 Myr, see e.g. Blitz et al., 2007; McKee & Ostriker, 2007; Miura et al., 2012) before they are destroyed by supernova explosions (SNe) and stellar winds.

Actively star-forming molecular clouds are known to be strongly supersonically turbulent with typical velocity dispersions around (e.g. Leroy et al., 2008; Stilp et al., 2013), or larger in interacting galaxies (e.g. Herrera et al., 2011). Regulation by supersonic turbulence is a good candidate to theoretically explain the observed properties of star formation, as it globally supports a molecular cloud against gravity, but locally produces over-dense filaments and knots that may collapse into stars. A variety of approaches have been developed in the past years to derive star formation efficiencies from statistical properties of gravo-turbulent fragmentation inside molecular clouds (Padoan & Nordlund, 2011; Krumholz & McKee, 2005; Hennebelle & Chabrier, 2011; Padoan et al., 2012; Federrath & Klessen, 2013, hereafter FK13). As supersonic turbulence decays on relatively short time-scales of the order of the sound crossing time, it has to be maintained by some production mechanism over the lifetime of a molecular cloud. Processes such as large scale shear and instabilities in galactic discs (e.g. Gómez & Cox, 2002; Wada et al., 2002; Kim et al., 2003; Kim & Ostriker, 2007; Agertz et al., 2009; Krumholz & Burkert, 2010), accretion of gas on to a galaxy (e.g. Hopkins et al., 2013; Genel et al., 2012; Elmegreen & Burkert, 2010; Klessen & Hennebelle, 2010), and merger events or other galactic interactions (e.g. Bournaud et al., 2011; Teyssier et al., 2010) come into question here, but also local processes like stellar winds (e.g. Wolf-Chase et al., 2000; Vink et al., 2000; Vink, 2011), radiation pressure (Krumholz & Thompson, 2012) and SNe (e.g. Ostriker & Shetty, 2011; Agertz et al., 2009; Vollmer & Beckert, 2003), or the effects of thermal instabilities (e.g. Wada & Norman, 2001; Kritsuk & Norman, 2002; Iwasaki & Inutsuka, 2014) are possible turbulence production mechanisms.

In order to numerically simulate a realistic galaxy as a whole, star formation and the entailing stellar feedback have to be taken into account. However, the resolution to properly follow the evolution inside star-forming clouds in a galactic scale simulation is far from being feasible with contemporary computational resources. Recent simulations of isolated disc galaxies (IDG) feature resolutions down to a few parsec or even a fraction of a parsec (e.g. Hopkins et al., 2013; Renaud et al., 2013; Dobbs & Pringle, 2013; Benincasa et al., 2013; Booth et al., 2013; Monaco et al., 2012), while simulations of galaxies from cosmological initial conditions reach resolutions of some 10 parsec (e.g. Agertz et al., 2009; Munshi et al., 2013; Kraljic et al., 2012). To tackle sub-resolution processes, an appropriate subgrid-scale (SGS) model has to be applied. In the last decade a wide range of different models have been devised to effectively describe star formation and stellar feedback (e.g. Agertz et al., 2013; Stinson et al., 2006, 2013; Wise et al., 2012) using resolved quantities and assumptions about the small-scale properties of the ISM.

In galaxy simulations, the star formation rate is usually modelled using a constant efficiency per free fall time, which locally enforces a KS relation

 ˙ρs=ϵffρτff∝ρ1.5, (1)

where is the local gas density and the local free fall time. To avoid spurious star formation, additional constraints are applied, for example, a threshold for the minimal density required for star formation and a maximal temperature. More sophisticated models distinguish between different gas components. Gnedin et al. (2009) suggest to relate the star formation rate to the density of molecular gas instead of the total gas density. Murante et al. (2010) use a simple multi-phase approach to determine the fraction of the gas density that is available for star formation. For the simulations presented in this article, we use a multi-phase model for the ISM, including an estimation of the amount of shielded molecular gas (Braun & Schmidt, 2012, hereafter BS12). The star formation efficiency in the molecular gas is dynamically computed from the numerically unresolved turbulence energy, which is determined by a SGS model for compressible turbulence (see Schmidt & Federrath, 2011, hereafter SF11). Since we incorporate the coupling between resolved and unresolved scales as turbulent stresses in the Euler equations, our galaxy simulations are large-eddy simulations (LES). Moreover, we apply the energy-conserving AMR techniques presented in Schmidt et al. (2014). The diagonal part of the turbulent stresses acts as non-thermal pressure that usually dominates over the thermal pressure in cold and dense environments. This allows us to apply both thermal and turbulent feedback by channeling a fraction of the SN energy into the production of SGS turbulence energy. As we will show, this has important consequences for the regulation of star formation. In a way, this is similar to kinetic feedback (see, e.g., Agertz et al., 2013), with the important difference that we assume that turbulent motions are mainly excited on length scales below the grid resolution. For the thermal feedback, a small portion of the SNe energy is stored in a non-cooling budget - decaying on a time-scale of - to mimic the effect of hot SNe bubbles on sub-resolution scales, while the rest of the gas is allowed to cool radiatively. Although we include only effects of ionizing radiation from massive stars and SNe II, we are able to reproduce several observational features of star formation and turbulence in quiescent, gas-rich (or high redshift ) disc galaxies and star-forming regions.

This paper is structured as follows. First we describe the numerical methods, the SGS model, and the setup of our isolated disc galaxy (IDG) simulations in Sections 2 and 3. In Section 4, we present results from a suite of four simulations with different treatments of feedback, followed by our conclusions in Section 5

2 Numerical implementation

We carried out IDG simulations using the cosmological hydrodynamics code Nyx (Almgren et al., 2013). Nyx, built on the BoxLib software framework, uses Adaptive Mesh Refinement (AMR) to provide higher numerical resolution in sub-volumes of particular interest. Nyx solves the standard Euler equations using an unsplit Piecewise Parabolic Method (PPM); additional source terms are treated via a predictor/corrector scheme. Nyx is capable of following the evolution of different collisionless massive particles in the -Body formalism using a Kick–Drift–Kick algorithm, and (self-)gravity is taken into account using a Particle Mesh scheme with multigrid solver. We extended Nyx to run adaptively refined LES of an IDG as described in the following.

To handle sub-resolution processes, such as star formation, stellar feedback, cooling, and thermal instability, we use a model based on BS12 with a few minor modifications, MIST (Multi-phase Interstellar medium, Star formation and Turbulence model). The key-features of MIST are the following.

• Atomic, metal line, and dust cooling, photoelectric heating on dust.

• Separation of gas into two phases due to thermal instabilities that exchange energy and material via different mixing, heating, or cooling processes. The two phases represent a diffuse warm component and a clumpy, cold component of the gas. Balance of effective (i.e. thermal plus turbulent) pressure between the phases is assumed to obtain their respective densities.

• Formation of stars from the molecular fraction of the cold phase that is shielded from dissociating radiation. The star formation rate is computed dynamically from the thermal and turbulent state of the gas.

• Depending on the age of a stellar population the stellar feedback is applied. Lyman-continuum radiation and SNe II are taken into account. The SNe not only enrich the gas with metals, but also deposit kinetic energy and thermal energy into the gas. The SNe ejecta are treated as an additional sub-phase of the warm phase that does not cool efficiently and is gradually mixed with the rest of the warm phase.

• Unresolved turbulence is coupled to almost all processes implemented in MIST. Besides the source terms that are related to the scale separation for LES, small-scale pressure gradients caused by phase separation and SNe are taken into account.

An overview of important variables and coefficients is given in Table 1.

2.1 Gas dynamics

As an extension to the standard compressible Euler equations we introduce a new degree of freedom in the form of the SGS turbulence energy density (see SF11) and its source terms in order to model the behaviour of unresolved turbulent fluctuations. Furthermore we include the source terms as needed for the MIST model. The set of conservation equations for the evolution of gas becomes

 ∂ρ∂t+∇⋅(ρu)=−˙ρs,SF+˙ρs,FB, (2)
 ∂(ρu)∂t+∇⋅[ρuu+(p+23ρK)−τ∗]=ρg−u˙ρs,SF+us,FB˙ρs,FB, (3)
 (4)
 (5)

Here is the gas density, the velocity vector, the total specific energy of the gas, the specific turbulent SGS-energy, the specific internal energy, the thermal gas pressure (where is the polytropic equation of state parameter), and the gravitational acceleration vector (see Section 2.1.1). The SGS turbulence model related quantities , , and are defined in Section 2.1.2. For a definition of the star formation rate see Section 2.2.2. The stellar feedback rate , the mass-weighted average velocity of SN-ejecta, and the specific supernova energy deposit are explained in Section 2.3.2. We apply a multiphase model to determine the specific thermal energy of the cold phase, and the net cooling rate (see Section 2.2.1) and the turbulence energy production via phase separation and via SN feedback energy deposit (see Section 2.2.2). To conserve the kinetic energy of SN-ejecta, is added to (see equation 34 in Section 2.3.2).

2.1.1 Gravity

The massive components in our IDG simulations are dark matter, baryonic gas, and stars, where the dark matter component is assumed to be a static halo, and the other two are dynamically evolved. The gravitational acceleration vector is computed as the sum of a static acceleration due to the dark matter halo and the negative gradient of the gravitational potential due to the dynamical components:

 g=gdm−∇Φdyn, (6)

with the static acceleration (see equation 37 in Section 3.1). represents the solution of Poisson’s equation

 △Φdyn=4πG(ρdyn−¯ρdyn), (7)

where , is the mean of , is the total stellar density, and is the gravitational constant.

2.1.2 Hydrodynamical turbulence model

The interaction between resolved and unresolved turbulent velocity fluctuations is modelled using the SGS turbulence stress tensor , which can be seen as an analogue to the viscous dissipation tensor in the Navier-Stokes equations. With its components following SF11 read

 τij=2Cτ1Δ(2ρK)1/2S∗ij−2Cτ2ρKui,kuj,kul,mul,m−23(1−Cτ2)ρKδij, (8)

where

 S∗ij=Sij−13δijd=12(ui,j+uj,i)−13δijuk,k (9)

is the trace-free rate of strain, and the grid scale. The trace free stress tensor , used in equations (3) and (4), is given by . The SGS turbulence energy production rate , the SGS turbulence dissipation rate (which does not appear in equation 4 as it is absorbed into ), and the SGS turbulence diffusivity are given by

 ΠSGS = τijSij, (10) (ρεSGS) = ρCεK3/2Δ, (11) κSGS = Cκ(2K)1/2. (12)

We use the closure coefficients , , , and as determined by SF11 for compressible turbulence.

2.2 Non-adiabatic physics

In the following we describe how unresolved physics such as heating, cooling, star formation and stellar feedback was implemented in the code Nyx. For a more detailed description of the underlying model we refer to BS12. As input for the computation of the non-adiabatic physics sources we need the hydrodynamical state, the source terms belonging only to the SGS turbulence model and the stellar feedback terms. Contrary to BS12, stellar feedback is considered an external source in the calculation, as it depends on the stellar population represented by -body particles but not on the hydrodynamical state. Given the SGS- and stellar feedback source terms, the actual sources are calculated by subcycling the BS12 model ODEs locally in a grid cell, to resolve all time-scales of relevant processes, particularly the cooling time-scale, and then averaging the rate of change over the hydro-step. To follow the metal enrichment, we calculate three hydrogen density , helium density and metal density . Their conservation equations are of the form

 ∂ρX∂t+∇⋅(ρXu)=−ρXρ˙ρs,SF−∂ρs,X∂t∣∣∣FB, (13)

where indicates one of the species H, He, or Z, and is the ejection rate of that species by SNe (see Section 2.3.2).

2.2.1 Cold and warm gas phases

To keep track of the multiphase state in a grid cell, we introduce an additional passively advected quantity, the cold-phase fractional density , from which we can easily reconstruct the warm-phase density . The warm phase thermal energy is given by with a constant specific thermal energy of the cold phase, corresponding to a temperature . The conservation equation of reads

 ∂ρc∂t+∇⋅(ρcu)=Λc+ΛwfTIew−ec−˙ρs,SF−ASN˙ρs,FB, (14)

where

 Λc=−ρcεSGS−ΓPAHc−ρcρΓLyc and Λw=Λradw−ρwεSGS−ΓPAHw−ρwρΓLyc

are the net cooling rates of the cool and warm phase, respectively. is effectively a heating rate. Material is removed from the cold phase and transferred to the warm phase instead of increasing .

 ASN=13826⎛⎜⎝ρw,pamHcm3⎞⎟⎠45(ℓcpc)−65(ρcρc,pa)35 (15)

is the SN cold-phase evaporation coefficient, and is the efficiency parameter for turbulence production by the thermal instability, if the indicator (see BS12). Here is the photoelectric heating rate, the heating rate due to Lyman continuum radiation from young, massive stars (see Section 2.3.2), and

 Λradw=^ρwρw,pa^Λradw(ρw,pa,Z,^Tw) (16)

the radiative cooling rate, which is interpolated from a cooling table. These tabled cooling rates were computed using the photo-ionization program package Cloudy (Ferland et al., 1998, version 08.00). and are the average densities of the warm and the cold phase, respectively. Those are computed from the fractional phase densities ( and ) and the thermal and turbulent energies (, , and ) by assuming balance of effective (thermal plus turbulent) pressure between the phases at cold clump scale , as explained in detail in BS12. is the metallicity of the gas. and are the temperature and the fractional density of the warm-phase gas, that is allowed to cool radiatively. Note that and may differ from and in areas affected by recent SNe feedback. The treatment of the third gaseous phase, the hot SNe ejecta, which is prevented from cooling, is described in Section 2.3.2. The total net cooling rate, used in equation (4), is then given by .

2.2.2 Star formation rate

Stars are assumed to form from the molecular fraction of the gas in the cold phase, at a rate (Krumholz et al., 2009)

 ˙ρs,SF=fH2ρcϵPNtc,ff, (17)

where is the formation rate of gravitationally bound cores per free fall time

 tc,ff=√3π32Gρc,pa. (18)

To calculate , we use the Padoan & Nordlund (2011, hereafter PN11) model in the single free-fall formulation of FK13

 ϵPN=(1−floss)r12crit2⎛⎜⎝1+erf⎡⎢⎣σ2c−2log(rcrit)(8σ2c)12⎤⎥⎦⎞⎟⎠. (19)

Here is the fraction of mass in gravitationally bound cores lost during prestellar collapse through winds etc., and the standard deviation of the assumed density probability density function (PDF) of log-normal shape. The broadening parameter is set depending on which of the three production terms of turbulence energy in equation (5) is locally the dominant one

 Πmax=max[ΠSN,ΠTI,ΠSGS], (20)

where

 ΠTI=ΛwfTIϵTI (21)

is the contribution due to thermal instability, and

 ΠSN=˙ρs,FBeSNϵSN (22)

describes turbulence production by SNe. is the fraction of the energy released by SN that is deposited in the form of turbulent energy. We define by

 b=⎧⎨⎩1/3 if ΠSGS=Πmax2/3 if ΠTI=Πmax1 if ΠSN=Πmax. (23)

Here we assume the large-scale driving to be mostly caused by shear, the SNe driving to be mostly compressive, and the thermal instability driving to be of intermediate type. The corresponding values follow from Federrath et al. (2010).

To obtain the turbulent Mach-number of the cold phase, the SGS energy has to be rescaled from the grid scale to the cold clump scale , assuming a Kolmogorov velocity scaling exponent : 111We assume to be the largest scale represented in the cold phase. Consequently, the scaling of the turbulent velocities is applied to those scales, on which only the warm phase exists. Turbulence in the warm gas is usually subsonic or transonic at most. In this regime the assumption of a Kolmogorov-type scaling behaviour with coefficient seems valid.

 M2c=2K(ℓcΔ)2ηγ(γ−1)ec. (24)

The critical over-density ratio , above which a bound object is formed, is given by (FK13)

 rcrit=0.0067×5×2K(ℓcΔ)2ηπGρc,paℓ2cM2c. (25)

The molecular fraction of cold gas is computed from the cold and warm phase fractional densities and energies, assuming effective pressure equilibrium, using a Stömgren-like approach. The penetration depth of impinging radiation into a spherical cold clump of diameter is determined by the balance between production and dissociation due to UV-photons. The dissociating radiation field is assumed to be homogenous and isotropic. However, in dense, cold environments, which are identified by and , is dimmed by a factor

 fIν=13+23max[1+TIν−T(e,Z)TIν,0], (26)

with , because of the assumed shielding from radiation by the environment.

The assumption of a log-normal shaped density PDF, which is an essential part of the theory of PN11, applies to turbulence in isothermal gas. For a consistent definition of the internal energy of the cold phase in MIST, an adiabatic exponent is required. However, we assume a constant average temperature of the cold phase because of processes which are not explicitly treated. Both observational and numerical studies on the density PDFs show that the density PDF of the cold phase of the ISM is well approximated by a log-normal PDF (e.g. Hughes et al., 2013; Schneider et al., 2014). Although a power-law tail is generally found for actively star-forming clouds in which dense cores undergo gravitational collapse, FK13 point out that this does not significantly affect the star formation efficiency because the log-normal turbulent density fluctuations feed the collapsing gas that populates the power-law tail at high densities. Despite of the underlying inconsistency, all currently available analytic models for the calculation of star formation efficiencies, including PN11, are based on this assumption. Substituting PDFs with power-law tails into these models does not amend the problem because this would lead to divergent integrals. As a consequence, the construction of consistent models of the star formation efficiency is an open problem.

2.3 Implementation of stellar particles

2.3.1 Stellar particle creation

A particle is characterized by its position , mass , velocity , and an arbitrary number of additional properties. To handle the dynamical evolution of stars, we implemented a particle type with three additional properties: the initial mass , creation time , and metallicity , which are needed for the application of stellar feedback. A stellar particle does not represent a single star, but a single stellar population with a normalized initial mass function (IMF) (where is the number of stars of individual initial mass per solar mass of stellar population) scaled by . We use the IMF of Chabrier (2001).
To avoid the repeated creation of particles in all cells where , we introduce a stellar density field , that acts as an intermediate buffer for the stellar mass. This is treated as an passively advected quantity with respect to the hydro-solver, but massive with respect to gravity. Its conservation equation reads

 ∂ρs,m∂t+∇⋅(ρs,mu)=˙ρs,SF−∂ρs,m∂t∣∣∣SP, (27)

where represents the mass transfer from into stellar particles. The total stellar mass density is given by .
A pair of stellar particles and is created in the cell centre, if the agglomerated exceeds the threshold (corresponding to a minimum particle pair mass for a cell size of ), the mass is removed from . The properties of the new particles are

 m1,2p=ρs,mΔ3/2,u1,2p=u±urnd,m1,2pi=ρs,mΔ3/2,t1,2pc=t+dt/2,Z1,2p=ρZ/ρ. (28)

Here is the hydro time-step, and is a random velocity (in opposite directions for and to conserve total momentum). This random velocity component is intended to reflect the unresolved motions of the cold clumps, which the stars originate from. Its absolute value is drawn from a Gaussian distribution with expectation value and variance

 σ2urnd=furnd2K(1−(ℓcΔ)2η). (29)

The fudge factor is designed to obtain as random velocity component at the end of the particle growth phase, when the final is reached. is the mean mass of all stellar particles and the upper threshold mass for particle growth.
A newly created stellar particle collects the stellar mass in along its path [using a Nearest Grid Point algorithm (NGP)]. Its properties are updated using

 ´mp=mp+ρs,mΔ3,´up=upmp+uρs,mΔ3´mp,´mpi=mpi+ρs,mΔ3,´tpc=tpcmp+(t+dt/2)ρs,mΔ3´mp,´Zp=Zpmp+ρZρρs,mΔ3´mp,´xp=xpmp+xρs,mΔ3´mp. (30)

The final mass is reached, when either or .

2.3.2 Feedback mechanism and hot phase

We consider the two physical stellar feedback processes, Lyman continuum heating and SNe explosions, and use the equations in BS12 to compute their contributions to the source terms for the update of the hydrodynamical state.
The stellar mass sink field acts on the gas only via Lyman continuum heating , where we assume zero age of the stellar population it represents. The amount of feedback (the mass of SNe ejecta , the mass in the different chemical species , and the heating rate due to Lyman continuum radiation ) during a hydro time-step is computed individually for every stellar particle according to its properties. To obtain , the feedback is mapped to the hydro-mesh using a cloud in cell (CIC) scheme and time-averaged over the time step . The SNe ejecta mass is computed by

 (31)

where is the initial mass of an individual star that goes SN at an age of (Raiteri et al., 1996). is removed from the particle mass: . The metal load of the ejecta is given by , where is the fraction of metals produced in the dying stars, the other species scale linearly with metallicity. The total Lyc heating rate is . A fraction of energy load of the ejecta is deposited into (via ) and , and , respectively. However, this does not account for the kinetic energy of the ejecta that must also be transferred to the hydro-mesh along with their mass:

 ∂(ρsEkin)∂t∣∣∣FB=∑pu2s,p2∂ρs∂t∣∣∣FB,p. (32)

Here we sum over all local contributions to kinetic energy from every individual particles ’’ The momentum of the ejecta transferred to bulk momentum of the gas changes the bulk kinetic energy of the gas by

 ∂(ρEkin)∂t∣∣∣FB=u2(2us−u)∂ρs∂t∣∣∣FB. (33)

To conserve energy, we add the difference to

 ∂(ρK)∂t∣∣∣exFB=∂(ρEkin)∂t∣∣∣FB−∂(ρEkin)∂t∣∣∣FB. (34)

The heating of the gas to very high temperatures by SNe is described by a hot phase density , obeying the conservation equation

 ∂ρh∂t+∇⋅(ρhu)+min[ρh∇⋅u,0]=∂ρs∂t∣∣∣FB−ρhexp(−t/τh), (35)

where describes the loss of thermal energy due to adiabatic expansion222Employing a ceiling prevents producing hot phase when gas is compressed., and the decay of due to successive mixing of the hot gas into the ISM. The half-lifetime-scale is defined such that a SNe bubble shell at typical expansion velocity (roughly the speed of sound in the hot phase) travels roughly during that period, which leads to . is the constant specific thermal energy of the hot phase gas.
For consistency, the input parameters and to the derivation of the radiative cooling rate (equation 16) are computed as follows:

 ^ρw=ρw−ρh,^ew=ρwew−ρheh^ρw. (36)

3 Simulations

3.1 Initial conditions

In the simulation domain with a volume of we initialize an isothermal gaseous disc with an exponential surface density profile residing in a static dark matter halo using the potential-method of Wang et al. (2010), which gives initial conditions similar to Agertz et al. (2009). The choice of this setup is advantageous compared to a setup using a constant vertical scale height of the disc alike that by Tasker & Tan (2009), because it is adiabatically stable.

In the absence of a stellar component, the exponential gaseous disc is defined by its mass , its orientation of the disc angular momentum assuming a radial scale length , an initially uniform metallicity and a temperature of the disc.

The dark matter is modelled by a static halo with a NFW-shaped density profile (Navarro et al., 1997). It only contributes to the dynamics via its gravitational acceleration

 gdm=−GMdmrlog(1+cdm)−cdm1+cdm(log(rs)r3−cdmrdmr2rs), (37)

at a given position with distance vector from the halo centre, its absolute value and the scaled dimensionless radius . The NFW profile used is fully characterized by the halo mass , the virial radius , and a concentration parameter .

3.2 Individual runs

We performed eight isolated galaxy runs with model parameters listed in Table 2. All runs were carried out with a root grid of cells. AMR levels with a factor of 2 spatial and temporal refinement were created using refinement criteria based on density; specifically, any cells with density above the minimum value of were tagged for refinement up unto a specified maximum number of levels. In runs ref, nE, nB, nEnB, sSF, and sSF2 a spatial resolution of was obtained using six levels of refinement. The effects of the feedback implementation are explored with the runs nE, nB, and nEnB in comparison with run ref, that features MIST with the reference parameters as given in Table 1.
The runs sSF and sSF2 run feature a simplified model for the ISM without phase separation and . A threshold controlled star formation recipe is applied here, according to which stars are allowed to form at a free fall time efficiency of , if the local density exceeds and the local temperature is lower than . While in sSF2 the whole SGS-turbulence model is switched off, the model is still active in sSF, but star formation is decoupled from and the SGS-energy production terms related to MIST are switched off in this case (i.e. and ). Stellar particles are not created as pairs (see Section 2.3.1) in sSF and sSF2, since is not defined in both cases.
In addition two runs with fewer refinement levels were performed in order to investigate the effects of numerical resolution on the results of our simulations. The runs lres5 and lres4 feature effective resolutions of and using five and four levels of refinement, respectively.
Up to densities around the resolution requirement by Truelove et al. (1997) is easily satisfied in all runs with MIST. This limit can be shifted towards much higher densities, if we consider the effective pressure instead of thermal pressure only. The Jeans length of the few densest cells may temporarily drop below though, but never below the size of a cell , before the dense region is disrupted again by feedback. The latter statements are not true in case of sSF and sSF2, in which the Jeans length criterion is always violated in the dense clusters.
The combined usage of LES and MIST in a simulation increases the amount of computational resources required by less than 10 per cent compared to runs without, and significantly less than an additional level of refinement (more than 100 per cent).
The isolated galaxies were evolved for at least one orbital time at 10 kpc radius from the centre ( Myr).

4 Results

4.1 Disk evolution in the ref run

Initially the disc is adiabatically stable, but as the gas is allowed to cool by radiation, it loses its thermal support in height and collapses into a thin cold disc. The disc becomes Toomre-unstable and fragments into clumps. In those clumps the gas eventually becomes dense enough to become molecular and consequently begins to form stars. The SNe feedback of the newly formed stars then eventually disperses the clumps. A fraction of the stars formed in some of those clumps may form a stellar cluster that survives over a much longer time than the lifetime of about 20 Myr of the gas clump from which it originated. The stellar clusters tend to move towards the centre of the disc as a result of dynamical friction, where they eventually merge into a central agglomeration of stars, if they are not disrupted before by the tidal forces in the disc. However, they do not build up a bulge as their velocity dispersion is too small. The majority of the stars form a rather smooth disc with a scale height of a few hundred parsec. The structure of the stellar disc is shown in the top left panel of Fig. 1. The stars stripped in the potential of the disc form tidal tails around their birth cluster.

The SNe also carve holes into the disc and launch a wind leaving the disc. The interplay between cooling, gravity, and SNe shapes the gaseous component into a fluffy disc, with holes and knots, as demonstrated in the top and bottom panels in the middle of Fig. 1. The wind mostly not only consists of metal enriched hot gas from SNe, but also carries a fraction of the original cold clump with it (see bottom left panel of Fig. 1). The wind originates from the cavities of the disc caused by the SNe at speeds ranging from 300 to , and pushes a shell of cold or warm gas outward. The ejected gas is either mixed into the wind, or falls back into the disc. Far away from its origin, the wind from all sources merges into a hot, but dilute, sub-sonically turbulent, and metal-rich () bubble that continues to expand. The bottom panels of Fig. 1 give an impression of the metallicity, mass, and velocity structures in the winds above the disc.

As seen in the top right panel of Fig. 1, active star formation occurs only in a few compact regions away from the centre. Because of the local metal enrichment due to SNe ejecta, the threshold density for star formation drops. This causes clumps to form stars earlier during their collapse, when the density is still relatively low, which lowers their star formation rate and makes them more prone to dispersal by SNe. The residual stellar clusters are fewer, lighter, and more easily disrupted in the galaxy’s potential.
The global star formation rate in the ref run is plotted as black line in Fig. 2. Initially the amount of star formation increases quickly, as the region of star formation grows. It reaches its peak around 300 Myr after start of simulation, and then gradually declines due the consumption and the metal enrichment of the gas reservoir in the inner disc. At this stage approximately 30 per cent of the initial gas mass has been converted into stars.

4.1.1 Differences in the nE run

In the nE run we set . Thus SNe feedback does not directly increase the unresolved turbulent energy (see equation 5). This has basically two effects on the overall evolution. On the one hand, star formation in a clump is active over a longer period of time because of the higher star formation efficiency in moderately turbulent clumps (), causing more stars to form before a clump of gas is dispersed. This leads to a slightly higher global star formation rate than in the ref run (see the blue line in Fig. 2). On the other hand, stars forming in a clump after the first stars have already produced SNe have a significantly reduced velocity dispersion compared to their analogues from the ref run. As a consequence, those stars are more likely to stay near the clump of their origin. In this case, stellar clusters tend to be more massive and more strongly gravitationally bound, and hence, their tidal tails are less prominent. Clusters are initially more abundant and have a longer lifetime. Eventually they merge into a few very massive clusters, which accrete gas, and host intermittent star formation, as the gas is driven apart due to feedback. They remain gravitationally bound and stable, as their stellar mass is sufficiently large and concentrated.
As a consequence, the resulting stellar and gaseous discs are more clumpy than in the ref run. This can be seen by comparing the projections of total gaseous density and stellar density in Fig. 3, where both stellar and gaseous clumps are fewer in numbers of appearance but more massive with higher central densities. During the simulation approximately 40 per cent of the initial gas mass was consumed by star formation after 1 Gyr.

4.1.2 Differences in the nB run

In the nB run we turned off the treatment of the hot SN ejecta phase, i.e. the hot SN gas is directly mixed into the warm phase. This allows the gas to radiate away the feedback energy more quickly. Thus star formation in a clump goes on for a longer time, consuming a larger fraction of the gaseous mass. Once almost all gas is depleted the feedback takes the lead, and drives stronger and faster winds than in the ref run. Like in the nE run, the resulting stellar clusters are more stable and massive, and subsequently merge to form very massive clusters, that cannot be disrupted in the tidal field of the disc. Those massive clusters are hosting continuous star formation. They keep accreting gas from the disc and converting it into stars, thereby sustaining their gaseous mass at the same level for a long time.
As in the nE run, the resulting gaseous and stellar discs are clumpy, but the stellar clusters are much more massive, and the gas is depleted on much shorter time-scales (see Fig. 3). The stellar density in the centre of the clusters reaches values much greater than . This is why we stopped this simulation after 0.5 Gyr. During this period of time about 30 per cent of gaseous mass was turned into stars.

4.1.3 Differences in the nEnB run

Setting and turning off the treatment of the hot SN ejecta phase combines the effects described above. This leads to an even more violent evolution with a few very massive stellar clusters in the inner regions, while the gas is depleted quickly (see Fig. 3). We stopped this run even earlier than the nB run. The amount of gas consumed by star formation was about 30 per cent in roughly 400 Myr.

4.2 Star formation

4.2.1 Time scales

The 20 Myr time-scale related to single star-forming regions reflects the period of time needed to produce enough stars, such that their combined feedback is able to quench star formation, and to evaporate the most dense, central region of the star-forming cloud. If the environment is still dense enough, star formation may continue in an compressed layer around the expanding bubble. Depending on the mass of the clump and its surroundings, star formation can continue in this mode. A time series depicting the evolution of one of those clumps is shown in Fig. 4. The impact of the most massive clumps is seen in the quick variations of the global star formation rate on time-scales around 10 to 40 Myr (see Fig. 5). Metal-enriched clouds tend to have a shorter lifetime compared to metal-poor clouds, as star formation can start at lower densities, because shielding from radiation is enhanced due to higher dust abundances.

By assuming different threshold densities to compute the reservoir of available gas, we derive global gas depletion time-scales

 τdep,i(t)=∑ρ(t)>ρthr,iρΔ3∑˙ρSFΔ3=Mgas,i˙MSF. (38)

To filter out quick variations, we apply a moving average

 ¯¯¯τdep,i(t)=1tav∫t+tav/2t−tav/2τdep,i(t)dt (39)

with . Plots of the resulting depletion times are shown in Fig. 5 for different choices of . The typical time-scales range from less than 100 Myr for the highest threshold of 1 to a few Gyr for the total gas content of the simulation domain. Gas above densities is most likely actively star-forming. Contrary to the other cases, for slightly shrinks with time. Because of the increasing metallicity the threshold density needed for star formation drops gradually. Star formation and subsequent feedback then prevent gas from becoming as dense as in metal-poorer environments. for remains almost constant at 0.25 Gyr after the initial transient phase. For lower , grows with time, and the growth rate tends to increase with decreasing . The growth of if is also related to the increase of metallicity. In metal-rich environments gas forms stars already at lower densities, implying greater , but the efficiency (see equation 17) remains about the same. For this reason is effectively lowered faster than gas supply shrinks. We expect that all will saturate, as the impact of increasing metallicity is further reduced in already enriched gas (see BS12). However, the magnitude of the depletion time-scales is well within range of observationally inferred ones (Daddi et al., 2010; Genzel et al., 2010).

4.2.2 Local efficiency

The star formation efficiency is defined by

 ˙ρSF=ϵffρτff, (40)

where a constant value for is commonly assumed. In our simulations is computed dynamically from the local turbulent hydrodynamical state. While the star formation efficiency in the cold phase, (see equation (17)), is almost constant – around – for all star-forming regions, the shielded molecular content and vary significantly, such that is boosted in high density regions. As a consequence, a large contribution to the global star formation rate comes from just a few temporarily very active spots in the disc, which explains the relatively large variations on short time-scales. As shown in Fig.6, we find a good correlation , which holds for all times in the ref run and the nE run. In the other runs without treatment of a hot phase, the power-law slope flattens towards large , which is clearly a sign that the effectiveness of thermal feedback plays a major role. The typical value of is, on average, reproduced in runs with hot phase treatment, while those without produce stars at significantly higher average efficiency.

4.2.3 H2 and HI

We find a very tight correlation between the star formation column density and the molecular column density , as shown on the left panel of Fig. 7. This correlation implies a robust power-law relation between and with an exponent slightly above one. This corresponds to an almost constant depletion time of the molecular gas in star-forming regions around 80 Myr. This correlation holds for all runs independent of the simulation time, as shown in the middle and right panels of Fig. 7.

The existence of the correlation itself is not surprising, as it is assumed in the MIST model equations that is proportional to (see equation 17). , however, is not imposed. Assuming a constant star formation efficiency , one would expect , but in our model and depend implicitly and non-linearly on , , and , which are in turn nonlinear functions of , , , , and . BS12 already showed that a correlation with can be found in the equilibrium solutions of the MIST model equations. This suggests that the interplay of the implicit dependencies effectively results in an almost linear relation between estimated -mass available to star formation and the star formation rate, which is an important feature of the self-regulating mechanisms implemented in MIST. The depletion time-scale of molecular gas from the equilibrium solutions (BS12) is about 1-2 Gyr, in agreement with that found by e.g. Bigiel et al. (2011) from CO-observations. We derive an average depletion time of  Myr in our simulations, or shorter ( Myr) in regions of higher molecular density. This is considerably, by a factor of around 25, shorter (see Fig. 7). But this time-scale matches the depletion time of molecular gas within star-forming molecular clouds, found by Murray (2011)333Murray (2011) find a GMC star-forming efficiency and related free fall times . A combination of both yields a depletion time-scale of ., and that for dense molecular gas inferred from HCN-observations by Gao & Solomon (2004)444Gao & Solomon (2004) find a linear correlation between the galactic star formation rate and the amount of what they call ’dense gas’ in a galaxy: . This corresponds to a dense gas depletion time of ., and the total gas depletion time-scale for . This discrepancy between the depletion time-scales of molecular gas is caused by the fact that we estimate the amount of gas involved in the processes of active star formation by calculating the molecular fraction of gas in the centres of cold cloud complexes in equilibrium with an external radiation field. However, there are also environments containing significant amounts of that are not in equilibrium with radiation or only partially shielded. Molecular material that is transported from shielded to not shielded regions by turbulent motions or SNe blast waves without being instantaneously dissociated is not tracked in our model for example. This material, as well as gas that is partially molecular but not dense enough to form stars, contributes to observed molecular hydrogen column densities, but is not present in our simulations. Our data on the -content of the ISM do reflect the amount of dense molecular gas that is actively forming stars, but do not reproduce the full amount of molecular gas seen in CO-observations. Coarse graining our -data changes neither the inferred depletion time-scales nor the slope significantly, which is consistent with Gao & Solomon (2004), who used galactic quantities.

Without following the chemical evolution of molecular and atomic gas, we cannot, unfortunately, distinguish between atomic gas and gas that would observationally be considered molecular. So we focus on the combined content. The general shape of over shown in the left panel of Fig. 8 is similar to observational relations(e.g. in Schruba et al., 2010), but our distribution is shifted towards higher densities and rates.
As demonstrated by Schruba et al. (2010), the distribution in - space changes from low to high between the vertical ’barrier’ behaviour in - space and the linear relation observed for . This transition happens around the critical density for atomic to molecular conversion, which is metallicity dependent. Increasing the density of gas only leads to an increase in the molecular density, if its density was already above the critical density. Because of that the atomic density saturates. This explains the ’vertical barrier’-shaped distribution in - space, if we take into account that the amount of surrounding atomic gas does not affect star formation in the molecular gas.
From the high-density tail in the left panel of Fig. 8 one would rather infer a power-law slope than a linear relation with , but also a power-law slope seems possible. At this point we cannot distinguish whether the reason for this is either that stellar feedback prevents the high density tail from being populated up to densities at which could be observed, or the super-linear relation also seen in the equilibrium solutions of BS12 is recovered. With increasing metallicity of the star-forming gas the critical density of the atomic to molecular transition drops gradually. This effect is demonstrated by the shift towards lower with increasing time in the middle panel of Fig. 8. A linear fit to the tail of the distribution in ref after 1 Gyr yields a gas depletion time , roughly consistent with (see Section 4.2.1). In the ref run roughly the same amount of stars and metals were produced after 1 Gyr as in the nEnB after 0.4 Gyr. Nevertheless the relations are very different. While the high-density tail indicating the transition from to becomes shallower and more prominent in the course of the disc evolution, the shape of over in the nEnB run at 0.4 Gyr is about the same as in the early stages of ref. This is a consequence of inefficient mixing of the hot, metal-rich material from SNe with cold dense, but still metal-poor material that possibly could be turned into stars. The shift between the relations for different runs in the right panel of Fig. 8 is mainly caused by gas consumption, as the produced metals have not been mixed into the star-forming material yet. So, the shift between our results and observational findings is, on the one hand, a consequence of the lower metallicity and higher gas contents, compared to observed local galaxies.

On the other hand the resolution scale , at which the relation between star formation and gas density is evaluated, influences the critical density of the atomic to molecular transition too, as demonstrated in Fig. 9. Once the averaging scale is considerably larger than a typical star-forming region, i.e. above  120 pc in the ref run, the transition range in space seems to be shifted towards lower densities. Observational data like in Bigiel et al. (2011) usually correspond to averages over much larger areas, as observations used for statistical analysis of star formation in disc galaxies have spatial resolutions about some , compared to our . The star formation rate of a given active region is put into relation with a larger volume that contains besides the dense star-forming gas large amounts of ambient, inactive gas. Since the tracers have a certain lifetime, observationally measured star formation rates are temporal averages as well, while our data follow the instantaneous star formation rate. However, the shift of the position of the knee in the space with coarsening is a combination of resolution effects and the already discussed effects of stellar to gaseous mass ratio and the metallicity in the dense gas of the discs. Within a resolution element the fraction of gas that is not directly involved in the star formation process increases with decreasing resolution due to averaging. The high density tail in - space indicating a correlation of star formation with gas density becomes therefore less pronounced in case of coarser resolution.