© ESO, 2013

1. Introduction

Numerical simulations of galaxy formation have now reached a high degree of sophistication and include increasingly detailed physics. This ranges from star-formation rates (SFRs) treated with sub-grid recipes based on the Kennicutt-Schmidt (KS) law to more direct processes triggered by Jeans instabilities (Stinson et al. 2006; Hopkins et al. 2011), feedback that is treated either through supernovae heating or through kinetic impulse, momentum-driven flows due to stellar winds (Sales et al. 2010; Ostriker & Shetty 2011), and in some cases, a multiphase gas (Maio et al. 2007; Gnedin et al. 2009). However, big unknowns remain as free parameters in the process: first, there is the nature of dark matter and its behaviour, leading to well-identified problems at the galactic scale, such as predicted cuspy profiles (while only cores are observed) or the prediction of a large number of dark satellites around each Milky-Way type galaxy. Second, a large unknown is also the nature and location of the missing baryons: observationally, only less than half of the baryons have been identified (e.g. Fukugita & Peebles 2004), and this missing baryonic component is certainly gaseous, given the results of microlensing (e.g. Wyrzykowski et al. 2009). This gas might be hot (with a temperature of the order of millions of kelvins) or cold, and most of it must be located in the intergalactic space to avoid overpredicting rotation curves.

A significant fraction of dark baryons should also exist in galaxies and could reside under the form of cold gas that is dense enough to be in the molecular phase (e.g. Pfenniger & Combes 1994; Grenier et al. 2005; Bournaud et al. 2007; Langer et al. 2010). A reservoir of cold gas may help to moderate star formation and explain almost stationary star-formation histories in late-type galaxy discs as observed today (e.g. Wyse 2009). When the thermal evolution of gas is studied during galaxy formation, an immediate conclusion is that almost all the baryons should have cooled and condensed in dark matter potential wells, where they are supposed to form stars, which are not seen. This overcooling problem can be solved through re-heating the gas, even before its inflow in the dark matter bound structures (Blanchard et al. 1992; Davé et al. 2001). However, the star-formation efficiency might be lower in globally low gas density environments, such as the outer parts of galaxies, and cold gas reservoirs are another solution to explore. The timescale spent by the gas in the cold phase is largely unknown, and this could change the baryon cycle. Cosmological simulations are unable to derive the abundance of cold gas by lack of spatial resolution, which leads to a reduction of density dynamics and underestimates the cooling (that increases with density) locally.

In recent years, the cold mode gas accretion to assemble mass onto galaxies has been revived versus the hot mode gas accretion. It has been realised that only part of the gas was heated to the virial temperature of a structure before cooling and inflowing, and most of the gas could be accreted cold (Kereš et al. 2005). The cold gas accretion along the cosmic filaments dominates for small haloes that are lighter than 10¹² M_⊙, while the hot quasi-spherical mode dominates for massive systems. There is also a redshift dependence on the importance of the two modes, where the cold mode is more dominant at high redshift (z larger than 2). Birnboim & Dekel (2003) and Dekel & Birnboim (2006) found similar conclusions, using one-dimensional simulations and analytic arguments. The limit between dominating cold and hot modes corresponds to a baryonic mass close to the scale that separates the red sequence and the blue cloud in galaxy populations (Kauffmann et al. 2003; Baldry et al. 2004), and the origin of bimodality in observed galaxies has been connected to these two different theoretical assembly modes. Even in massive structures, cosmological simulations have shown that cold filaments of gas may penetrate their hot gaseous halo (Birnboim et al. 2007; Dekel et al. 2009).

Another well-known problem in the baryon cycle and galaxy formation is to correctly reproduce the galaxy stellar mass function. If the large number of predicted dwarf galaxies is suppressed by supernovae feedback that is assumed to expel gas out of the small potential wells, gas is still accreted in massive systems and forms growing blue star-forming galaxies at the massive end of the mass function. Solutions have been proposed through AGN feedback (either quasar mode through heating and outflows, or radio mode through powerful radio jets) with limited success. In particular, the gas expulsion is more efficient again in low-mass galaxies. Even when ejected with a speed larger than the escape velocity, gas is re-accreted onto galaxies, slowed down not only through gravity but also by hydrodynamical interactions with the halo gas. The re-accretion has been called a halo fountain by Oppenheimer & Davé (2008). According to the environment, gas that outflows from low or intermediate mass haloes can be re-accreted by more massive companion galaxies – a process that has been dubbed the intergalactic fountain (Kereš et al. 2009). Globally, the frequent re-accretion of the ejected gas limits the efficiency of the AGN feedback, maintaining some overcooling problem and keeping massive galaxies active in star formation until late times with no definite quenching. It is possible to introduce a cold gas reservoir in the baryon cycle with a sufficiently long timescale in this phase to account for the observed properties of galaxies. In the present paper, we want to explore this possibility, beginning by isolated galaxy simulations. We consider a cold dense molecular phase in the gas component, which can cool down to temperatures much lower than the 10⁴ K that is usually assumed for “cold gas” in cosmological simulations.

Thanks to the always increasing computing capacities, there has been a growing number of simulations in recent years considering the cold phase. These works come from two different domains with different goals:

• 1.

  Cosmologically oriented work, tending to higher spatial resolution, of the order of 100 pc. The main interest here is to determine whether the small scale structure impacts the large-scale transfer of angular momentum, the stability of discs, formation of bulges, concentration of the galaxy, transformation of galaxies from late to early type, and whether feedback can destroy dark matter cusps, etc. (Maio et al. 2007; Gnedin et al. 2009; Schaye et al. 2010; Murante et al. 2010; Bournaud et al. 2010; Hopkins et al. 2011).

• 2.

  Star-formation oriented work, tending to enlarge the field of view from molecular clouds to the whole galaxy structure. These simulations are in general multi-phase with the goal of answering questions about the cloud formation and the star-formation efficiency. The galaxy models are often not self-consistent with no transfer of angular momentum from gas to stars or dark matter halo. (Slyz et al. 2005; Tasker & Bryan 2006, 2008; Dobbs & Bonnell 2008; Dobbs et al. 2011; Shetty & Ostriker 2008; Wada et al. 2011).

Some works are at the transition between the two scales. When the simulations are unable to resolve the cloud fragmentation scale, these works have tried to represent the multiphase aspect of the gas using sticky particles. These particles move ballistically in the potential and are able to collide and coagulate to form a whole mass spectrum of clouds (e.g. Semelin & Combes 2002; Booth et al. 2007; Revaz et al. 2009).

In the present work, our goal is to deal with the multi-phase of the gas in a self-consistent way, allowing the gas to cool down to 100 K, which is a reasonable temperature for the quiescent dense molecular phase, and to follow its interaction with the other phases through cooling/heating, star formation, and feedback through a Tree-SPH code. The H₂ molecule formation is simplified by a recipe related to the density, metallicity, and self-shielding of the gas. The star formation recipe is based on the Kennicutt-Schmidt law, and core-collapse supernovae feedback is considered.

The numerical techniques and the initial conditions of the simulations carried out are described in Sect. 2 with the details of the baryonic physics adopted to deal with cold molecular hydrogen, star formation, and feedback. Section 3 describes the results of the simulations. The influence of the feedback and inclusion of molecular hydrogen cooling is detailed. Our conclusions are drawn in Sect. 4.

2. Numerical techniques

We use the Gadget-2 code (Springel 2005) that computes gravitational and hydrodynamical forces. Gravity is computed by a Tree algorithm with a Barnes Hut opening angle θ of 0.7, considering only the monopole moments of the gravitational field for computational reasons (see Springel 2005 for details). We take the same constant softening length ϵ for all particle types (stars, gas, and dark matter). Hydrodynamics is treated with a smooth particle hydrodynamics (SPH) algorithm with an individual smoothing length h, that is computed such that a constant mass is contained in a sphere of radius h and that is allowed to be as low as 0.1 ϵ. We add baryonic physics to this code, as described in this section.

2.1. Cooling

Gas in the interstellar medium (ISM) loses internal energy by radiation due to physical processes occurring at microscopic scales. This cooling can be encapsulated in a cooling function that represents the corresponding volume rate of energy loss and depends on the chemical composition of the gas, its density, and temperature. In simulations of galaxies, cooling is often considered only down to a temperature of ~10⁴ K, since the most efficient cooling processes are due to ionised hydrogen and helium, which are almost only present in the atomic form below 10⁴ K. This usually implies the coldest dense gas is at an equilibrium temperature of about 10⁴ K in these simulations, because dense gas experiences heating pressure forces and shocks if it reaches lower temperatures: the gas is prevented from collapsing further. We include cooling by metals, molecular hydrogen, and HD down to 100 K, which is a more realistic temperature minimum for the ISM.

In dense media, cooling by metals might dominate at solar metallicity. However, metals are abundant only in the central parts of galaxies, and radial abundance gradients are observed in giant spiral galaxies with an exponential decline of about 0.6 dex in 10 kpc on average (e.g. Henry & Howard 1995; van Zee et al. 1998). Therefore, the outer parts of galaxies are reminiscent of the primordial galaxies with low metal abundance. Cooling through the H₂ and HD molecules is then important and could change the physics of star formation in these regions. Ultraviolet observations by GALEX have shown that star formation can be active at large radii, much farther than the optical radius R₂₅, in regions where H_α observations are not able to reveal moderate-age populations of stars. UV-bright discs extend up to 3–4 times the optical radius in about 30% of spiral galaxies (Thilker et al. 2005; Gil de Paz et al. 2005, 2007). The presence of molecular hydrogen and star formation in outermost discs of spirals is of prime importance to study cold gas accretion, which is considered one key factor in galaxy evolution (e.g. Kereš et al. 2005; Dekel et al. 2009).

2.1.1. H, He, and metals cooling above 10⁴ K

We use the cooling functions computed by Sutherland & Dopita (1993) for a plasma in collisional ionisation equilibrium above 10⁴ K. These cooling functions include cooling due to ionised H, He, and metals and are tabulated for different metallicities.

2.1.2. Metal-line cooling below 10⁴ K

Following Maio et al. (2007), we compute cooling functions of a few metals below 10⁴ K, CII, OI, SiII, and FeII, as they are the most abundant heavy elements released by stars in the ISM and are thus the main ingredients for cooling. We consider their collisions with H atoms and electrons. For low densities obtained in galaxy simulations, the populations of energy levels needed to compute the cooling rates differ from the Boltzmann populations of the local thermal equilibrium (LTE). The populations here depend on the abundances of species that can collide with the metal and allow it to change energy level, which makes the rate of energy loss per metal dependent on these abundances (in contrast to the LTE case, where the rate depends only on temperature). The resolution of the equations governing the populations is included in the code using the quantum data given in Maio et al. (2007), and the cooling functions are computed using the obtained populations and the transition probabilities, as in Maio et al. (2007).

For the metals we consider here we choose the same abundances as in Sutherland & Dopita (1993; solar abundances and primordial ratios) and a low electronic fraction , assuming the cold gas is quasi-neutral.

The cooling curves for the hydrogen (atomic or ionised) number density n_H = 1 cm^-3 and different metallicities Z = [Fe/H] are plotted in Fig. 1.

2.1.3. H₂ cooling

We are interested in studying the influence of molecular hydrogen cooling on galaxies, especially in regions where metals are not abundant. We tabulate the LTE cooling function as a function of temperature from the Boltzmann energy level populations and transition probabilities. For low densities where the cooling function is lower than the LTE cooling function, we use data from Glover & Abel (2008) for a 3:1 ortho-para ratio of H₂. We assume H₂ cooling below 10⁴ K takes place in an almost neutral medium and thus consider collisions with other H₂ molecules, H atoms, and He atoms. Collisions with H atoms are often the only one considered. However, as we are interested in describing media with a high fraction of H₂, considering the collisions with H₂ and He is necessary to compute a cooling that would otherwise vanish in regions poor in atomic hydrogen.

The final expression for the volume cooling rate (in erg cm^-3 s^-1) we use in the code is: (1)Λ_H2,LTE is the LTE cooling rate per H₂ molecule in erg s^-1, and Λ_H2,H, Λ_H2,H2, and Λ_H2,He, in erg cm³ s^-1, account for collisions with respectively H, H₂ and He, for low densities. This formula, similar to what Hollenbach & McKee (1979) used, interpolates between the low density and high density regimes. In Fig. 2 we plot this cooling rate for a fixed hydrogen-nuclei number density n_Hnucl = 1 cm^-3 (such that n_Hnucl = n_H + 2n_H2) and different H₂ mass fractions f_H2, with . We do not consider cooling due to collisions of H₂ molecules with metals, as collisional coefficients are not well determined. We also neglect the cooling due to H, because although its abundance can be equal to the H₂ abundance around shocks, it then drops by orders of magnitude below. H cooling might play a role in strong starburst galaxies or in primordial galaxies at high redshift (e.g. Ricotti et al. 2001; Ahn & Shapiro 2007; Petkova & Maio 2012).

In Fig. 3 we show the volume cooling rate due to H₂ and metals for a fixed hydrogen-nuclei number density n_Hnucl = 1 cm^-3. We show the volume cooling rates for different metallicities and H₂ mass fractions. What is plotted is , where Λ_H2 is the volume cooling rate of Eq. (1), and Λ_N,met is the volume cooling rate of Fig. 1 divided by (1 cm^-3)² (the scaling being such for these low densities of hydrogen atoms). We see that H₂ can bring a significant contribution to the cooling even for solar metallicity gas.

2.1.4. HD cooling

We also consider HD cooling. Despite its lower abundance than molecular hydrogen, the molecule HD can help cool the gas because of its permanent electric dipole. The quadrupole transition probabilities are similar to the ones of H₂, but the dipole transition probabilities are orders of magnitude greater, which makes the cooling per molecule significantly greater. We use quantum data to compute the LTE cooling function in a similar way that was used for H₂. For low densities, we took the function of Lipovka et al. (2005). This cooling function considers the collisions of HD with H atoms and only dipole transitions, which is justified by their much bigger contribution to the cooling than quadrupole transitions. We followed Glover & Abel (2008) by assuming the cooling function is simply proportional to the H-atom density for low densities, so that the cooling function per molecule is n_HΛ_HD,H(n_H = 1 cm^-3) for low densities.

We write the volume cooling rate, like for H₂, as: (2)where Λ_HD,LTE is the LTE cooling rate per HD molecule in erg s^-1, and Λ_HD,H, in erg cm³ s^-1, accounts for collisions with H atoms.

In the simulations, we assume that n_HD = 10^-5 n_H2 because the fractionation of HD with respect to H₂ occurs only when the gas temperature becomes as low as 100 K. The abundance can then reach 200 times the normal abundance. The molecule HD can help the cooling of the gas down to 30 K in a primordial gas without metals (e.g. Yoshida 2006). However, these temperatures are below our minimum considered. With the abundance of n_HD = 10^-5 n_H2, the influence of HD in cooling the gas is little in the simulations of this paper.

2.1.5. Cooling algorithm

The interstellar gas component is modelled as an ideal gas with an adiabatic index γ that is set to . Gadget-2 integrates the evolution of the specific entropy rather than the internal energy for accuracy reasons developed in Springel & Hernquist (2002). The evolution of the specific entropy of a particle i is governed by (3)where Λ is the volume cooling rate in erg cm³ s^-1 (the sum of the contributions described in above), Π_ij is a factor accounting for the artificial viscosity, and is the symmetrised SPH kernel.

As the cooling time can be lower than the time-steps equal to the minimum of (where η is an accuracy parameter that we keep at 0.025 as in the default Gadget-2 parametrisation, ϵ is the gravitational softening length, and a is the acceleration) and the Courant time limitation for gas particles, we use an implicit cooling scheme to stabilise the resolution. We solve iteratively the implicit equation (4)where   dt is the considered time-step. This implicit scheme is however not critical for our types of simulations, which have little shocks and have no atomic/molecular heating terms.

The relation between temperature and specific entropy is (5)where is the mean molecular weight, with x_i the number fraction of the species i and m_p the proton mass. For a gas composed of atomic and molecular hydrogen and helium and if the metals contribution is neglected,

(6)where X is the hydrogen-nuclei mass fraction that we set to 0.76. We do not compute the ionisation fraction in the simulations. We assume the gas is quasi-neutral below 10⁴ K. Above that temperature, we assume there is no molecular hydrogen and the gas is fully ionised (H⁺ and He²⁺). Because the mean molecular weight is higher for a neutral gas than for an ionised one, there is a specific entropy range for which the temperature obtained from Eq. (5) is above 10⁴ K by assuming a neutral gas while it is below 10⁴ K by assuming an ionised gas. We set the temperature of the gas to 10⁴ K in this specific entropy range.

2.2. Resolution

Distinguishing physical effects from numerical ones can be difficult if the resolution of a simulation is not well adapted. The number of particles setting the mass resolution, the gravitational softening length, and the hydrodynamical smoothing length need to be chosen such that they do not introduce unwanted artefacts while allowing for the simulation to run at a reasonable pace.

We want to describe the gas down to low temperatures. However if the temperature is too low, the Jeans length and Jeans mass of the gas are more poorly resolved, which can lead to numerical artefacts (Bate & Burkert 1997). Schaye & Dalla Vecchia (2008) use an effective equation of state at high densities that causes the Jeans mass to be independent of density, while Hopkins et al. (2011) take a density-dependent pressure floor to ensure the Jeans length is resolved. Bournaud et al. (2010) uses an equation of state that mimics the effect of cooling at all temperatures and a temperature floor at high densities. We take a temperature floor of 100 K for all densities. We note that the gas is stabilised by the rotation of the disc since we simulate disc galaxies.

The number of particles should ideally be the largest possible to have a good mass resolution, which allows the Jeans mass to be well resolved. We run simulations with 1 200 000 initial particles: a third are gas particles, a third are stellar particles, and a third are dark matter particles. Table 1 shows the particle masses that we have in the particular case of our galaxy model that is detailed in Sect. 2.6. New stellar particles have a mass that depends on the parameter N_g→⋆, the number of stellar particles produced out of one gas particle, as described in Sect. 2.4.

The choice of gravitational softening length ϵ depends on the number of particles: a value that is too small for a given number of particles introduces unphysical two-body relaxation in collisionless media (stellar components of galaxies and dark matter haloes), while one that is too large decreases the spatial resolution by “blurring” density features and does not allow the Jeans length to be gravitationally resolved. Gravitational softening lengths ϵ are usually taken as scaling with the mean inter-particle distance, which is considered to be of the order of the number of particles to the power for a 3-dimensional simulation. We take a softening length of ϵ = 100 pc for all particle types. This particular value is derived from the GalMer simulations (e.g. Di Matteo et al. 2007), which took ten times less particles for a softening of 280 pc. The GalMer simulations were isothermal at 10⁴ K, and since the temperature reaches down to 100 K here, we take a softening length that is a little inferior to the cubic root, while still allowing for an efficient computation on a few tens of computing cores.

The value of the gas smoothing length h should be small enough to allow for a good density resolution. Gadget-2 uses variable smoothing lengths: the densities and smoothing lengths are computed together, so that the mass contained in a sphere of radius h is a constant. We take a minimum h equal to a tenth of the gravitational softening. This minimum can have subtle consequences on the structure of the gas. In the presence of strong cooling and with no sufficient star formation or efficient feedback, a number of particles can be stuck at this minimum, making the corresponding regions increasingly denser with no possibility of getting more diffuse. This can be computationally very demanding and lead to a large work-load imbalance between processors.

2.3. H₂ fraction

Our goal is not to determine the molecular abundance through a detailed chemical scheme, which would be too sophisticated for our present approximated treatment (which does not consider radiative transfer). Krumholz et al. (2008, 2009), McKee & Krumholz (2010) have derived an analytic expression of the mass fraction of molecular hydrogen f_H2 in an idealised spherical cloud submitted to a uniform and isotropic radiation field. Including radiative transfer and the formation and destruction of H₂ through dust, assuming a steady state, they obtain (7)with (8)The parameter τ_C is the dust optical length and χ is a scaled UV flux, divided by the hydrogen-nuclei number density. The mass fraction f_H2 is set to zero when s > 2.

The dust optical length is . The dust cross-section to the Lyman-Werner radiation per H nucleus σ_d is set, by using the reference Milky Way value, to σ_d = Z′10^-21 cm² (e.g. Krumholz & Gnedin 2011) where Z′ is the metallicity normalised by the solar metallicity. The variable μ_H is the mean mass per H nucleus, and Σ is a column density Σ = ρL obtained from a local scale height: . We compute (∇ρ)_i, the gradient of the density of the particle i by: ρ_i(∇ρ)_i = ∑ _jm_j(ρ_j − ρ_i)∇_iW(r_ij,h_i). The scale height considers the variation in the density: it increases with density but is inversely proportional to its gradient, so it is lower in the case of large gradients encountered on the outer parts of density features like clumps or spiral arms. It is thus well adapted in computations of the column density used to determine the shielding from radiation.

We compute the χ parameter using the new stars’ radiation field G. We assume that the stars whose age is inferior to 10 Myr radiate in the Lyman-Werner band and can dissociate the H₂ molecules. We do not consider that older stars radiate in this band, so they do not have any effect on the fraction of H₂ in our simulations. We take with a scaling factor that depends on the resolution of the simulations and tune it so that the obtained H₂ fractions are consistent with observations. We do not follow the radiative transfer of the photons from young stars. The radiation flux produced by stars decreases with the distance r to a star in , as the gravitational field (we assume an ISM escape fraction equal to unity). We insert the computation of the flux received from stars younger than 10 Myr in the Gadget-2 gravitational tree functions. It is similar to computing the gravitational force due to the young stars only and can be easily done by adding a range of variables representing this contribution only. We thus use the gravitational tree to compute a flux proportional to the mass of new stars over the squared distance.

2.4. Star formation

We implement star formation in a stochastic way that reproduces a Schmidt law (Katz 1992): (9)where ρ_∗ is the volume stellar density, ρ_g the volume gas density, t_ff is the free fall time , and c_⋆ is the star-formation efficiency per free-fall time.

The Schmidt law can be enforced by giving a gas particle a probability to spawn a star at each time-step Δt: (10)This implementation means that a gas particle of mass m_g can spawn a star particle of mass m_∗ with this probability at each time-step. The mass of the gas particle is then reduced by the amount of m_∗ until there is no more mass in the gas particle. The number of N_g→⋆ of star particles created by a gas particle is a compromise between a good mass and time resolution of star formation and CPU cost. If several stellar particles are created from a single gas particle, star formation is better resolved but the total number of particles can increase significantly, which slows down the execution of the code. In the case where N_g→⋆ is greater than one, we note that this spawning scheme implies that gas particles have different masses, of either the initial one m_g, or a smaller multiple of . The density and smoothing length are computed by Gadget-2, so that the mass contained in a sphere of radius h_i is fixed with h_i being the smoothing length of the particle i. We were careful about modifying the algorithm to consider the different masses, so that this mass condition is still satisfied.

We use common selection rules for the particles that are allowed to spawn stars: they must have a density higher than a threshold n_Hmin, have a temperature lower than a maximum temperature T_max, and be in a converging flow (∇.v < 0).

2.5. Feedback

Our simulations include core-collapse supernovae kinetic feedback. We consider that a supernova explosion releases the canonical value E_SN = 10⁵¹ erg. A fraction α_fb of this energy is given to the ISM, while the rest is radiated away. Each new stellar particle of mass m_⋆ inputs an energy E_input = α_fbϵ_SNm_⋆ where ϵ_SN is the number of formed supernovae per formed stellar mass multiplied by the canonical energy E_SN = 10⁵¹ erg. The parameter ϵ_SN is thus the supernovae energy released per formed stellar mass. We consider the number of formed supernovae per formed stellar mass to be 0.01, which is of the order of the fractions obtained from commonly used initial mass functions (IMF) such as a Salpeter IMF with a slope = −1.35 and lower and upper limits of 0.1 and 40 M_⊙, respectively, where stars more massive than 8 M_⊙ are considered to be supernovae, so we have .

Each neighbour i of a new star particle 0 receives an energy weighted by its distance to the new stellar particle: (11)so that the sum of the energies given to the neighbours is E_input. If the newly created stellar particle has left a remnant gas particle, no feedback energy is given to the remnant. We recompute the smoothing length and list of neighbours at the position of a new stellar particle, considering only the neighbouring gas that has not been turned into stars (and not considering a possible gas particle remnant at the exact position of the new stellar particle). Therefore, the sum of the fractions involving kernels is indeed unity, and the mass of the gas affected by each supernova explosion is a constant.

The neighbours are given a velocity kick , directed along the line joining the stellar particle and the neighbour and away from the new stellar particle. The feedback is only kinetic. This input kinetic energy can however be converted into internal energy through the SPH viscosity.

The number N_g→⋆ of stellar particles created out of one gas particle has an influence on the distribution of feedback energy. Feedback is more gradually input for a larger N_g→⋆.

2.6. Initial conditions

We consider the case of a giant Sb galaxy.

The gaseous and stellar discs follow Miyamoto-Nagai density profiles. The density of the gas disc is in cylindrical coordinates: (12)and the stellar disc density is: (13)The Miyamoto-Nagai profiles are chosen because the associated potential is analytic, and the initial velocity dispersions can thus be easily computed, but relaxation makes the profiles quickly exponential as observed.

We use the gravitation routines of Gadget-2 to obtain the forces acting on particles. We compute the circular velocity of disc particles using these gravitational accelerations and add an analytic asymmetric drift correction to have a more realistic velocity profile.

The radial velocity dispersion is derived from , where Q is the Toomre parameter that we set to 1 for both discs, Σ is the surface density, and κ is the epicyclic frequency derived from the potential. The azimuthal velocity dispersion is obtained from where Ω is the angular speed. The vertical dispersion is set by the isothermal equilibrium of the disc.

The spherical stellar bulge and dark matter halo have Plummer density profiles. The bulge density is thus given by (14)and the halo density is (15)All the masses and characteristic lengths are specified in Table 2.

The velocity dispersion of the spherical components is chosen to be isotropic and is derived from the second moment of the Jeans equation. The radial dispersion is thus (16)The velocity curve for this analytic model with the contributions of the different components is shown in Fig. 4. The rotation is due mainly to the stellar components near the centre of the galaxy and to dark matter at large radii.

We prepare the initial conditions by letting the galaxy evolve for 300 Myr with only gravitational forces for any type of particles. This allows us to start the simulations with dynamically relaxed discs that especially exhibit no annuli instabilities. The initial surface densities of gas and stars are shown in Fig. 5, and a snapshot of the gas is shown on Fig. 6.

3. Simulations and results

We assume that there is a metallicity gradient in the gas. We take a central metallicity Z = [Fe/H], such that Z(R = 0) = Z_⊙ + 0.5 and assume that the metallicity decreases by 1 dex per 10 kpc: (17)where R is the cylindrical radius. The influence of the cooling by metals thus decreases with the distance to the centre of the galaxy. For simulations with no molecular hydrogen, we expect the outer regions to have fewer density features and to be warmer than the central ones. If, however, we add some molecular hydrogen, depending on its fraction and on the details of star formation and feedback, the gas may be also able to form clumps in the outer regions and form more stars there.

The gas temperature is initially 100 K. The temperature floor we apply is more exactly a specific energy floor, corresponding to a temperature of 100 K for purely atomic gas and going up to 186 K for gas with hydrogen present only in the molecular form.

The simulations presented here all have a number of stars spawned by a gas particle N_g→⋆ = 4. Unless otherwise specified, the star-formation efficiency per free-fall time is set to c_⋆ = 0.1, and the number density threshold for star formation is n_Hmin = 10^-1 cm^-3. The number density n_H is hereafter the total hydrogen number density that was written as n_Hnucl in Sect. 2. We use a maximal temperature for star formation T_max = 30 000 K, which happens not to be selective, because only a few diffuse particles can reach this temperature in our simulations. We show some consequences of a variation in some of the star-formation parameters in Sect. 3.2.

3.1. Simulations without H₂

We first perform simulations with no molecular hydrogen with varying feedback efficiencies α_fb. We have four different feedback efficiencies: either no feedback (α_fb = 0), α_fb = 1%, α_fb = 10%, or α_fb = 40%.

Density maps of the gas discs for these runs are shown in Fig. 7 at three simulation times: 0.5 Gyr, 1 Gyr, and 3 Gyr. The gas can reach smaller values of local volume density when feedback is included, as star formation in dense regions inputs some energy in the ISM and prevents it from getting denser, which slows down the star formation. For the run without feedback, the central part of the galaxy exhibits a strong density contrast on the snapshots at 0.5 Gyr and 1 Gyr, many small clumps and thin spiralling filaments can be seen. If the feedback is switched on, we can see only a few clumps for a feedback efficiency of α_fb = 1% and smoother spiral features, and we see no clumps for α_fb = 10% and α_fb = 40%. A central thin bar is formed in all the cases. After 3 Gyr, the bar is still clearly visible for the higher feedback efficiencies but is less obvious otherwise, because the central parts have been depleted from gas by star formation.

Figure 8 shows the temperature-number density histograms after 0.5 Gyr of evolution, a time at which gaseous discs are actively forming stars in all the simulations. On the top and right of each plot, the marginal probability density functions (PDFs) of hydrogen-nuclei number density and temperature, respectively, are shown. The temperature PDFs show that the fraction of gas at T ≃ 10⁴ K increases with feedback, while the cold dense gas fraction decreases. On the density PDFs, we can see that gas reaches smaller maximum densities for higher feedback efficiencies, and there is an increasingly high fraction of diffuse gas. In the low feedback runs, a significant fraction of the dense gas has a temperature close to the minimum: this is the dense gas of the central parts, subject to metal-line cooling. The fraction decreases for higher feedback efficiencies because of the dissipation of energy by feedback, making the densest gas of the simulations warmer. If the gas is heated by pressure forces, viscous shocks, or feedback, its temperature does not reach much beyond 10⁴ K because of the stronger cooling that it undergoes above this temperature. The diagonal branches observed on the left of each plot account for the gas that is subject to very little cooling away from the centre of the disc because of the metallicity gradient and thus follows adiabatic paths.

We plot in Fig. 9 the time evolution of the total mass of gas present in the simulations. All the gas is originally in the disc but can leave it under the effect of gravitational heating or stellar feedback. The characteristic time of consumption of the gas increases with feedback efficiency because of the moderating effect of feedback on gas density, which regulates star formation. The curves have increasing horizontal asymptotes y-values: this is due to gas leaving the disc middle plane and also to the inability of the gas to reach the star-formation threshold. The gas in the outer parts of the disc remains diffuse with almost no star formation in all cases, as the lower abundance of metals does not allow the gas to cool down enough to form stars. More details concerning the star formation in these simulations are presented in the following section when compared to simulations including molecular hydrogen.

3.2. Simulations with H₂

We now study the impact of the inclusion of H₂ on the gas physical state, star formation, and structure of the discs.

3.2.1. H₂ fraction

We first set the feedback efficiency to α_fb = 1% and tune the χ factor in the H₂ mass fraction f_H2 so that we find a global H₂ mass fraction (global ratio of molecular to total hydrogen mass) that is consistent with observations while having enough molecular hydrogen to study its effects. Our metallicity gradient implies that H₂ tends to be less present in the outer parts of the galaxy, where the metallicity and dust abundance are low and the gas is less shielded from Lyman-Werner radiation. However, this can be compensated by the lower SFR in these regions, which decreases the ambient Lyman-Werner luminosity. We plot the evolution of the global H₂ mass fraction in Fig. 10 for four different UV flux scaling factors. The global mass fraction first decreases with time, because newly formed stars contribute to the dissociating radiation field, and becomes stable when the SFR becomes almost null (see the red solid lines of Fig. 21 for the evolution of the gas mass and the SFR in the χ × 50 simulation). The exact exposition of the gas to UV flux is not well known, since it depends on the physics at a very small scale, which is below our resolution, and on complex radiation transfer through gas clumps and associated dust. Figure 11 shows the surface density of H₂ and atomic hydrogen gas HI after 0.5 Gyr of evolution for the case of χ × 50 that we have chosen. Such a surface density profile is similar to some observed profiles of local galaxies described in Young & Scoville (1991) and more recently in Bigiel et al. (2008).

We plot the mass fraction of H₂ versus the hydrogen-nuclei number density at t = 0.5 Gyr in Fig. 12. The H₂ mass fraction is null for diffuse gas and equal to unity for the densest gas, while the transition presents some scatter. This is due to the variation in both the structure of the gas and the amount of nearby recent star formation. The transition density is low compared to similar work by Gnedin et al. (2009) or Christensen et al. (2012b), because of our calibration to obtain a realistic global H₂ mass fraction and of our generally lower densities compared to Gnedin et al. (2009) and Christensen et al. (2012b). These lower densities arise partly from our choice of a lower threshold density for star formation, which prevents the gas from becoming as dense. For simulations with a higher threshold density for star formation, the calibration factor (that is varied on Fig. 10) should be higher.

This method can be sensitive to the resolution of the simulation through the dependence of the molecular mass fraction on the column density proportional to and the χ factor proportional to . We have run simulations with a feedback efficiency of 10% and a varying number of particles from our reference N_ref = 1 200 000 to simulations with , , N_ref × 2, and N_ref × 4. The softening length ϵ has been modified respecting ϵ ∝ N^−1/3 with N being the number of particles, so that it is proportional to the average inter-particle distance. The other parameters, including c_∗, have been kept constant. The comparison between the simulations is delicate, because simulations of different resolutions have a different density structure evolution, impacting the dust optical depth, and a different star-formation history, impacting the χ factor. We represent the dust optical depth as a function of hydrogen-nuclei number density for the different resolutions at the same simulation time on Fig. 13. It shows practically no dependence on the resolution. The slight difference in the transition density between the mostly atomic and molecular hydrogen when the resolution is changed is due to the higher average UV field. This can be seen on the plot representing the distribution of the factor ln (1 + 0.6χ + 0.01χ²) entering the computation of the molecular mass fraction versus density.

3.2.2. Gas physical state

Figure 14 shows the aspect of the disc for α_fb = 1% at the same simulation epochs as in Fig. 7, whose second row is the simulation for the same feedback efficiency but without H₂. Cooling by collisions of atomic hydrogen with metals or by H₂ make the central parts similar in both simulations. The difference is in the outer parts of the galaxy, where metal-line cooling is poorly efficient because of our assumed metallicity gradient. In this case, the gas remains diffuse with no other cooling processes, but the inclusion of H₂ allows the gas to be clumpier: we see density features that were absent in the purely atomic simulation. The surface density of H₂ is also plotted. It indeed follows the density features of the gas: the clumps and filamentary structures can be seen in H₂. We represent the corresponding power spectra of the gas surface density on Fig. 15 for total gas masses normalised to the same value. It can be seen that the disc with molecular hydrogen has more small and intermediate size structures. We also plot the clumping factor evolution with time on Fig. 16. We define this clumping factor, as in Springel & Hernquist (2003), by (18)Figure 16 shows molecular hydrogen makes the gas clumpier.

The temperature-number density histograms of Fig. 17 and the marginal PDFs show higher fractions of gas in the coldest dense phase (with exact temperatures depending on the molecular fractions) than for the same feedback efficiencies without molecular hydrogen (see Fig. 8 for comparison). There is no clear diagonal branch, because there is some efficient cooling in the whole disc. We separate the gas in two ranges of temperature, below and above 1000 K, and study the evolution of the fractions of cold and warm gas (gas lying below or above this threshold). The majority of star formation happens in the first Gyr in all the simulations (especially when the feedback efficiency is low and stars form quickly), so we focus on this period and plot the cold gas fraction as a function of time on Fig. 18 for two feedback efficiencies and with or without the inclusion of cooling by molecular hydrogen. All the gas is initially in the cold phase. Feedback decreases the fraction of cold gas: the kinetic energy given to particles in dense star-forming regions is transformed into thermal energy by pressure forces and viscosity, even more so as the feedback efficiency is high. For a given feedback efficiency, the cold gas phase represents a much higher fraction of the gas if H₂ cooling is considered and stars are formed from colder and denser gas, increasing the star-formation efficiency.

As we are especially interested in studying the state of the gas as a function of galactocentric radius, we plot the histograms of gas density versus radius in Fig. 19 for the atomic and molecular simulations having α_fb = 1% after 0.5 Gyr of evolution. The effect on density is clear: the gas exhibits density peaks at larger radii and has a mean higher density with the inclusion of H₂. Some gas is denser than the star-formation threshold even at large radii. Star formation is detailed in the next section.

The interstellar gas we obtain is heterogeneous but could be even more so with different parameters for star formation. We also run simulations with a lower star-formation efficiency per free-fall time c_∗ = 0.01 and the same threshold density n_Hmin = 10^-1 cm^-3 and simulations with the same c_∗ = 0.1 but a higher n_Hmin = 10¹cm^-3. Reducing c_∗ or increasing n_Hmin both allow for higher density peaks, as seen on Fig. 20. The effect of including H₂ is also visible in these sets of simulations, but the H₂ abundance is overestimated here as the density is on average higher than in our other simulations (we run these new simulations without changing the calibration of the χ factor in the molecular mass fraction). When c_∗ is lowered or n_Hmin is heightened, the gas forms stars more difficultly and can thus collapse more efficiently, resulting in a more fractioned gas.

3.2.3. Star formation

The effect of H₂ cooling on the star-formation efficiency is visible in Fig. 21 that shows the total gas mass evolution with time and the SFRs for the different simulations. The depletion time of the gas decreases when molecular hydrogen cooling is included. We also show the SFRs as a function of time for our alternative star-formation parameters presented at the end of Sect. 3.2.2 in Fig. 22. The reduced ability of the gas to form stars results in a delayed peak of star formation after some time when the gas has got denser and denser through cooling before being consumed by star formation.

Having inserted the formation time of stars in our simulation outputs, we can track star formation spatially. Our outputs are temporally spaced by 10 Myr. We define the SFR as being the mass of stars formed during 10 Myr divided by this time, which is similar to the SFRs obtained from observations in H_α. We plot the cumulative SFR averaged on the first Gyr as a function of radius in Fig. 23, and we see that about the same amount of star formation occurs in the central parts, but molecular hydrogen starts playing a role at large radii. The difference occurs at a larger radius for the simulation with no feedback, which can be explained by the already high clumping in central parts of the disc, making star formation very efficient even without H₂. Figure 24 shows the projected density of stars formed since the beginning of the simulation for α_fb = 1%: the disc of new stars is more extended if we include H₂. There is a few stellar clumps and also a very clear stellar bar that is maintained in the stellar component after a few Gyrs in both cases. Clumps of young stars can be seen on the density maps: they follow the gas clumps. These stellar clumps gradually lose energy and are eventually absorbed by the central bar.

We study further the star formation by drawing KS diagrams representing the surface density of SFR as a function of the gas surface density. Since we are limited in mass resolution for star formation (stellar particles of a fixed mass of ~10⁴ M_⊙ are created stochastically, which is completely different from some smooth star formation), we combine data corresponding to 50 snapshots from t = 200 Myr to t = 700 Myr in our study. We use a polar grid with a given number n_R of bins in cylindrical radius R and a given number n_θ of bins in azimuthal angle θ. Using this kind of grid allows for a more uniform signal/noise repartition than with an orthogonal grid and optimises the number of new stars per cell. In Fig. 25, we have plotted KS diagrams for simulations without feedback and with a feedback efficiency α_fb = 10%. These can be compared to Agertz et al. (2013)’s plots of azimuthally averaged KS diagrams of disc galaxies for different feedback intensities. Very similarly, we observe a global diminution of the SFR surface density when we increase the feedback strength. The figure quantifies how, on average, for the same gas surface densities, the SFR is lower with higher feedback. The feedback makes the gas more diffuse and destroys clumps. Two cells containing the same total amount of gas but different fractions of diffuse gas (cells are larger than the clumps size) have different star-formation efficiency. This explains the smaller scatter in the simulations with feedback: as the gas is more homogeneous, the relation between surface densities of SFR and gas is better determined. In Fig. 25, lines of constant gas depletion time are indicated. The gas depletion time is defined as . It can be seen that the high SFR and gas surface density regions of the galaxies have depletions times that are as low as a few hundred Myr in the simulations with no feedback, while the low SFR and density regions have depletion times up to 10 Gyr. The outer parts of the disc with a low surface density form stars much less efficiently than the central parts. We also see a difference for low surface densities between purely atomic and atomic + molecular simulations: the low surface density regions form stars more efficiently in atomic+molecular simulations. This is because H₂ cooling allows the gas to be locally denser and especially allows it to be more concentrated in the disc plane, as shown in the following section. Then the gas is denser in volume density and forms stars more efficiently for a given surface density.

Figure 25 also shows simulations with an alternative lower star-formation efficiency per free-fall time c_∗ = 0.01. In this case, the KS diagrams are shifted towards the right because the higher clumping allow densities to be higher and star formation is less efficient for a given volume density, which roughly translates in a given surface density. The scatter for α_fb = 10% is a little higher for c_∗ = 0.01 than c_∗ = 0.1, because the gas is clumpier for a lower c_∗ and the relation between volume and surface density departs there again from a one-to-one relation.

For simulations including H₂, the SFR is shown separately as a function of atomic and molecular components on the bottom of Fig. 25. The atomic hydrogen surface density is confined to low values for our galaxies, and the SFR-HI diagrams show a large scatter, because HI is too diffuse to track the star-forming gas. H₂ is determined as a better tracer of star formation in observations of nearby galaxies (e.g. Bigiel et al. 2008), and our results agree with a lot of theoretical work (e.g. Maio et al. 2010; Glover & Abel 2008; Greif et al. 2010; Wada et al. 2009; Christensen et al. 2012a; Gnedin et al. 2009).

3.2.4. Vertical structure of the disc

The inclusion of molecular hydrogen has a significant impact on the vertical structure of the disc, because the cold and dense gas concentrates in the disc middle plane and because of the impact of gas clumping on the distribution of feedback energy. Christensen et al. (2012a,b) include non-equilibrium formation of H₂, self-shielding, and dust shielding of both HI and H₂ in galaxies extracted from cosmological simulations and explore the influence of including H₂ formation for a fixed feedback efficiency. Like them, we find that the introduction of H₂ makes the outer parts of discs denser, while allowing the gas to go further away from the disc plane. The higher concentration in the disc plane, which is due to H₂ cooling, increases the SFR and thus the supernovae feedback. The feedback regulates star formation, but the SFR is still higher when H₂ can cool the gas and make it denser. The gas clumping also makes the feedback more efficient than in purely atomic simulations. In our feedback scheme, particles get velocity kicks weighted by the SPH kernel, so that the kicks are larger for particles closer to the new stellar particles. We have performed simulations with various feedback efficiencies, helping to check this effect. Figure 26 shows the fraction of gas that is further than 1 kpc from the disc (or radially further than 60 kpc from the centre of the galaxy): without feedback, the gas is only gravitationally heated and the effect of denser gas and a higher clumping factor due to H₂ makes the fraction of gas leaving the disc smaller than for a purely atomic hydrogen gas. With feedback, however, there is a higher fraction of gas outside the disc when H₂ is included. In our simulations, the effect is especially striking in the outer parts of discs when H₂ formation is considered, as the vertical restoring force is lower there.

We checked that the vertical density distribution is consistent with both a more concentrated disc and a higher fraction of gas leaving the disc. Figure 27 shows the vertical mass profile of the gas at large radii for α_fb = 40%: the gas is more concentrated in the disc plane, but the distribution has higher density “tails”. Figure 28 shows the characteristic height z_1/2 of the gas, which is the distance from the disc plane for which the density equals half of the central density, as a function of radius for the various simulations. Especially at large radii, the characteristic height is lower for simulations with H₂ for all feedbacks as H₂ concentrates the gas in the middle plane. The fraction of gas beyond this height however increases with the feedback efficiency at large radii when H₂ is considered, as gas is then efficiently sent away from the disc plane by feedback. Without feedback, some difference remains, probably due to gravitational heating produced by a higher clumping. We see indeed a slight difference in vertical velocity dispersion.

3.2.5. Gas density profile

The galaxies we have considered until now have a low gas surface density. We also run simulations with smaller characteristic radii, hence higher surface densities. Figure 29 displays the surface density for a Miyamoto-Nagai gas radius r_g of 3.6 kpc. The transition radius between H₂-dominated and HI-dominated regions, which is the radius at which Σ_H2 = Σ_HI, has then a value similar to the average observed by Bigiel & Blitz (2012; their average value for nearby spiral galaxies observations is 14 M_⊙ pc^-2). Previous galaxies belong to the lower surface density group observed by Bigiel & Blitz (2012). The effect of H₂ cooling is reduced in galaxies with higher surface density, since a higher fraction of the gas belongs to the central regions with much more efficient star formation than the outer parts. The cooling by H₂ is more important when there is more gas in the metal-poor outer regions.

4. Conclusions

To explore the influence of molecular hydrogen in the physics of spiral discs, star formation, and gas reservoirs in galaxy evolution, we have implemented detailed cooling by metals for temperatures as low as 100 K and cooling by H₂ due to collisions with H, He, and other H₂ molecules in Gadget-2. The determination of the H₂ density is inspired by the KMT recipe (Krumholz et al. 2008, 2009; McKee & Krumholz 2010) and uses the stellar UV flux from young stars. We study the influence of the cold, dense molecular phase on star formation. The simple method used to determine the H₂ mass fraction requires some calibration, when the resolution is changed. Using the instantaneous UV flux from young stars also means that some calibration must be done depending on the average SFR. We have also implemented a stochastic star formation and a kinetic supernovae feedback whose efficiency was varied in simulations with or without cooling by molecular hydrogen. The evolution of the ISM of a galaxy depends on the parameters chosen for star formation and feedback. While we have focused on a star-formation efficiency per free-fall time c_∗ = 0.1 and a minimum density for star formation n_Hmin = 10^-1 cm^-3, changing these parameters can give a more heterogeneous ISM, but the resolution required to attempt to resolve the Jeans mass is more important.

The influence of H₂ in the formation of dense gas and star formation is very important in the outer extended disc. Molecular hydrogen influences the vertical structure of the discs, especially when there is some stellar feedback: first, H₂ makes the gas more concentrated in the middle layer of the disc plane. Second, the gas is also more susceptible of being ejected far from the disc due to the higher efficiency of feedback in high density regions. Correlating SFR and gas surface density, it is found that molecular gas is a much better tracer of star formation than atomic gas, as what is observed in nearby galaxies. We find that including molecular hydrogen allows some slow star formation to occur in the low metallicity outer parts of galaxies. If gas is accreted by the discs, it may help store some cold gas with a slow star formation.

Acknowledgments

We thank Yves Revaz, Martin Stringer, Maxime Bois and Benjamin L’Huillier for stimulating discussions. Computations have been done partly on the GENCI/TGCC machines and partly on the cluster funded by the European Research Council under the Advanced Grant Program 267399-Momentum.

References

All Tables

All Figures