© The Authors 2024

1. Introduction

Despite the substantial progress achieved in the last decades, both observationally and theoretically driven, a complete understanding of the various processes involved in galaxy formation has yet to be reached. Among them, modelling the pivotal phenomenon of star formation (SF) faces formidable challenges arising from the intricate interplay of complex physical mechanisms acting across a wide range of scales, from the large-scale structure (of the order of megaparsec) to the sub-stellar ones (≲10⁻⁸ pc). The ultimate reservoir for SF is the densest and coldest phase of the interstellar medium (ISM), dominated in mass by molecular hydrogen (H₂) and organized in self-gravitating dusty molecular clouds (MC) featuring an approximate mass range from about 10⁴ to 10⁶ M_⊙ (e.g. Kennicutt & Evans 2012; Miville-Deschênes et al. 2017). The gravitational collapse of these clouds initiates the cooling and SF processes within them. Indeed, while the observed star formation rate (SFR) surface density of galaxies correlates with the neutral gas surface density (the Schimdt-Kennicutt law e.g. Schmidt 1959; Kennicutt 1998), the correlation becomes, unsurprisingly, much stronger when only the molecular gas surface density is considered (e.g. Wong & Blitz 2002; Bigiel et al. 2008). While the most widespread interpretation of the latter correlation is that the formation of molecular gas controls the SFR, it has also been proposed that both of these processes, the MC formation and the onset of SF, might share a common underlying cause (e.g. Krumholz et al. 2011; Mac Low & Glover 2012).

Numerical simulations of galaxy formation in a cosmological context are far from the resolution required to resolve the structure of MC, the minimal requirement to model SF from first principles. Moreover, they do not incorporate many of the physical mechanisms that determine the MC evolution to form stars. Instead, they rely on various proposed sub-resolution prescriptions that connect the SFR with the estimated availability of gas eligible for SF. In most cases, this simply means considering gas that, at a resolved level, turns out to be cold and dense enough (T ≲ a few × 10⁴ K; n ≳ 0.1–10 cm⁻³), and that the local SFR is assumed to be proportional to its density divided by an SF timescale usually proportional to the local free-fall time. The Schimdt-Kennicutt law loosely inspires the latter relationship. This kind of approach is usually adopted when simulating large cosmological volumes. Examples include Illustris (Vogelsberger et al. 2014), Magneticum (Dolag 2015), NIHAO (Wang et al. 2015), IllustrisTNG (Pillepich et al. 2018), the Massive-Black (Khandai et al. 2015) and Fable (Henden et al. 2018) simulations, the Horizon-AGN suite (Dubois et al. 2016, 2021) and the DIANOGA simulations (Ragone-Figueroa et al. 2018; Bassini et al. 2020).

Somewhat more elaborated treatments compute the local SFR by including a sub-resolution prescription to estimate the fraction of cold and dense gas that is in molecular form. The prescription can be derived by theoretical modelling (e.g. Krumholz et al. 2008, 2009; Krumholz & Gnedin 2011), as in Kuhlen et al. (2012), the Mufasa and Simba simulations (Davé et al. 2017, 2019), and the FIRE suite (Hopkins et al. 2014, 2018), or observations (e.g. Blitz & Rosolowsky 2006), as in Murante et al. (2015)¹.

At a higher level of complexity, some simulations, thus far dealing with smaller volumes and/or single halos, predict the amount of molecular gas following the main processes of H₂ formation and destruction and take advantage of this prediction to compute the SFR (e.g. Pelupessy et al. 2006; Gnedin et al. 2009; Christensen et al. 2012; Tomassetti et al. 2015). The dominant H₂ formation mechanisms to consider in metal-enriched environments are catalytic reactions occurring on the surfaces of dust grains (for a review see Wakelam et al. 2017). Direct formation in the gas phase does indeed turn out to be very inefficient due to the difficulty of the axis-symmetric-forming molecule in radiating away the excess energy to remain stable. This process becomes relevant only when the ISM is extremely dust-poor. (e.g. Gould & Salpeter 1963; Krumholz et al. 2009; Bekki 2013; Bron et al. 2014; Tomassetti et al. 2015; Valdivia et al. 2016; Maio et al. 2022). However, the aforementioned simulation works lack comprehensive dust creation and destruction modelling. Instead, they rely on a fixed dust-to-gas (D/G) ratio linearly scaling with metallicity to estimate the dust content roughly. A further complication arises from the dust-catalyzed H₂ formation rate being proportional to the grain’s surface rather than its mass. Therefore, it would be relevant to have information on the dust size distribution to estimate it. The primary H₂ destruction mechanism is believed to be photo-dissociation via the two-step Solomon process caused by photons in the Lyman–Werner (LW) band 912 Å < λ < 1108 Å. It is worth pointing out that dust also substantially contributes to shield regions where H₂ forms efficiently on grains from the destructive effect of LW radiation emitted by OB stars.

Clearly, dust plays a pivotal role in determining the H₂ content of the ISM and, consequently, the reservoir available for SF. However, despite its significance, current cosmological simulations have yet to integrate the simultaneous evolution of dust content and H₂ promoted by dust. Prior studies in this field, including H₂ formation on grain surfaces, have adopted a simplistic method, merely scaling the dust content with gas metallicity while overlooking the impact of grain size distribution on dust surfaces. Notably, Romano et al. (2022) have conducted a simulation which accounts for both dust content and H₂ concurrently. However, it focuses solely on one isolated disc galaxy and spans a relatively short period of 2 billion years of evolution, during which the SFR is not linked to the molecular content.

In this work, we present an update of the MUPPI (MUlti Phase Particle Integrator) sub-resolution model (Murante et al. 2015), in which the evolution of dust content, molecular (H₂) gas and SF are mutually linked. Previous versions of MUPPI featured an advanced description of a multi-phase ISM and proved successful in zoomed-in simulations of late-type galaxies (e.g. Valentini et al. 2020; Granato et al. 2021; Valentini et al. 2023, and references therein). The molecular gas content, and thus the SFR, was estimated adopting the Blitz & Rosolowsky (2006) or Krumholz et al. (2009) laws. Here, we take instead advantage of the treatment of dust evolution introduced by Granato et al. (2021) and Parente et al. (2022) to predict the evolution of H₂, which in turn is used to estimate the local SFR.

In Granato et al. (2021), we included in MUPPI a treatment of dust formation and evolution and studied the dust properties of a Milky Way (MW)-like galaxy. This model allows us to predict dust grains’ chemical composition and size distribution, utilizing the Hirashita (2015) two-size approximation. Furthermore, we have accounted for hot gas cooling due to collisions of ions with dust particles, following the Dwek & Werner (1981) model. More recently, in Parente et al. (2022), we used MUPPI and its treatment of dust evolution for the first time to simulate a cosmological box and test its performance by comparing the simulated galaxy properties with observational constraints concerning a broader population of galaxies. However, in that work, we deliberately adopted the MUPPI modelling, as previously calibrated, by using zoom-in simulations of a MW-like galaxy without modifications. Substantial tensions regarding low dust content at high redshift and low star formation rate density (SFRD) were reported in that work. The dust evolution modelling, in that case, originated only minor modifications to the behaviour of MUPPI related to its effects on cooling. Conversely, introducing a treatment of molecular gas formation promoted by dust and its use as a reservoir for SF makes the treatment of dust deeply intertwined with the core of MUPPI. Therefore, in the present work, we have revised some aspects of its modelling and parameters, obtaining a simulation set that, at the same time, essentially solves the tensions pointed out by Parente et al. (2022) in their cosmological volumes (Section 4.1) while preserving the good results discussed in previous papers on zoom-in simulations of MW-like galaxies.

The structure of the paper unfolds as follows: In Section 2, we provide a concise overview of the numerical simulations, the sub-resolution model MUPPI, and the adjustments implemented in this iteration of the code to simulate a cosmological box reflecting galaxy properties consistent with observations. Section 3 describes the modelling for the formation of H₂ on dust grain surfaces and the approximation employed to simulate the molecule destruction resulting from LW radiation. Our findings, showcasing a satisfactory replication of observed properties through our model galaxies, including the H₂ mass function and the cosmic evolution of the H₂ mass density, are presented in Section 4. In Sections 5, 6 and 7, we delve into the analysis of three applications of our model: the integrated molecular Kennicutt-Schmidt (mKS) relation, the determination of H₂ mass from dust, and the transition from atomic to molecular gas. Finally, our conclusions are presented in Section 8.

2. Numerical simulations

The sub-resolution model that we employ to account for the unresolved SF, stellar feedback processes and treatment of AGN feedback is MUPPI (MUlti Phase Particle Integrator, Murante et al. 2010, 2015; Valentini et al. 2017, 2019, 2020, 2023), with some modifications detailed in Section 2.2. MUPPI is implemented in the GADGET3 code, a customized version of GADGET2 (Springel 2005). Our version includes the improved SPH (smoothed particle hydrodynamics) formulation proposed by Beck et al. (2016). MUPPI assumes that every SPH particle capable of performing SF embodies a multi-phase (MP) ISM, characterized by a hot, a cold, and a virtual stellar phase. The mass and energy flows among the different components, as depicted in Fig. 1, are governed by a set of ordinary differential equations. A gas particle becomes MP when its number density exceeds 0.005 cm⁻³ and its temperature falls below 5 × 10⁴ K. Entering MP does not imply the immediate onset of SF. Instead, it indicates that the gas particle is considered to model the ISM. Once in MP, the hot phase coexists in pressure equilibrium with the cold phase, assumed to maintain a constant temperature of T_c = 300 K.

Fig. 1. Schematic view of the mass fluxes in a multi-phase gas particle. Gas cooling moves mass from the hot to the cold phase at a rate Ṁcool. Within the cold phase, an H2 reservoir begins to grow. ${\dot{M}}_{\mathrm{form}{\mathrm{H}}_2}$ is the H2 formation rate determined by the particle dust content and the number density of atomic hydrogen, which in turn is computed from the cold phase mass. ${\dot{M}}_{\mathrm{dest}{\mathrm{H}}_2}$ is the H2 destruction rate by LW photons, which transports mass from the molecular reservoir back to the cold phase. Ṁsf is the rate at which molecular gas creates a virtual stellar phase, which is contained within the gas particle until a stellar particle is stochastically spawned. Further destruction of the molecular mass occurs due to the action of massive stars and AGN, at rates Ṁev and ${\mathrm{f}}_{\mathrm{H}2}{\dot{M}}_{\mathrm{AGN}}$, respectively, where fH2 represents the fraction of H2 in the cold phase. AGN also moves material from the cold to the hot phase at a rate ṀAGN

2.1. Initial conditions and resolution

To ensure robustness in our results and a sampling of galaxy population sufficient for our purposes, we simulate five cosmological boxes measuring 26 cMpc per side, from distinct initial conditions². We checked that the z = 0 halo mass function is satisfactorily sampled up to 2 × 10¹³ M_⊙ by our total simulated volume. The cosmological parameters are taken from Planck Collaboration XIII (2016): Ω_m = Ω_DM + Ω_b = 0.3089, Ω_b = 0.0486 and Ω_Λ = 0.6911, h = 0.6774. The power spectrum has a primordial index n_s = 0.9667 and normalization σ₈ = 0.8159.

We evolve 256³ DM particles and 256³ gas particles. The mass resolution of DM particles is ≃3.75 × 10⁷ M_⊙, while gas particles have an initial mass ≃7 × 10⁶ M_⊙. For the computation of the gravitational force, we use a Plummer-equivalent softening length of 960 pc, constant in comoving units for z > 6 and constant in physical units at lower redshift.

Our simulations achieve a resolution similar to that employed in the largest boxes of state-of-the-art suites such as Eagle (Schaye et al. 2015), IllustrisTNG (Pillepich et al. 2018), and Simba (Davé et al. 2019), which is relatively modest. In this work, we compare our results with globally integrated properties of galaxies derived from masses in the various components and overall sizes. However, Murante et al. (2015) showed that the advanced sub-resolution model MUPPI allows to capture to some level the structure and dynamics of MW-size galaxies, even at the resolution of their AqC6 and GA1 initial conditions (readers can refer to their Table 1 and Figures 17 and 18), which is very similar to that adopted in this work.

2.2. Changes with respect to previous MUPPI models

In this work, we introduced a series of modifications to the MUPPI setup described by Valentini et al. (2023), which was calibrated using zoom-in simulations of a MW-type galaxy. We remark that the latter work and most recent papers using MUPPI are based on simulations with a mass resolution higher by more than one order of magnitude than that in our cosmological boxes. Resolution dependence of the ‘best’ MUPPI setup has to be expected at some level, as in any sub-resolution modelling, meaning that the changes discussed below may be at least in part demanded by the differences in resolutions. Together with our treatment of the molecular mass fuel for SF, these modifications significantly improved the agreement with the basic properties of the galaxy population in the cosmological boxes, which are presented in Section 4. These adjustments essentially preserve the key achievements noted in earlier zoom-in simulations. Figure 2 presents face-on and edge-on maps of the various baryonic components of two spiral galaxies extracted from the simulations. One galaxy has a baryonic mass of approximately 6 × 10¹⁰ M_⊙, comparable to that of the Milky Way (e.g. McMillan 2017), while the other is about twice as massive.

Fig. 2. Two examples of disc galaxies from our simulations at z = 0. The four rows depict the mass surface densities of stars, gas, dust, and H2 on pixels of ∼0.5 kpc per side. Columns 1 and 2 illustrate face-on and edge-on views of a galaxy with M200 = 1.5 × 1012 M⊙, Mstar = 5.3 × 1010 M⊙, Mgas = 5.0 × 109 M⊙, Mdust = 7.0 × 107 M⊙, and MH2 = 2.2 × 109 M⊙. Columns 3 and 4 show the corresponding maps for a galaxy with M200 = 3.6 × 1012 M⊙, Mstar = 1.2 × 1011 M⊙, Mgas = 4.6 × 109 M⊙, Mdust = 6.3 × 107 M⊙, and MH2 = 2.0 × 109 M⊙.

The mentioned modifications to MUPPI are listed below, and Table 1 provides a summary of them.

• Among the different IMF explored in Valentini et al. (2019), we adopt that by Chabrier (2003), which is likely the most popular when dealing with observations.

• Following Valentini et al. (2019) (readers can refer to their Table 1 and Section 3), we adopt here a lower value for the kinetic wind efficiency f_FB, kin = 0.08, more appropriate for an IMF enriched with massive stars.

• We introduce a dependence of the lifetime of particles in the wind, t_wind, on the velocity dispersion σ_DM of neighbouring DM particles. Wind particles are designed to sample galactic outflows. A gas particle exits its MP stage after a maximum allowed time given by the free-fall time of the cold gas. When a gas particle exits a MP stage, it has a probability P_kin (a model parameter) of being kicked and becoming a wind particle for a time interval t_wind. Wind particles are decoupled from hydro forces during the aforementioned interval t_wind. Still, they are affected by radiative cooling and receive kinetic energy from neighbouring star-forming gas particles. The latter is used to increase their velocity. For further details, readers can refer to Appendix A and references therein. The dependence of t_wind on σ_DM introduced in the present work reads as follows:

  $$\begin{array}{c}\hfill t_{\mathrm{wind}}={\left(\frac{{\sigma }_0}{{\sigma }_{\mathrm{DM}}}\right)}^4\frac{t_{\mathrm{ff},\mathrm{c}}}{2},\end{array}$$(1)

  where ${\sigma }_0=50{\mathrm{km}\mathrm{s}}^{-1}$ and t_ff, c is the free-fall time of the sub-grid cold phase. For particles with σ_DM greater (smaller) than σ₀, their time spent in the wind is shortened (prolonged) compared to half their free-fall time. As a result, in less massive haloes, where the typical σ_DM is smaller, wind particles require more time to return and become eligible for MP once again. This adjustment is crucial to prevent the formation of low-mass galaxies. The value of σ₀ was adjusted to reproduce the Stellar-Halo mass relation at z = 0. In Appendix A, we illustrate how this choice impacts the SFRD, the stellar-halo mass relation (SHMR), and the stellar mass function (SMF) (Fig. A.1).

• We adopt the low metallicity feedback (LMF) introduced in Valentini et al. (2023) to mimic feedback by an early population of stars. We assume that star-forming particles with low metallicity contribute to the ISM with kinetic energy with an efficiency larger by a factor of 20 with respect to higher metallicity star-forming particles. In Valentini et al. (2023), this factor was 10. The impact of this choice is shown in Appendix A (Fig. A.2).

• The density threshold for entering MP is halved with respect to Valentini et al. (2023), so that now it is set to 0.005 ${\mathrm{cm}}^{-3}$. This adjustment aids in shifting the peak of the SFRD towards lower redshifts, aligning it with observed determinations. The SF efficiency is also halved and set at f_* = 0.01.

• The black hole (BH) radiative and feedback efficiencies are set to 0.1. The halo mass for BH seeding is $1.3\times {10}^{11}{\mathrm{M}}_{\odot }$, with the seeded BH having a mass of $1.2\times {10}^5{\mathrm{M}}_{\odot }$. Our simulations do not incorporate radio-mode feedback. As for the pinning of BHs at the centre of the hosting galaxy, our model uses the method detailed in Ragone-Figueroa et al. (2018), which avoids non-physical behaviours that can affect less careful BH repositioning (see also Damiano et al. 2024).

• As for the AGN feedback energy partition between the hot and cold phases of MP particles, we abandon the fiducial solution by Valentini et al. (2020), based on an evolving, thought parameterized estimate of the covering factor of cold clouds. We adopt instead their simpler prescription in which the energy is evenly distributed between the two phases, dubbed both in that paper. The latter prescription yielded qualitatively similar outcomes in their zoom-in MW-like galaxy simulations. However, it lowers the SFRD at low redshifts in our cosmological boxes.

2.3. The dust model

In this work, the treatment of dust pollution from evolved stars (AGBs and SNae) and subsequent evolution of the grain population follows the model illustrated in our recent works (Gjergo et al. 2018; Granato et al. 2021; Parente et al. 2022), including the effect of hot gas cooling due to collision with grains. We adopted the choice of parameters calibrated in the last paper. The treatment of dust size distribution and the size dependence of related processes adopts the computationally cheap two-size approximation devised by Hirashita (2015). It was derived by comparing with results of one-zone computation, and its satisfactory performances have also been tested against more detailed treatments also by other groups (e.g. Aoyama et al. 2020).

3. Treatment of H₂ evolution in MUPPI

Building on the groundwork of the already implemented dust formation and evolution modelling in MUPPI (see Granato et al. 2021), the present study takes a step forward by integrating a model for the formation of H₂ on dust grains. This prompts us to introduce an additional equation to the existing set of five differential equations. These equations delineate the evolution of the cold, hot, and “virtual” stellar component, the thermal energy of the hot phase, and the AGN feedback energy acting on the particle (Eqs. (28) to (31) in Valentini et al. 2020). The new equation describes the mass evolution of the H₂ portion (see below), fully encompassed within the cold phase.

MUPPI assumes that the molecular gas is the fuel for SF. Therefore, the SFR can be written as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{sf}}=f_{\ast }\frac{f_{\mathrm{mol}}{M}_{\mathrm{c}}}{t_{\mathrm{ff}}},\end{array}$$(2)

where f_* is an efficiency factor, the local free-fall time³t_ff refers to the cold gas and f_mol is the fraction of the latter in molecular clouds. In previous works, we computed the latter using external prescriptions, either theoretical (Krumholz et al. 2009) or phenomenological (Blitz & Rosolowsky 2006) in nature (see Valentini et al. 2023, and references therein). With our new implementation of molecular hydrogen evolution we compute the virtual SFR of MUPPI gas particles as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{sf}}=f_{\ast }\frac{M_{\mathrm{H}2}}{(1-Y)t_{\mathrm{ff}}},\end{array}$$(3)

where Y = 0.24 is the Helium fraction, and M_H2 is obtained by adding to MUPPI a new differential equation that describes the physical processes involved in the H₂ evolution:

$$\begin{array}{cc}\hfill {\dot{M}}_{{\mathrm{H}}_2}& ={\dot{M}}_{{\mathrm{formH}}_2}-{\dot{M}}_{{\mathrm{destH}}_2}-{\dot{M}}_{\mathrm{sf}}-{\dot{M}}_{\mathrm{ev}}-{\mathrm{f}}_{\mathrm{H}2}{\dot{M}}_{\mathrm{AGN}}.\hfill \end{array}$$(4)

Here the ${\dot{M}}_{{\mathrm{formH}}_2}$ and ${\dot{M}}_{{\mathrm{destH}}_2}$ terms account for the formation of H₂ on dust grains and its destruction by dissociating UV radiation, respectively. The terms ${\dot{M}}_{\mathrm{ev}}$ and Ṁ_AGN were already present in MUPPI, transferring mass from the cold to the hot phase (Eq. (29) in Valentini et al. 2020). The former represents the “evaporation” of the molecular phase due to stellar winds and is set to $0.1{\dot{M}}_{\mathrm{sf}}$. This destruction term is fully applied to the molecular phase in our treatment. The latter, Ṁ_AGN, represents the destruction due to the share of AGN feedback energy that couples to the cold phase. A fraction f_H2 = M_H2/M_c of this destruction rate is applied to the molecular phase. A schematic representation of the mass fluxes just described is shown in Fig. 1, depicting a MP particle.

3.1. H₂ Formation on dust grains

Various theoretical studies, starting with Gould & Salpeter (1963)⁴, have shown that the H₂ formation in the gas-phase within the ISM cannot account for the inferred abundance of H₂ in space. According to their estimates, the most important mechanism for forming molecular hydrogen is the encounter of hydrogen atoms on the surface of interstellar grains. This conclusion is currently widely accepted (see Wakelam et al. 2017, for a recent review) and adopted in theoretical works to estimate the H₂ formation rate in dust enriched ISM (e.g. Krumholz et al. 2009; Bekki 2013; Bron et al. 2014; Tomassetti et al. 2015; Valdivia et al. 2016).

Hollenbach et al. (1971) computed the rate per unit volume of molecular formation on dust grains as

$$\begin{array}{c}\hfill R_{\mathrm{G}}=\frac{1}{2}\gamma ⟨v_{\mathrm{H}}⟩n_{\mathrm{H}}{n}_{\mathrm{gr}}⟨{\sigma }_{\mathrm{gr}}⟩;\end{array}$$(5)

the product of the last four factors gives the rate per unit volume at which HI atoms, with number density n_H and mean velocities ⟨v_H⟩, collide with dust grains, with number density n_gr and mean geometric cross-section ⟨σ_gr⟩. The efficiency factor γ represents the fraction of HI atoms that combine on the surface of dust grains and detach from it as H₂ molecules, referred to as the recombination coefficient.

Using Eq. (5) and the two-size approximation adopted by our dust evolution model, the total rate of H₂ formation, including both small and large grains, reads as follows:

$$\begin{array}{c}\hfill R_{\mathrm{G}}=\frac{1}{2}\gamma ⟨v_{\mathrm{H}}⟩n_{\mathrm{H}}(n_{\mathrm{gr},\mathrm{s}}⟨{\sigma }_{\mathrm{gr},\mathrm{s}}⟩+n_{\mathrm{gr},\mathrm{l}}⟨{\sigma }_{\mathrm{gr},\mathrm{l}}⟩)\end{array}$$(6)

where subscripts gr,s and gr,1 stand for small and large grains. Following our previous works, we adopt for them the reference sizes of $0.005\mathrm{\mu }\mathrm{m}$ and $0.05\mathrm{\mu }\mathrm{m}$, respectively^5,6. Finally, by setting $\gamma ⟨v_{\mathrm{H}}⟩=5\times {10}^4\mathrm{cm}{\sec }^{-1}$ (Hollenbach et al. 1971), we can estimate the term for H₂ formation in Eq. (4) as

$$\begin{array}{c}\hfill {\dot{M}}_{{\mathrm{formH}}_2}=C_{\rho }2m_{\mathrm{p}}{R}_{\mathrm{G}}{V}_{\mathrm{c}}\end{array}$$(7)

where m_p is the proton mass, V_c is the cold phase volume⁷ and the clumping factor C_ρ = 10 is intended to boost the density of the cold phase to account for the unresolved molecular cloud densities in our model.

Since gas particles initiate their evolution without dust, the formation of H₂, and consequently SF, would never begin without some initial “seeding” of H₂. To address this issue, whenever a particle enters MP, and its D/G ratio is lower than 0.03 (D/G)_⊙, with (D/G)_⊙ = 0.01, SF is conducted using the empirical prescription from Blitz & Rosolowsky (2006), consistent with the previous version of MUPPI. We verified that our results are stable against changes in these thresholds by a factor of at least two.

3.2. The photo-dissociation rate

Photons from recently formed stars in the LW region ($\sim 940-1070\mathrm{\AA }$) are particularly harmful to H₂. The rate at which H₂ dissociates is estimated following Abel et al. (1997):

$$\begin{array}{c}\hfill k_{\mathrm{LW}}=C_{\mathrm{LW}}4\pi J_{\nu ,\mathrm{LW}}\end{array}$$(8)

where J_ν, LW denotes the specific radiative intensity averaged over angles at $h\nu =12.87\mathrm{eV}$, and C_LW = 1.1 × 10⁸ erg⁻¹cm² sr Hz.

However, we point out that several model computations (e.g. Gnedin et al. 2009; Mac Low & Glover 2012; Krumholz et al. 2009; Christensen et al. 2012) suggest that this process may be subdominant in setting the molecular gas fraction at galactic scales in situations where H₂ production is dominated by dust surface catalysis, that is in the conditions to which our modelling applies. This conclusion derives from the high efficiency of H₂ self-shielding and dust shielding, effectively protecting most of the volume occupied by ${\mathrm{H}}_2$. Nevertheless, we implemented an approximate treatment of LW H₂ dissociation to have a first-order mean to estimate its possible importance.

Rigorously accounting for this H₂ destruction effect in cosmological simulations is currently unfeasible. It would require on-the-fly radiative transfer computation in a dusty medium, including the H₂ self-absorption. However, besides the prohibitive numerical cost, the distributions of the relevant emitting sources (stars younger than 5–10 Myr) and the dense absorbing medium are unresolved in such simulations. Therefore, previous works dealing with this aspect in a similar context (e.g. Gnedin et al. 2009; Christensen et al. 2012) adopted simplified treatments. Our approach, described in the following, focuses on modelling the phenomena at unresolved scales (namely within MP gas particles), where we believe the bulk of the effect arises. We establish a simple connection between the dissociation rate of H₂ and the recent SFR within an MP particle.

The starting point is to notice that most LW photons come from stars still embedded or relatively close to (Giant) MC that are optically thick at these short wavelengths. We pre-compute the specific intrinsic (meaning before any absorption effect) luminosity for a constant SFR with GRASIL (Silva et al. 1998; Granato et al. 2000), at the LW frequency, as a function of age t. For the Chabrier IMF adopted in this work, the result is well described by the following polynomial:

$$\begin{array}{c}\hfill log\frac{L_{\nu ,\mathrm{LW}}[{\mathrm{erg}\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}]}{\mathrm{SFR}[{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}]}=0.126x^3-0.774x^2+1.573x+27.03\end{array}$$(9)

where x = log(t[Myr]). This expression reproduces almost exactly the population synthesis result for t ranging from 0.3 to 300 Myr and stellar metallicity close to 0.3 Z_⊙. After 100 Myr the LW luminosity saturates to the 100 Myr value, with a contribution from stars older than 20–30 Myr lower than 20%. The metallicity dependence is weak enough to ensure accuracy within 20% between 0.001 Z_⊙ to 3 Z_⊙.

We use Eq. (9) to compute, at a given time the LW radiation emitted by stars (usually virtual stars before the spawning of a new stellar particle) born in the MP particles during the previous timesteps of the MP cycle, by finite difference. However, dust attenuates LW radiation efficiently, substantially protecting H₂ from LW photo-dissociation. Moreover, also H₂ self-shielding has a decisive protective effect. We refer the reader to the Appendix B for details on our estimate of the shielding factor. Here, we point out the primary outcome of this estimate. If the physical properties of the MC in the Milky Way (Miville-Deschênes et al. 2017) are representative of the general population of MC, we expect the term ${\dot{M}}_{{\mathrm{destH}}_2}$ in Eq. (4) to play a negligible role. This conclusion is further strengthened if, as observations suggest (Dessauges-Zavadsky et al. 2019), MC in galaxies at higher redshift and/or with higher star formation activity exhibit higher column densities. Even by artificially enhancing the LW radiation field given by Eq. (9) by a factor of 10, the results are almost indistinguishable from those obtained by completely neglecting the LW destruction, indicating that the conclusion is robust.

4. Model comparison with observational data

To validate our model, this section contrasts key statistical characteristics of simulated galaxies across our five boxes, introduced in Section 2.1, with existing observational data. We adjusted the model parameters to find the best match with the observed cosmic SFRD, the stellar mass function, and the stellar mass-halo mass function at z = 0. Consequently, the remaining results presented here can be considered predictions of our model. A visual rendition of the outcome of our model at redshift z = 0 can be seen in Fig. 3, which illustrates the density maps for the various components within one of our simulated volumes.

Fig. 3. Stellar, gas, dust, and H2 mass density maps in one of the (26 Mpc)3 simulation box computed in pixels of 200 × 200 kpc2.

To identify galaxies, we employ SUBFIND (Springel et al. 2001; Dolag et al. 2009), considering objects with stellar masses larger than 3 × 10⁷ M_⊙, which corresponds to an average of 25 stellar particles. This low number of minimum resolved stellar elements should be sufficient to reasonably define the basic properties of galaxies considered in the present work, i.e. masses of the various components. The integrated masses and SFR in this study are evaluated within twice the half-stellar mass radius. These computations focus solely on particles that are linked to a specific subhalo as identified by SUBFIND.

We begin by examining properties related to the stellar and dust components. These comparisons were feasible even before integrating the H₂ evolution into our modelling (see Parente et al. 2022). However, we have now achieved notable improvements in aligning with the data, thanks to the modifications we introduced in MUPPI and, to some extent, the inclusion of the H₂ modelling.

4.1. Stellar properties of the galaxy population

In this subsection, we test our simulations against the main observed statistical properties of the stellar component of galaxies, namely the cosmic SFRD, the SMF, the main sequence of star-forming galaxies (MS), and the SHMR, the latter three at z = 0, 1, and 2. We reiterate that the SFRD, SMF, and the SHMR at z = 0 have been used to calibrate the model parameters, and therefore they are not model predictions. Nevertheless, it is not trivial that the model can reproduce those observables with a proper choice of free parameters that are fixed for the rest of our study.

4.1.1. The star formation rate density

The green curve in Fig. 4 shows the average SFRD from the five simulated boxes and the full dispersion of the five individual curves. Our model utilizes a Chabrier (2003) IMF, and the displayed data have been reported to this same IMF whenever necessary. The agreement with observational determinations is quite good, comfortably fitting within the existing range of variations among them. This improvement is pronounced compared to the previous version of MUPPI model (e.g. Parente et al. 2022), which relied on estimating molecular mass through the empirical law proposed by Blitz & Rosolowsky (2006), alongside a setup calibrated only on zoom-in simulations of MW-like galaxies (see Section 2.2).

Fig. 4. Average cosmic SFRD calculated from five ${(26\mathrm{Mpc})}^3$ simulation boxes, with the shaded region indicating the full dispersion of the curves. We compare with some determinations from literature: solid and dashed black lines represent the fits to the cosmic SFRD obtained by Madau & Dickinson (2014) and Davies et al. (2016) respectively, while black circles show the determination of Driver et al. (2018). Following their recommendation, the data from Madau & Dickinson (2014) have been scaled by a factor of 0.63 to convert from a Salpeter (1955) to a Chabrier (2003) IMF. The remaining data points require no further conversion.

4.1.2. The stellar mass function

The observational galaxy SMF at redshift z = 0.1,  1.0, and 2, shown in Fig. 5, is also well recovered by our model, up to the stellar mass statistically accessible by the simulated volume. We compare our model at z = 0.1 with the best fit obtained in Rodríguez-Puebla et al. (2020) for their compilation of data from different works. The shaded grey region represents systematic errors, encompassing variations arising from differences in mass-to-light ratios and mass definitions (photometry) within the observed comparison samples utilized in that study. A similar good agreement is found when contrasting our findings with data from Ilbert et al. (2010) and Mortlock et al. (2011) for z > 0. However, our SMF does not encompass the higher mass ranges observed in their data. This limitation stems from this study’s relatively small simulation volumes, (26 Mpc)³.

Fig. 5. Stellar mass function is shown by the green curves at z = 0.1, z = 1, and z = 2, based on the five simulated boxes in our fiducial model. The error bars represent Poissonian errors for each stellar mass bin, typically smaller than the size of the point. The shaded green region encompasses the full dispersion of the individual SMFs corresponding to each of the five simulation boxes. We compare our model at z = 0.1 with the best-fit model of Rodríguez-Puebla et al. (2020); the grey area is an estimate of systematic errors which captures variations due to differences in the mass-to-light ratios and mass definitions (photometry) of the observational samples in that work. Simulations at z = 1 are confronted with the data from Ilbert et al. (2010) at two redshift ranges and at z = 2 with the data in Mortlock et al. (2011).

4.1.3. The main sequence

The Mstar vs. SFR relation is illustrated in Fig. 6, with green contours depicting our complete sample of galaxies. We contrast our predicted MS of star-forming galaxies, represented by the green curve connecting the data points, with the observational findings outlined in Popesso et al. (2023), which summarizes an extensive homogenized compilation of studies published since 2014. Our results refer to a sub-sample of simulated galaxies selected using the same criteria as in the latter paper, namely 8.5 < log(Mstar/M_⊙) < 11.5 and 0.01 < SFR/(M_⊙ yr⁻¹) < 500. The normalization is generally reasonably well reproduced, with a tendency to slightly under-predict the observed MS at a given stellar mass.

Fig. 6. Distribution of our simulated galaxies in the stellar mass-SFR plane is shown by the contours. Solid black lines are the MS determinations for observed galaxies presented in Popesso et al. (2023), resulting from an extensive homogenized compilation of studies published since 2014 (their Eq. (10)). The shaded grey band visualize the typical dispersion of determinations used to get their fit. The green curves connecting the points represent the median SFR per mass bin for a sub-sample of simulated galaxies selected using the same criteria as in Popesso et al. (2023), namely 8.5 < log(Mstar/M⊙) < 11.5 and 0.01 < SFR/( M⊙ yr−1) < 500. Bars delineate the 25–75% percentiles.

Across the various determinations included in the Popesso et al. (2023) compilation, there is a variation in the redshift and stellar mass at which the MS bends. Nevertheless, after converting all observations to a common calibration, they find that the main sequence displays a curvature towards higher stellar masses at all redshifts. Notably, our findings indicate the continuous curvature of this relationship in the log-log plane, consistent with the observational estimate. In contrast, other cosmological simulations, including Illustris, Eagle, Mufasa, IllustrisTNG, and Simba, predict a linear correlation, at least up to high masses (e.g., Donnari et al. 2019; Hahn et al. 2019; Davé et al. 2019, and references therein). In our model, a similar bending is also present in the MH2-Mstar relation (Section 4.3). Considering that H₂ is the fuel for SF, the latter implies that the bending of the MS is mainly driven by a relative lack of H₂ in the most massive galaxies.

Finally, we note that our model also predicts a population of passive galaxies, which is increasingly clear with decreasing redshift. At z = 0.1 the population is located at log Mstar/M_⊙ ∼ 10.5 and SFR ∼ 0.1 M_⊙/yr, broadly consistent with observations (e.g., Renzini & Peng 2015).

4.1.4. The stellar-halo mass relation

Figure 7 shows the z = 0 relationship between the M₂₀₀ halo mass and the stellar mass of central galaxies, delineated by the baryon conversion efficiency ϵ = Mstar/(f_b M₂₀₀), with f_b = Ω_b/Ω_m. We compare our results with the abundance matching estimate by Moster et al. (2013) and with the empirical model by Moster et al. (2018). The two latter determinations are strongly data-driven, particularly the abundance matching determination, which agrees with our result. Although our simulation volume does not allow us to extensively test our model in massive (≳10¹³ M_⊙) halos, simulated galaxies feature a decreasing ϵ at M₂₀₀ ≳ 10¹² M_⊙, in agreement with the reported data-based determinations.

Fig. 7. Stellar mass-halo mass relation at z = 0 is assessed by considering galaxies across all our five simulated boxes, fb = 0.157 is the baryon fraction Ωb/Ωm. The stellar mass is computed within 0.1R200. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 and 0.2 dex-wide bins along the x and y axes, respectively. We compare our results with abundance matching estimate by Moster et al. (2013), and with the empirical model by Moster et al. (2018).

4.2. Dust mass function

Figure 8 shows the predicted dust mass function (DMF) at different redshifts and compares it with observational data (black symbols) and other simulations (grey lines). In general, previous hydrodynamic simulations (McKinnon et al. 2017; Hou et al. 2019; Li et al. 2019) had difficulties in reproducing the dust mass function simultaneously at low and high redshift. Our simulation outcomes are in good agreement with local observations. However, accurately probing redshifts of z = 1 and z = 2 presents some challenges. Observations lack coverage of the range of small dust masses that our simulations sample. Conversely, the limited volume of our simulations could lead to an inability to capture the high masses observed in real data, consequently resulting in an under-prediction of the DMF at the highest dust mass bin.

Fig. 8. Dust mass function at z = 0.1, z = 1, and z = 2 considering galaxies in the five simulated boxes are represented by green dots connected by solid lines. The error bars represent the Poissonian error for each dust mass bin. The shaded green region encompasses the full dispersion of the individual DMF corresponding to each of the five simulation boxes. The decline of the DMF at low masses is due to the imposed stellar mass resolution limit. We compare our results with observational data, represented by symbols and simulations delineated by lines. Observations at z = 0.1 correspond to Pozzi et al. (2020), Beeston et al. (2018), Dunne et al. (2011), and Vlahakis et al. (2005), at z = 1 (Pozzi et al. 2020) and at z = 2 (Pozzi et al. 2020 and Dunne et al. 2003). Simulations correspond to the works of Parente et al. (2022), Hou et al. (2019), Li et al. (2019), and McKinnon et al. (2017).

We reiterate the significant improvement over our previous model, as discussed in Parente et al. (2022), represented by the solid grey line in the figure. Importantly, this enhancement was achieved by keeping the dust model unchanged while modifying aspects of the galaxy evolution model (see Section 2.2) and including the H₂ modelling.

4.3. H₂ in model galaxies

In the present sub-section, we discuss quantities related to the H₂ content of galaxies, namely the H₂ cosmic density, the H₂ mass function, and the stellar-H₂ mass relation at various redshifts.

Figure 9 illustrates the evolution of the cosmic H₂ content, showcasing the molecular gas density ρ_H2 for our model alongside constraints derived from observations. These constraints comprise an update of the compilation presented in Maio et al. (2022). Despite being by far the most abundant molecule, the H₂ content of galaxies cannot be estimated directly due to its lack of emission at the cold temperatures at which most of its mass resides in MC. The reported estimates are mainly based on CO emission, assuming some conversion factor from CO to H₂, or on IR emission from dust, assuming a D/G ratio. Both methods suffer from large systematic uncertainties, ∼2, widely debated in the literature and likely increasing with redshift (Tacconi et al. 2020). Nevertheless, CO-based, dust-based methods and our model predict a molecular gas cosmic density peak at z ∼ 1.5, and align with the broad trend of decreasing over the last ∼9 billion years. The evolution of ρ_H2 in our simulated box is somewhat milder than the observational data suggested. However, the choice of the CO-to-H₂ conversion factor influences the latter. Indeed, Decarli et al. (2019) find, using the ASPECS-Pilot survey, that the decline of ρ_H2 by a factor of approximately 6 from its peak at z ∼ 1.5 to the present would reduce to a factor of 3 if they had used a conversion factor of 2${\mathrm{M}}_{\odot }{(\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^{-2})}^{-1}$ instead of their favoured choice 3.6 ${\mathrm{M}}_{\odot }{(\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^{-2})}^{-1}$ for galaxies at z > 1.

Fig. 9. Mean cosmic evolution of the molecular hydrogen density considering the five simulated boxes. The shaded region encloses the full dispersion of the five curves. The right-hand axis shows ΩH2, the H2 mass density normalised to the present-day critical density. Observational constraints are ASPECS (ALMA Spectroscopic Survey in the Hubble Ultra Deep Field, Decarli et al. 2020); VLASPECS (which uses the NSF’s Karl G. Jansky Very Large Array, VLA, covering part of the ASPECS footprint, Riechers et al. 2020); COLDz (VLA CO Luminosity Density at High Redshift in the COSMOS and GOODS-North fields combined; Riechers et al. 2019); PHIBBS2 (Plateau de Bure High-z Blue Sequence Survey 2, which uses the NOrthern Extended Millimeter Array – NOEMA – to observe field sources, Lenkić et al. 2020); xCOLD GASS (eXtended CO Legacy Database for the Galaxy All Sky Survey, Fletcher et al. 2021); ALMACAL-abs (blind search of CO absorbers against the ALMA Calibrator archive sources, Klitsch et al. 2019); ALMACAL-CO (CO emission-line survey Hamanowicz et al. 2023); Herschel PEP (PACS Evolutionary Probes Herschel survey in the GOODS-N -S fields, Berta et al. 2013); and UKIDS-UDS (UK InfraRed Telescope – UKIRT – Infrared Deep Sky Survey in the Ultra-Deep Survey, Garratt et al. 2021). Data listed in the top right legend do not use direct measurement of CO emission to estimate molecular masses.

While taking the various conversion factors adopted by the data in Fig. 9 at face value, our prediction, though exhibiting the same qualitative behaviour, underestimates the measurements at z ∼ 1–2 by a factor ∼2 and lies in the lower limit of observational constraints for higher redshift. We cannot exclude that at very high redshifts when the dust cosmic abundance was likely much lower than it is now, alternative H₂ formation channels, other than dust catalysis, could improve the match with the data. These processes are not included in our computations. However, this is very unlikely to occur in the redshift range 1–2, where the tension with current observational estimates is more evident. In fact, in this redshift regime, cosmic dust abundance is generally estimated to be even higher than in the present-day Universe (see e.g. Fig. 10 in Parente et al. 2023).

We also analyse the H₂ mass function (H2MF) at different redshifts. In Fig. 10, the H2MF in the Local Universe is compared with the determination of Rodríguez-Puebla et al. (2020), Fletcher et al. (2021), and Andreani et al. (2020). Our predictions closely reproduce the results of Andreani et al. (2020), particularly those adopting the CO-to-H₂ conversion factor depending on luminosity. At z = 1 and z = 2, the observational validation of the H2MF is limited to the highest masses. An inspection of the figure reveals that at z = 1 the model H2MF falls, in the most massive bin, within the observed constraints outlined in Decarli et al. (2016) whose determination corresponds to the redshift bin z = [0.695 − 1.174] with ⟨z⟩ = 0.95 Similarly, at redshift z = 2, our H2MF at the highest mass bin falls within the constraints by Decarli et al. (2016), despite the latter corresponding to a higher redshift bin. The other two observational constraints depicted in the plot are derived also using galaxy samples at higher redshifts than our determination: Decarli et al. (2019), which mirrors the methodology of Decarli et al. (2016) but in a larger volume, and Riechers et al. (2019).

Fig. 10. H2 mass function at z = 0.1, 1.0, and z = 2.0. The error bars represent Poissonian errors for each stellar mass bin, typically smaller than the size of the point. The shaded green region encompasses the full dispersion of the individual H2MFs corresponding to each of the five simulation boxes. In the local Universe, we report for comparison the H2 mass function obtained by Andreani et al. (2020), who use a fixed (COfix) or a luminosity dependent (CO-L) CO-to-H2 conversion factor and the best-fit model of Rodríguez-Puebla et al. (2020) and Fletcher et al. (2021). At higher redshifts, the boxes show constraints from Decarli et al. (2016) (ALMA ASPECS-Pilot survey), Decarli et al. (2019) (ALMA ASPECS-LP survey), and Riechers et al. (2019) (COLDz survey). The H2MF at redshift z = 0 is also shown in the middle and right panels using a thin dotted line to facilitate comparison.

We obtain little evolution in the H2MF from z ∼ 2 to the present. This becomes apparent when comparing the determinations of the H2MF at redshifts z = 1 and z = 2 with the thin dotted line in Fig. 10, representing the H2MF at redshift z = 0. This finding contrasts with, for example, Riechers et al. (2019) and Decarli et al. (2019), who report a shift of the characteristic CO luminosity towards lower values of one order of magnitude from z ∼ 2 to 0 (an equal shift in the H2MF for a constant CO-to-H₂ conversion factor).

In Fig. 11, we depict the relationship between the stellar mass and the H₂ mass content of galaxies at z = 0, 1, and 2. We also plot the z = 0 relationship in the central and right panels to facilitate comparison across these redshifts. It can be observed that galaxies of equal mass exhibit higher H₂ content as redshift increases from z = 0 to z = 2, particularly at higher stellar masses (a factor ∼4). At redshift z = 0.1, we compare our predicted H₂ masses with the determinations by Saintonge et al. (2017). In their study, galaxies without CO detections are assigned upper limits, represented by downward arrows in the figure. The contours representing our full population of galaxies align well with the observational data. The median H₂ mass per stellar mass bin falls below observational detections at the highest stellar masses. However, considering the non-detections, our median might be a fair representation of the real population, within which galaxies undergoing quenching are more likely to go undetected.

Fig. 11. Molecular hydrogen mass as a function of stellar mass for galaxies in our simulation. In addition to contours, solid circles display the median $\mathrm{MH}2$ within bins of stellar mass, along with the 25%–75% percentiles for model galaxies with MH2 > 0. To facilitate comparison, the dotted line in the central and right panels corresponds to the z = 0 relation. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes. Empty circles in the left panel represent the xCOLD GASS sample reported in Saintonge et al. (2017); non-detection in this sample is represented by downward arrows (0.01 < redshift < 0.05).

Finally, we notice that the MH2 vs. Mstar relationship bends in the log-log plane at high stellar masses, driving the bending of the MS discussed in Section 4.1.3 (see Fig. 6). The latter bending is somewhat less pronounced than the former. Indeed, the depletion time MH2/SFR tends to decrease with increasing mass in our simulations, as shown in Fig. 12. This occurs because more massive galaxies feature, on average, a shorter free fall time t_ff in the cold phase of MP particles (see Eq. (3)). These galaxies, which tend to be spheroid-dominated, are quenched because AGN feedback decreases their SF fuel and not because morphological quenching (e.g. Martig et al. 2009; Gensior et al. 2020) stabilises the cold gas against gravitational collapse.

Fig. 12. Depletion time as a function of stellar mass for model galaxies. In addition to contours, solid circles display the median depletion time within bins of stellar mass, along with the 25%–75% percentiles for model galaxies with $\mathrm{MH}2>0$. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

5. The Kennicutt-Schmidt relation

Over the decades, a well-defined correlation between the disc-averaged gas and SFR surface densities in galaxies has been identified and extensively studied (Schmidt 1959; Kennicutt 1989). This relationship is now commonly referred to as the (in its disc-averaged form, either global or integrated) Kennicutt-Schmidt (KS) law. This law holds over several orders of magnitude and has inspired most of the sub-resolution prescriptions adopted in galaxy formation models, including ours. More recently, observational studies have concentrated on analysing this correlation at resolved scales and distinguishing between the atomic and molecular gas surface densities. The general result is that the KS relation holds even within galaxies, and it is primarily driven by the molecular component (see de los Reyes & Kennicutt 2019, and references therein). The latter point confirms that molecular clouds are closely linked to the fuel for star formation.

The integrated KS relation for our model galaxies at various redshifts is shown in Fig. 13. For each galaxy, an arrow indicates the shift from its position in the plot when considering only molecular gas (the molecular Kennicutt-Schmidt, mKS) to its position when both atomic and molecular gas are included (the total Kennicutt-Schmidt, KS). The arrow colour encodes the H₂ mass fraction of the galaxy. Both the colour and the length of the arrows indicate a clear increase of this fraction with the surface density. Moreover, at low densities, the arrowheads display a somewhat larger horizontal dispersion than the tails: the molecular gas alone is more closely related to SFR than the total gas. Conversely, at high surface density, HI is sub-dominant, and the total mKS approaches the total KS relation. The black solid line corresponds to a fit to Local Universe data compiled by de los Reyes & Kennicutt (2019), consolidating information from 114 publications on the KS. Its slope is n = 1.03 and is reported in all panels for reference. To approximate the observational sample, the simulated sample includes galaxies with stellar masses greater than 10⁸ M_⊙ and bulge-to-total ratios B/T < 0.5. We apply the same selection criteria for z = 1 and z = 2 galaxies, shown in the central and right panels.

Fig. 13. Integrated KS relation. The tail of the arrows represents ΣH2, while the head indicates ΣH2 + HI. We show the results for galaxies with stellar masses Mstar > 108 M⊙ and bulge-to-total ratios B/T < 0.5. The colour bar encodes the H2 fraction of each galaxy. Red/blue dashed line segments are linear fits in the log-log plane for galaxies with H2 fraction lower/greater than 0.5. The corresponding slopes are labelled in each panel for H2 and H2 + HI gas phases. Solid line represents one of the fits (slope n = 1.03) from de los Reyes & Kennicutt (2019) to the ΣH2-ΣSFR relation obtained from a compilation of observational samples of local spirals.

The simulation captures the integrated mKS relation (the arrow tails) at redshift z = 0, exhibiting a reasonable alignment with the n = 1 slope at ${\mathrm{\Sigma }}_{\mathrm{H}2}\lesssim 4{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}$. Above this limit, a steepening of the relationship is evident. At higher redshift z = 1 and 2 (central and right panels of Fig. 13), we find that the two-slope shape of the relation is preserved, but the normalization increases by a factor ∼2 and 4, respectively.

Recently, Teng et al. (2024) analysed the central 1.5 kpc regions of 80 galaxies from the PHANGS-ALMA survey (Leroy et al. 2021), by adopting a CO-to-H₂ conversion factor depending on the local CO velocity dispersion at a 150 pc scale. They concluded that at the high Σ_H2 levels reached in galaxy centres, the mKS relation becomes steeper, a phenomenon that goes undetected when a fixed CO-to-H₂ conversion factor is used. While our estimate of the mKS is integrated over the entire galaxy, it is suggestive that we get a clear steepening to a similar slope at the high Σ_H2 end.

In our model, the steepening of the mKS relation at higher Σ_H2 has the following interpretation. The SFR is computed in MP particles using Eq. (3). Thus, we expect that Σ_SFR ∝ Σ_H2/t_ff, where the free-fall time is computed using the density of the cold phase $t_{\mathrm{ff}}\propto {\rho }_{\mathrm{c}}^{-0.5}$. The molecular fraction tends to be low at low Σ_H2 values, as illustrated by the colour coding in the figure. Then, the free-fall time is substantially unaffected by the H₂ density, being the latter a sub-dominant component of the cold gas (see scheme in Fig. 1). Conversely, as Σ_H2 increases, so does the H₂ fraction, and the free-fall time tends to scale progressively as $t_{\mathrm{ff}}\propto {\rho }_{\mathrm{H}2}^{-0.5}$. In the former regime, we then expect a linear proportionality between the H₂ and the SFR densities, steepening towards ${\mathrm{\Sigma }}_{\mathrm{SFR}}\propto {\mathrm{\Sigma }}_{\mathrm{H}2}^{1.5}$ in the latter regime. Moreover, the correlation shows a larger dispersion in the former regime. When the H₂ fraction is low, the density and the free-fall time are determined mainly by HI. As a result, in this regime Σ_SFR depends on the density of two species, HI and H2, with the former exhibiting dispersion at any given value of the latter. Conversely, at high Σ_H2, where almost all the cold mass is made of H₂, SFR computation is based solely on one mass species, reducing its dispersion.

6. H₂ mass determination from dust

As mentioned before, the emission of H₂ molecules is not directly measured. Therefore, any evaluation of the molecular gas content in galaxies is limited to indirect methods. As molecules predominantly form on the surface of grains, which also act as a shield against dissociating radiation, one such method involves utilizing dust as a tracer for molecular gas.

In a recent work, by using star-forming galaxies with stellar masses > 10¹⁰ M_⊙ and z ≲ 0.2, Bertemes et al. (2018) establish a formula to derive molecular masses from dust masses (their Eq. (15)): log(M_mol) = log(M_dust)+1.83 + 0.12[(12 + log(O/H)) + 8.67]. This involved computing CO-based molecular masses, which introduces a dependency on metallicity through the conversion factor CO-to-H₂. In their work, the coefficient of 0.12 multiplying the metallicity term is obtained by adopting a geometric mean of the CO-to-H₂ conversion factors derived in Genzel et al. (2012) (G12) and Bolatto et al. (2013) (B13). This coefficient changes to −0.37 or 0.61 if the geometric mean is replaced by the CO-to-H₂ conversion factors derived in G12 and B13, respectively. Their formulas yield similar values around the solar metallicity but diverge significantly at lower and higher Z. In G12, the CO-to-H₂ conversion factor is derived using a mixed sample of main-sequence star-forming galaxies with z ≳ 1 plus five Local Group galaxies. The whole sample metallicity range is −0.9 ≲ log(Z/Z_⊙)≲0.3. On the other hand, B13 formula is derived from theoretical grounds.

In Figure 14, contours show H₂-to-dust ratio as a function of metallicity for the simulated galaxies. We over-plot the two lines corresponding to the Bertemes et al. (2018)H₂-to-dust equation (without He fraction) but using the mentioned G12 and B13 conversion factors. The circles connected by solid lines illustrate the H₂-to-dust ratios plotted against metallicity for galaxies within the MS, specifically those with a deviation $\mid \log (\mathrm{SFR})-\log ({\mathrm{SFR}}_{\mathrm{MS}})\mid <0.5$. The MS at each redshift were determined by fitting a functional form as the one proposed in Popesso et al. (2023), which is suitable for our model (see Subsection 4.1.3 and Fig. 6). The general trend across all redshifts suggests that galaxies with higher metallicities exhibit smaller H₂-to-dust ratios, with a higher normalization associated with higher redshifts. This trend is at odd with Eq. (15) in B18, based on the geometric mean of G12 and B13 prescription and featuring a positive slope. However, it aligns with their findings when only the G12 prescription is considered. The match is particularly good at z = 0 and z = 1, where the curve for model galaxies approaches G12 line in the range of metallicities used by that study.

The relationship between the H₂-to-dust ratio and stellar mass is illustrated in Fig. 15. As in the previous figure, the H₂-to-dust fraction for the global population is smaller now than in the past. In our simulated galaxy sample, the most massive galaxies ($10.75\lesssim \log (\mathrm{Mstar}/{\mathrm{M}}_{\odot })\lesssim 11.25$) exhibit an H₂-to-dust ratio of approximately 65 at z = 2 and 30 at z = 0. While at the lowest mass end, the typical H₂-to-dust ratio is about 160 at z = 2 and 50 at z = 0. However, the dispersion in the H₂-to-dust values at these lower masses is significantly higher compared to the high mass end. Consistently, this feature is also present at low metallicities in Fig. 14 .

Fig. 14. Contours depict the H2-to-dust ratio as a function of gas metallicity for all model galaxies. MS galaxies ($\mid \log (\mathrm{SFR})-\log ({\mathrm{SFR}}_{\mathrm{MS}})\mid <0.5$) at the given redshift are represented by the circles, which depict the median H2-to-dust ratio within stellar mass bins, accompanied by the 25%–75% percentiles. The two lines labelled as G12 and B13 correspond to fits obtained in Bertemes et al. (2018) when using the CO-to-H2 conversion factors derived in Genzel et al. (2012) or Bolatto et al. (2013), respectively. The colour contours enclose grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

Fig. 15. Contours depict the H2-to-dust ratio as a function of stellar mass for all model galaxies. MS galaxies ($\mid \log (\mathrm{SFR})-\log ({\mathrm{SFR}}_{\mathrm{MS}})\mid <0.5$) at the given redshift are represented by the circles, which depict the median H2-to-dust ratio within stellar mass bins, accompanied by the 25%–75% percentiles. The colour contours enclose grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

7. The atomic to molecular transition

Figure 16 shows the dependence of MP gas particle H₂ fraction on the cold component sub-grid number density at redshift 0.1, 1, and 2. Since the cold gas in MP is assumed to have a given temperature of 300 K, its number density translates univocally into pressure, as shown on the top axis scale. The dotted black line is the relationship used in previous MUPPI works to estimate the molecular fraction from the gas pressure, inspired by the Blitz & Rosolowsky (2006) findings⁸. The colour coding in the top panel corresponds to the D/G ratio of particles. The fraction of H₂ increases with particle density and D/G ratio, as can be understood based on the linear dependence of the H₂ formation rate on atomic hydrogen and grain densities (see Eq. (6)). To illustrate this trend, the median H₂ fraction and number density are presented for particles grouped into six bins of D/G ratios (each spanning 0.5 dex), represented by coloured coded circles. These circles show that the typical H₂ fraction decreases as redshift decreases for particles within a specific D/G bin. For instance, the H₂ fraction of the most “dusty” particles in our simulations, characterized by $-2.5<\log (\mathrm{D}/\mathrm{G})$, exhibit a median value of approximately 1, 0.95, and 0.7 at redshifts z = 2, 1, and 0, respectively.

Fig. 16. Transition from atomic to molecular hydrogen in MP particles illustrated for the three analysed redshifts as a function of the sub-grid number density of the cold phase. The latter is assumed to have a temperature 300 K hence the pressure is proportional to the number density and shown in the top axis (k is the Boltzmann constant). The colour code shows the mean D/G mass ratio (top panel) and the gas phase metallicity (metals-to-gas mass ratio, bottom panel) of the gas particles falling within each pixel (150 pixels per axis). colour-coded circles represent the median H2 fraction and number density for particles within 0.5 dex bins of log(D/G) (top panel) and $\log (\mathrm{Zgas}/{\mathrm{Z}}_{\odot })$ (bottom panel). Meanwhile, stars indicate the same, but exclusively for particles in the transition region with H2 fraction between 0.4 and 0.6.

The transition from low to high molecular fractions at a given value of D/G is sharp and tends to shift to slightly higher values of n with increasing redshift. The latter is evidenced by the coloured stars, which show the median number density for particles of a given D/G having H₂ fraction between 0.45 and 0.55. This effect can be attributed to the decreasing importance of small grains in the dust budget at high redshift. While large grains dominate the dust mass, small grains dominate the surface area, and the H₂ formation process is a surface process. By taking into account these considerations, a practical application could involve developing prescriptions to estimate the typical H₂ fraction within a specific region of the ISM for models where H₂ evolution is not explicitly computed. This estimation would consider factors such as gas density and the budgets of small and large grains.

Finally, in the bottom row of Fig. 16, the pixels are colour-coded according to gas metallicity. It is evident that gas metallicity exhibits a much weaker correlation with the molecular fraction at a given density than D/G (top panel). This weak correlation is not unexpected, given that the relationship between dust and metal content is not straightforward: Metal abundance does not directly control the formation rate of H₂. Instead, it is primarily governed by dust abundance and gas density. Additionally, when considering the content of small grains, which, as mentioned earlier, significantly impact the rate of H₂ formation, the complexity of this relationship becomes even more apparent. Nevertheless, it can still be envisaged that, as expected, the molecular transition occurs typically at higher densities when the metallicity decreases.

8. Summary and conclusions

From the perspective of galaxy formation models, one of the most necessary tasks for the coming years is integrating cosmological simulations with models capable of reproducing the creation and evolution of the gaseous molecular phase. In this work, we have taken a significant step towards this direction within the MUPPI sub-resolution model (see Valentini et al. 2023, and references therein) to describe star formation from a multi-phase ISM, taking advantage of our recent incorporation of dust evolution in cosmological simulations (Granato et al. 2021). By establishing a connection between SF, ISM metal enrichment, dust production, and dust-promoted H₂ formation, we effectively model the star formation-dust-H₂ loop present in nature, resulting in a more physically grounded model for non-primordial ISM.

The main findings of this study can be summarized as follows:

• Our model, which is calibrated based on the cosmic star formation rate density, the stellar mass function at z = 0, and the stellar mass-halo mass relation at z = 0, satisfactorily matches several other observational data. This includes the stellar mass function at z = 1 and z = 2, the dust and H₂ mass functions at z = 0, z = 1, and z = 2, the cosmic evolution of H₂ density, and the relation between molecular and stellar mass in galaxies (see Figs. 4 to 11).

• The integrated molecular Kennicutt-Schmidt relation emerges in our model galaxies as early as redshift z ∼ 2, although with a higher normalization than at redshift 0 (Fig. 13). The slope remains close to unity across the considered redshifts for surface densities below ∼5 M_⊙/pc². However, at higher surface densities, we observe a slope more consistent with 1.5, which is reminiscent of a recent observational finding concerning the dense central 1.5 kpc regions of galaxies (Teng et al. 2024). The normalization of the simulated Kennicutt-Schmidt relation at redshift zero aligns well with determinations from the Local Universe. The emergence of a molecular KS relation and its steepening at high molecular density in our simulations are natural consequences of the adopted SF law (Eq. (3)). However, the normalization of the curve results from the calibration process, and its evolution can be considered a model prediction.

• The H₂-to-dust ratio diminishes with increasing metallicity across all examined redshifts (Fig. 14). Our model favours the G12 CO-to-H₂ conversion factor. A decrease of the H₂-to-dust ratio is also observed with increasing stellar masses (Fig. 15), particularly at the highest masses. The normalization of these curves decreases from z = 2 to z = 0, implying that the H₂-to-dust fraction for the global population is now smaller than it was in the past. Galaxies in our most massive bin ($10.75\lesssim \log (\mathrm{Mstar}/{\mathrm{M}}_{\odot })\lesssim 11.25$) have an H₂-to-dust ratio of ∼65,  30 at z = 2,  0 respectively. At the lowest mass end, the typical H₂-to-dust varies from ∼160 at z = 2 to ∼50 at z = 0, though the dispersion in the H₂-to-dust values at these masses (or metallicities) is much higher than at the high mass end (or metallicities).

• On a particle-by-particle basis, the transition from atomic to molecular hydrogen at a specific D/G value is abrupt (Fig. 16). It tends to shift towards higher number densities with increasing redshift. This behaviour is primarily influenced by the crucial role of small grains, which dominate the surface per unit mass and whose relative importance diminishes with increasing redshift compared to large grains.

• The typical H₂ fraction decreases as redshift decreases for particles with a given D/G (Fig. 16).

In upcoming studies, we plan to devise methods for estimating the typical H₂ fraction within specific regions of the ISM to be incorporated in models that do not explicitly compute H₂ or dust evolution. To comprehensively describe galaxy evolution, especially at high redshifts, we also aim to model the direct formation of H₂ in the gas phase, a process relevant in dust-free conditions.

Acknowledgments

We thank the anonymous referee for carefully reading our manuscript and the interesting comments provided. This project has received funding from the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) (PIP-2021-11220200102832CO), from the Agencia Nacional de Promoción de la Investigación, el Desarrollo Tecnológico y la Innovación de la República Argentina (PICT-2020-03690), and from the European Union’s HORIZON-MSCA-2021-SE-01 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement number 101086388 – Project (LACEGAL). MV is supported by the Italian Research Center on High Performance Computing, Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2, within the activities of Spoke 3, Astrophysics and Cosmos Observations, and by the INFN Indark Grant. Simulations have been carried out at the computing centre of INAF (Italy). We acknowledge the computing centre of INAF-Osservatorio Astronomico di Trieste, under the coordination of the CHIPP project (Bertocco et al. 2019; Taffoni et al. 2020), for the availability of computing resources and support.

References

Appendix A: Refinements to MUPPI

This section illustrates the effect of the two most relevant modifications made to MUPPI to achieve a cosmological simulation box exhibiting galaxy statistical properties consistent with observational data, as discussed in Section 2.2.

Wind particles are designed to sample galactic outflows. The details about their behaviour can be found in Murante et al. (2015) and Valentini et al. (2017). A gas particle exits its multi-phase stage after a maximum allowed time given by the free-fall time of the cold gas. When a gas particle exits a multi-phase stage, it has a probability P_kin (a model parameter) of being kicked and becoming a wind particle for a time interval t_wind. Wind particles are decoupled from hydrodynamic forces. While they are affected by radiative cooling, they cannot contribute to the artificial thermal conduction implemented in the P-Gadget SPH. Moreover, they receive kinetic energy from neighbouring star-forming gas particles, which increases their velocity along their least resistance path since they are kicked against their own density gradient. The wind stage can conclude before t_wind whenever the particle density drops below a chosen density threshold.

Figure A.1 shows the importance of scaling a particle duration in the wind phase t_wind with the velocity dispersion of its neighbouring DM particles, as in Equation 1. In the top panel, we present the SFRD for our reference model alongside two alternative runs for comparison. The latter two do not scale t_wind with σ_DM, instead adopting t_wind = 45 Myr − t_ff, c as in previous MUPPI runs, and in particular as in Parente et al. (2022). They differ in the adopted SF efficiency f_*. In the run with f_* = 0.01, the efficiency is the same as in the fiducial version of the model. However, the very low SFRD results in this solution being completely unsatisfactory. Therefore, we also show a simulation with f_* increased to 0.02, as in previous versions of MUPPI. The SFRD for this latter case approaches the observed SFRD, still without satisfactorily reproducing it. Further increasing f_* is not viable because, as evident in the middle and bottom panels, this assumption for computing t_wind results in a misrepresentation of how halos are populated by galaxies and an inaccurate count of galaxies of a given mass, issues that are not solved by a further increase in f_*.

Another notable improvement, aligning our model galaxies more closely with observational data, was achieved by increasing the feedback efficiency from low-metallicity stars (LMF). As outlined in Valentini et al. (2023) (and also in other numerical works as Maio et al. 2010; Schaye et al. 2015; Maio et al. 2016; Ma et al. 2017; Pillepich et al. 2018, and references therein), this boosted feedback was introduced to simulate the increased feedback energy expected from SNe generated by Population III stars. Specifically, the model star-forming particles with metallicities $\mathrm{Z}<0.05{\mathrm{Z}}_{\odot }$ produce feedback energy that is boosted with respect to the usual ${\mathrm{E}}_{\mathrm{SN}}={10}^{51}\mathrm{erg}$ by a factor determined by a model parameter, which we have doubled in this version of the code and set to 20. Fig. A.2 illustrates the impact of this choice in one of the simulated boxes. As expected, the SFRD with 10×LMF is higher than in our fiducial (20×) case, particularly at higher redshift. Additionally, it peaks at higher redshifts compared to the case with 20×. The middle panel shows that the excess of SF mainly occurs in galaxies populating small halos, which is also reflected in the excess of low-mass galaxy counts in stellar mass function in the bottom panel.

Fig. A.1. Cosmic SFRD from our reference simulation is compared with two alternative runs for one of the simulated boxes in the top panel. Both these runs do not scale the particle wind phase duration with σDM, as in the fiducial model, instead adopting ${\mathrm{t}}_{\mathrm{wind}}=45\mathrm{Myr}-{\mathrm{t}}_{\mathrm{ff},\mathrm{c}}$. They differ for the adopted SF efficiency f*, which is either 0.01 or 0.02. Observational determinations are the same as in Fig. 4. In the middle and bottom panels we show the SHMR and the SMF for the f* = 0.02 case.

Fig. A.2. Same as in Fig. A.1 but comparing the fiducial model with a run where the low metallicity feedback is reduced by a factor 2.

Appendix B: Estimate of shielding

To estimate both dust shielding and H₂ self-shielding, we adopt the schematic geometry depicted in Fig. B.1. The molecular mass is divided into giant molecular clouds (GMC), characterized by a typical mass M_MC and radius R_MC. At the resolution adopted in this work, even GMC are about one order of magnitude less massive than individual gas particles. The GMC are randomly distributed in the particle volume, which, consistently with the MUPPI model, is estimated as the ratio between the particle mass and the SPH density. The corresponding radius is R_p. Moreover, we assume that stars younger than an escape timescale t_esc ∼ a few Myr (henceforth young stars; in our computations, we set t_esc = 3 Myr ) are still embedded in the parent GMC whose medium practically completely absorbs their UV radiation; older stars (old stars) are distributed in the particle volume outside the GMC ⁹. Therefore, H₂ is affected in a different way by the radiation emitted by stars younger and older than t_esc.

Fig. B.1. Representation of the idealized molecular mass distribution inside a multi-phase gas particle. The destruction of H2 is carried out following this schematic. The basic idea is that the H2 mass (+helium) is organized into Giant Molecular Clouds (GMC), represented by the smaller spheres, with mass MMC and radius RMC. These GMC are considered randomly distributed inside the gas particle, the volume of which is estimated as the ratio between the particle mass and the SPH density, with a corresponding size Rp. Stars younger than a given tesc are considered to be still within the GMC, while older stars are considered to be outside.

The term accounting for destruction in Eq. 4 is given by the following integral over the volume V_aMC of all the N_MC molecular clouds:

$$\begin{array}{c}\hfill {\dot{M}}_{{\mathrm{destH}}_2}={\int }_{V_{\mathit{aMC}}}{\rho }_{{\mathrm{H}}_2}(r_g)C_{\mathrm{LW}}4\pi J_{\nu ,\mathrm{LW}}(r_g)d^3{r}_g,\end{array}$$(B.1)

where we use the subscript g (gas) for the dumb integration position, for reasons that will be clear in the following. We estimate separately the contribution to the LW radiation field J_ν, LW arising from the young stars inside the same cloud J_{ν, LW, y} and from old stars outside any cloud J_{ν, LW, o}. The former is given by the following integral over the volume V_MC of one of the identical MC:

$$\begin{array}{cc}\hfill J_{\nu ,\mathrm{LW},\mathrm{y}}& (r_g)=\hfill \\ \hfill & {\int }_{V_{1MC}}\frac{L_{\nu ,\mathrm{LW},1\mathrm{y}}}{V_{1MC}}\frac{1}{4\pi r_{\mathit{yg}}^2}exp(-{\tau }_{d,yg})S_{H_2}(N_{H_2,yg})d^3{r}_y.\hfill \end{array}$$(B.2)

Here, L_{ν, LW, 1y} is the total LW specific luminosity of all young stars inside a given MC. The first fraction represents their constant volume emissivity, estimated assuming a continuous smooth distribution of sources ¹⁰, r_yg ≡ ⃗r_y − ⃗r_g, τ_d, yg ≡ τ_d(|⃗r_y − ⃗r_g|) is the dust optical depth between the positions ⃗r_y and ⃗r_g, and S_H2(N_H2, yg) is the H₂ self shielding factor that depends on the H₂ column density N_H2, yg ≡ N_H2(|⃗r_y − ⃗r_g|). Using Eq. B.2 into Eq. B.1, we obtain the contribution of young stars to H₂ destruction:

$$\begin{array}{c}\hfill {\dot{M}}_{{\mathrm{destH}}_2,\mathrm{y}}=M_{{\mathrm{H}}_2}{C}_{\mathit{LW}}{L}_{\nu ,\mathrm{LW},1\mathrm{y}}\times \mathcal{A},\end{array}$$(B.3)

where

$$\begin{array}{c}\hfill \mathcal{A}=\frac{1}{V_{1MC}^2}\underset{V_{1MC}}{\int }\underset{V_{1MC}}{\int }\frac{exp(-{\tau }_{d,yg})S_{H_2}(N_{H_2,yg})}{r_{\mathit{yg}}^2}{d}^3{r}_g{d}^3{r}_y,\end{array}$$(B.4)

is the mean over all pairs of points in the MC of the integrand. To derive the former expression we used ρ_H2(⃗r_g) = M_H2/(N_MCV_1MC), and the fact that with our assumptions the integral appearing in Eq. B.1 over the volume of all MC is N_MC times the integral over the volume of one MC.

As for the H₂ self shielding factor S_H2, we adopt the approximations by Draine & Bertoldi (1996):

$$\begin{array}{c}\hfill S_{H_2}(N_{H_2})=min[1,{\left(\frac{N_{H_2}}{{10}^{14}{\mathrm{cm}}^{-2}}\right)}^{-0.75}].\end{array}$$(B.5)

We assume the H₂ column density of a GMC to be approximately N_H2 ≃ 1.15 × 10²¹ cm⁻² from the surface to the center. The column density is proportional to M_MC/R_MC². For example, this value is consistent with a GMC having a representative mass of M_MC = 10⁶ M⊙ and a radius of R_MC = 100 pc, assuming a helium fraction of 0.24. Indeed Miville-Deschênes et al. (2017) found that half of the MW molecular mass is contained in MC more massive than 8.4 × 10⁵ M_⊙, and a well-defined correlation between the mass and the radius of MC in the MW, yielding a corresponding radius of 96 pc and N_H2 ≃ 1.04 × 10²¹ cm⁻². We note, however, that the universality of MC physical properties in galaxies is debated due to observations suggesting a significant increase of N_H2 at higher redshift and SF activity Dessauges-Zavadsky et al. (2019). Were this the case, our conclusion that LW destruction is practically negligible for the purposes of the present paper would become even stronger.

The corresponding dust optical depth can be evaluated by recalling that our model of dust evolution approximates the size distribution with just two representative sizes: large grains (L), featuring a radius a_L ≃ 0.05 μm and small ones (S) with a_S ≃ 0.005 μm. For this section, we can adopt for both silicate and carbonaceous grains an intermediate value of the material density s ≃ 3gr cm⁻² (Granato et al. 2021; Parente et al. 2022). We approximate the grain cross-section at the LW wavelengths with the geometrical one πa². The column density N_d, a of grains of a given size a is related to N_H2 by

$$\begin{array}{c}\hfill \frac{4}{3}\pi a^{3}s{N}_{d,a}=\frac{1}{X}2m_H{N}_{H_2}{D}_a,\end{array}$$(B.6)

where D_a is the D/G mass ratio, and X is the H mass fraction. Therefore, the corresponding dust optical depth is

$$\begin{array}{c}\hfill {\tau }_{d,a}\simeq \frac{3m_H{N}_{H_2}{D}_a}{2sXa}\end{array}$$(B.7)

For our two size dust population this yields

$$\begin{array}{c}\hfill {\tau }_d\simeq 2.2(\frac{D_L}{0.01}+\frac{10D_S}{0.01})\left(\frac{N_{H_2}}{{10}^{21}{\mathrm{cm}}^{-2}}\right).\end{array}$$(B.8)

Using Eqs. B.5 and B.8, and the assumed N_H2 along the radius, we can then estimate by a Montecarlo procedure the mean 𝒜 (Eq. B.4), as a function of the dust to gas ratio(s) in each MP particle. In practice, we use a table of pre-computed values, well approximated by

$$\begin{array}{c}\hfill log(\mathcal{A}{R}_{\mathit{MC}}^2)=-0.03x^3-0.12x^2-0.17x-4.50,\end{array}$$(B.9)

where x = max(log τ_d, −2).

H₂ can also be affected by the destruction from LW radiation coming from outside the MC, generated by stars older than t_esc (but as discussed above, younger than ∼20 Myr; older stars do not contribute significantly to the LW region) spread in the whole particle volume. If L_{ν, LW, o} is the total LW luminosity of all old stars in the particle, and V_P is the particle volume outside MC, we can write an expression analogous to Eq. B.2 to account for the contribution to the radiation field inside a molecular cloud due to all the stars:

$$\begin{array}{cc}\hfill J_{\nu ,\mathrm{LW},\mathrm{o}}& (r_g)=\hfill \\ \hfill & {\int }_{V_P}\frac{L_{\nu ,\mathrm{LW},\mathrm{o}}}{V_P}\frac{1}{4\pi r_{\mathit{og}}^2}exp(-{\tau }_{d,og})S_{H_2}(N_{H_2,og})d^3{r}_o,\hfill \end{array}$$(B.10)

which inserted into Eq. B.1 yields:

$$\begin{array}{c}\hfill {\dot{M}}_{{\mathrm{destH}}_2,\mathrm{o}}=M_{{\mathrm{H}}_2}{C}_{\mathit{LW}}{L}_{\nu ,\mathrm{LW},\mathrm{o}}\times \mathcal{B},\end{array}$$(B.11)

where

$$\begin{array}{c}\hfill \mathcal{B}=\frac{1}{V_{\mathit{aMC}}^2{V}_P}{\int }_{V_{\mathit{aMC}}}{\int }_{V_P}\frac{exp(-{\tau }_{d,og})S_{H_2}(N_{H_2,og})}{r_{\mathit{og}}^2}{d}^3{r}_g{d}^3{r}_o.\end{array}$$(B.12)

This expression is similar to Eq. B.4, but now the mean refers to all pairs of points "og", the first representing an old star outside MC and the second one representing H₂ gas inside the specific MC. To estimate by Montecarlo method this average, we accept the simplifying assumption that the segment connecting the two points does not intersect other MC, and we neglect the contribution to τ_d, og from less dense dusty ISM outside the MC. Again, we notice that giving up to one or both of these simplifications would produce a greater shielding factor, reinforcing our conclusion that LW destruction is negligible for our purposes. We found that ℬ (Eq. B.12) is well approximated by

$$\begin{array}{c}\hfill log(\mathcal{B}{R}_{\mathit{MC}}^2)=-0.28x^2+(0.10y-0.69)x-1.71y-5.24,\end{array}$$(B.13)

where x = max(log τ_d, −1) and y = log(R_p/R_MC). This approximation holds for x < 2.5 and 0.5 < y < 2.5, ranges more than sufficient to cover the requirements of our simulations.

All Tables

All Figures

Fig. 1. Schematic view of the mass fluxes in a multi-phase gas particle. Gas cooling moves mass from the hot to the cold phase at a rate Ṁcool. Within the cold phase, an H2 reservoir begins to grow. ${\dot{M}}_{\mathrm{form}{\mathrm{H}}_2}$ is the H2 formation rate determined by the particle dust content and the number density of atomic hydrogen, which in turn is computed from the cold phase mass. ${\dot{M}}_{\mathrm{dest}{\mathrm{H}}_2}$ is the H2 destruction rate by LW photons, which transports mass from the molecular reservoir back to the cold phase. Ṁsf is the rate at which molecular gas creates a virtual stellar phase, which is contained within the gas particle until a stellar particle is stochastically spawned. Further destruction of the molecular mass occurs due to the action of massive stars and AGN, at rates Ṁev and ${\mathrm{f}}_{\mathrm{H}2}{\dot{M}}_{\mathrm{AGN}}$, respectively, where fH2 represents the fraction of H2 in the cold phase. AGN also moves material from the cold to the hot phase at a rate ṀAGN

In the text

Fig. 2. Two examples of disc galaxies from our simulations at z = 0. The four rows depict the mass surface densities of stars, gas, dust, and H2 on pixels of ∼0.5 kpc per side. Columns 1 and 2 illustrate face-on and edge-on views of a galaxy with M200 = 1.5 × 1012 M⊙, Mstar = 5.3 × 1010 M⊙, Mgas = 5.0 × 109 M⊙, Mdust = 7.0 × 107 M⊙, and MH2 = 2.2 × 109 M⊙. Columns 3 and 4 show the corresponding maps for a galaxy with M200 = 3.6 × 1012 M⊙, Mstar = 1.2 × 1011 M⊙, Mgas = 4.6 × 109 M⊙, Mdust = 6.3 × 107 M⊙, and MH2 = 2.0 × 109 M⊙.

In the text

	Fig. 3. Stellar, gas, dust, and H2 mass density maps in one of the (26 Mpc)3 simulation box computed in pixels of 200 × 200 kpc2.
In the text

Fig. 4. Average cosmic SFRD calculated from five ${(26\mathrm{Mpc})}^3$ simulation boxes, with the shaded region indicating the full dispersion of the curves. We compare with some determinations from literature: solid and dashed black lines represent the fits to the cosmic SFRD obtained by Madau & Dickinson (2014) and Davies et al. (2016) respectively, while black circles show the determination of Driver et al. (2018). Following their recommendation, the data from Madau & Dickinson (2014) have been scaled by a factor of 0.63 to convert from a Salpeter (1955) to a Chabrier (2003) IMF. The remaining data points require no further conversion.

In the text

Fig. 5. Stellar mass function is shown by the green curves at z = 0.1, z = 1, and z = 2, based on the five simulated boxes in our fiducial model. The error bars represent Poissonian errors for each stellar mass bin, typically smaller than the size of the point. The shaded green region encompasses the full dispersion of the individual SMFs corresponding to each of the five simulation boxes. We compare our model at z = 0.1 with the best-fit model of Rodríguez-Puebla et al. (2020); the grey area is an estimate of systematic errors which captures variations due to differences in the mass-to-light ratios and mass definitions (photometry) of the observational samples in that work. Simulations at z = 1 are confronted with the data from Ilbert et al. (2010) at two redshift ranges and at z = 2 with the data in Mortlock et al. (2011).

In the text

Fig. 6. Distribution of our simulated galaxies in the stellar mass-SFR plane is shown by the contours. Solid black lines are the MS determinations for observed galaxies presented in Popesso et al. (2023), resulting from an extensive homogenized compilation of studies published since 2014 (their Eq. (10)). The shaded grey band visualize the typical dispersion of determinations used to get their fit. The green curves connecting the points represent the median SFR per mass bin for a sub-sample of simulated galaxies selected using the same criteria as in Popesso et al. (2023), namely 8.5 < log(Mstar/M⊙) < 11.5 and 0.01 < SFR/( M⊙ yr−1) < 500. Bars delineate the 25–75% percentiles.

In the text

Fig. 7. Stellar mass-halo mass relation at z = 0 is assessed by considering galaxies across all our five simulated boxes, fb = 0.157 is the baryon fraction Ωb/Ωm. The stellar mass is computed within 0.1R200. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 and 0.2 dex-wide bins along the x and y axes, respectively. We compare our results with abundance matching estimate by Moster et al. (2013), and with the empirical model by Moster et al. (2018).

In the text

Fig. 8. Dust mass function at z = 0.1, z = 1, and z = 2 considering galaxies in the five simulated boxes are represented by green dots connected by solid lines. The error bars represent the Poissonian error for each dust mass bin. The shaded green region encompasses the full dispersion of the individual DMF corresponding to each of the five simulation boxes. The decline of the DMF at low masses is due to the imposed stellar mass resolution limit. We compare our results with observational data, represented by symbols and simulations delineated by lines. Observations at z = 0.1 correspond to Pozzi et al. (2020), Beeston et al. (2018), Dunne et al. (2011), and Vlahakis et al. (2005), at z = 1 (Pozzi et al. 2020) and at z = 2 (Pozzi et al. 2020 and Dunne et al. 2003). Simulations correspond to the works of Parente et al. (2022), Hou et al. (2019), Li et al. (2019), and McKinnon et al. (2017).

In the text

Fig. 9. Mean cosmic evolution of the molecular hydrogen density considering the five simulated boxes. The shaded region encloses the full dispersion of the five curves. The right-hand axis shows ΩH2, the H2 mass density normalised to the present-day critical density. Observational constraints are ASPECS (ALMA Spectroscopic Survey in the Hubble Ultra Deep Field, Decarli et al. 2020); VLASPECS (which uses the NSF’s Karl G. Jansky Very Large Array, VLA, covering part of the ASPECS footprint, Riechers et al. 2020); COLDz (VLA CO Luminosity Density at High Redshift in the COSMOS and GOODS-North fields combined; Riechers et al. 2019); PHIBBS2 (Plateau de Bure High-z Blue Sequence Survey 2, which uses the NOrthern Extended Millimeter Array – NOEMA – to observe field sources, Lenkić et al. 2020); xCOLD GASS (eXtended CO Legacy Database for the Galaxy All Sky Survey, Fletcher et al. 2021); ALMACAL-abs (blind search of CO absorbers against the ALMA Calibrator archive sources, Klitsch et al. 2019); ALMACAL-CO (CO emission-line survey Hamanowicz et al. 2023); Herschel PEP (PACS Evolutionary Probes Herschel survey in the GOODS-N -S fields, Berta et al. 2013); and UKIDS-UDS (UK InfraRed Telescope – UKIRT – Infrared Deep Sky Survey in the Ultra-Deep Survey, Garratt et al. 2021). Data listed in the top right legend do not use direct measurement of CO emission to estimate molecular masses.

In the text

Fig. 10. H2 mass function at z = 0.1, 1.0, and z = 2.0. The error bars represent Poissonian errors for each stellar mass bin, typically smaller than the size of the point. The shaded green region encompasses the full dispersion of the individual H2MFs corresponding to each of the five simulation boxes. In the local Universe, we report for comparison the H2 mass function obtained by Andreani et al. (2020), who use a fixed (COfix) or a luminosity dependent (CO-L) CO-to-H2 conversion factor and the best-fit model of Rodríguez-Puebla et al. (2020) and Fletcher et al. (2021). At higher redshifts, the boxes show constraints from Decarli et al. (2016) (ALMA ASPECS-Pilot survey), Decarli et al. (2019) (ALMA ASPECS-LP survey), and Riechers et al. (2019) (COLDz survey). The H2MF at redshift z = 0 is also shown in the middle and right panels using a thin dotted line to facilitate comparison.

In the text

Fig. 11. Molecular hydrogen mass as a function of stellar mass for galaxies in our simulation. In addition to contours, solid circles display the median $\mathrm{MH}2$ within bins of stellar mass, along with the 25%–75% percentiles for model galaxies with MH2 > 0. To facilitate comparison, the dotted line in the central and right panels corresponds to the z = 0 relation. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes. Empty circles in the left panel represent the xCOLD GASS sample reported in Saintonge et al. (2017); non-detection in this sample is represented by downward arrows (0.01 < redshift < 0.05).

In the text

Fig. 12. Depletion time as a function of stellar mass for model galaxies. In addition to contours, solid circles display the median depletion time within bins of stellar mass, along with the 25%–75% percentiles for model galaxies with $\mathrm{MH}2>0$. The colour contours delineate grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

In the text

Fig. 13. Integrated KS relation. The tail of the arrows represents ΣH2, while the head indicates ΣH2 + HI. We show the results for galaxies with stellar masses Mstar > 108 M⊙ and bulge-to-total ratios B/T < 0.5. The colour bar encodes the H2 fraction of each galaxy. Red/blue dashed line segments are linear fits in the log-log plane for galaxies with H2 fraction lower/greater than 0.5. The corresponding slopes are labelled in each panel for H2 and H2 + HI gas phases. Solid line represents one of the fits (slope n = 1.03) from de los Reyes & Kennicutt (2019) to the ΣH2-ΣSFR relation obtained from a compilation of observational samples of local spirals.

In the text

Fig. 14. Contours depict the H2-to-dust ratio as a function of gas metallicity for all model galaxies. MS galaxies ($\mid \log (\mathrm{SFR})-\log ({\mathrm{SFR}}_{\mathrm{MS}})\mid <0.5$) at the given redshift are represented by the circles, which depict the median H2-to-dust ratio within stellar mass bins, accompanied by the 25%–75% percentiles. The two lines labelled as G12 and B13 correspond to fits obtained in Bertemes et al. (2018) when using the CO-to-H2 conversion factors derived in Genzel et al. (2012) or Bolatto et al. (2013), respectively. The colour contours enclose grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

In the text

Fig. 15. Contours depict the H2-to-dust ratio as a function of stellar mass for all model galaxies. MS galaxies ($\mid \log (\mathrm{SFR})-\log ({\mathrm{SFR}}_{\mathrm{MS}})\mid <0.5$) at the given redshift are represented by the circles, which depict the median H2-to-dust ratio within stellar mass bins, accompanied by the 25%–75% percentiles. The colour contours enclose grids containing a specific number of galaxies as indicated in the colour bar, utilizing 0.2 dex-wide bins along both the x and y axes.

In the text

Fig. 16. Transition from atomic to molecular hydrogen in MP particles illustrated for the three analysed redshifts as a function of the sub-grid number density of the cold phase. The latter is assumed to have a temperature 300 K hence the pressure is proportional to the number density and shown in the top axis (k is the Boltzmann constant). The colour code shows the mean D/G mass ratio (top panel) and the gas phase metallicity (metals-to-gas mass ratio, bottom panel) of the gas particles falling within each pixel (150 pixels per axis). colour-coded circles represent the median H2 fraction and number density for particles within 0.5 dex bins of log(D/G) (top panel) and $\log (\mathrm{Zgas}/{\mathrm{Z}}_{\odot })$ (bottom panel). Meanwhile, stars indicate the same, but exclusively for particles in the transition region with H2 fraction between 0.4 and 0.6.

In the text

Fig. A.1. Cosmic SFRD from our reference simulation is compared with two alternative runs for one of the simulated boxes in the top panel. Both these runs do not scale the particle wind phase duration with σDM, as in the fiducial model, instead adopting ${\mathrm{t}}_{\mathrm{wind}}=45\mathrm{Myr}-{\mathrm{t}}_{\mathrm{ff},\mathrm{c}}$. They differ for the adopted SF efficiency f*, which is either 0.01 or 0.02. Observational determinations are the same as in Fig. 4. In the middle and bottom panels we show the SHMR and the SMF for the f* = 0.02 case.

In the text

	Fig. A.2. Same as in Fig. A.1 but comparing the fiducial model with a run where the low metallicity feedback is reduced by a factor 2.
In the text

Fig. B.1. Representation of the idealized molecular mass distribution inside a multi-phase gas particle. The destruction of H2 is carried out following this schematic. The basic idea is that the H2 mass (+helium) is organized into Giant Molecular Clouds (GMC), represented by the smaller spheres, with mass MMC and radius RMC. These GMC are considered randomly distributed inside the gas particle, the volume of which is estimated as the ratio between the particle mass and the SPH density, with a corresponding size Rp. Stars younger than a given tesc are considered to be still within the GMC, while older stars are considered to be outside.

In the text