EDP Sciences
Free Access
Issue
A&A
Volume 559, November 2013
Article Number A55
Number of page(s) 19
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201220952
Published online 13 November 2013

© 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 1012 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 104 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:

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 H2 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 ~104 K, since the most efficient cooling processes are due to ionised hydrogen and helium, which are almost only present in the atomic form below 104 K. This usually implies the coldest dense gas is at an equilibrium temperature of about 104 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 H2 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 R25, 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 104 K

We use the cooling functions computed by Sutherland & Dopita (1993) for a plasma in collisional ionisation equilibrium above 104 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 104 K

Following Maio et al. (2007), we compute cooling functions of a few metals below 104 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 nH = 1 cm-3 and different metallicities Z = [Fe/H] are plotted in Fig. 1.

thumbnail Fig. 1

Volume cooling rates for nH = 1 cm-3, following Sutherland & Dopita (1993) for high temperatures (T > 104 K), and Maio et al. (2007) for low temperatures (T < 104 K), as a function of temperature and metallicity Z = [Fe/H] (Z = 0 is solar). In the code, we use exact density-dependent cooling functions for the low temperature part, cf Sect. 2.1.

Open with DEXTER

2.1.3. H2 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 H2. We assume H2 cooling below 104 K takes place in an almost neutral medium and thus consider collisions with other H2 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 H2, considering the collisions with H2 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 H2 molecule in erg s-1, and ΛH2,H, ΛH2,H2, and ΛH2,He, in erg cm3 s-1, account for collisions with respectively H, H2 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 nHnucl = 1 cm-3 (such that nHnucl = nH + 2nH2) and different H2 mass fractions fH2, with . We do not consider cooling due to collisions of H2 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 H2 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).

thumbnail Fig. 2

Volume cooling rate defined in Eq. (1) for a hydrogen-nuclei number density of 1 cm-3 as a function of temperature for different H2 to total hydrogen mass ratios fH2. For high fractions of H2 and at such densities, the cooling can become slightly more important just below 104 K when the H2 mass fraction is reduced, because the contribution of collisions with atomic hydrogen dominates in this range of temperature.

Open with DEXTER

thumbnail Fig. 3

Volume cooling rates due to H2 and metals for a fixed hydrogen-nuclei number density of 1 cm-3 and different metallicities Z = [Fe/H] (Z = 0 is solar) and H2 mass fractions.

Open with DEXTER

In Fig. 3 we show the volume cooling rate due to H2 and metals for a fixed hydrogen-nuclei number density nHnucl = 1 cm-3. We show the volume cooling rates for different metallicities and H2 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)2 (the scaling being such for these low densities of hydrogen atoms). We see that H2 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 H2, 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 H2. 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 nHΛHD,H(nH = 1 cm-3) for low densities.

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

In the simulations, we assume that nHD = 10-5 nH2 because the fractionation of HD with respect to H2 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 nHD = 10-5 nH2, 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 cm3 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 xi the number fraction of the species i and mp 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 104 K. Above that temperature, we assume there is no molecular hydrogen and the gas is fully ionised (H+ and He2+). 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 104 K by assuming a neutral gas while it is below 104 K by assuming an ionised gas. We set the temperature of the gas to 104 K in this specific entropy range.

2.2. Resolution

Table 1

Resolution: Gravitational softening and particle masses.

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 Ng→, 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 104 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. H2 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 fH2 in an idealised spherical cloud submitted to a uniform and isotropic radiation field. Including radiative transfer and the formation and destruction of H2 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 fH2 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 cm2 (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 = ∑ jmj(ρj − ρi)∇iW(rij,hi). 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 H2 molecules. We do not consider that older stars radiate in this band, so they do not have any effect on the fraction of H2 in our simulations. We take with a scaling factor that depends on the resolution of the simulations and tune it so that the obtained H2 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, tff 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 mg 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 Ng→ 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 Ng→ is greater than one, we note that this spawning scheme implies that gas particles have different masses, of either the initial one mg, 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 hi is fixed with hi 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 nHmin, have a temperature lower than a maximum temperature Tmax, 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 ESN = 1051 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 Einput = αfbϵSNm where ϵSN is the number of formed supernovae per formed stellar mass multiplied by the canonical energy ESN = 1051 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 Einput. 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 Ng→ 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 Ng→.

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.

Table 2

Sb galaxy parameters.

thumbnail Fig. 4

Initial rotation curve. Top solid (blue) thick line: total rotation curve. Dashed (purple) line: bulge. Dotted (green) line: stellar disc. Dash-dot (green) line: gas. Solid (red) line: dark matter halo.

Open with DEXTER

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.

thumbnail Fig. 5

Initial surface densities after 300 Myr of initial conditions preparation. Solid blue line: stellar disc. Dashed green line: gas disc. The black lines (solid for stars and dashed for gas) represent the analytic profiles taken at the beginning of the initial conditions preparation.

Open with DEXTER

thumbnail Fig. 6

Initial density map of the gas disc after 300 Myr of initial conditions preparation. The box size is [80 kpc × 80 kpc].

Open with DEXTER

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 Ng→ = 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 nHmin = 10-1 cm-3. The number density nH is hereafter the total hydrogen number density that was written as nHnucl in Sect. 2. We use a maximal temperature for star formation Tmax = 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 H2

thumbnail Fig. 7

Projections of the gas density after 0.5 Gyr, 1 Gyr, and 3 Gyr of evolution. Box sizes are [60 kpc × 60 kpc] for face-on views and [60 kpc × 10 kpc] for edge-on views. Each row corresponds to a feedback efficiency indicated on the left. The column integrated density scale is the same for all plots. This is done with the visualisation software SPLASH (Price 2007).

Open with DEXTER

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 ≃ 104 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 104 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.

thumbnail Fig. 8

Temperature-hydrogen number density histogram of the gas after 0.5 Gyr of evolution. The 2D histograms and marginal 1D histograms are all mass weighted and normalised. From top to bottom: runs with an increasing feedback efficiency.

Open with DEXTER

thumbnail Fig. 9

Time evolution of the mass of the gas component for different feedback efficiencies.

Open with DEXTER

3.2. Simulations with H2

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

3.2.1. H2 fraction

We first set the feedback efficiency to αfb = 1% and tune the χ factor in the H2 mass fraction fH2 so that we find a global H2 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 H2 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 H2 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 H2 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).

thumbnail Fig. 10

Global H2 mass fraction versus time for αfb = 1%. The UV flux increases from top to bottom.

Open with DEXTER

thumbnail Fig. 11

Radial distribution of the surface density of H2, atomic gas HI, and total hydrogen gas after 0.5 Gyr for αfb = 1% and the selected UV scaling factor.

Open with DEXTER

We plot the mass fraction of H2 versus the hydrogen-nuclei number density at t = 0.5 Gyr in Fig. 12. The H2 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 H2 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.

thumbnail Fig. 12

Mass fraction of H2 versus hydrogen-nuclei number density after 0.5 Gyr for αfb = 1%. The mass weighted 2D histogram is in grey scale, while the mass-weighted 1D histogram is represented as the blue solid line.

Open with DEXTER

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 Nref = 1 200 000 to simulations with , , Nref × 2, and Nref × 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χ2) entering the computation of the molecular mass fraction versus density.

thumbnail Fig. 13

Different resolutions with αfb = 10% after 0.5 Gyr of evolution. Top: H2 mass fraction as a function of hydrogen-nuclei number density. Middle: dust optical depth as a function of hydrogen-nuclei number density. Bottom: factor involving the UV flux entering the computation of the H2 mass fraction (see Eqs. (78)) as a function of hydrogen-nuclei number density.

Open with DEXTER

3.2.2. Gas physical state

thumbnail Fig. 14

Projections of the gas density (top row) and H2 density (bottom row) after 0.5 Gyr, 1 Gyr, and 3 Gyr of evolution. Box sizes are [60 kpc × 60 kpc] for face-on views and [60 kpc × 10 kpc] for edge-on views. The column integrated density scale is the same for all plots. This is done with the visualisation software SPLASH (Price 2007).

Open with DEXTER

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 H2. Cooling by collisions of atomic hydrogen with metals or by H2 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 H2 allows the gas to be clumpier: we see density features that were absent in the purely atomic simulation. The surface density of H2 is also plotted. It indeed follows the density features of the gas: the clumps and filamentary structures can be seen in H2. 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.

thumbnail Fig. 15

Gas surface density power spectra after 0.5 Gyr of evolution for αfb = 1%.

Open with DEXTER

thumbnail Fig. 16

Clumping factor (defined in Eq. (18)) evolution with time αfb = 1%.

Open with DEXTER

thumbnail Fig. 17

Temperature-hydrogen number density histogram of the gas after 0.5 Gyr of evolution. The 2D histograms and marginal 1D histograms are all mass weighted and normalised. From top to bottom: runs with an increasing feedback efficiency.

Open with DEXTER

thumbnail Fig. 18

Fraction of cold gas as a function of time for two different feedback efficiencies αfb = 1% and αfb = 10%. Solid lines: with H2. Dashed lines: without H2.

Open with DEXTER

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 H2 cooling is considered and stars are formed from colder and denser gas, increasing the star-formation efficiency.

thumbnail Fig. 19

PDFs of hydrogen number density as a function of galactocentric radius after 0.5 Gyr of evolution with or without including H2 cooling in the simulations for αfb = 1%. The mass weight per pixel of density and radius increases with darkness. Horizontal black line: threshold density for star formation.

Open with DEXTER

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 H2. 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 nHmin = 10-1 cm-3 and simulations with the same c = 0.1 but a higher nHmin = 101cm-3. Reducing c or increasing nHmin both allow for higher density peaks, as seen on Fig. 20. The effect of including H2 is also visible in these sets of simulations, but the H2 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 nHmin is heightened, the gas forms stars more difficultly and can thus collapse more efficiently, resulting in a more fractioned gas.

thumbnail Fig. 20

PDFs of hydrogen number density as a function of galactocentric radius after 0.5 Gyr of evolution with or without including H2 cooling in the simulations for αfb = 1% and alternative star-formation parameters. Left: runs with nHmin = 10-1 cm-3 and c = 0.01. Right: runs with nHmin = 101 cm-3 and c = 0.1. Horizontal black line: threshold density for star formation.

Open with DEXTER

thumbnail Fig. 21

Top: evolution of the gaseous mass of the galaxy. Bottom: SFR evolution. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER

3.2.3. Star formation

The effect of H2 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.

thumbnail Fig. 22

SFR evolution for αfb = 1% and alternative star-formation parameters. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER

thumbnail Fig. 23

Cumulative SFR as a function of galactocentric radius, averaged on the first Gyr. Solid lines: run with H2. Dashed lines: run without H2.

Open with DEXTER

thumbnail Fig. 24

Projection of the density of stars formed during the simulations after 0.5 Gyr, 1 Gyr and 3 Gyr of evolution. Top row: feedback efficiency αfb = 1% and no H2. Bottom row: feedback efficiency αfb = 1% and H2. Box sizes are [80 kpc × 80 kpc] for face-on views and [80 kpc × 10 kpc] for edge-on views.

Open with DEXTER

thumbnail Fig. 25

KS diagrams. Two left columns: no feedback. Two right columns: αfb = 10%. Star-formation efficiencies per free fall time c are indicated on the top. Top row: simulations without H2. Second row: simulations with H2. Two bottom rows: SFR versus the HI only and H2 only surface densities for simulations with H2.

Open with DEXTER

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 H2. 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 H2. 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 ~104 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 nR 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 H2 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 H2, 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. H2 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 H2, self-shielding, and dust shielding of both HI and H2 in galaxies extracted from cosmological simulations and explore the influence of including H2 formation for a fixed feedback efficiency. Like them, we find that the introduction of H2 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 H2 cooling, increases the SFR and thus the supernovae feedback. The feedback regulates star formation, but the SFR is still higher when H2 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 H2 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 H2 is included. In our simulations, the effect is especially striking in the outer parts of discs when H2 formation is considered, as the vertical restoring force is lower there.

thumbnail Fig. 26

Fraction of gas beyond 1 kpc from the disc plane. Solid lines are runs with H2 and dashed ones are runs without H2.

Open with DEXTER

thumbnail Fig. 27

Vertical mass profile of the disc for radii R > 15 kpc for runs with αfb = 40%. Solid line: run with H2. Dashed line: run without H2. The two vertical lines mark the characteristic vertical heights z1/2.

Open with DEXTER

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 z1/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 H2 for all feedbacks as H2 concentrates the gas in the middle plane. The fraction of gas beyond this height however increases with the feedback efficiency at large radii when H2 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 rg of 3.6 kpc. The transition radius between H2-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 H2 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 H2 is more important when there is more gas in the metal-poor outer regions.

thumbnail Fig. 28

Vertical structure of the gas. First row: characteristic gas height z1/2 as a function of galactocentric radius after 0.5 Gyr of evolution. Second row: fraction of gas beyond the characteristic height. Third row: vertical velocity dispersion. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER

thumbnail Fig. 29

Radial distribution of the surface density of H2, atomic gas HI and total hydrogen gas after 100 Myr for αfb = 40% and a χ × 500 UV scaling factor.

Open with DEXTER

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 H2 due to collisions with H, He, and other H2 molecules in Gadget-2. The determination of the H2 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 H2 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 nHmin = 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 H2 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, H2 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

Table 1

Resolution: Gravitational softening and particle masses.

Table 2

Sb galaxy parameters.

All Figures

thumbnail Fig. 1

Volume cooling rates for nH = 1 cm-3, following Sutherland & Dopita (1993) for high temperatures (T > 104 K), and Maio et al. (2007) for low temperatures (T < 104 K), as a function of temperature and metallicity Z = [Fe/H] (Z = 0 is solar). In the code, we use exact density-dependent cooling functions for the low temperature part, cf Sect. 2.1.

Open with DEXTER
In the text
thumbnail Fig. 2

Volume cooling rate defined in Eq. (1) for a hydrogen-nuclei number density of 1 cm-3 as a function of temperature for different H2 to total hydrogen mass ratios fH2. For high fractions of H2 and at such densities, the cooling can become slightly more important just below 104 K when the H2 mass fraction is reduced, because the contribution of collisions with atomic hydrogen dominates in this range of temperature.

Open with DEXTER
In the text
thumbnail Fig. 3

Volume cooling rates due to H2 and metals for a fixed hydrogen-nuclei number density of 1 cm-3 and different metallicities Z = [Fe/H] (Z = 0 is solar) and H2 mass fractions.

Open with DEXTER
In the text
thumbnail Fig. 4

Initial rotation curve. Top solid (blue) thick line: total rotation curve. Dashed (purple) line: bulge. Dotted (green) line: stellar disc. Dash-dot (green) line: gas. Solid (red) line: dark matter halo.

Open with DEXTER
In the text
thumbnail Fig. 5

Initial surface densities after 300 Myr of initial conditions preparation. Solid blue line: stellar disc. Dashed green line: gas disc. The black lines (solid for stars and dashed for gas) represent the analytic profiles taken at the beginning of the initial conditions preparation.

Open with DEXTER
In the text
thumbnail Fig. 6

Initial density map of the gas disc after 300 Myr of initial conditions preparation. The box size is [80 kpc × 80 kpc].

Open with DEXTER
In the text
thumbnail Fig. 7

Projections of the gas density after 0.5 Gyr, 1 Gyr, and 3 Gyr of evolution. Box sizes are [60 kpc × 60 kpc] for face-on views and [60 kpc × 10 kpc] for edge-on views. Each row corresponds to a feedback efficiency indicated on the left. The column integrated density scale is the same for all plots. This is done with the visualisation software SPLASH (Price 2007).

Open with DEXTER
In the text
thumbnail Fig. 8

Temperature-hydrogen number density histogram of the gas after 0.5 Gyr of evolution. The 2D histograms and marginal 1D histograms are all mass weighted and normalised. From top to bottom: runs with an increasing feedback efficiency.

Open with DEXTER
In the text
thumbnail Fig. 9

Time evolution of the mass of the gas component for different feedback efficiencies.

Open with DEXTER
In the text
thumbnail Fig. 10

Global H2 mass fraction versus time for αfb = 1%. The UV flux increases from top to bottom.

Open with DEXTER
In the text
thumbnail Fig. 11

Radial distribution of the surface density of H2, atomic gas HI, and total hydrogen gas after 0.5 Gyr for αfb = 1% and the selected UV scaling factor.

Open with DEXTER
In the text
thumbnail Fig. 12

Mass fraction of H2 versus hydrogen-nuclei number density after 0.5 Gyr for αfb = 1%. The mass weighted 2D histogram is in grey scale, while the mass-weighted 1D histogram is represented as the blue solid line.

Open with DEXTER
In the text
thumbnail Fig. 13

Different resolutions with αfb = 10% after 0.5 Gyr of evolution. Top: H2 mass fraction as a function of hydrogen-nuclei number density. Middle: dust optical depth as a function of hydrogen-nuclei number density. Bottom: factor involving the UV flux entering the computation of the H2 mass fraction (see Eqs. (78)) as a function of hydrogen-nuclei number density.

Open with DEXTER
In the text
thumbnail Fig. 14

Projections of the gas density (top row) and H2 density (bottom row) after 0.5 Gyr, 1 Gyr, and 3 Gyr of evolution. Box sizes are [60 kpc × 60 kpc] for face-on views and [60 kpc × 10 kpc] for edge-on views. The column integrated density scale is the same for all plots. This is done with the visualisation software SPLASH (Price 2007).

Open with DEXTER
In the text
thumbnail Fig. 15

Gas surface density power spectra after 0.5 Gyr of evolution for αfb = 1%.

Open with DEXTER
In the text
thumbnail Fig. 16

Clumping factor (defined in Eq. (18)) evolution with time αfb = 1%.

Open with DEXTER
In the text
thumbnail Fig. 17

Temperature-hydrogen number density histogram of the gas after 0.5 Gyr of evolution. The 2D histograms and marginal 1D histograms are all mass weighted and normalised. From top to bottom: runs with an increasing feedback efficiency.

Open with DEXTER
In the text
thumbnail Fig. 18

Fraction of cold gas as a function of time for two different feedback efficiencies αfb = 1% and αfb = 10%. Solid lines: with H2. Dashed lines: without H2.

Open with DEXTER
In the text
thumbnail Fig. 19

PDFs of hydrogen number density as a function of galactocentric radius after 0.5 Gyr of evolution with or without including H2 cooling in the simulations for αfb = 1%. The mass weight per pixel of density and radius increases with darkness. Horizontal black line: threshold density for star formation.

Open with DEXTER
In the text
thumbnail Fig. 20

PDFs of hydrogen number density as a function of galactocentric radius after 0.5 Gyr of evolution with or without including H2 cooling in the simulations for αfb = 1% and alternative star-formation parameters. Left: runs with nHmin = 10-1 cm-3 and c = 0.01. Right: runs with nHmin = 101 cm-3 and c = 0.1. Horizontal black line: threshold density for star formation.

Open with DEXTER
In the text
thumbnail Fig. 21

Top: evolution of the gaseous mass of the galaxy. Bottom: SFR evolution. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER
In the text
thumbnail Fig. 22

SFR evolution for αfb = 1% and alternative star-formation parameters. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER
In the text
thumbnail Fig. 23

Cumulative SFR as a function of galactocentric radius, averaged on the first Gyr. Solid lines: run with H2. Dashed lines: run without H2.

Open with DEXTER
In the text
thumbnail Fig. 24

Projection of the density of stars formed during the simulations after 0.5 Gyr, 1 Gyr and 3 Gyr of evolution. Top row: feedback efficiency αfb = 1% and no H2. Bottom row: feedback efficiency αfb = 1% and H2. Box sizes are [80 kpc × 80 kpc] for face-on views and [80 kpc × 10 kpc] for edge-on views.

Open with DEXTER
In the text
thumbnail Fig. 25

KS diagrams. Two left columns: no feedback. Two right columns: αfb = 10%. Star-formation efficiencies per free fall time c are indicated on the top. Top row: simulations without H2. Second row: simulations with H2. Two bottom rows: SFR versus the HI only and H2 only surface densities for simulations with H2.

Open with DEXTER
In the text
thumbnail Fig. 26

Fraction of gas beyond 1 kpc from the disc plane. Solid lines are runs with H2 and dashed ones are runs without H2.

Open with DEXTER
In the text
thumbnail Fig. 27

Vertical mass profile of the disc for radii R > 15 kpc for runs with αfb = 40%. Solid line: run with H2. Dashed line: run without H2. The two vertical lines mark the characteristic vertical heights z1/2.

Open with DEXTER
In the text
thumbnail Fig. 28

Vertical structure of the gas. First row: characteristic gas height z1/2 as a function of galactocentric radius after 0.5 Gyr of evolution. Second row: fraction of gas beyond the characteristic height. Third row: vertical velocity dispersion. Solid lines: runs with H2. Dashed lines: runs without H2.

Open with DEXTER
In the text
thumbnail Fig. 29

Radial distribution of the surface density of H2, atomic gas HI and total hydrogen gas after 100 Myr for αfb = 40% and a χ × 500 UV scaling factor.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.