© The Authors 2025

1 Introduction

In the context of the standard cosmological model, the formation of galaxies results from the condensation of gas within dark matter haloes, which gravitationally collapse from small density perturbations and later grow through the aggregation of matter and smaller substructures (Rees & Ostriker 1977; Silk 1977; White & Rees 1978; Blumenthal et al. 1984). Since the early stages of galaxy formation, several physical processes have been shaping the properties of galactic systems, acting as fundamental components of the intricate network of mechanisms unfolding across cosmic time. Gravitational collapse, gas cooling, mass accretion, star formation (SF), chemical contamination, and various forms of feedback act together, in an interconnected way, determining the properties of galaxies as they evolve (Dalgarno & McCray 1972; Efstathiou 1992; White & Frenk 1991). The occurrence of mergers, interactions, and mass accretion might also trigger morphological transformations and variations in the dynamical, chemical, and structural properties of galactic systems.

While disentangling the complex interplay of these processes and their individual effects on the evolution of galaxies remains elusive, numerical simulations developed over the last decades have proven instrumental in understanding how galaxies form and evolve within a cosmological context. Current simulations can produce systems that resemble the overall properties of real galaxies, including our Galaxy (Grand et al. 2017; Springel et al. 2018; Libeskind et al. 2020).

One of the crucial ingredients in simulations of the growth and evolution of galaxies is the proper modelling of SF and the corresponding feedback. Stars form as the result of the collapse and fragmentation of dense gas clouds in the centre of dark matter haloes along the discs of spiral galaxies. Initially triggered by gas cooling – which determines the amount and distribution of dense gas – SF works as a source of chemical elements and energy for the interstellar medium (ISM). The balance between cooling and the heating resulting from supernova explosions is key in determining the star formation rates (SFRs) over time and space. Additionally, factors such as the circulation of gas due to feedback, accretion of fresh gas from the intergalactic medium and galactic fountains, internal instabilities, and interactions and collisions influence the ISM properties in non-trivial ways. One of the most important advantages of cosmological simulations is that all these processes participating in the evolution of galaxies are naturally accounted for. However, neither SF nor stellar feedback can be included from first principles, as the scales at which they occur are not resolved. For this reason, sub-grid recipes are included to account for the impact of unresolved physics on larger scales, thereby allowing galaxy formation simulations to be run on cosmological volumes at high resolutions.

In the case of SF, models assuming a simple correlation between the SFR density and the gas density have been used since the pioneering studies of Katz (1992) and Steinmetz & Mueller (1994). In such models, cold dense gas transforms into star particles at a rate that is proportional to the gas density and inversely proportional to a typical timescale. This approach, if accompanied by an efficient modelling of feedback, predicts galaxies that reproduce the slope of the observed relation between the SFR and the gas density (Schmidt 1959; Kennicutt 1998), as well as the SFR levels, as long as a tuneable proportionality factor is introduced, known as the star formation efficiency (SFE), and the typical SF timescale is the dynamical or free-fall time of the gas (Scannapieco et al. 2006).

During the last decade, more recent observations made it evident that the correlation between the gas density and the SFR is even stronger when molecular gas is considered (Wong & Blitz 2002; Leroy et al. 2008; Bolatto et al. 2011; Schruba et al. 2011 and references therein), both at resolved scales (Baker et al. 2021) and at integrated (i.e. galaxy-wide) scales across different redshifts (Baker et al. 2022). Additionally, the typical timescale of SF may vary depending on the properties of the ISM where SF occurs (Bigiel et al. 2008, 2010; Bolatto et al. 2011; Leroy et al. 2013). Motivated by these results, a number of studies have developed more sophisticated SF models within simulations and semi-analytic models. These models rely on a better characterisation of the ISM (either increasing numerical resolution or developing better sub-grid models), opening up the possibility to follow the content of molecular hydrogen as a function of time and to introduce H₂-based SF laws (Murante et al. 2010; Mollá et al. 2015; Millán-Irigoyen et al. 2020).

In contrast to SF, modelling the effects of feedback in simulations has proven extremely complex. In part, the reason for this is that a detailed understanding of stellar feedback and its interaction with the surrounding medium has not been fully achieved, and discussion on numerical aspects of the modelling still remain (see e.g. Agertz et al. 2013). Moreover, different ways of implementing stellar feedback show significant variations in predictions for the stellar masses, morphologies, and other galaxy properties, even for a given initial condition (IC) (Scannapieco et al. 2012), although significant progress has been made in this direction (Kim et al. 2013, 2016). The modelling of feedback is naturally linked to the modelling of SF, and thus the various numerical parameters that are included to regulate SF and feedback levels are unique to a given model, lacking any corresponding physical parameter whose value can be estimated in experiments or observations. Finally, additional forms of feedback including black holes and cosmic rays have also been explored, as they also participate in the process of galaxy formation but may not be as important as stellar feedback for galaxies up to the Milky Way (MW) mass scale (Naab & Ostriker 2017).

While feedback is what ultimately regulates the SF activity in galaxies, it is clear that an adequate description of the ISM is key to determine the conditions for SF. Various sub-resolution models, which describe the multi-phase structure of the ISM, have been proposed, with the aim of providing a better modelling of SF in the context of simulations. Yepes et al. (1997) and Hultman & Pharasyn (1999) incorporated in Eulerian and Lagrangian codes the McKee & Ostriker (1977) model, assuming that each gas element was composed of a hot and cold phase in pressure equilibrium. Cold clouds form in the hot medium due to thermal instability and produce stars, which transfer mass back to the hot phase and can evaporate the cold clouds. Springel & Hernquist (2003) improved and extended these models in the GADGET code, also adding galactic winds driven by SF as a form of feedback.

As for numerical implementations including molecular hydrogen, in Pelupessy (2006) a sub-grid model for the formation of molecular gas was implemented, using a time-dependent model for the HI → H₂ transition. The model was used for an integrated study of SF and H₂ gas on a dwarf galaxy, in the setting of non-cosmological galaxy models. In the isolated galaxy formation simulations of Robertson & Kravtsov (2008), H₂ was tracked assuming equilibrium abundances based on the gas metallicity. Gnedin et al. (2009) implemented a model for tracking the non-equilibrium abundances of H₂ in their adaptive mesh refinement (AMR) code and ran cosmological simulations of low-metallicity galaxies, exploring the interactions between atomic gas, molecular gas, and dust. Christensen et al. (2012) adapted this model to incorporate it into their smoothed particle hydrodynamics (SPH) code and ran simulations of isolated dwarf galaxies in a cosmological context up to z=0. Murante et al. (2014) implemented a molecular-based SF law in cosmological simulations of a MW-type halo. In their model, SF is linked to the molecular abundance within gas particles computed using the phenomenological relation of Blitz & Rosolowsky (2006), which only depends on gas pressure. Valentini et al. (2022) used cosmological simulations of a MW-size halo, focusing on the low-metallicity regime, to compare the predictions of two different models for computing the molecular fraction: the theoretical model of Krumholz et al. (2009) and the phenomenological approach of Blitz & Rosolowsky (2006). They found that the former model produces a clumpier ISM and a more complex H₂ distribution in better agreement with observations of nearby disc galaxies. Hopkins et al. (2014) also assumed an H₂-regulated SF prescription for their galaxy formation simulation, computing the molecular hydrogen density with the local column density and metallicity as in Krumholz & Gnedin (2011).

In this work, we further explore the influence of assuming a molecular-based SF law on the formation of galaxies. We base our work on the semi-analytical model of Millán-Irigoyen et al. (2020), which assumes that the ISM has a multi-phase structure that includes atomic, ionised, and molecular gas, as well as metals and dust (which is especially important for molecular hydrogen formation in low metallicity environments at high redshift). Our model is grafted onto the moving-mesh, magnetohydrodynamical code AREPO (Springel 2010), which discretises the gas component into cells. We assume that each of these cells represents a multi-phase structure with ionised, atomic, and molecular hydrogen phases, as well as a stellar phase (note that metals are followed as well, using the standard chemical evolution model of AREPO). The evolution of each of these phases is governed by a set of differential equations, which quantify the amount of mass that is exchanged between them because of different physical mechanisms. The SFR associated with the stellar phase is linked to the fraction of molecular hydrogen and is regulated by a density-dependent SF timescale.

While our implementation of SF is explicitly tied to the molecular hydrogen fraction, we acknowledge that the causal role of molecular gas in SF remains debated. Numerical simulations by Glover & Clark (2012) have shown that SF can proceed efficiently in cold atomic gas, and that a cloud's ability to shield itself from the interstellar radiation field – rather than its chemical composition – is key to determine whether it can cool and collapse to form stars. Observationally, SF has also been detected in H i-dominated environments, including dwarf irregular galaxies with little or no detectable molecular gas (Roychowdhury et al. 2015). These findings suggest that H₂ may function more as a tracer of cold, dense gas than as a strict prerequisite for SF. Therefore, the assumption of an H₂-based SF prescription in our model should be interpreted with the caveat that it may not be universally valid, particularly in low-metallicity or poorly shielded environments.

This work is organised as follows. A description of the simulation code, the included physics, and the new SF model is provided in Section 2. In Section 3, we apply our model to the formation of a MW-mass galaxy in a cosmological context and discuss our predictions for the properties of the simulated system, focusing on the molecular, atomic, and ionised phases, and on the SF activity. In Section 4, we present the predictions of our model for the molecular Kennicutt-Schmidt (KS) law and for the SFE per free-fall time associated with gas cells, and in Section 5 we discuss the implications of assuming a molecular-based SF law. Finally, in Section 6, we give our conclusions.

2 The simulation code and the star formation model

In this section, we describe the implementation of our SF subgrid model, which follows the ionised, atomic, and molecular phases of hydrogen and links the SFR to the molecular phase. We first summarise the main characteristics of the base code where we implemented our model, and then, the model itself.

2.1 The AREPO code

We used the N-body, magnetohydrodynamics (MHD) code AREPO to perform the simulations, which we briefly describe here. For a more comprehensive explanation, we refer to Springel (2010). AREPO is a moving-mesh code that tracks MHD and collisionless dynamics in a cosmological setting. Gravitational forces are computed using a conventional TreePM method (Springel 2005), which employs a Fast Fourier Transform technique for long-range forces and a hierarchical oct-tree algorithm (Barnes & Hut 1986) for short-range forces, in conjunction with adaptive time-stepping. AREPO employs a dynamic unstructured mesh, constructed from a Voronoi tessellation of the simulation box. This allows for a finite-volume discretisation of the MHD equations. The MHD equations are solved using a second-order Runge-Kutta integration scheme with high-accuracy least-squares spatial gradient estimators of primitive variables (Pakmor et al. 2015).

A distinctive feature of AREPO is the ability to transform the mesh at any time step through a mesh reconstruction, a capability not found in standard grid-based methods. The mesh reconstruction ensures that each cell contains a specific target mass (within a certain tolerance), meaning that areas of high density are resolved with more cells than areas of low density. Additionally, the mesh generating points move with the fluid flow, allowing AREPO to overcome the Galilean non-invariance problem that standard Eulerian mesh codes have and significantly reduce the advection errors that appear in complex supersonic flows. The quasi-Lagrangian nature of the method makes it comparable to other Lagrangian methods such as SPH, although it eliminates several of its limitations, such as the need for artificial viscosity.

AREPO provides several modules to model different physical phenomena. In particular, we apply the galaxy formation model used in the Auriga Project, which includes primordial and metal-line cooling, a uniform ultraviolet background field

for reionisation, magnetic fields (Pakmor et al. 2014, 2017, 2018), energetic and chemical feedback from Type II supernovae, and mass loss and metal return owing to Type Ia supernovae and asymptotic giant branch stars (Vogelsberger et al. 2013; Marinacci et al. 2013; Grand et al. 2017). We did not include the active galactic nuclei feedback module, given that its effects are subdominant for a MW-mass galaxy.

Our code, including the SF, feedback and cooling routines, is identical to the standard AREPO implementation, except for the calculation of the SF probability associated with each gas cell (see next section). Using this probability, star particles are formed stochastically from gas cells that are eligible to do so. As we only vary the SF probability leaving the rest of the implementation unchanged, any difference in the results when our code is used will only be due to the changes induced in the SF and subsequent feedback levels, and their effects on the spatial distribution and properties of the ISM.

2.2 The star formation model

In our SF model, star particles are formed from gas cells following a stochastic approach. In this scheme, described in Springel & Hernquist (2003) (see also Grand et al. 2017), gas elements that are denser than a density threshold (ρ_th) are eligible to form stars, with a probability given by $$p=1-exp(-\frac{{\dot{m}}_{*}\Delta t}{M_{\text{cell}}}),$$(1) where ṁ_* is the SFR associated with the gas cell of mass M_cell and Δ t is the integration time step. The exponent in Eq. (1), ṁ_* Δ t/M_cell, thus represents the fraction of stellar mass formed in the cell during the time step.

In the standard formulation of AREPO, the SFR of a gas cell with density ρ and cold mass fraction x is given by $${{\dot{m}}_{*}|}_{\text{STD}}=\frac{x{M}_{\text{cell}}}{t_{*}\left({\rho }_{\text{cell}}\right)},$$(2) where t_*, the SF timescale, is $$t_{*}\left({\rho }_{\text{cell}}\right)=t_0^{*}{\left(\frac{{\rho }_{\text{cell}}}{{\rho }_{\text{th}}}\right)}^{-1/2},$$(3) and $t_0^{*}$, the maximum SF timescale, is a parameter of the model.

In our model, instead, the probability of forming a star is $$p=1-exp(-f_s),$$(4) where f_s is the fraction of stellar mass in a gas cell. This fraction is obtained by solving a system of differential equations, which follows four components: an ionised phase with a temperature of ∼ 10⁴ K, an atomic phase with a temperature of ∼ 100 K, a molecular component with a temperature of ∼ 10 K, and stars¹. The four components are parametrised as mass fractions relative to the total cell mass: f_X=M_X/M_cell where X can be i: ionised, a: atomic, m: molecular, or s: stellar. The mass of the cell at each time t is therefore $$M_{\text{cell}}\left(t\right)=M_i\left(t\right)+M_a\left(t\right)+M_m\left(t\right)+M_s\left(t\right),$$(5) and the sum of the corresponding fractions add up to unity $$f_i+f_a+f_m+f_s=1.$$(6)

It is important to emphasise that these fractions represent the internal composition of each gas cell within our sub-grid model, rather than physical particle masses in the simulation. In particular, the stellar mass fraction (f_s) does not correspond to actual star particles; instead, it is used solely to compute the instantaneous SFR of a gas cell, which in turn governs whether a new star particle is created. Moreover, in practice we find that f_a+f_m+f_i=1−f_s ≍ 1, since f_s remains small in all gas cells, typically of the order of 3 × 10⁻⁴. As illustrated in Fig. 1, the four phases exchange mass via the following physical processes: recombination of ionised atoms, condensation of atomic hydrogen, photodissociation of molecular gas, photoionisation of atoms, SF, and mass return from stars to the interstellar gas.

In terms of the mass fractions, the equations governing the evolution of the phases are $$\frac{\mathrm{d}}{\mathrm{d}t}{f}_i\left(t\right)=-\frac{f_i\left(t\right)}{{\tau }_{\text{rec}}\left(t\right)}+{\eta }_{\text{ion}}\psi \left(t\right)+R\psi \left(t\right),$$(7a) $$\frac{\mathrm{d}}{\mathrm{d}t}{f}_a\left(t\right)=\frac{f_i\left(t\right)}{{\tau }_{\text{rec}}\left(t\right)}-{\eta }_{\text{ion}}\psi \left(t\right)-\frac{f_a\left(t\right)}{{\tau }_{\text{cond}}\left(t\right)}+{\eta }_{\text{diss}}\psi \left(t\right),$$(7b) $$\frac{\mathrm{d}}{\mathrm{d}t}{f}_m\left(t\right)=\frac{f_a\left(t\right)}{{\tau }_{\text{cond}}\left(t\right)}-{\eta }_{\text{diss}}\psi \left(t\right)-\psi \left(t\right),$$(7c) $$\frac{\mathrm{d}}{\mathrm{d}t}{f}_s\left(t\right)=\psi \left(t\right)-R\psi \left(t\right),$$(7d) where τ_rec and τ_cond are the time-scales for recombination and condensation, respectively, η_ion and η_diss represent the efficiencies of ionisation and dissociation of hydrogen, respectively, R is the mass return fraction from stars, and ψ(t) is the fractional instantaneous SFR, $${\frac{\mathrm{d}}{\mathrm{d}t}{f}_s\left(t\right)|}_{\mathrm{SFR}}=:\psi \left(t\right)=\frac{f_m}{{\tau }_{\mathrm{star}}}.$$(8)

We assume it depends only on the amount of molecular gas and is regulated by a SF timescale (τ_star), as in Millán-Irigoyen et al. (2020).

The exchange of mass between phases due to recombination, condensation, and SF are modelled using timescales (τ_rec, τ_cond, and τ_star, respectively), which depend on the properties of the gas and are therefore time- and cell-dependent. In contrast, the terms corresponding to photoionisation, photodissociation, and mass return are proportional to ψ(t), with constant proportionality factors (η_ion, η_diss, and R, respectively). A detailed description and computation of all these factors is provided in Section 2.3.

Our system of equations was solved for each gas cell that is eligible to form stars, and integrated during the corresponding time step, with the assumption that the mass, density, and metallicity of each gas cell was constant during integration. Gas cells enter our routine with initial values for the ionised and atomic fractions given by the standard treatment of ionised and neutral phases in AREPO, provided that they have neither molecular nor stellar mass. It is worth noting that because AREPO is not purely Lagrangian, there might be mass fluxes between cells. If inflowing gas into a cell does not have the same relative mass fractions than the present cell, the f_X values should in principle be affected. For now, we ignored such fluxes and treated every cell as a 'closed box' during the integration of our set of equations.

It is important to note that the density threshold for SF, ρ_th, which determines whether a gas cell enters our SF routine at each time step, does not have a large impact on our results. While we adopted a fiducial value of 0.19 cm⁻³, as used in cosmological simulations of similar resolution, gas cells with densities below ρ_th would have a very low stellar fractions (as we show in Appendix A) and therefore a very low probability of producing star particles. We tested this with a simulation using 0.019 cm⁻³ for the ρ_th, and getting very similar results compare to our fiducial simulation. In other words, the minimum density of a gas cell needed to produce stars is encoded in our SF law via the link between SF and the amount of molecular gas in the cell.

Fig. 1 Physical processes responsible for transformations between the different phases. Arrows going inwards indicate that the named process increases the amount of the corresponding phase, and vice versa for arrows going outward.

2.3 Parameters of the model

2.3.1 Star formation timescale

For the SF timescale we use the free-fall time of the gas cell, as in Krumholz & McKee (2005) $${\tau }_{\mathrm{star}}=\sqrt{\frac{3\pi }{32G{\rho }_{\mathrm{g}}}},$$(9) where G is the gravitational constant and ρ_g, in the context of our model, is the density of the corresponding cell, ρ_cell. It is worth noting that we adopt the free-fall time of the gas cell as the star formation timescale, following the definition in Krumholz & McKee (2005). However, unlike their formulation, we do not include a separate efficiency parameter ε_ff; instead, SF is self-regulated in our model via the coupled evolution of the gas phases described by Eq. (7), which naturally determines the effective SFR.

2.3.2 Recombination timescale

The recombination rate of hydrogen atoms can be written as (Osterbrock & Ferland 2006; Gnedin et al. 2009; Christensen et al. 2012; Millán-Irigoyen et al. 2020) $${\frac{\mathrm{d}}{\mathrm{d}t}{n}_a|}_{\mathrm{rec}}=\frac{n_{\mathrm{i}}}{{\tau }_{\mathrm{rec}}}={\langle \sigma v\rangle }_{\mathrm{B}}{n}_e{n}_i,$$(10) where τ_rec is the recombination timescale and 〈σ v〉_B is the thermally averaged Case-B recombination cross-section (we do not consider transitions to the ground state). We take a value of 〈σ v〉_B=2.6 × 10⁻¹³ cm³ s⁻¹, which is appropriate for a ionised phase with a temperature of 10⁴ K (Nebrin 2023). n_e is the electron number density, which, by mass balance and within our model, is equal to the ionised hydrogen number density n_i. The number densities are calculated from the fractions as $$n_X=\frac{N_X}{V_{\text{cell}}}=\frac{f_X}{m_X}{\rho }_{\text{cell}},$$(11) where ρ_cell is the density of the gas cell, N_X is the number of particles of component X, and m_X is the mass of a single element of the X component.

In terms of the fraction of atomic hydrogen, f_a, the rate at which atomic hydrogen is created can be written as $${\frac{\mathrm{d}}{\mathrm{d}t}{f}_{\mathrm{a}}|}_{\text{rec}}={\langle \sigma v\rangle }_{\mathrm{B}}\frac{{\rho }_{\text{cell}}}{m_P}{f}_i^2=\frac{f_i}{{\tau }_{\text{rec}}},$$(12) where m_P is the proton mass and we defined the timescale τ_rec as $${\tau }_{\text{rec}}=\frac{m_P}{{\langle \sigma v\rangle }_{\mathrm{B}}{\rho }_{\text{cell}}{f}_i}.$$(13)

It is important to note that this is not a constant, as it depends not only on the density of the cell, but also on the ionised fraction f_i.

2.3.3 Condensation timescale

The rate of molecular hydrogen formation due to the condensation of atomic gas in the surface of dust grains can be written as (Hollenbach & Salpeter 1971; Hollenbach et al. 1971) $${\frac{\mathrm{d}}{\mathrm{d}t}{n}_m|}_{\text{cond}}=R_d{n}_H{n}_a,$$(14) where R_d is the dust grain H₂ formation rate coefficient, n_H is the hydrogen number density, and n_a is the atomic hydrogen number density.

In the literature n_H is generally taken as n_H=n_a+2 n_m, because most studies consider cold gas clouds (T ∼ 100 K) dominated by molecular and atomic hydrogen. In contrast, we have atomic, molecular, and ionised gas within our gas cell; so, we use n_H=n_a+2 n_m+n_i.

In terms of our adopted notation, we can write Eq. (14) as $${\frac{\mathrm{d}}{\mathrm{d}t}{f}_m|}_{\text{cond}}=2R_d(f_a+f_m+f_i)\frac{{\rho }_{\text{cell}}}{m_P}{f}_a.$$(15)

Similarly to our approach of the previous section, we can rewrite this equation as $${\frac{\mathrm{d}}{\mathrm{d}t}{f}_m|}_{\text{cond}}=\frac{f_a}{{\tau }_{\text{cond}}},$$(16) and define a condensation timescale, τ_cond, as $${\tau }_{\text{cond}}=\frac{m_P}{2R_d{\rho }_{\text{cell}}(f_a+f_m+f_i)},$$(17) which depends on the fractions of atomic, molecular, and ionised hydrogen.

Theoretical and experimental work has shown that R_d ∝ T^1/2 (Black & van Dishoeck 1987), and R_d ∝ n_dust, the number density of dust grains (Hollenbach et al. 1971). Assuming the simplest dust model, where the amount of dust is proportional to the metallicity Z, we have R_d ∝ T^1/2 Z (Pelupessy 2006; Wolfire et al. 2008). However, following previous prescriptions (Pelupessy 2006; Gnedin et al. 2009; Christensen et al. 2012; Mollá et al. 2017; Millán-Irigoyen et al. 2020), we assume that R_d only scales with the metallicity, using ∼ 100 K for the temperature of the cold neutral gas. We also normalise R_d to the solar metallicity (Wolfire et al. 2008) $$R_{\odot }=3.5\times {10}^{-17}{\mathrm{cm}}^3{\mathrm{s}}^{-1},$$(18) and thus we have $$R_d=\frac{R_{\odot }}{Z_{\odot }}Z,$$(19) where Z_☉=0.0127, the same value used within AREPO.

It is important to note that, in the case of uncontaminated gas (Z=0), R_d goes to 0 and thus, τ_cond diverges. This is solved by introducing an effective metallicity Z_eff, which is added to Z, resulting in the following expression for τ_cond $${\tau }_{\text{cond}}=\frac{m_P{Z}_{\odot }}{2R_{\odot }(Z+Z_{\text{eff}}){\rho }_{\text{cell}}(f_a+f_m+f_i)}.$$(20)

Adding an effective metallicity has also a physical meaning. If R_d ∝ Z, the underlying hypothesis is that the only conversion channel for HI → H₂ is the condensation on the surface of dust grains. Assuming instead R_d ∝ Z+Z_eff, we also consider molecular formation at very low metallicities. Following Glover & Jappsen (2007), we use Z_eff=10⁻³ Z_☉.

Finally, we use an additional adjustable parameter, the clumping factor C_ρ=〈ρ²〉 /〈ρ〉² (Gnedin et al. 2009; Christensen et al. 2012) with a fiducial value of C_ρ=100, in order to phenomenologically account for any potential uncertainties in τ_cond. So, we end up with $${\tau }_{\text{cond}}=\frac{m_P{Z}_{\odot }}{2R_{\odot }{C}_{\rho }(Z+Z_{\text{eff}}){\rho }_{\text{cell}}(f_a+f_m+f_i)}.$$(21)

It is worth noting again that τ_cond is not a constant value, as depends not only on the density, but on the atomic, molecular, and ionised fractions (or, equivalently, on the stellar fraction, as these four fractions add up to unity) and on the metallicity.

2.3.4 Photodissociation and photoionisation efficiencies

The disassociated and ionised mass rates per unit of created stellar mass are $${\eta }_{\text{diss}}=\frac{{\dot{M}}_{\text{diss}}}{\mathrm{SFR}}=2m_P{c}_{\text{diss}}\frac{{\dot{N}}_{\text{diss}}}{\mathrm{SFR}},$$(22) $${\eta }_{\text{ion}}=\frac{{\dot{M}}_{\text{ion}}}{\mathrm{SFR}}=m_P\frac{{\dot{N}}_{\text{ion}}}{\mathrm{SFR}},$$(23) where Ṅ_diss is the number of photodissociating photons produced per unit time (in the Lyman-Werner band, 912 Å to 1107 Å), Ṅ_ion the number of ionising photons produced per unit time (between 0 and 912 Å), and c_diss is an efficiency factor. For the molecular dissociation reaction, we consider the efficiency given by Draine & Bertoldi (1996), where it is shown that dust grains may absorb up to ∼ 60 percent of the photons capable of dissociating hydrogen molecules and a large fraction of the remaining photons excite different rotational and vibrational states, reducing their dissociation probability to ∼ 15 percent. We thus have c_diss=0.4 × 0.15=0.06. For the ionisation reaction, we assume 100% efficiency.

Using the number of photons produced per unit time, Q(t^′, Z^′), by a single stellar population of 1 M_☉, age t^′, and metallicity Z^′, we can calculate Ṅ as $$\dot{N}\left(t\right)={\int }_0^{t}SFR(t-t^{\text{'}})Q(t^{\text{'}},Z^{\text{'}})\mathrm{d}{t}^{\text{'}}$$(24) where SFR(t−t^′) is the instantaneous SFR at the moment of birth of the stellar population of age t^′, and Z^′=Z(t−t^′) is defined as the metallicity at that moment. Notice how the integral is in stellar age t^′, from 0 to the present time t. Because most of the contributions to the integral comes from young blue stars that die in the first 10 to 100 Myr (Q is up to 6 orders of magnitude larger for young stars than for older ones) it is possible to do the approximations $$SFR(t-t^{\text{'}})\approx SFR\left(t\right),$$(25) $$Z(t-t^{\text{'}})\approx Z\left(t\right),$$(26) if one assumes that the SFR and Z vary slowly in those timescales.

Then, the number of photons per unit mass and unit time can be written as $$Q(t^{\text{'}},Z)={\int }_{{\lambda }_1}^{{\lambda }_2}\frac{\lambda L_{\lambda }(t^{\text{'}},Z)}{hc}\mathrm{d}\lambda ,$$(27) where L_λ(t^′, Z) is the luminosity per unit of wavelength of a stellar population of 1 M_☉, age t^′, and metallicity Z, and therefore we have $${\eta }_{\mathrm{diss}}=2m_{\mathrm{P}}{c}_{\mathrm{diss}}{\int }_0^t{Q}_{\mathrm{LW}}(t^{\text{'}},Z)\mathrm{d}{t}^{\text{'}},$$(28) $${\eta }_{\text{ion}}=m_{\mathrm{P}}{\int }_0^t{Q}_{\mathrm{H}}(t^{\text{'}},Z)\mathrm{d}{t}^{\text{'}},$$(29) where the subscript LW indicates that we consider the Lyman-Werner band for dissociation (λ₁=921 Å and λ₂=1107 Å) and the subscript H indicates that the integration used for photoionisation is between λ₁=0 and λ₂=912 Å.

The photoionisation and photodissociation efficiencies (η_ion, η_diss) are evaluated locally within each gas cell as a function of its metallicity, rather than being derived from the star particles that may eventually form stochastically. Since our method does not include explicit radiative transfer, it does not capture photon propagation or escape, and should therefore be interpreted with this limitation in mind. Although the underlying physical processes are the ionisation and dissociation of hydrogen, the photon production rate is tied to recent SF. Accordingly, our ionisation rates are expressed per unit SFR and must be multiplied by ψ(t).

Using the values from PopStar (Mollá et al. 2009), we compute Q for the Chabrier (2003) IMF at the six metallicities for which we have yield tables −0.0001, 0.0004, 0.004, 0.008, 0.02, and 0.05 – and for ages between 0.1 Myr and 15 Gyr. These values are used as interpolation tables to compute η_diss and η_ion at a given metallicity and integration time, which we take as constants during the integration of Eqs. (7). The asymptotic value of η_diss ranges from 1353 for the lowest metallicity to 816 for the highest one, and η_ion ranges from 8489 for the lowest metallicity to 2360 for the highest one.

We note that our formulation of the dissociation and ionisation processes implicitly assumes that every photon available contributes directly to altering the local chemical state of the gas, i.e. that the local reservoir of molecular and atomic gas can absorb the full photon budget. This is clearly an approximation: in reality, not all photons result in ionisations or dissociations, and some may escape or be absorbed by dust. In fact, the constant efficiency factor c_diss=0.06 defined above partially accounts for this effect, as it represents the attenuation due to dust shielding. We verified that our assumption is valid in most relevant cases: in approximately 90% of the star-forming gas cells, the mass of atomic gas exceeds the amount required to be ionised given the computed value of η_ion. Thus, while the current formulation may overestimate feedback coupling in certain edge cases, it remains a reasonable approximation for the bulk of the star-forming regions in our simulations.

2.3.5 Mass recycling

We parametrise the mass returned from the stellar to the gas component using the parameter R, defined as the mass fraction of a stellar population that is returned to the ISM under the instantaneous recycling approximation (IRA), first explicitly used by Schmidt (1963). Note that this hypothesis can be relaxed considering the lifetimes of stars, as done in Millán-Irigoyen et al. (2020), which would make R time-dependent. In our simulations, this would produce an additional computational cost, which is not justified, given the short integration timescales over which we solve our system of differential equations.

The returned mass fraction R of a stellar population with masses between m+d m and Initial Mass Function (IMF) φ(m) is (Pagel 2009; Matteucci 2003) $$R=\frac{{\int }_{m_{\mathrm{ir}}}^{m_{max}}(m-m_{\mathrm{rem}}\left(m\right))\varphi \left(m\right)\mathrm{d}m}{{\int }_{m_{min}}^{m_{max}}m\varphi \left(m\right)\mathrm{d}m}$$(30) where m_min and m_max are the minimum and maximum masses of the IMF, m_ir is the mass limit for the IRA, and m_rem(m) is the remnant stellar mass given by the yield model. We use the stellar yields model of Portinari et al. (1998) compiled by Mollá et al. (2015), where the parameters have the following values: m_ir=8 M_☉, m_min=0.08 M_☉, and m_max=100 M_☉. We calculated R for the six metallicities for which we have yield tables (see previous section). R ranges from 0.16 for the lowest metallicity to 0.18 for the highest one. We use these values as a interpolation table to compute R, which remain constant during the integration of Eq. (7).

2.3.6 A comparison of the model timescales regulating mass exchange

The amount of mass exchanged between the different phases, which ultimately determines the SF levels, is regulated by the various timescales involved in the Eq. (7) and their dependencies on the properties of the gas. Fig. 2 compares the range of possible values for the timescales linked to SF (τ_star), condensation (τ_cond), and recombination (τ_rec). While all three decrease with increasing gas density (the colour scale), τ_cond has a secondary dependency on the gas metallicity, and τ_rec with the ionised fraction (note that this makes the set of Eq. (7) non-linear). As can be seen from Fig. 2, at a given density, τ_star is the largest of the three, τ_rec the smallest, and τ_cond lies in between. Even if we consider the dependencies other than the density, the conversion from ionised to the atomic material is almost always the fastest, followed by the transition from atomic to molecular gas, and finally from the molecular to the stellar phase. This behaviour is what determines the predictions of the model on individual gas cells, as we discussed in more detail in Appendix A.

Fig. 2 Comparison of the different timescales (τstar, τcond, and τrec) in our model. The colour scale indicates the cell density, and the y-axis the secondary dependency of each timescale when appropriate. τcond depends on the sum fa+fm+fi too (see Section 2.3), but for this plot we take fa+fm+fi=1.0.

2.4 Simulation setup and initial conditions

We use for our simulations the ICs of the halo denoted as Au6 in the Auriga Project (Grand et al. 2017). These correspond to zoom-in ICs of a halo that has, at z=0, a virial mass of M₂₀₀= 1.05 × 10¹² M_☉ and is relatively isolated. The mass resolution of the simulation (level 4 in Auriga) is ∼ 3 × 10⁵ M_☉ for dark matter particles and ∼ 5 × 10⁴ M_☉ for baryons.

The gravitational softening length for dark matter and star particles is fixed in comoving coordinates at 500 h⁻¹ pc up to z= 1, and fixed in physical units at 369 pc thereafter. The softening length of gas cells scales with the mean radius of the cell; the minimum and maximum softening lengths in physical units are 369 pc and 1.85 kpc, respectively. The cosmological parameters adopted are taken from the Planck Collaboration XVI (2014): Ω_m=0.307, Ω_b=0.048, Ω_Λ=0.693, and a Hubble constant of H₀=100 h km s⁻¹ Mpc⁻¹, where h=0.6777. The rest of the input parameters of the code are identical to those used in the Auriga Project (Grand et al. 2017).

As we show in the next section, our simulation, referred to as Au6_MOL, produces a galaxy that, at z=0, has a spiral structure, with an extended, rotationally supported disc and a bulge, in line with the expectations for the formation of a bulge-disc system in a galaxy, which has evolved quietly over several gigayears. This is consistent with Au6 being a galaxy with a quiet merger history and no significant merger events in its recent history (see Grand et al. 2017).

3 Formation of molecular hydrogen and evolution of the gas and stellar phases

The final stellar distribution of the simulated galaxy can be observed in Fig. 3, where we show stellar density maps for the face-on and edge-on projections (the z-axis is aligned with the stellar angular momentum). Au6_MOL presents a disc-like morphology, where ordered rotation (evident from the velocity field included in this figure) is observed up to a radius of about 25 kpc. The final morphology of the galaxy is the result of its SF history which, as shown in Fig. 4, grows rapidly during the first Gyr of evolution, reaching a maximum at about ∼ 6−8 Gyr, and declining smoothly up to z=0. The blue and red lines show separately the contributions of stars that, at z=0, are in the bulge and disc regions, respectively (separated here only by the stellar radius). At z=0, the SFR of the galaxy is 5.4 M_☉ yr⁻¹, and its final stellar mass is 6.5 × 10¹⁰ M_☉.

In our model, we assume that molecular gas is required to form stars, and thus the balance between H₂ formation and destruction at early epochs is what sets the time at which stars start to form. Fig. 4 shows that this occurs at ∼ 1 Gyr, indicating that the formation of molecular hydrogen has already started. The very first stars form from pristine gas, which (as we show in Appendix A) is highly inefficient at forming molecular hydrogen, even at high densities. However, after the collapse of the gas within the halo, high densities are achieved and a large amount of gas cells are eligible to form stars, entering our routine. Due to the stochastic nature of SF, some of these metal-free gas cells, even having a low molecular fraction, will be able to produce stars. The upper panel of Fig. 5 shows a histogram of the birth time of stars formed during the first 1.6 Gyr of evolution. We can observe that only 20 star particles form in the first Gyr, indicating that molecular hydrogen has not yet been formed in a significant amount.

After the first Gyr, a transition to significantly higher SF is found, with more than 3000 star particles forming in the following 600 Myr. This transition is driven by the chemical enrichment of the ISM produced by the death of the first stars, which increases the production efficiency of star particles as molecular hydrogen forms more efficiently at higher metallicities (due to the inverse dependence of the condensation timescale on metallicity, see Eq. (20)). The distribution of stellar metallicities for stars formed in the first 1.6 Gyr of evolution is shown in the lower panel of Fig. 5 (the colours correspond to the same time intervals of the upper panel). From this plot we can observe that stars formed in the first Gyr are either metal-free (10%) or with metallicities in the range 10⁻⁵−10^−2.5 Z_☉. In the case of the other time bins, we find that none of the newly formed stars have Z=0; and while some stars are still relatively metal-poor, even before 1.2 Gyr, many stars have already reached metallicities of at least 0.1 Z_☉, which is the typical metallicity where the transition to forming molecular hydrogen starts to be more efficient (see Appendix A).

As molecular hydrogen forms and SF is triggered, the simulated galaxy grows in a feedback-regulated regime, where the amount of gas in the different phases remains approximately constant, as we show in Fig. 6². The fact that neither the molecular mass nor the mass in the atomic/ionised components vary significantly after the first ∼ 3 Gyr of evolution, while SF continues, indicates an approximate equilibrium of mass exchange between the gas phases, consumption into stars, and mass return from SF. After the first 3 Gyr, the gas is distributed in the atomic, ionised, and molecular phases, respectively, with fractions of ∼ 63%, ∼ 33%, and <3%.

The growth of the galaxy during the early phases of its formation, where our model is more relevant, is illustrated in Fig. 7, where we show the stellar and the molecular, atomic, and ionised gas distributions at 2, 3, 4, and 5 Gyr of evolution. While molecular gas is quite extended already at 2 Gyr, SF is only efficient in the central, most dense regions of the galaxy. A well-defined, extended disc structure in the gas is already present at 4 Gyr. This structure is maintained up to the present time, allowing the formation of the stellar disc. The three gas phases have a similar spatial structure at all times, in the innermost regions, although the atomic gas is more extended than the molecular and ionised phases.

After 5 Gyr, the galaxy grows in a smooth fashion with a slowly declining SFR and, by z=0, the gas distributes itself in a disc-like configuration, as shown in Fig. 8 (left-hand panels). The three gas phases (subsequent panels) are spatially correlated in the disc region (note that we kept the same colour scale in all cases, for a better comparison). Molecular gas is only present in the densest regions of the disc, while the atomic and ionised phases are more extended, both radially and vertically. The stellar mass and the mass gas in the different components, both in the inner 40 kpc and up to the virial radius, are given in Table 1. In the table we also include the corresponding values for the R₉₅ parameter, defined as the radius, which encloses 95% of the total mass of any given component (see Iza et al. 2022)³. As evident when we compare Figs. 8 and 3, the gas disc has approximately twice the size of the stellar disc, which is confirmed by the values for R₉₅: 15.6 kpc for the stars and ∼ 35 kpc for the ionised gas phase. The R₉₅ value of the molecular gas (∼ 20 kpc) is in line with the stellar size.

The way in which the gas is distributed spatially in the different phases can be better seen from the surface density profiles shown in Fig. 9 and the profiles of the relative fractions of the three gas phases. From this plot it is clear that the atomic phase largely dominates the gas mass in the disc region and that the transition to the halo, where ionised gas starts dominating, occurs at ∼ 32 kpc. The molecular gas is always a small fraction of the gas, although it is worth noting that this component is constantly being turned into stars. As shown in Table 1 the relative proportion between atomic and ionised gas in the disc region is approximately 2.3:1, and the molecular gas comprises less than 3% of the total gas mass.

Fig. 3 Projected stellar mass distributions, in face-on and edge-on views, for the simulation Au6_MOL at z=0. The white arrows indicate the corresponding velocity field.

Fig. 4 Evolution of the SFR for the simulation Au6_MOL, computed at z=0. The red and blue lines show the contribution of stars with r ≤ 2 kpc (mostly bulge stars) and 2 kpc<r ≤ 40 kpc (disc stars), respectively.

Fig. 5 Upper panel: birth time distribution of the stars in the main sub-halo of simulation Au6_MOL at 1.6 Gyr. The colours indicate different birth time ranges. Lower panel: corresponding distribution of stellar metallicities. Metal-free stars are included artificially in the left-most bin. Metal-free stars are only formed during the first time bin (0.0 Gyr ≤ t<1.0 Gyr).

Fig. 6 Evolution of the mass and relative fractions of the different gas phases for the main sub-halo of simulation Au6_MOL up to a radius of 40 ckpc. We include the stellar mass and total gas mass for reference.

Fig. 7 Projected mass density of the stars and the molecular, atomic, and ionised gas, at 2, 3, 4, and 5 Gyr for the simulation Au6_MOL, in face-on views. The colour scales are fixed for each row to highlight the similarities and differences of the various distributions.

Fig. 8 Projected mass density of the total, molecular, atomic, and ionised gas for the simulation Au6_MOL at z=0, in face-on (upper panels) and edge-on (lower panels) views. The colour scales are fixed to highlight the similarities and differences of the various distributions.

4 Predictions for the SF law and SF efficiency

To test our SF model, in this Section we investigate its predictions for the molecular Kennicutt-Schmidt (mKS) relation, and for the SFEs. As shown in Fig. 10, the predicted mKS relation for our simulation Au6_MOL at z=0 is obtained using young stellar particles (up to stellar ages of 200 Myr). Even though in our model the KS law is not imposed (see Section 2), we find a well-defined correlation between the H₂ and SFR surface densities. Our data points cover more than 2 orders of magnitude in both axes, with the largest H₂ densities corresponding to the inner regions of the galaxy. The best-fit parameters associated with a linear fit to the data are also included in the figure.

In the figure we also include the resolved observations of Sun et al. (2023), who analysed the mKS relation across 80 nearby star-forming spiral galaxies at a resolution of 1.5 kpc. In both our model and the observations, the SFR surface density, Σ_SFR, increases with molecular gas surface density, Σ_H2, following a power-law trend. While the general trend obtained in our simulation is similar to the observational findings, we find a higher slope, of 1.16 ± 0.03, which lies within the upper end of the observational range reported by Sun et al. (2023) (0.88 to 1.21, depending on methodology). In particular, the simulation predicts higher values of Σ_SFR at Σ_H2 ≳ 6 M_☉ pc⁻². It is possible that our simplified treatment of gas shielding is partly responsible for this discrepancy; we will explore the effects of self-shielding on the mKS relation in a forthcoming paper.

The predicted distribution of SF efficiencies (SFE) for gas cells can also be used as a test for the model, to quantify the degree of agreement with available measurements. We adopt the usual parametrisation of the SFR of gas elements, i.e. $${\dot{\rho }}_{\star }={ϵ}_{\mathrm{ff}}\frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}},$$(31) where ρ_g is the gas density and t_ff is the free-fall time of the cloud, to obtain ε_ff, a quantity that can be physically interpreted as the cloud SFE per free-fall time. In the simulations, we can compute ε_ff for each cell, which represents the average SFE per free-fall time at the resolved scale. In fact, we cannot only obtain the predicted ε_ff distributions at all times, but we can also study possible dependencies on properties such as the gas density and metallicity, which are the two most important parameters of our SF prescription that correlate with the SFR (see Appendix A).

The distribution of ε_ff in our simulation Au6_MOL, at different times, is shown in Fig. 11. As discussed in the previous section, our SF model has a higher impact at early times, in the low-metallicity regime of the gas. The upper and lower rows in the figure include, respectively, the results obtained for three different ranges in metallicity and density. Several important results are revealed by this figure. First, the values of the SFE cover a wide range between ∼ 0.001% to ∼ 10%. As expected, the lowest values of ε_ff and the largest spread is found at t=2 Gyr. Second, there is no variation in the maximum possible values of ε_ff, which are always of the order of 10%. Third, the median values of the distributions (indicated in the figure with the blue arrows) are 0.4−0.6% for all times except for t=2 Gyr. Finally, the SFE increases with increasing metallicity and density (coloured lines). As discussed previously, in our model the SFR per cell depends on various parameters in a complex way, including the density and the metallicity (see also next section).

We note that the increasing trend of SFE with metallicity seen in Fig. 11 is not a prediction that emerges independently from the model's dynamics, but rather a natural consequence of our assumptions regarding H₂ formation. Specifically, the model assumes that molecular hydrogen formation is more efficient in higher-metallicity environments due to more abundant catalytic sites, leading to higher molecular fractions and thus higher SFE. While the absolute values of SFE are a prediction of the model, the qualitative dependence on metallicity reflects this built-in assumption.

We can also calculate the values of ε_ff for gas cells that produced star particles, right before this occurs. In this case, we find that, as expected, they belong to the high tail of the distributions as they have higher probability of forming stars. The median values of the distributions for gas cells that produced stars are indicated by the red arrows in Fig. 11, and are between 3−4% at all times, except for t=2 Gyr where the median value is 0.9%.

These results for the SFE per free-fall time in our simulation agree very well with observations from giant molecular clouds, which suggest that ε_ff varies from 0.1% to 10% (e.g. Krumholz et al. 2019; Schinnerer & Leroy 2024), and resolved measurements in nearby star-forming spiral galaxies and in the MW, with values of the order of ∼ 0.7% (Utomo et al. 2018; Sun et al. 2023). Moreover, both the variations and values obtained with our model are consistent with those found in SF simulations in a turbulent ISM (Semenov et al. 2016), when averaged over the resolved scales.

Fig. 9 Upper panel: projected mass density profiles for the gas and stellar components of the simulation Au6_MOL at z=0. We show separately the contributions of the ionised, atomic, and molecular phases, and include the stellar component, for comparison. Lower panel: radial profiles of the global mass fraction of the ionised, atomic, and molecular gas components, normalised to the total mass of gas.

Fig. 10 Molecular Kennicutt-Schmidt relation, at z=0, for our simulations Au6_MOL (black dots). Each point corresponds to a 1.5 kpc × 1.5 kpc parcel of a 30 kpc square region centred on the main sub-halo. We also include a linear fit (blue line) for our data points and the observational values of Sun et al. (2023) (orange).

Fig. 11 Distribution of the SFE per free fall time of the gas in the disc (r<40 ckpc) of simulation Au6_MOL, at 2, 6, 10, and 13.8 Gyr. The first row shows the distribution separated by metallicity ranges (the colour lines show the ranges while the black lines show the full distribution), and the second row shows the distribution separated by density ranges (ρth in the legend is the density threshold for SF). The blue arrows indicate the median of the total distribution (black line), while the red arrows indicate the median of the distribution for gas that has been transformed into stars (this distribution is not shown).

5 Discussion

Due to the complex cycle between gas, SF, and feedback in galaxies, the development of a physically motivated model for SF that works at the resolved scales (ideally capturing also our knowledge on smaller, unresolved scales) is not an easy task. Moreover, the tuning of the input parameters associated with SF cannot be done separately from other parameters of the model, most notably those related to feedback, which provides SF regulation and sets the SF levels. Our model has input parameters within the system of equations, but most of them can be obtained from observations or theoretical considerations (see Section 2.3). As a first test to our SF model, we compared the molecular KS law and the distribution of SFE per free-fall time with several observations.

In this section, we compare the results of our model to those of two other runs that use alternative SF laws, commonly adopted in zoom-in simulations of galaxy formation. The differences between the three models are in the estimation of the SFR of gas cells, which alters the probability of forming stars, i.e., $$p=1-exp(-F_s),$$(32) where F_s represents the fraction of stellar mass in the gas cell that was formed during time step Δ t (see Section 2.2). The assumptions for F_s in the different models are as follows:

• Au6_MOL assumes that F_s=f_s where f_s is the stellar mass fraction obtained using our system of equations.

• Au6_STD uses the standard AREPO formulation (as adopted in the Auriga Project) where F_s depends on the cold mass fraction (x), and the assumed SF timescale (t_*, with a maximum SF timescale of $t_0^{*}=2.2\mathrm{Gyr}$), as described in Section 2.2: $$F_s=\frac{x}{t_{*}}\Delta t$$(33)

• Au6_BLT uses a molecular-based SF probability, which depends on an estimated molecular fraction f_m per cell, obtained using the observationally motivated relation of Blitz & Rosolowsky (2006), i.e.: $$f_m=\frac{1}{1+P_0/P},$$(34) where P is the gas pressure of the cell and P₀/k=3.5 × 10⁴ K cm⁻³. The F_s value in this model is $$F_s=f_{*}\frac{f_m}{t_{\mathrm{ff}}}\Delta t$$(35) where f_* is an efficiency parameter and t_ff is the free-fall time of the cell. Following Murante et al. (2010, 2014), who adopt a similar approach, we take f_*=0.02 (see also Valentini et al. 2022).

The comparison between the three simulations are summarised in Fig. 12 where we show the evolution of the SFR and of the gas mass in the different components and in Fig. 13 where we show the z=0 distributions of the stars and the ionised, atomic, and molecular gas (note that there is no estimation of the molecular mass for Au6_STD). Looking at the top panel of Fig. 12, the SFRs of the Au6_MOL and Au6_STD models are very similar, while Au6_BLT lies below. More important is the fact that the two models using H₂-dependent SF laws show a delay in the onset of SF of about 0.5 Gyr with respect to Au6_STD, which results from the need to form molecular gas first before forming stars. In our model, this is a consequence of the low enrichment levels of the ISM at these early times. The evolution of the mass in the ionised, neutral (i.e. atomic plus molecular in the Au6_MOL and Au6_BLT cases), and molecular phases also presents some differences between the models, even though the z=0 values are very similar in all cases. We also observe that Au6_BLT produces a lower amount of molecular gas in comparison with our model, which explains the observed differences in the SFRs.

As it can be seen in Fig. 13, the galaxies end up with similar morphologies but somewhat different disc sizes for the stellar and gas components. Interestingly the molecular gas is more extended in Au6_BLT than in Au6_MOL, whereas the corresponding stellar discs show the opposite trend. In contrast to our model, where the amount of H₂ depends on the density and metallicity, in Au6_BLT the H₂ fraction depends only on the gas pressure. The structure of the gas in the different components also shows some differences that are more evident in the ionised phase (which is consistent with the results of Valentini et al. 2022, but see also Christensen et al. 2012).

In the left-hand panel of Fig. 14, we compare the evolution of the median SFE per free-fall time (ε_ff) for gas cells, also including the shaded range from the 25% to the 75% percentiles of the data. In the case of Au6_STD, ε_ff depends on the cold mass fraction and the assumed SF timescale (Section 2.2), while in Au6_BLT it depends on the fraction of molecular hydrogen and the assumed SF efficiency parameter (f_* in Eq. (35)). In Au6_STD, ε_ff stays approximately constant during the whole evolution, at a value of 3%. This results from the fact that there are no significant variations in the cold gas fraction of star-forming gas cells (note that, although unnoticeable at the scale of the plot, there is a non-zero spread in ε_ff). In contrast to Au6_STD, the two models where the SF law is based on H₂ show a noticeably evolution of the typical SFE, as well as a large spread, most notably in the case of our Au6_MOL model.

In the middle- and right-hand panels of Fig. 14 we further explore how the SFEs depend on the properties of the gas, showing ε_ff as a function of gas density and metallicity, at z=0 (similar results are found for other redshifts). Both Au6_MOL and Au6_BLT present a positive correlation between the SFE and density/metallicity, and these correlations are significantly stronger in Au6_MOL compared to Au6_BLT, even though both simulations show similar SFEs, particularly in the last 6 Gyr of evolution (left-hand panel of this figure). This is explained by the fact that in Au6_BLT, there is a direct relation between ε_ff and gas density, because in this model the molecular fraction only depends on the pressure; and the dependence on metallicity results from correlations between the gas density and the metallicity.

We note that, additionally to the use of the molecular KS law, the distribution, evolution, and correlations of the SFE per freefall time predicted by different models could be useful as a test for SF routines in simulations where the structure of the ISM is not resolved. This might be particularly useful in the case of our model, where the SF routine can be easily adapted to link SF to the molecular gas, the atomic gas (likely precursor of SF at high redshift), or a combination of both.

Fig. 12 Comparison of the SFR and evolution of the gas mass in the ionised, neutral, and molecular components between the three simulations with different SF laws: Au6_MOL (our model), Au6_BL (an alternative H2-based SF law), and Au6_STD (the standard AREPO implementation, where SF is linked to the total gas density). Only gas cells within a sphere with r=40 ckpc are considered.

6 Conclusions

In this work, we presented a model to describe the various components of hydrogen – ionised, atomic, and molecular – in the context of numerical simulations of galaxy formation and analysed the effects of linking SF to the molecular phase. Our implementation aims at providing a more realistic model for SF in simulations, particularly to allow a varying SFE of the gas cells, which can be linked to the local properties of the gas. The model is based on the formulation of Millán-Irigoyen et al. (2020) and grafted onto the cosmological, magnetohydrodynamical moving-mesh code AREPO.

In our sub-grid model, we assumed that each gas cell is composed of ionised, atomic, and molecular gas, as well as stars. These components exchange mass through recombination (ionised to atomic), condensation (atomic to molecular), SF (molecular to stellar), photodissociation (molecular to atomic), photoionisation (atomic to ionised), and mass return (stellar to ionised). While the modelling is relatively simple, it allows us to account for variations in the corresponding efficiencies and typical timescales of the various processes, according to the properties of the local ISM. In our model, the SFR per gas cell is assumed to be directly proportional to the amount of H₂, with a typical timescale that depends on gas density, and with no additional factors (i.e. we do not introduce any particular SFR efficiency parameter). Although we change the SF prescription to link the SFR to the amount of H₂, we adopt the standard AREPO modelling for primordial and metal-line cooling, reionisation from a uniform ultraviolet background field, magnetic fields, energy feedback from Type II supernovae, and chemical enrichment via supernovae II and Ia and AGB stars. However, in this work, we switched off the effects of active galactic nuclei, as their effects are subdominant for MW-mass galaxies.

We applied our model to a cosmological zoom-in IC of the formation of a MW-mass galaxy, taken from the Auriga Project (Grand et al. 2017). In order to compare our model with other prescriptions, we also ran two additional simulations, using the standard SF law of AREPO, and an H₂-based SF model where the abundance of molecular hydrogen is obtained from an observationally motivated relation with the local gas pressure (Blitz & Rosolowsky 2006).

Our model predicts the formation of a spiral-like galaxy, with a young, rotationally supported disc component, and a small bulge. The gas disc is twice as large as the stellar one and composed mainly of atomic (∼ 63%) and ionised (∼ 33%) material. The molecular gas is always a small fraction of the total gas (≲3%), which results from a rapid consumption of H₂ into stars. The relative fractions of atomic and ionised material in the disc remain approximately constant during the whole evolution, even at times between 1 Gyr and 3 Gyr where the SFR changes rapidly. In our model, the formation of molecular hydrogen is not efficient for metal-free gas, even for timescales larger than 1 Gyr. However, the efficiency increases rapidly as the ISM gets chemically enriched at levels of the order of 0.1 Z_☉, following an increase in the amount of molecular hydrogen. This allows the formation of an extended gas disc of molecular material.

Our model predicts a correlation between the molecular gas and SFR density, similar to the molecular KS law. It is worth noting that this relation is not imposed in our model but results naturally from the system of equations of our sub-grid treatment. The correlation between the Σ_H2 and Σ_SFR in Au6_MOL is somewhat steeper than the observational estimates obtained from the PHANGS-ALMA sample (Sun et al. 2023).

An advantage of our model is that gas particles can have varying SFEs that depend on the local ISM properties rather than being an (imposed) fixed value. We find that the SFE per free fall time varies in the range ∼ 0.1% to ∼ 10% with typical values of 0.4−0.6%, except for very early times with a median of 0.08% (measured at t=2 Gyr). The SFE per free-fall time is an increasing function of metallicity and density and shows a mild evolution for times ≳ 3−4 Gyr. The obtained values for the SFEs are consistent with those derived from observations of molecular gas and from simulations of a turbulent ISM.

We compared our model to the standard implementation linking SF to total gas density and to an alternative approach tying SF to H₂ abundance, assuming a scaling with local gas pressure. The two simulations with H₂-based SF laws show a delay in the SF onset of about 0.5 Gyr with respect to the standard formulation. We find different distributions for the SFE per free-fall time of gas in the three runs, as well as different dependencies (if at all) on the properties of the ISM that might potentially be tested against observations.

In particular, our model does not impose any correlation between the SFR and the total or molecular gas surface densities and does not require an ad hoc efficiency parameter. In the current implementation, however, the molecular hydrogen is considered a precursor for SF, although this can be easily changed in the models where SF is linked to the atomic gas or a combination of the two components. Moreover, the dust model included here is extremely simple, and it will be extended in future work allowing for a better treatment of dust and its relation to molecular hydrogen and SF. Additional effects such as photodissociation of H₂ by metagalactic Lyman-Werner radiation also need to be considered. Overall, our sub-grid model allows us to link the SF to the local ISM properties at the resolved scales, while being simple enough for simulating galaxy formation in a cosmological context, ignoring the high complexity of the ISM at ≲100 pc scales.

Fig. 13 Projected mass densities of the stellar, ionised, neutral, and molecular gas components at z=0 for the three simulations with different SF laws: Au6_MOL (our model), Au6_BL (an alternative H2-based SF law), and Au6_STD (the standard AREPO implementation, where SF is linked to the total gas density).

Fig. 14 Left: time evolution of the median of εff for the three simulations, together with the median absolute deviation (shades). Middle: correlation between the SFE parameter and the gas density at z=0. Right: correlation between the SFE parameter and the gas metallicity at z=0.

Acknowledgements

We thank the anonymous referee for their helpful comments and suggestions, which improved our paper. EL acknowledges funding from the HORIZON-MSCA-2021-SE-01 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No. 101086388 (LACEGAL-III: Latin American-Chinese-European Galaxy Formation Network). CS acknowledges stimulating and interesting discussions with Andi Burkert, Vadim Semenov, Andrey Kravtsov, and Oscar Agertz. CS and SEN are members of the Carrera del Investigador Científico of CONICET. They acknowledge funding from Agencia Nacional de Promoción Científica y Tecnológica (PICT-201-0667).

Appendix A Predictions of the model for fixed initial conditions

In this section, we describe the behaviour of Eqs. (7) for fixed ICs, which cover the range of parameters (density, metallicity, and relative fraction of ionised and atomic gas) typical of our galaxy formation simulations and the adopted resolution (which sets the maximum gas densities achieved). A first approach to understanding the behaviour of the model equations which govern the evolution of the gas and stellar phases was given in Section 2.3.6, where we compared the typical timescales involved in the transformation of mass between the different phases. As explained in the text, the three timescales of our model – i.e. recombination, condensation, and SF – have a similar dependency on the density (τ ∝ ρ⁻¹ for recombination and condensation, and τ_star ∝ ρ^−1/2 for SF); the recombination timescale also depends on the ionised fraction, and the condensation timescale has additional dependencies on the gas metallicity and the stellar fraction (which sets the value of f_i+f_a+f_m, see Eq. (20)). Furthermore, η_diss, η_ion, and R depend on the gas metallicity.

It is important to note that, in a simulation, as we integrate our system of equations in time, the timescales vary. In the results shown in this section, the density and metallicity are kept fixed during integration, unlike τ_rec and τ_cond, which depend on the variables, ionised and stellar fractions, respectively. In Fig. A.1 we show the time-evolution of the various fractions, for different initial metallicities and densities (see Fig. A.2 for understanding the selected density/metallicity values shown here). We show results up to 1 Gyr, to understand the predictions of the model up to times where star particles start to form at various physical conditions. It is worth noting that in our cosmological simulation, all the gas cells that have been transformed into stars were integrated for a total time <10 Myr.

Let us first focus on the solid lines, which show a case where the initial atomic and ionised fractions are equal (with the stellar and molecular phases set to zero). The Fig. A.1 shows a first transition, from ionised to atomic material, at times determined by the recombination timescale τ_rec (see Fig. 2). This timescale is inversely proportional to the density, explaining why the transition occurs earlier in the right-hand panels (i.e. the highest density ones) and later for the left-hand ones. A second transition occurs from the atomic to the molecular phase, regulated by τ_cond. In this case, the transition occurs earlier as we go to higher densities, where the transition also becomes more efficient. This results in a significant production of stellar mass. For the metallicity dependence, as the gas gets enriched (from the top to the bottom rows) only τ_cond is affected, such that higher molecular/stellar fractions are obtained. Once a stellar phase is formed, the additional dependencies regulating the processes of photodissociation, photoionisation, and mass recycling will make our system of equations even more complex.

Finally, the dashed lines in this figure show what happens if we start from a different IC where almost all material is in the atomic phase (f_a(0)=0.99. and f_i(0)=0.01). In this case, we see that results converge quite fast to the same values as the case where the initial atomic and ionised phases are 0.5. This means that our system of equations has low sensitivity to the ICs (not so much with the parameters ρ_cell and Z), and in all cases the stellar fraction starts growing only after this convergence takes place.

In Fig. A.1 we showed the predictions for our model for a set of fixed ICs, and discussed that the initial ionised/atomic fractions are unimportant considering the relevant timescales of the model. In Fig. A.2 we better show how predictions depend on the gas density and metallicity. The range of densities has been selected between the density threshold for SF (which determines the minimum density a gas cells must have to form stars, and thus to enter our routine) and the typical maximum density seen in our simulation (which depend on numerical resolution). In the case of metallicity, we cover from metal-free to highly enriched gas. We present in Fig. A.2 the fractions of ionised, atomic, molecular, and stellar phases, as a function of density and for various metallicities (note the different y-ranges of the various rows, which were chosen for the sake of clarity). The three columns correspond to different integration times of 1, 10, and 100 Myr.

The behaviour of the different phases as a function of density is determined by the dependence of all timescales on density: the larger the density, the smaller the timescales; depending on the time observed, different transitions have already taken place. In general, the ionised and atomic fractions are more efficiently transformed into the molecular/stellar phases for higher densities, and this effect is enhanced as the system evolves for longer times. The dependence on metallicity is such that metal-free gas is unable to produce a significant amount of stellar mass; however, as the metallicity reaches about 10⁻¹ Z_☉, the fraction of molecular gas increases significantly, and so does the fraction of stellar mass (which in turn, consumes the molecular gas).

The results of this appendix show that, depending on the time-evolution of our equations, equilibrium solutions can be valid for a given period of time; however, as the ICs of each gas cell vary with time, as well as the timescales involved in the equations (which depend on the properties at the time a particle enters the routine), it is not trivial to predict the impact of our SF prescription in a cosmological simulation.

Fig. A.1 Evolution of the gas phases when integrating Eqs. (7) up to t=1 Gyr. Different panels corresponds to different values of gas metallicity and density, as indicated above each panel. Solid lines correspond to the ICs fi=fa=0.5, and dashed lines with fi=0.01 and fa=0.99. In both cases the initial values for fm and fs are 0.

Fig. A.2 Final value of the different phases (rows) when integrating Eqs. (7) for 1, 10, and 100 Myr (columns). Each colour corresponds to a different value of metallicity.

References

All Tables

All Figures

	Fig. 1 Physical processes responsible for transformations between the different phases. Arrows going inwards indicate that the named process increases the amount of the corresponding phase, and vice versa for arrows going outward.
In the text

	Fig. 2 Comparison of the different timescales (τstar, τcond, and τrec) in our model. The colour scale indicates the cell density, and the y-axis the secondary dependency of each timescale when appropriate. τcond depends on the sum fa+fm+fi too (see Section 2.3), but for this plot we take fa+fm+fi=1.0.
In the text

	Fig. 3 Projected stellar mass distributions, in face-on and edge-on views, for the simulation Au6_MOL at z=0. The white arrows indicate the corresponding velocity field.
In the text

	Fig. 4 Evolution of the SFR for the simulation Au6_MOL, computed at z=0. The red and blue lines show the contribution of stars with r ≤ 2 kpc (mostly bulge stars) and 2 kpc<r ≤ 40 kpc (disc stars), respectively.
In the text

	Fig. 5 Upper panel: birth time distribution of the stars in the main sub-halo of simulation Au6_MOL at 1.6 Gyr. The colours indicate different birth time ranges. Lower panel: corresponding distribution of stellar metallicities. Metal-free stars are included artificially in the left-most bin. Metal-free stars are only formed during the first time bin (0.0 Gyr ≤ t<1.0 Gyr).
In the text

	Fig. 6 Evolution of the mass and relative fractions of the different gas phases for the main sub-halo of simulation Au6_MOL up to a radius of 40 ckpc. We include the stellar mass and total gas mass for reference.
In the text

	Fig. 7 Projected mass density of the stars and the molecular, atomic, and ionised gas, at 2, 3, 4, and 5 Gyr for the simulation Au6_MOL, in face-on views. The colour scales are fixed for each row to highlight the similarities and differences of the various distributions.
In the text

	Fig. 8 Projected mass density of the total, molecular, atomic, and ionised gas for the simulation Au6_MOL at z=0, in face-on (upper panels) and edge-on (lower panels) views. The colour scales are fixed to highlight the similarities and differences of the various distributions.
In the text

Fig. 9 Upper panel: projected mass density profiles for the gas and stellar components of the simulation Au6_MOL at z=0. We show separately the contributions of the ionised, atomic, and molecular phases, and include the stellar component, for comparison. Lower panel: radial profiles of the global mass fraction of the ionised, atomic, and molecular gas components, normalised to the total mass of gas.

In the text

	Fig. 10 Molecular Kennicutt-Schmidt relation, at z=0, for our simulations Au6_MOL (black dots). Each point corresponds to a 1.5 kpc × 1.5 kpc parcel of a 30 kpc square region centred on the main sub-halo. We also include a linear fit (blue line) for our data points and the observational values of Sun et al. (2023) (orange).
In the text

Fig. 11 Distribution of the SFE per free fall time of the gas in the disc (r<40 ckpc) of simulation Au6_MOL, at 2, 6, 10, and 13.8 Gyr. The first row shows the distribution separated by metallicity ranges (the colour lines show the ranges while the black lines show the full distribution), and the second row shows the distribution separated by density ranges (ρth in the legend is the density threshold for SF). The blue arrows indicate the median of the total distribution (black line), while the red arrows indicate the median of the distribution for gas that has been transformed into stars (this distribution is not shown).

In the text

	Fig. 12 Comparison of the SFR and evolution of the gas mass in the ionised, neutral, and molecular components between the three simulations with different SF laws: Au6_MOL (our model), Au6_BL (an alternative H2-based SF law), and Au6_STD (the standard AREPO implementation, where SF is linked to the total gas density). Only gas cells within a sphere with r=40 ckpc are considered.
In the text

	Fig. 13 Projected mass densities of the stellar, ionised, neutral, and molecular gas components at z=0 for the three simulations with different SF laws: Au6_MOL (our model), Au6_BL (an alternative H2-based SF law), and Au6_STD (the standard AREPO implementation, where SF is linked to the total gas density).
In the text

	Fig. 14 Left: time evolution of the median of εff for the three simulations, together with the median absolute deviation (shades). Middle: correlation between the SFE parameter and the gas density at z=0. Right: correlation between the SFE parameter and the gas metallicity at z=0.
In the text

	Fig. A.1 Evolution of the gas phases when integrating Eqs. (7) up to t=1 Gyr. Different panels corresponds to different values of gas metallicity and density, as indicated above each panel. Solid lines correspond to the ICs fi=fa=0.5, and dashed lines with fi=0.01 and fa=0.99. In both cases the initial values for fm and fs are 0.
In the text

	Fig. A.2 Final value of the different phases (rows) when integrating Eqs. (7) for 1, 10, and 100 Myr (columns). Each colour corresponds to a different value of metallicity.
In the text