Star formation and molecular hydrogen in dwarf galaxies: a non-equilibrium view

Star formation and molecular hydrogen in dwarf galaxies: a non-equilibrium view

Chia-Yu Hu, Thorsten Naab, Stefanie Walch, Simon C. O. Glover, Paul C. Clark
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild Strasse 1, D-85740 Garching, Germany
Physikalisches Institut der Universität zu Köln, Zülpicher Strasse 77, D-50937 Köln, Germany
Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
School of Physics and Astronomy, Cardiff University, 5 The Parade, Cardiff CF24 3AA, Wales, UK

We study the connection of star formation to atomic (HI) and molecular hydrogen (H) in isolated, low metallicity dwarf galaxies with high-resolution ( = 4 M, = 100) SPH simulations. The model includes self-gravity, non-equilibrium cooling, shielding from a uniform and constant interstellar radiation field, the chemistry of H formation, H-independent star formation, supernova feedback and metal enrichment. We find that the H mass fraction is sensitive to the adopted dust-to-gas ratio and the strength of the interstellar radiation field, while the star formation rate is not. Star formation is regulated by stellar feedback, keeping the gas out of thermal equilibrium for densities 1 cm. Because of the long chemical timescales, the H mass remains out of chemical equilibrium throughout the simulation. Star formation is well-correlated with cold ( T 100 K ) gas, but this dense and cold gas - the reservoir for star formation - is dominated by HI, not H. In addition, a significant fraction of H resides in a diffuse, warm phase, which is not star-forming. The ISM is dominated by warm gas (100 K T K) both in mass and in volume. The scale height of the gaseous disc increases with radius while the cold gas is always confined to a thin layer in the mid-plane. The cold gas fraction is regulated by feedback at small radii and by the assumed radiation field at large radii. The decreasing cold gas fractions result in a rapid increase in depletion time (up to 100 Gyrs) for total gas surface densities 10 Mpc, in agreement with observations of dwarf galaxies in the Kennicutt-Schmidt plane.

methods: numerical, galaxies: ISM, galaxies: evolution

1 Introduction

Dwarf galaxies are thought to be the building blocks of larger galaxies in the hierarchical picture of galaxy formation. While they contribute little to the total mass budget of the galaxy population, they are the most numerous type of galaxy in the local Universe. Dwarf galaxies tend to contain fewer heavy elements compared to more massive ones (Tremonti et al., 2004; Hunter et al., 2012). As such, they are ideal laboratories for studying physical processes of the interstellar medium (ISM) under chemically simple conditions, which can be quite different from those found in normal spiral galaxies such as the Milky Way.

On kpc-scales, spatially resolved observations of nearby star-forming spiral galaxies have demonstrated that the surface density of the star formation rate has a clear correlation with the molecular hydrogen (H) surface density and little correlation with the atomic hydrogen (HI) surface density (Bigiel et al., 2008, 2011). In light of the observational evidence, H-dependent sub-resolution recipes for star formation have been widely adopted in hydrodynamical simulations as well as semi-analytic models, where the H abundances are calculated either by directly incorporating non-equilibrium chemistry models (e.g. Gnedin, Tassis & Kravtsov, 2009; Christensen et al., 2012) or by analytical ansatz assuming chemical equilibrium (e.g. Fu et al., 2010; Kuhlen et al., 2012; Thompson et al., 2014; Hopkins et al., 2014; Somerville, Popping & Trager, 2015). The implicit assumption is that star formation only takes place in H-dominated clouds, and the HI-to-H transition has been thought to be responsible for the so-called star formation threshold, the surface density of gas below which star formation becomes extremely inefficient.

However, recent theoretical studies have cast doubt on the causal connection between H and star formation (Krumholz, Leroy & McKee, 2011; Glover & Clark, 2012b), based on the insensitivity of the gas thermal properties to the H abundances. Radiative cooling from H is negligible when the gas temperature drops below a few hundred Kelvin. The formation of carbon monoxide (CO), which is indeed an efficient coolant at very low temperatures, does require the presence of H. Nevertheless, singly ionized carbon (C) provides cooling that is almost equally efficient as CO at all but the lowest temperatures through fine structure line emission. It is therefore possible for gas devoid of both H and CO to cool down to low enough temperatures and form stars by gravitational collapse. In this picture, the correlation between H and star formation originates from the fact that both H formation and star formation take place in regions well-shielded from the interstellar radiation field (ISRF) instead of the former being a necessary condition for the latter. Such a correlation is expected to break down in low-metallicity environments where the H formation timescales are much longer. Indeed, Krumholz (2012) predicts that star formation would proceed in HI-dominated clouds if the gas metallicity is below a few percent of solar metallicity (Z) based on timescale arguments. Glover & Clark (2012c) also examine this issue, using simulations of isolated clouds, and show that the star formation rate of the clouds is insensitive to their molecular content and that for metallicities and below, the star-forming clouds are dominated by atomic gas.

Observationally, the detection of H in star-forming dwarf galaxies has proven to be very challenging. CO emission, as the standard approach to derive the H abundance, tends to be very faint in these galaxies. Schruba et al. (2012) observed 16 nearby star-forming dwarf galaxies and only detected CO successfully in five galaxies with oxygen abundance 12 + log (O/H) 8.0 with very low CO luminosities per star formation rate compared to those found in massive spiral galaxies. CO was not detected in the other 11 galaxies with 12 + log (O/H) 8.0 even if stacking techniques were used. The interpretation was that the CO-to-H conversion factor is significantly higher in low metallicity environment, under the assumption that there should be much more H given their star formation rate (i.e., assuming constant H depletion time 2 Gyr). The dependence of the CO-to-H conversion factor on metallicity is already observed in Local Group galaxies using dust modeling to estimate the H mass (Leroy et al., 2011) and has considerable theoretical support (Wolfire, Hollenbach & McKee, 2010; Glover & Mac Low, 2011; Shetty et al., 2011a, b; Bolatto, Wolfire & Leroy, 2013). However, an alternative interpretation is that H is also rare and star formation in these systems is taking place in regions dominated by atomic hydrogen (see Michałowski et al., 2015 who adopted such an interpretation for their observed galaxies).

In this paper we conduct numerical simulations to study the ISM properties in an isolated star-forming dwarf galaxy. Our model includes gravity, hydrodynamics, non-equilibrium cooling and chemistry, shielding from the ISRF, an H-independent star formation recipe, stellar feedback and metal enrichment in a self-consistent manner. We investigate the relationship between H and star formation and explore the effects of varying the strength of ISRF and the dust-to-gas mass ratio (DGR). In Section 2 we describe the details of our numerical method. In Section 3 we show the ISM properties when it is in thermal and chemical equilibrium. In Section 4 we present the results of our numerical simulations. In Section 5 we discuss the implications of our results and the potential caveats. In Section 6 we summarize our work.

2 Numerical Method

2.1 Gravity and Hydrodynamics

We use the gadget-3 code (Springel, 2005) where the collisionless dynamics of gravity is solved by a tree-based method (Barnes & Hut, 1986), while the hydrodynamics is solved by the smoothed particle hydrodynamics (SPH) method (Lucy, 1977; Gingold & Monaghan, 1977). We have implemented a modified version of SPH, called SPHGal, in gadget-3 which shows a significantly improved numerical accuracy in several aspects (Hu et al., 2014). More specifically, we adopt the pressure-energy formulation of SPH111 We use the pressure-energy SPH instead of the pressure-entropy SPH as used in (Hu et al., 2014) for reasons described in Appendix A. (Read, Hayfield & Agertz, 2010; Saitoh & Makino, 2013) which is able to properly follow fluid instabilities and mixing without developing severe numerical artifacts commonly found in traditional SPH (e.g. Agertz et al., 2007). The so-called ‘grad-h’ correction term is included following Hopkins (2013) to ensure the conservation properties when the smoothing length varies significantly (e.g. at strong shocks). The smoothing length is set such that there are = 100 particles within a smoothing kernel. We use a Wendland kernel function that has been shown to be stable against pairing instability for large (Dehnen & Aly, 2012), which is necessary for reducing the ‘E error’ (Read, Hayfield & Agertz, 2010) and therefore improving the numerical convergence rate (Read & Hayfield, 2012; Dehnen & Aly, 2012). We adopt artificial viscosity to properly model shocks with an efficient switch that only operates at converging flows, similar to the prescriptions presented in Morris & Monaghan (1997) and Cullen & Dehnen (2010). We also include artificial thermal conduction (Price, 2008) but only in converging flows to smooth the thermal energy discontinuities, which can lead to severe noise at strong shocks when the pressure-energy formulation is used. The viscosity and conduction coefficients are varied in the range of [0.1,1] and [0,1] respectively. Our SPH scheme shows satisfactory behaviors and accuracies in various idealized numerical tests presented in Hu et al. (2014).

2.2 Chemistry Model

Our model of the chemistry and cooling follows closely the implementation in the SILCC project (Walch et al., 2015; Girichidis et al., 2015), based on earlier work by Nelson & Langer (1997), Glover & Mac Low (2007) and Glover & Clark (2012a). We track six chemical species: H, H, CO, H, C, O and free electrons. Only the first three species are followed explicitly, i.e., their abundances are directly integrated based on the rate equations in our chemistry network. The fractional abundance of neutral hydrogen is given by


where denotes the fractional abundance of species ; note that all fractional abundances quoted here are relative to the number of H nuclei. Silicon is assumed to be present in singly ionized form (i.e. as Si) throughout the simulation, while carbon and oxygen may either be present as C and O, or in the form of CO, which leads to


where and are the abundances of the total carbon and oxygen respectively. Finally, the abundance of free electron is given by


A list of our chemical reactions for H and H is summarized in Table 1 of Micic et al. (2012). The H is formed via collisional ionization of hydrogen with free electrons and cosmic rays, and is depleted through electron recombination in both the gas phase and on the surfaces of dust grains. The H is formed on the surfaces of dust grains and is destroyed via ISRF photo-dissociation, cosmic ray ionization, and collisional dissociation (with H, H and free electrons).

Our chemical model also includes a treatment of carbon chemistry, following the scheme introduced in Nelson & Langer (1997). It assumes that the rate-limiting step of CO formation is the process . The CH may either react with atomic oxygen and form CO or be destroyed via ISRF photo-dissociation. However, we will not go into detailed investigation about CO in this work, as a proper modeling for CO formation is beyond our resolution limit. This is especially so in low metallicity environments where CO only resides in regions of very high density (Glover & Clark, 2012c).

2.3 Cooling/Heating Processes

We include a set of important non-equilibrium cooling and heating processes. The term ’non-equilibrium’ refers to the fact that the processes depend not only on the local density and temperature of the gas but also on its chemical abundances of species, which may not be in chemical equilibrium. Cooling processes include fine-structure lines of C, O and Si, the rotational and vibrational lines of H and CO, the hydrogen Lyman-alpha line, the collisional dissociation of H, the collisional ionization of H, and the recombination of H in both the gas phase and on the surfaces of dust grains. Heating processes include photo-electric effects from dust grains and polycyclic aromatic hydrocarbons (PAHs), ionization by cosmic rays, the photo-dissociation of H, the UV pumping of H and the formation of H.

We do not follow the non-equilibrium cooling and heating processes in high temperature regimes. For T K we adopt a cooling function presented in Wiersma, Schaye & Smith (2009), which assumes that the ISM is optically thin and is in ionization equilibrium with a cosmic UV background from Haardt & Madau (2001). The total cooling rate depends on the temperature and density of the gas as well as the abundance of heavy elements. We trace eleven individual elements (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) for both gas and star particles based on the implementation in Aumer et al. (2013).

2.4 Shielding of the Interstellar Radiation Field

Shielding from the ISRF affects both the chemistry and the cooling/heating processes. For the hydrogen chemistry, the H ISRF photo-dissociation rate, , is attenuated by both the dust shielding and the H self-shielding:


where s is the unattenuated photo-dissociation rate, is the strength of the ISRF relative to the Solar neighborhood value estimated by Habing (1968). The total attenuation factor is where and are the attenuation factors by dust extinction and by H self-shielding, respectively. We adopt


where is the DGR relative to the Milky Way value (1 %) and cm is the averaged cross section of dust extinction. The H self-shielding is related to the H column density using the relation given by Draine & Bertoldi (1996). A similar treatment is used for the carbon chemistry. The CO photo-dissociation rate is attenuated by the dust extinction, H shielding and CO self-shielding. A more detailed description can be found in Walch et al. (2015).

Dust extinction also reduces the photo-electric heating rate by blocking the radiation in the energy range between 6 and 13.6 eV. We adopt the heating rate given by Bakes & Tielens (1994) and Bergin et al. (2004) of:


where is the attenuated radiation strength, is the number density of the gas and is the heating efficiency defined as


where and is the electron number density.

To calculate the column densities relevant for shielding we have incorporated the TreeCol algorithm (Clark, Glover & Klessen, 2012) into our version of Gadget-3. Unlike the extragalactic UV background that is external to the simulated galaxy, the sources of the ISRF are the young stars embedded in the ISM. This means that the column densities should not be integrated over the entire disc, but have to be truncated at certain length scales. Ideally, one should integrate the column densities up to individual stars that contribute to the local radiation field separately, but this would entail performing a computationally expensive radiative transfer calculation on every timestep, which is impractical. Instead, we make the simplifying assumption that the material beyond a predefined shielding length from a local gas particle is not relevant for the shielding (see e.g. Dobbs et al., 2008 or Smith et al., 2014, who make similar approximations in their galactic simulations). We justify this assumption by noting that for gas particles in dense clouds – the only ones which are significantly shielded – the dominant contribution to the shielding generally comes from gas and dust in the cloud itself or in its immediate vicinity, rather than from the diffuse gas between the clouds. Therefore, provided we adopt a value for that is larger than the typical size of the dense clouds, our results should be insensitive to the precise value chosen for 222 Along a few lines of sight which intercept other dense clouds, taking a small causes us to potentially underestimate the shielding. However, the local radiation field will always be dominated by other lines of sight with low column density. .

For each gas particle, TreeCol defines equal-area pixels using the healpix algorithm (Gorski & Hivon 2011) and computes the column density within each pixel out to . We set and  pc throughout this work. However, we explore in Appendix C.2 the effect of varying within the range 20–200 pc and show that as expected our results are insensitive to the precise value of within this range.

2.5 Star Formation Model

Unlike cloud-scale simulations, our mass resolution is not sufficient to follow the gravitational collapse of the gas to densities where it will inevitably end up in stars. Instead, we define an instantaneous star formation rate for each gas particle to estimate how much gas is expected to be converted into stars: , where is the local free-fall time, is the star formation efficiency, and is the volumetric density of gas. We set = 0.02 to account for the inefficient star formation which might originate from the sub-resolution turbulent motions (Krumholz & McKee, 2005). We assume the gas is ‘star-forming’ () only if the gas has , and negative velocity divergence. We choose = 100 cm as this is the typical densities of the giant molecular clouds in our Galaxy, the reservoir gas for star formation for which = 0.02 is defined. We also set = 100 K to ensure that we do not attempt to form stars in hot dense gas which has a high Jeans mass. In practice, most gas with would actually be colder than 100 K (cf. Fig. 9). Our definition of ‘star-forming gas’ is very simplistic. In reality, stars may still form out of a gas cloud with if the sub-resolution density structure is very clumpy. Note that the appropriate choice of and also depends on resolution. With higher resolutions, one should be able to follow the gravitational collapse to smaller scales and denser environments (e.g. individual molecular cores), and therefore both and should be set higher.

We adopt the commonly used stochastic approach of star formation: in each timestep of size t, a star-forming gas particle has a probability of to be converted into a star particle of the same mass. Since almost always holds, the conversion test in each timestep is a rare event Bernoulli trial and the number of ’success’ events during a time period of for a given follows a Poisson distribution with the parameter = 1 (which is the expectation value), i.e., one gas particle would be converted into a star particle in a time period of on average. Note that the ratio of the standard deviation to the expectation value for a Poisson distribution is , so the actual star formation timescale can deviate from the input by 100% if we look at a single particle. Only when we look at the averaged star formation rate of a group of particles would the random fluctuation be reduced to a satisfying level.

The instantaneous star formation rate of gas particles, however, is not an observable. In fact, it is merely an estimate for the star formation that will happen in the next timesteps. Therefore, we measure the star formation rate in a given region by the total mass of newly formed star particles with age less than in that region divided by , which is set to be 5 Myr in this work. Such a definition is more compatible with what is measured in observations than the instantaneous star formation rate assigned to the gas particles based on the adopted star formation model. These two definitions of star formation rate give almost identical results if we sum over a large enough region (e.g. the total star formation rate of the galaxy), although locally (both in space and in time) they can be quite different (cf. Section 4.7.1). Throughout this work we will adopt this definition when we present the star formation rate (except for Fig. 16 where both definitions are shown).

2.6 Stellar Feedback and Metal Enrichment

2.6.1 Supernova type II (SNII)

We assume each star particle represents a stellar population of mass and calculate the corresponding mass that will end up in SNII. For our adopted Kroupa initial mass function (Kroupa, 2001), we have . When the age of a star particle reaches 3 Myr, we return of mass to the ISM with enriched metal abundances according to the metallicity dependent yields given by Woosley & Weaver (1995). The returned mass is added to the nearest 100 neighboring gas particles weighted by the smoothing kernel. The energy budget for a stellar population is where erg is the typical energy for a single SN event and is the number of SN events.

Physically, a supernova remnant (SNR) should first go through a non-radiative phase where momentum is generated (the Sedov phase) until radiative cooling kicks in and the total momentum gradually saturates (e.g. Ostriker & McKee, 1988; Blondin et al., 1998). However, numerically, the resolution requirement for modeling the correct evolution of a SNR is very demanding (see e.g. Walch & Naab, 2015). It has long been recognized that insufficient resolution leads to numerical over-cooling: most of the injected energy is radiated away before it can significantly affect the ISM. As shown in Dalla Vecchia & Schaye (2012), this occurs when the cooling time is much shorter than the response time of the gas for a given resolution. For usual implementations in SPH simulations where , this means that the returned mass is always much smaller than and the SNR would be poorly resolved. The situation would not be alleviated by reducing as long as is assumed. Dalla Vecchia & Schaye (2012) circumvent this issue by injecting the energy to fewer neighboring particles and, if necessary, stochastic injection which groups several SN events into a single energetic one. By doing so they guarantee that the gas would always be heated to the local minimum of the cooling function and has more time to develop the blast wave.

Although grouping several SN events into one energetic explosion enhances the dynamical impact of SN feedback on the ISM, it also coarsens the granularity (both spatial and temporary) of the SN events. As shown in Kim & Ostriker (2015), a single energetic explosion over-produces both the momentum and energy compared to a series of spatially coherent explosions with the same amount of total energy. The difference would probably be even more severe if the explosions occur at different locations. Therefore, for a given energy budget, coarser SN sampling can over-estimate the impact of SN feedback. With the SN sampling as a free parameter, the effect of feedback becomes tunable or even arbitrary. Indeed, in large-scale cosmological simulations with necessarily compromised resolutions, the feedback has to be calibrated by fitting to observations (Schaye et al., 2015) and thus can only be regarded as a phenomenological model. However, if one’s goal is to directly resolve individual blast waves without tunable parameters, as we try to do in this work, then the SN sampling should be set to the physical one by making sure that each SN event has the canonical energy of erg.

As the star particle mass reaches M, the energy budget for a star particle becomes smaller than (i.e. ), which is also unphysical. In this work we adopt an ansatz similar to the stochastic injection in Dalla Vecchia & Schaye (2012). We let star particles with the appropriate age explode as SNII with a probability of /100 M and with an energy of . In this work, = 4 M and therefore each star particle has a 4% chance of producing a type II SN. Our resolution is close to the = 1 M requirement for a reasonably converged energy and momentum evolution of an SNR as shown in the resolution study in Appendix B.

Note that despite the fact that our particle mass (4 M) is comparable to the mass of a single star, the star particles should still be considered as stellar populations instead of individual stars. A star particle, once formed, is only a representative point of the collisionless distribution function that describes the stellar system of the galaxy. With this interpretation, there is no conceptual issue even when the particle mass reaches sub-solar scales, though the system might be over-resolved. The collisional dynamics of star clusters is suppressed by our gravitational softening (2 pc), and thus could only be included in a sub-resolution model, which is not considered in this work.

2.6.2 Supernova type Ia (SNIa) and asymptotic-giant-branch (AGB) stars

We include feedback by SNIa and AGB stars based on the implementation presented in Aumer et al. (2013). For SNIa, we adopt a delay time distribution (DTD), the SN rate as a function of time for a given stellar population formed in a single burst. The DTD has a power-law shape where is the stellar age, with the normalization of 2 SNIa events per 1000 M (Maoz & Mannucci, 2012). The amount of mass returned to the ISM is calculated by sampling the DTD with a 50 Myr time bin, and the metal yields based on Iwamoto et al. (1999). Similarly, the mass returned by the AGB stars is calculated from the metal yields presented in Karakas (2010) with the same time bin as SNIa. Assuming an outflow velocity = 3000 km/s and 10 km/s for SNIa and AGB stars, respectively, we return energy of 0.5 into the ISM where is the returned mass in a given time bin, though in our simulations their effect is sub-dominant compared to the SNII feedback.

2.7 Timestep Limiter

We include a timestep limiter similar to Saitoh & Makino (2009) and Durier & Dalla Vecchia (2012) to correctly model the strong shocks. For each active particle , we identify any neighboring particles within its smoothing kernel whose timesteps are four times longer than the timestep of particle and force them to be active at the next synchronization point. In addition, we re-calculate the hydrodynamical properties for the particles which we inject the feedback energy into and make them active such that their timesteps will be assigned correctly at the next synchronization point. Note that by modifying a particle’s timestep before it completes its full kick-drift-kick procedure we necessarily violate energy conservation, but to a much lesser extent than a numerical scheme without the timestep limiter.

2.8 Numerical Resolution

The mass resolution of an SPH simulation depends not only on the particle mass () but also the number of particles within a smoothing kernel (). For a given kernel function and , using more particles in a kernel means smoothing over more mass and hence worse resolution. However, because of the low-order nature of SPH, a relatively large is required to reduce the so-called ’E-error’, a zeroth order numerical noise induced by particle disorder (see e.g., Read, Hayfield & Agertz, 2010; Dehnen & Aly, 2012). We adopt = 100 as a compromise between suffering too much from the E-error ( too small) and over-smoothing ( too large). It seems reasonable to regard the kernel mass as the resolution of SPH simulations. However, different kernel functions entail different extent of smoothing and different scales of compact support (). A Gaussian kernel is an extreme example with infinite and while its smoothing scale is obviously finite. The same (and hence ) can therefore mean very different resolutions depending on the adopted kernel. A more physical meaningful way is required to define a length scale which reflects the true extent of smoothing. Dehnen & Aly (2012) proposed to define the smoothing scale as , where is the second moment of the kernel function (or ’standard deviation’). Following such a definition, for the commonly used cubic spline kernel and for our adopted Wendland kernel. Therefore, the mass resolution would be . For in our simulations (see Section 4.1), this means that is the mass scale for which we define the local density of a particle. Everything below is blurred by smoothing.

2.8.1 Jeans mass criterion

In hydrodynamical simulations that include self-gravity (as we do in this work), one important scale to resolve is the Jeans mass such that


where is the Jeans mass and is a prescribed number to ensure is well-resolved (Bate & Burkert, 1997; Robertson & Kravtsov, 2008). When the Jeans mass is not resolved (), perturbations on all scales down to the resolution limit () should collapse physically. However, perturbations can be created by numerical noise, triggering gravitational collapse, which tends to occur when the perturbations are only marginally resolved. Eq. 8 makes sure that perturbations near the resolution limit are physically stable and all perturbations that collapse are well-resolved and are of physical origin rather than numerical noise. The choice of is somewhat arbitrary, as it is difficult to judge whether a collapse is physical or numerical. Commonly suggested values in the literature are in the range of = 4 15 (e.g. Robertson & Kravtsov, 2008; Hopkins, Quataert & Murray, 2011; Richings & Schaye, 2015), which seems to be quite stringent considering the resolution is about 0.1 . This may be related to the smaller commonly used which is more prone to noise. As such, we take = 1, which means that we require one kernel mass to resolve the Jeans mass. According to the gas distribution in the phase diagram (see Fig. 9), our maximum number density satisfying is about 200 cm, which corresponds to the smoothing length 2 pc. This motivates us to choose our gravitational softening length to be also 2 pc.

A commonly adopted approach to ensure that the Jeans mass is formally resolved throughout the simulations is to introduce a ’pressure floor’ which artificially stabilizes the gas (by adding pressure by hand) when it becomes Jeans-unresolved (violating Eq. 8) (e.g. Robertson & Kravtsov, 2008; Hopkins, Quataert & Murray, 2011; Renaud et al., 2013; Vogelsberger et al., 2014; Schaye et al., 2015). However, physically, a Jeans-unresolved perturbation should collapse. In fact, the only way to properly follow the evolution of the Jeans-unresolved gas is to increase the resolution so that it becomes Jeans-resolved. When a gas cloud becomes Jeans-unresolved, the credibility of gravitational effect is lost irrespective of whether a pressure floor (or any other similar approach that makes the equation of state artificially stiffer) is used or not. Besides, the pressure floor makes the gas artificially stiff and thus also sabotages the accuracy of hydrodynamical effects. We therefore refrain from using a pressure floor and simply try to keep the majority of gas Jeans-resolved by using high enough resolution.

3 The ISM in Equilibrium

Before delving into the non-equilibrium evolution of the ISM and the complicated interplay of all the physical processes, in this section we investigate the ISM properties in both thermal and chemical equilibrium under typical conditions in dwarf galaxies. To see how a system evolves towards chemical and thermal equilibrium in a uniform and static medium, we run the code with only the cooling and chemistry modules while turning the gravity and SPH solvers off. The density and column density are directly specified as input parameters for each particle rather than calculated with SPHGal or TreeCol. The initial temperature is set to be K. We set the metallicity Z = 0.1 Z and the cosmic ray ionization rate s.

The parameters that directly affect H formation and destruction are and (Eq. 10 and 11). The strength of the ISRF in the diffuse ISM is expected to correlate with the star formation rate density. As will be shown in Section 4, the total star formation rate of our simulated galaxy is M yr and the radius of the star-forming region is 2 kpc. Thus, the star formation rate surface density is M yr kpc, which is about ten times smaller than the value in the solar neighborhood. Assuming the correlation between star formation rate surface density and is linear (e.g. Ostriker, McKee & Leroy, 2010), = 0.17 would be a plausible choice. On the other hand, the emergent radiation from the less dusty star-forming regions in low-metallicity environments seems to be stronger than in metal-rich ones due to the higher escape fraction of UV photons (see e.g. Bolatto et al., 2011; Cormier et al., 2015). Therefore, the resulting in dwarf galaxies might be as strong as that in the solar neighborhood. The observed DGR of a galaxy has been shown to scale linearly with its metallicity (though with significant scatter). This implies = 0.1 for our adopted metallicity Z = 0.1 Z. However, below Z Z, the DGR starts to fall below the linear DGR-Z relation (Rémy-Ruyer et al., 2014).

In this work we explore three different combinations of and . We will use the following naming convention:

  • G1D01: = 1.7, = 0.1

  • G1D001: = 1.7, = 0.01

  • G01D01: = 0.17, = 0.1

3.1 Chemical Equilibrium

3.1.1 The H formation/destruction timescale

We adopt the H formation rate:


where is the rate coefficient (Hollenbach & McKee, 1979), and , and are the number densities of molecular hydrogen, atomic hydrogen, and hydrogen nuclei, respectively. is the DGR normalized to the Milky Way value such that = 1 means the DGR = 1%. At low temperatures ( 100 K), the rate coefficient cm s is relatively insensitive to temperature variations, and the solution to Eq. 9 is with the formation timescale:


where .

Meanwhile, the dominant process of H destruction is photo-dissociation by the ISRF in the Lyman and Werner bands (11.2 13.6 eV). From Eq. 4, the destruction timescale is


Therefore, H clouds exposed to an unattenuated radiation field will be quickly destroyed, independent of the gas density. However, the photo-dissociation itself is also an absorption process, which attenuates the UV radiation efficiently and thus lengthens as the H column density accumulates (H self-shielding).

To illustrate the formation and destruction of H for different ISM parameters, we set up hydrogen column densities from cm to cm sampled by 32768 points equally spaced in and evolve the chemistry network including cooling and heating. This corresponds to a one-dimensional absorbing slab in a homogeneous medium. The H column density is obtained by direct integration:


We define the H mass fraction and therefore = 1 means that the hydrogen is fully molecular.

In the upper three rows in Fig. 1 we show the time evolution of in a uniform medium of = 100 cm in G1D01, G1D001 and G01D01 from top to bottom. Initially, the medium is totally atomic ( = 0) in the left panels and totally molecular ( = 1) in the right panels. The formation times are = 100 Myr, 1000 Myr and 100 Myr for G1D01, G1D001 and G01D01 respectively. If the medium is initially atomic, in a few , gas at sufficiently high becomes totally molecular and the system reaches chemical equilibrium, which is consistent with Eq. 10. If the initial medium is molecular, the gas would be photo-dissociated starting from the surface (low ). In the unattenuated region, Eq. 11 suggests rapid destruction on time scales of kyr, much shorter than . However, only the thin surface directly exposed to the UV radiation would be destroyed very promptly. Due to H self-shielding, the dissociating front propagates slowly into high and the system reaches chemical equilibrium on time scales much longer than 1 kyr. On the other hand, the free-fall time at = 100 cm is 5 Myr. Therefore, in the ISM of dwarf galaxies, the timescale to reach equilibrium, either staring from an atomic or molecular medium, is quite long compared to the local free-fall time.

3.1.2 The H fraction in chemical equilibrium

An unperturbed cloud will eventually reach a steady state where the H formation rate equals its destruction rate:


Together with and (assuming the ionization fraction of hydrogen is zero), the ratio of H to HI at chemical equilibrium can be written as


In unshielded regions, tends to be very low for typical conditions. Only in well-shielded environments where very few H-dissociating photons are present can grow to a significant value. The H profile is dictated by the value of , which can be viewed as a dimensionless parameter representing the capability of H formation relative to its destruction. In the bottom row of Fig. 1 we show the equilibrium H profile ( vs. ) with different and . For a given value of , increasing (decreasing) by a factor of 10 has the same effect as decreasing (increasing) by the same factor. The equilibrium H profile is thus a self-similar solution depending on the value of (Sternberg, 1988; McKee & Krumholz, 2010; Sternberg et al., 2014). Note that although the equilibrium profile is self-similar, the time required for reaching the equilibrium ( or ) is not.

Figure 1: Upper three rows: time evolution of in a uniform medium of = 100 cm in G1D01 (1st row), G1D001 (2nd row) and G01D01 (3rd row). The medium is initially totally atomic ( = 0) in the left panels and totally molecular ( = 1) in the right panels. Bottom row: the equilibrium H profile ( vs. ) with different (left panel) and (right panel). For a given , the equilibrium H profile is self-similar depending on the dimensionless parameter .

3.2 Thermal Equilibrium

We set up the hydrogen number density from to sampled by 32768 points equally spaced in and let the system cool from K until it reaches thermal equilibrium. Here we assume the column density to be = where is the SPH smoothing length. The H column density is obtained by the approximation = . This is the minimum column density that would be assigned to a gas particle for the given density in our simulations if there were no other contributing material integrated by TreeCol. The smoothing length is calculated assuming the same mass resolution ( = 4 M, = 100) as in our simulations in Section 4.

In the upper three panels in Fig. 2 we show the individual cooling and heating processes in thermal equilibrium in G1D01, G1D001 and G01D01. The 4th panel shows the equilibrium temperature vs. density, and the 5th panel shows the equilibrium pressure vs. density. For almost the entire range of density, the photo-electric effect is the main heating mechanism. In G1D001, the photo-electric effect becomes less efficient (due to the low DGR), and thus cosmic ray ionization dominates at cm. The dust-gas collision rate rises with but only becomes an important cooling mechanism at the highest densities which rarely occur in our simulations. The most dominant coolant in the range of - cm is the C fine structure line emission. Below cm, as the temperature increases, the OI fine structure emission and the hydrogen Lyman-alpha emission start to dominate over C. Above cm the gas becomes optically thick to the UV radiation and so the photo-electric heating rate decreases. Meanwhile, the C cooling is taken over by the CO cooling as most carbon is in the form of CO in this regime. The H cooling is unimportant in all cases as it is typically only abundant at densities where the temperature is too low to excite H. In G1D001 and G01D01, the photo-electric heating is less efficient than in G1D01 due to lower DGR and weaker ISRF respectively (cf. Eq. 6). Therefore, the equilibrium temperatures are the highest in G1D01 (bottom panel of Fig. 2). Note that in low metallicity and low density regimes, the time required to reach thermal equilibrium can become quite long (Krumholz, 2012; Glover & Clark, 2014). Therefore, the gas is expected to be out of thermal equilibrium at low densities (cf. Fig. 9).

The maximum density with an equilibrium temperature about K is = 0.25, 0.03 and 0.1 cm in G1D01, G1D001 and G01D01, respectively. As will be shown in Fig. 13, determines the maximum radius for star formation as it marks the onset of thermal-gravitational instability (Schaye, 2004).

Figure 2: The cooling and heating rates for gas in thermal equilibrium in G1D01 (1st panel), G1D001 (2nd panel) and G01D01 (3rd panel). The 4th panel shows the equilibrium temperature vs. density, and the 5th panel shows the equilibrium pressure vs. density. The maximum density with an equilibrium temperature about K is = 0.25, 0.03 and 0.1 cm in G1D01, G1D001 and G01D01, respectively.

4 Simulations

4.1 Initial Conditions

We set up the initial conditions using the method developed in Springel (2005). The dark matter halo has a virial radius = 44 kpc and a virial mass = , and follows a Hernquist profile with an NFW-equivalent (Navarro, Frenk & White, 1997) concentration parameter = 10. The spin parameter for the dark matter halo is = 0.03. A disc comprised of both gas and stellar components is embedded in the dark matter halo. The total mass of the disc is with a gas fraction of 66 %. The small baryonic mass fraction (0.3 %) is motivated by the results of abundance matching (Moster et al., 2010; Moster, Naab & White, 2013). The disc follows an exponential profile with scale-length of 0.73 kpc. The scale-length is determined by assuming the disc is rotationally supported and the total angular momentum of the disc is 0.3 % of that of the dark matter halo, which leads to a one-to-one relation between and the scale-length. Since dwarf galaxies are usually expected to have a relatively thick disc (e.g. Elmegreen & Hunter, 2015), we set the scale-height of the disc to be 0.35 kpc. The initial metallicity (gas and stars) is 0.1 Z uniformly throughout the disc with the relative abundances of the various metals the same as in solar metallicity gas. The cosmic ray ionization rate is s. The particle mass is for the dark matter and for the baryons (both stars and gas). The gravitational softening length is 62 pc for dark matter and 2 pc for baryons. As discussed in Section 3, we explore three different combinations of and as shown in Table 1. Run G1D01_noFB has the same and as G1D01 while the stellar feedback is switched off.

We use cylindrical coordinates and to describe the simulations, where is the galactocentric radius and is the rotation axis of the disc. The origin is chosen at the center of mass of the stellar disc, so that = 0 is the mid-plane of the disc. We also use the spherical coordinate to indicate the distance from the origin to a certain radius.

Name feedback
G1D01 1.7 0.1 yes
G1D001 1.7 0.01 yes
G01D01 0.17 0.1 yes
G1D01_noFB 1.7 0.1 no
Table 1: Simulation runs and the corresponding setup

4.2 Morphology

In Fig. 3 we show, from left to right, the column density maps of HI (1st column), H (2nd column) and H (3rd column) and temperature maps (slices, 4th column) at = 500 Myr, where is the simulation time. The upper two rows are the face-on and edge-on views for run G1D01, while the lower two rows are for run G1D001. Fig. 4 shows the same maps as Fig. 3 but for run G01D01 (upper two rows) and run G1D01_noFB (lower two rows). The temperature maps are slices across the origin so that the face-on slices are at the mid-plane . We use different ranges of the color bars for different chemical species for display purposes.

The ISM is strongly dominated by HI in all cases. It has clumpy density structures with many SN-driven bubbles. Most H resides in a thin layer in the mid-plane with little H at large . The H approximately traces dense gas. The H traces the shells of the SN-driven bubbles and is not confined to the mid-plane. The warm (green) gas occupies most of the volume in almost all cases. The hot (red) gas also fills up a non-negligible amount of volume (mainly in the SN bubbles) while the cold (blue) gas occupies very little volume.

Comparing the runs with feedback, H is most abundant in G01D01 while almost non-existent in G1D001. This is expected since the strong ISRF and low DGR in G1D001 is a hostile environment for H formation. If we turn off feedback, H becomes more abundant than in G01D01 (cf. Fig. 6, although it occupies little volume because the gas collapses into massive clumps.

From the face-on images, the region where the density structures can be found is the largest in G1D001 because the radius beyond which gas can not cool and collapse is the largest. This will be shown more quantitatively in Fig. 13.

Figure 3: Images of G1D01 and G1D001 at = 500 Myr. From left to right: the column density maps of HI (1st column), H (2nd column) and H (3rd column) and temperature maps (slices) (4th column). The top two rows are the face-on and edge-on views for run G1D01 while the bottom two rows are for run G1D001. Note that we use different ranges of the color bars for different chemical species for display purpose, as the system is strongly dominated by HI.
Figure 4: Same as Fig. 3 but for run G01D01 (upper two rows) and run G1D01_noFB (lower two rows).

In Fig. 5 we show the face-on column density for the G1D001 run at = 500 Myr with zoom-in’s at different scales. The central large panel shows the disc within R 2 kpc and the other four panels are zoom-in views of the central one. The top left panel shows a filamentary structure about 300 pc long. The bottom left panel is an example of a SN bubble with a size of about 200 pc. The top right panel shows a region with plenty of dense clouds and the bottom right panel is a further zoom-in from the top right. Note that the spatial resolution we have in dense gas is about 2 pc and therefore the clumps in the bottom right panel are expected to be still well-resolved. The ISM is highly inhomogeneous with a complex density structure.

Figure 5: The face-on column density maps of HI for the G1D001 run at = 500 Myr at different scales. The color scale is the same as in Fig. 3 and 4. Central panel: the entire star-forming region (box size = 4 kpc). Top left: a filamentary structure that is about 300 pc long. Bottom left: a 200-pc scale bubble driven by SN feedback. Top right: a group of dense clouds. Bottom right: further zoom-in of the top right. The effective spatial resolution is about 2 pc (see Section 2.8). The ISM is highly inhomogeneous with complex density structures.

4.3 Time Evolution of Global Properties

4.3.1 Star formation rate and mass fraction

In the upper panel of Fig. 6 we show the time evolution of the total star formation rate of the galaxy (SFR) for D1G01, D1G001, D01G01 and D1G01_noFB. In the first 50 Myr the gas collapses onto a thin disc as it cools down and starts forming stars. This phase of global collapse is quickly terminated by SN feedback and the system gradually settles to a steady state, with SFR M yr. The gas depletion time of the galaxy /SFR 40 Gyr is so long that the gas reservoir is able to sustain the SFR throughout the simulation time. This is true even if galactic outflows are considered, as we will show in Section 4.3.3. The evolution of the SFRs is very similar in all runs with feedback, suggesting that the thermal properties of gas are not very sensitive to the strength of the ISRF or the DGR. The only run that shows significant difference is the no-feedback one G1D1_noFB. Here the SFR first gets as high as 0.07 M yr at = 200 Myr, more than an order of magnitude higher than all the other runs, and quickly declines afterwards. This is due to a relatively short gas depletion time 500 Myr. The gas reservoir is rapidly consumed by star formation.

In the lower panel of Fig. 6 we show the H mass fraction of the galaxy (solid lines), defined as the total H mass divided by the total gas mass in the ISM (, where and are the mass and the H mass fraction of individual particle ). The ISM is defined as gas within the region 2 kpc (which is roughly the edge of star formation activities, see Fig. 13) and 1 kpc (to exclude the halo region). The dashed lines show the corresponding (chemical-)equilibrium H mass fraction , where is calculated by Eq. 13 (i.e. assuming all particles are in chemical equilibrium). Note that the self-shielding factor is still obtained from the non-equilibrium from the simulations. Unlike the SFR, differs significantly between different runs. In G1D01 is about 0.05%. Lowering by a factor of ten (G1D001) decreases by more than an order of magnitude, as the H formation rate is directly proportional to . With ten times smaller (G01D01) the difference is slightly smaller: increases by about a factor of three. The no-feedback run (G1D1_noFB) shows the highest (up to 4%) as the global gravitational collapse forms massive dense clumps which provides an ideal environment for forming H. After = 200 Myr, declines because most of the massive clumps have turned into stars.

Comparing and in each run indicates that the equilibrium prediction systematically overestimates the H mass fraction. This difference is expected due to the slow H formation rate: = Gyr. In the ISM constantly stirred by the SN-driven turbulence, it is unlikely for a gas cloud to stay unperturbed for a few and form H. On the other hand, the destruction rate can also be slow in well-shielded regions, as shown in Section 3.1.1, which also leads to some over-abundant H gas that should have been destroyed if it were in chemical equilibrium. The over-abundant and under-abundant H gas compensate each other and therefore the global H fraction is only slightly lower than , although locally they are far out of chemical equilibrium (as we will show in Section 4.4.2).

In the no-feedback run (G1D01_noFB) and agree pretty well. This is partly due to the relatively ’quiescent’ ISM (due to the lack of SN feedback) allowing for H formation in a less disturbed environment. More importantly, the high-density clumps result in much shorter , driving the system towards chemical equilibrium much faster.

Figure 6: Time evolution of global quantities of the galaxy. Upper panel: total star formation rate (SFR). The SFR is relatively insensitive to the choice of and but differs strongly if stellar feedback is switched off. Lower panel: H mass fraction in the ISM. The ISM is defined as gas within the region 2 kpc and 1 kpc. Solid lines () are from the simulations while dashed lines () are calculated assuming chemical equilibrium. The H mass fraction is very sensitive to variations of and . The equilibrium prediction systematically over-estimates the H mass fraction.

4.3.2 Mass and volume fraction of different phases

In Fig. 7 we show the time evolution of the mass fraction (, solid lines) and volume fraction (, dashed lines) of the ISM ( 2 kpc and 1 kpc) in hot ( K), warm (100 K K) and cold ( K) phases. The volume fraction is obtained by integrating the volume (estimated by the smoothing kernel size, )-weighted temperature histogram for different phases. The warm gas dominates the ISM both in mass () and also in volume (). As visualized in Fig. 3 and 4, the hot phase occupies a non-negligible volume fraction (), though it contributes only in mass. On the other hand, the volume occupied by the cold gas is negligibly small, though its mass fraction can be as high as . The cold-gas fraction is a bit higher in G1D001 and G01D01 due to less efficient photo-electric heating in the low-DGR and weak-ISRF conditions (Eq. 6). In principle, a lower DGR also results in less dust-shielding, which would potentially result in less cold gas. However, dust-shielding only operates at , which is rare in our simulations except for the densest gas. There is also a slight increase in the cold-gas fraction over time as a result of metal enrichment, which will be shown in Fig. 8. In the no-feedback run G1D01_noFB, there is almost no hot gas as the hot gas is mostly produced by SN blastwaves, and the cold-gas fraction decreases over time due to the high SFR depleting the reservoir.

Figure 7: Time evolution of the mass fraction (, solid lines) and volume fraction (, dashed lines) of the ISM ( 2 kpc and 1 kpc) in hot ( K), warm (100 K K) and cold ( K) phases. The warm gas dominates the ISM both in mass and in volume. In all runs except for G1D01_noFB, the hot gas occupies a non-negligible volume () but contains little mass, while the cold gas shows the opposite behavior and slightly increases over time due to metal enrichment. In G1D01_noFB, the hot gas is absent and the cold gas decreases over time due to the rapid depletion of gas by star formation.

4.3.3 Feedback-driven galactic outflows

Supernova explosions inject energy into the ISM and may push the gas out of the disc and drive galactic outflows. We define the outflow rate as the mass flux integrated over a chosen surface area :


where is the gas density, is the gas velocity and is the normal direction of the surface (in the obvious sense of outwards from the disc). In SPH simulations, can be estimated as where and are the mass and velocity of particle , and is a predefined thickness of the measuring surface. The summation is over particles within the shell of the surface with 0. The choice of is a compromise between being too noisy ( too small) and over-smoothing ( too large). The inflow rate is defined in the same way but with the direction of reversed.

As a proxy for how much gas is leaving or entering the disc, we measure the outflow rate and inflow rate at the planes kpc, truncated at kpc. The choice of kpc is to avoid counting the thick disc material as outflows/inflows, especially at larger radii where the gas disc flares, (though the distinction between disc material and outflows/inflows is somewhat arbitrary). The thickness of the planes is set to 0.2 kpc but the results are not sensitive to the chosen value. To assess how much gas is escaping the dark matter halo, we also measure the outflow/inflow rate at the virial radius = 44 kpc (denoted as ). Here the thickness of the sphere is set to 2 kpc due to the much lower density around the virial radius. The mass loading factor is defined as . There should be a time difference between when the star formation occurs and when the outflowing material causally linked to that star formation event reaches the measuring surface. Our definition of the mass loading factor does not account for this effect and thus may seem inappropriate. This is especially true for as it takes about 200 Myr for the gas to reach assuming the outflowing velocity of 200 km/s. However, as the SFR does not change dramatically throughout the simulations, can still be a useful proxy for how efficient feedback is at driving outflows.

In Fig. 8 we show the time evolution of in panel (a), and in panel (b), and and in panel (c), respectively. First we focus on the kpc planes. Comparing with Fig. 6, rises sharply right after the first star formation occurs. Since there is no gas above the disc in the initial conditions, the initial outflows gush out freely into vacuum and thus the mass loading is very high, ( 10). This initial phase of violent outflow lasts about 100 Myr. A significant fraction of gas falls back onto the disc due to the gravitational pull, which gives rise to the inflows starting from 150 Myr. From Myr, the system reaches a quasi-steady state with a small but nonzero net outflow . The large fluctuations of indicate that the outflows are intermittent on timescales of about 50 Myr. On the timescale of 1 Gyr, however, is rather constant and is slightly higher than the star formation rate ().

Now we turn to the spherical surface at = 44 kpc. Due to the time delay for gas to travel from the disc to , the first nonzero appears at 200 Myr. The peak value of is much lower than the peak of at 80 Myr. While this is partly due to the gravitational pull that slows down the outflowing gas, the more important reason is that not all of the gas which passed through the 2kpc planes eventually made it to . Some would later fall back onto the disc (galactic fountain) and some would fill the space in the halo region. No inflows at are detected ( = 0 at all times) so can be regarded as a net outflow. The mass loading first reaches about the order of unity and then decreases as the gas accumulates in the halo and hinders the subsequent outflows.

In Fig. 8 panel (d) we show the averaged velocity at and at 2kpc, respectively. The averaged velocity is mass-weighted over the particles within the measuring shells. It may seem inconsistent that is much larger than . This is not because the outflowing particles are accelerated, but is again due to the fact that not all particles crossing the 2kpc planes have enough kinetic energy to eventually escape the halo, and the averaged velocities are thus much smaller. The initial high is due to the vacuum initial conditions in the halo region, and over time decreases as the halo gas hinders the subsequent outflows.

In Fig. 8 panel (e) we show the total gas mass in different spatial regions as a function of time. The disc region is defined as 6 kpc and 2 kpc, and the halo is defined as excluding the disc region. We denote , , and as the total gas mass in the disc, in the halo, and outside of the halo, respectively. The total mass of stars that have formed in the simulations is denoted as and is well confined within the disc region. Except for the initial phase of strong outflow, there is roughly the same amount of gas ejected into the halo as the amount of stars formed in the simulations (), both of which are more than one order of magnitude lower than . Therefore, hardly changes throughout the simulations, and as a consequence the system is able to settle into a quasi-steady state. The escaped mass is about 10-30% of (and ), which is consistent with 1 shown in panel (c).

In Fig. 8 panel (f) we show the mass weighted mean metallicity of the gas in the disc, Z, in the halo, Z, and for gas that has escaped the halo, Z. Starting from 0.1 Z, the disc metallicity slowly increases due to metal enrichment from the stellar feedback. At the end of simulations ( = 1 Gyr), Z only increased to about 0.15 Z in G1D001 and G01D01 and about 0.12 Z in G1D01. The slightly slower enrichment in G1D01 is due to its lower , a consequence of more efficient photo-electric heating that suppresses gas cooling and star formation. The metallicity in the halo Z is only slightly higher (about 20%) than Z, as not only the highly enriched SN-ejecta but also the low-metallicity gas in the ISM is pushed out of the disc, and Z is thus diluted by the low-metallicity gas. The metallicity of the escaped gas Z shows a sharp peak as high as 0.3 Z initially, but also drops to about the same level as Z very quickly (less than 50 Myr). In general, the outflows in our simulations are enriched winds, moderately.

Figure 8: Time evolution of outflow-related physical quantities in G1D01 (black), G1D001 (blue) and G01D01 (red). Panel (a): the inflow rate at = 2 kpc (). Panel (b): the outflow rate at = 2 kpc () and at = (). Panel (c): the mass loading factor at = 2 kpc () and at = (). Panel (d): the mass-weighted outflow velocity at = 2 kpc () and at = (), and inflow velocity at = 2 kpc (). Panel (e): the total gas mass in the disc (), in the halo (), outside of the halo (), and the total mass of stars formed in the simulations (). Panel (f): the mass-weighted metallicity in the disc (Z), in the halo (Z) and outside of the halo (Z).

4.4 ISM Properties

Having described the global properties of the simulated galaxies, we now turn to the local properties of the ISM and its multi-phase structure.

4.4.1 Thermal properties of the ISM

In Fig. 9 we show the phase diagrams (temperature vs. density, left panels) and the corresponding cooling and heating rates of different processes (median value, right panels) in the ISM region ( 2 kpc and 1 kpc) at = 300 Myr in run G1D01, G1D001, G01D01 and G1D01_noFB (from top to bottom). The black curves in the left panels show the (thermal-)equilibrium temperature-density relation as shown in Fig. 2. The dashed lines show the contours where . Particles lying below the dashed lines are unresolved in terms of self-gravity. Except for G1D1_noFB where the gas has undergone runaway collapse, the majority of gas lies in the Jeans-resolved region up to 200 cm. In the right panels, the thick dashed black curves represent the shock-heating rate, i.e., the viscous heating. Note that this shock heating rate does not include the thermal injection of the SN feedback, which instantaneously heats up the nearest 100 gas particles as described in Section 2.6.1. The thick solid blue and dashed brown lines show respectively the cooling and heating rate caused by adiabatic expansion and compression where is the specific internal energy and is the velocity divergence.

For 100 cm, the gas is approximately in thermal equilibrium at around few tens of Kelvin as the cooling time is much shorter than its local dynamical time (can be estimated by ). A small bump at the highest densities is due to H formation heating. Below = 100 cm, the scatter starts to increase due to feedback-driven turbulent stirring. For 1 cm, the adiabatic heating/cooling dominates over other radiative processes, and the gas spends several dynamical times before it can cool. Once the gas cools below K, it is likely to be shock heated to K in a dynamical time (gas being shock heated to K is possible but at such temperatures the gas will quickly cool down to K). Therefore, turbulence drives the gas temperatures away from the equilibrium curve and keeps them at K at 1 cm (see also Walch et al., 2011 and Glover & Clark, 2012c).

The distribution of gas in the phase diagram is not very sensitive to the choice of and . This is because SN feedback, which enhances turbulence and drives the gas out of thermal equilibrium in the range of cm where the equilibrium curve is most sensitive to and (Fig. 2), does not depend explicitly on and . In G1D01 the gas distribution follows the equilibrium curve better as the photo-electric heating is more efficient and dominates over shock heating. In G1D1_noFB the gas also follows the equilibrium curve and shows the least scatter, as the turbulent stirring is much weaker without SN feedback. There is no hot gas in G1D1_noFB because the hot gas is only generated by the SN feedback.

Figure 9: Phase diagrams (temperature vs. density, left panels) and the cooling and heating rates of different processes (median value, right panels) in the ISM region ( 2 kpc and 1 kpc) at = 300 Myr in run G1D01, G1D001, G01D01 and G1D01_noFB. The black solid lines shows the (thermal-)equilibrium temperature-density relation as shown in Fig. 2. The dashed line shows the contour where . SN feedback keeps the gas temperatures at K in the range of cm even if the equilibrium temperatures are much lower. As a result, the distribution is not sensitive to and , both of which affect photo-electric heating but not SN feedback.

4.4.2 Molecular hydrogen

In the top row of Fig. 10, we show the H mass fraction vs. the number density . The dashed lines represent the star formation threshold = 100 cm. Since most gas particles with are also cooler than 100 K (cf. Fig. 9), we can regard those located to the right of dashed lines as star-forming gas while those to the left are not star forming. In general, increases with as higher density leads to a higher H formation rate and usually implies more shielding. Except for G1D01_noFB, a large fraction of star-forming gas has low . The H mass fraction of gas with 100 cm is = 10%, 0.7%, 18% and 70% in G1D01, G1D001, G01D01 and G1D01_noFB, respectively, though the precise values depend on our definition of star-forming gas (). This implies that most star formation actually occurs in regions dominated by atomic hydrogen (except for G1D01_noFB), in agreement with previous theoretical predictions (Glover & Clark, 2012c).

On the other hand, there are also particles whose densities are too low to be star forming but that have 1. The fraction of the total H mass that resides in 100 cm is 42%, 25%, 70% and 4% in G1D01, G1D001, G01D01 and G1D01_noFB, respectively. This demonstrates that, just as in our own Galaxy, H-dominated regions are not necessarily star-forming333 We note that if we were able to follow the collapse to much smaller scales and much higher densities (e.g. prestellar cores), eventually all star formation will occur in H. Our definition of star-forming gas corresponds to the reservoir gas that correlates (both spatially and temporary) with star formation.. This is in line with Elmegreen & Hunter (2015) who also concluded that there should be a significant amount of diffuse H in dwarf galaxies.

As discussed in Section 3.1.1, the formation and destruction timescales of H are long compared to the free-fall time. In the presence of feedback, turbulence further reduces the local dynamical time. Therefore, the gas can easily be out of chemical equilibrium. In the bottom row of Fig. 10, we show vs. the dimensionless quantity (see Section 3.1 for the definition). The black curves indicate the - relation in chemical equilibrium (Eq. 14). Gas particles with under/over-abundant lie below/above the black curves. The gas would follow the equilibrium curves if the chemical timescales were much shorter than the local dynamical timescales. The large scatter of the distribution indicates the opposite: the gas is far from chemical equilibrium locally. However, since there are both significant amount of particles under-abundant and over-abundant in H, the total H mass of the galaxy is not very different from the equilibrium predictions (cf. Fig. 6).

As shown in Fig. 6, the global value of is sensitive to the assumed and . This is also true locally: the particle distribution in Fig. 10 differs dramatically when adopting different and . This is in contrast to the particle distribution in the phase diagram (Fig. 9) which is rather insensitive to variations of and . The largest scatter is in G01D01, as the weak ISRF leads to a lower H destruction rate. A significant fraction of over-abundant H can therefore survive for a longer time. In G1D001, the low DGR makes it difficult to have high even for the densest gas. Low means little self-shielding (dust-shielding is almost absent given its DGR) and thus shorter H destruction time. In addition, although its H formation time is the longest among all runs, it takes much less time to reach the low equilibrium mass fraction . The scatter in G1D001 is therefore the smallest. Despite the absence of feedback, the G1D1_noFB run also shows some scatter as a result of its turbulent motions due to thermal-gravitational instabilities. Nevertheless, there is a large fraction of H-dominated ( 1) gas that follows the equilibrium curve quite well. This explains the excellent agreement between the global and the equilibrium prediction in Fig. 6.

Figure 10: Top panels: the H mass fraction vs. the number density at = 300 Myr. The dotted line indicates the star formation threshold density = 100 cm. Most star formation actually occurs in regions dominated by atomic hydrogen, while H-dominated regions are not necessarily star-forming. Bottom panels: vs. the dimensionless quantity at = 300 Myr. The solid line shows the - relation in chemical equilibrium (cf. Eq. 14). Particles with under/over-abundant lie below/above the black curves. The gas is far from chemical equilibrium locally.

To see in which phase most of the H can be found, we show in Fig. 11 the phase diagram weighted by the H mass. The projected histograms of density and temperature (also weighted by ) are shown on the upper and right-hand sides of each panel, respectively. Most of H can be found in the cold and dense phase of the ISM. In G1D01 and G01D01 there is some H in the warmer (T 1000 K) and less dense phase. This is the H-over-abundant, high gas shown in Fig. 10, which would not be predicted by the equilibrium approach. Note that the H distribution does not necessarily reside in the star forming regimes (T 100 K and 100 cm). Again, this implies that non-star-forming regions can also be H-dominated.

Figure 11: Phase diagrams weighted by the H mass () at = 300 Myr. The projected histograms of density and temperature (also weighted by ) are shown on the upper and right-hand sides of each panel, respectively. Most H is located in the cold (T 100 K) gas. A significant fraction of H is located in the warm and diffuse phase in G1D01 and G01D01 due to non-equilibrium effects.

4.5 Radial Variations

4.5.1 Scale height of the gas disc

We define the scale height as the altitude where 75% of the gas mass between is enclosed between in each radial bin. The choice of = 2 kpc is to discard the particles that clearly do not belong to the disc in any reasonable definition. Setting = 1 kpc gives essentially identical results. The definition of 75% enclosed mass is motivated by the distribution for a self-gravitating disc (). Similarly, we define in the same way as but for 95% of enclosed mass (as ). We also define the scale height of the cold ( 100 K) gas and the H in the same fashion, denoted as and , respectively.

In Fig. 12 we show (green), (blue) and (red) in the solid/dashed curves as a function of , (arithmetic-)averaged over time from = 250 Myr to = 1 Gyr with a time interval of 1 Myr. The filled areas show the bands where is the standard deviation of in each radial bin. The size of each radial bin is 0.1 kpc. The scale height of the total gas and both increase with , especially at larger radii. On the other hand, 75%(95%) of the cold gas and H is found within a thin layer 0.1(0.2) kpc without flaring, in line with what Elmegreen & Hunter (2015) inferred. There is some apparent flaring of and for the band (but not the mean value) at 2 kpc in both G1D001 and G01D01. However, very little cold gas and H can be found at such radii (cf. Fig. 13). Star formation mainly takes place in the mid-plane 0.1 kpc, as this is where the gas density is high enough to cool efficiently and trigger gravitational collapse. The cold and dense environment in the mid-plane also favors H formation. Once the H is pushed above the disc by SN feedback, it would be destroyed by the unshielded ISRF very quickly. Therefore, the H is also confined to the mid-plane.

Figure 12: Scale height of the disc as a function of . is defined as the altitude where 75% (95%) of the gas mass between is enclosed between . The scale height of the cold ( 100 K) gas and the H are defined in a similar way. Results are averaged over time from = 250 Myr to = 1 Gyr with a time interval of 1 Myr. The solid lines represent the average values in each radial bin and the filled areas show the bands. The disc thickens at large radii. On the contrary, the cold atomic and molecular gas are confined within a thin layer in the mid-plane without flaring.

4.5.2 Radial profile

In Fig. 13 we show the radial profile of the surface density of the total gas (), star-forming gas (), H (), star formation rate (), as well as the mid-plane number density () and mid-plane pressure (), averaged over time from = 250 Myr to = 1 Gyr with a time interval of 1 Myr. The solid lines are the arithmetic average and the filled areas are bands in the given radial bins ( bands become negative in some bins and thus can not be shown in log-scale). The mid-plane is defined as a thin layer within 0.1 kpc.

An interesting feature in Fig. 13 is that drops much faster than as increases, i.e., the gas depletion time