# The natural emergence of the correlation between *H*\_2 and star formation rate surface densities in galaxy simulations

Alessandro Lupi, Stefano Bovino, Pedro R. Capelo, Marta Volonteri, Joseph Silk

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alessandro Lupi, Stefano Bovino, Pedro R. Capelo, Marta Volonteri, Joseph Silk. The natural emergence of the correlation between *H*\_2 and star formation rate surface densities in galaxy simulations. Monthly Notices of the Royal Astronomical Society, 2018, 474 (3), pp.2884-2903. 10.1093/mn-ras/stx2874. hal-01703758

HAL Id: hal-01703758

https://hal.science/hal-01703758v1

Submitted on 9 Jan 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The natural emergence of the correlation between H₂ and star formation rate surface densities in galaxy simulations

Alessandro Lupi,<fnref n="1" />★ Stefano Bovino,<fnref n="2" />,<fnref n="3" /> Pedro R. Capelo,<fnref n="4" /> Marta Volonteri<fnref n="1" /> and Joseph Silk<fnref n="1" />,<fnref n="5" />,<fnref n="6" />,<fnref n="7" />

Accepted 2017 November 3. Received 2017 October 3; in original form 2017 July 28

## ABSTRACT

In this study, we present a suite of high-resolution numerical simulations of an isolated galaxy to test a sub-grid framework to consistently follow the formation and dissociation of H₂ with non-equilibrium chemistry. The latter is solved via the package KROME, coupled to the meshless hydrodynamic code GIZMO. We include the effect of star formation (SF), modelled with a physically motivated prescription independent of H₂, supernova feedback and mass-losses from low-mass stars, extragalactic and local stellar radiation, and dust and H₂ shielding, to investigate the emergence of the observed correlation between H₂ and SF rate surface densities. We present two different sub-grid models and compare them with on-the-fly radiative transfer (RT) calculations, to assess the main differences and limits of the different approaches. We also discuss a sub-grid clumping factor model to enhance the H₂ formation, consistent with our SF prescription, which is crucial, at the achieved resolution, to reproduce the correlation with H₂. We find that both sub-grid models perform very well relative to the RT simulation, giving comparable results, with moderate differences, but at much lower computational cost. We also find that, while the Kennicutt–Schmidt relation for the total gas is not strongly affected by the different ingredients included in the simulations, the H₂-based counterpart is much more sensitive, because of the crucial role played by the dissociating radiative flux and the gas shielding.

**Key words:** ISM: molecules – galaxies: evolution – galaxies: formation – galaxies: ISM.

## 1 INTRODUCTION

After the pioneering work by Schmidt (1959), several observations showed that the rapidity with which galaxies form stars correlates with the surface density of the total available gas. However, because of the time-consuming and challenging effort to measure this relation with very high resolution, only recently we have been able to collect large amounts of data to better constrain it. In particular, recent observations by Bigiel et al. (2010) have shown that this relation seems to deviate from the canonical single power law with slope 1.4 (Kennicutt 1998) for low surface densities.

Observations of nearby star-forming galaxies have also evidenced a linear correlation between the star formation rate (SFR) surface density and the molecular hydrogen (H₂) surface density on sub-kpc scales (Bigiel et al. 2008, 2010). On the other hand, the correlation with atomic hydrogen seems to be very weak.

While many numerical simulations still use collisional ionization equilibrium chemistry and do not follow the formation and evolution of molecular hydrogen, several groups have started to use molecular hydrogen-based star formation (SF) prescriptions, either following self-consistently the formation and dissociation of H₂ with non-equilibrium chemistry (Gnedin, Tassis & Kravtsov 2009; Feldmann, Gnedin & Kravtsov 2011; Christensen et al. 2012; Tomassetti et al. 2015) or via sub-grid prescriptions based on the assumption of chemical equilibrium (Pelupessy, Papadopoulos & van der Werf 2006; Hopkins et al. 2014; Thompson et al. 2014;

<footnote n="1">¹Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="2">²Departamento de Astronomía, Facultad Ciencias Físicas y Matemáticas, Universidad de Concepción, Av. Esteban Iturra s/n Barrio Universitario, Casilla 160, Concepción, Chile</footnote>

<footnote n="3">³Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="4">⁴Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="5">⁵AIM-Paris-Saclay, CEA/DSM/IRFU, CNRS, Univ Paris 7, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="6">⁶Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="7">⁷BIPAC, University of Oxford, 1 Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>*E-mail: lupi@iap.fr</footnote>
Somerville, Popping & Trager 2015; Hopkins et al. 2017; Orr et al. 2017; Pallottini et al. 2017).

However, recent theoretical studies have revealed a lack of causal connection between H₂ and SF (Krumholz, Leroy & McKee 2011; Clark et al. 2012), and that the low temperatures of star-forming regions can be reached via molecular cooling as well as metal line cooling (see Capelo et al. 2017, for a detailed study of equilibrium and non-equilibrium metal cooling). They have also pointed out that dust-shielding from ultraviolet (UV) radiation is the main process driving the formation of H₂ and suggested that, more likely, the latter is controlled by SF, or, in general, by the gravitational collapse of atomic gas, not vice versa (Mac Low 2016).

Moreover, theoretical models of turbulent, magnetized clouds have tackled the problem of how SF should proceed inside giant molecular clouds (GMCs) (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011) without any assumption about the correlation with H₂, and found very good agreement with simulations and observations (Federrath & Klessen 2012).

Because of these uncertainties, a clear consensus is still missing, and the observational limits for measuring H₂, currently based on the abundance of CO, in particular for low-metallicity dwarf galaxies where CO appears to be underabundant (Schruba et al. 2012), prevent us from fully disentangling causes and consequences.

In view of this, a key role seems to be played by gas shielding from UV radiation. However, most of the simulations up to date do not consistently take into account local radiation from stellar sources and its effect on the interstellar medium, but only consider a uniform extragalactic background (e.g. Fiacconi, Feldmann & Mayer 2015; Kim et al. 2016). However, a correct inclusion of this effect requires very expensive full radiative transfer (RT) calculations. Several studies have started to incorporate it into numerical simulations, some with on-the-fly RT (Gnedin et al. 2009), and most of the others with sub-grid prescriptions (Christensen et al. 2012; Tomassetti et al. 2015; Forbes et al. 2016; Hu et al. 2016, 2017; Richings & Schaye 2016).

In this study, we present and test a model to self-consistently account for the formation and dissociation of H₂, via the chemistry package KROME (Grassi et al. 2014), including SF, supernova (SN) feedback, and extragalactic and local stellar radiation. We consider the attenuation of the radiation via dust and H₂ shielding (in KROME), and shielding from gas via two different sub-grid models and on-the-fly RT calculations, to investigate the emergence of the correlation between H₂ and SF.

The paper is organized as follows: in Section 2, we describe the model components and the simulation setup; in Section 3, we describe the basic tests of the chemical network with and without photochemistry; in Section 4, we present the results of our simulations including all the aforementioned processes; in Section 6, we discuss the limits of the model and we draw our conclusions.

## 2 SIMULATION SETUP

Here, we present and describe all the components of our model, which we implement in the code GIZMO (Hopkins 2015). Except for the mesh-less finite mass hydrodynamic scheme and the SN feedback prescription, which are the same as in Lupi, Volonteri & Silk (2017), all the other sub-grid prescriptions have been significantly revised with more physically motivated models. We also describe the initial conditions of our simulation suite, which are those of an isolated galaxy with structural parameters typical of z = 3 galaxies.

### 2.1 The hydrodynamic code GIZMO

GIZMO (Hopkins 2015), descendant from GADGET3, itself descendant from GADGET2 (Springel 2005), implements a new method to solve hydrodynamic equations, aimed at capturing the advantages of both usual techniques used so far, i.e. the Lagrangian nature of smoothed particle hydrodynamics codes, and the excellent shock-capturing properties of Eulerian mesh-based codes, and at avoiding their intrinsic limits. The code uses a volume partition scheme to sample the volume, which is discretized amongst a set of tracer 'particles' which correspond to unstructured cells. Unlike moving-mesh codes (e.g. AREPO; Springel 2010), the effective volume associated with each cell is not defined via a Voronoi tessellation, but is computed in a kernel-weighted fashion. Hydrodynamic equations are then solved across the 'effective' faces amongst the cells using a Godunov-like method, like in standard mesh-based codes. In these simulations, we employ the mesh-less finite-mass method, where particle masses are preserved, and the standard cubic-spline kernel, setting the desired number of neighbours to 32. Gravity is based on a Barnes-Hut tree, as in GADGET3 and GADGET2. While the gravitational softening is kept fixed at 80 and 4 pc for dark matter (DM) and stars, respectively, for gas we use the fully adaptive gravitational softening, which ensures that both hydrodynamics and gravity are treated assuming the same mass distribution within the kernel, avoiding the main issues arising when two different resolutions are used. The minimum allowed softening, which sets the effective resolution of the simulation, is set to 1 pc, corresponding to an inter-particle spacing of ~1.4 pc, but this value is reached only in the very high-density regions before a burst of SF occurs.

### 2.2 Star formation

SF is implemented using a stochastic prescription, where the SFR is defined as

$$ \dot{\rho}_{\text{SF}} = \varepsilon \frac{\rho_{\text{g}}}{t_{\text{ff}}} \qquad (1) $$

with $\varepsilon$ the SF efficiency parameter, $\rho_g$ the local gas density, $t_{\text{ff}} = \sqrt{3\pi/(32G\rho_g)}$ the free-fall time and G the gravitational constant. The commonly used prescription for SF assumes a constant value for $\varepsilon$, usually a few per cent, which is calibrated to reproduce the local Schmidt-Kennicutt law (hereafter KS; Schmidt 1959; Kennicutt 1998).

However, such a simplistic assumption is unable to capture the clustered SF in GMCs, resulting in a more diffuse SF across the entire galaxy.

Simulations of star-forming regions (Federrath, Klessen & Schmidt 2009; Federrath et al. 2010) have demonstrated that the usual assumption of a constant SF efficiency used in numerical simulations does not correspond to what really occurs in GMCs, where turbulence can locally trigger very efficient SF while keeping the global GMC stable to collapse. Federrath & Klessen (2012) compared different theoretical models of SF in GMCs with both numerical simulations of turbulent clouds and observations and showed that the multi-free-fall version of the model by Padoan & Nordlund (2011; hereafter PN) exhibits the best agreement with both observations and simulations. In this study, we implement the PN thermo-turbulent model in GIZMO, assuming that our gas particles represent an entire GMC. The SF efficiency is defined by three parameters describing the degree of turbulence in a GMC:
$\alpha_{vir} = 2E_{kin}/|E_{grav}|$,<fnref n="1" /> where $E_{kin}$ is the turbulent kinetic energy and $E_{grav}$ the gravitational energy of the GMC; the Mach number $M$; and the turbulent parameter $b$ describing the ratio between compressive and solenoidal turbulence. For a statistical mixture of solenoidal and compressive turbulence, the typical value is $b = 0.4$ (Federrath et al. 2010). The PN model assumes that the density distribution inside a GMC follows a log-normal probability density function (PDF) with mean density $\rho_0$ and variance $\sigma_{ws}^2 = \ln(1 + b^2 M^2)$ and that SF only occurs above a density $\rho_{crit}/\rho_0 = 0.067\theta^{-2}\alpha_{vir}M^2$, where $\theta = 0.97 \pm 0.10$ corresponds to the thickness of the post-shock layer in the turbulent GMC in units of the cloud size. With these parameters, we can compute the mean SF efficiency of the cloud by integrating the actual SF efficiency over the entire PDF, obtaining

$$ \varepsilon = \frac{\varepsilon_*}{2\phi_t} \exp\left(\frac{3}{8}\sigma_s^2\right) \left[1 + \text{erf}\left(\frac{\sigma_s^2 - s_{\text{crit}}}{\sqrt{2\sigma_s^2}}\right)\right], \quad (2) $$

where $\varepsilon_* = 0.5$ is the local SF efficiency to match observations (Heiderman et al. 2010), $1/\phi_t = 0.49 \pm 0.06$ is a fudge factor to take into account the uncertainty in the free-fall time-scale and $s_{crit} = \ln(\rho_{crit}/\rho_0)$.

Compared to the standard prescription used, which needs to be calibrated, this model is more tied to the simulation, since all the information required to compute $\varepsilon$ can be obtained from the gas properties themselves, using the kernel-weighting scheme available in the code. However, it cannot be used at arbitrarily low resolution, but only when single GMCs can be resolved with at least a few resolution elements.

In this study, we allow SF when gas matches two criteria: (i) $\rho_g > \rho_{SF, thr}$, where we assume that $\rho_g = \rho_0$ in the SF model and (ii) $M > 2$. The first criterion is not strictly necessary, but it is useful to prevent SF in very low density regions where the kernel weighting would extend to very large scales, totally uncorrelated with the typical sizes of GMCs. The second criterion, instead, is fundamental to limit SF in regions dominated by supersonic turbulence and avoid SF in trans-sonic and sub-sonic ones, where the assumption of a log-normal PDF is not valid anymore. Note that, unlike other recipes, we do not have an explicit criterion on the gas temperature. Although the criterion on the Mach number, for a fixed velocity dispersion, gives a temperature criterion, we will show in Appendix A that a natural temperature barrier arises for SF, without the Mach number threshold playing any crucial role.

When the two criteria are matched, we stochastically spawn a new stellar particle with the same mass of the progenitor gas cell mass.

## 2.3 Chemistry and radiative cooling

Chemistry and radiative cooling are computed with the KROME chemistry package (Grassi et al. 2014), solving the non-equilibrium rate equations for nine different species (H, H+, He, He+, He++, H−,

$H_2$, $H_2^+$ and $e^−$). We use model 1a in Bovino et al. (2016),<fnref n="2" /> which includes photoheating, $H_2$ UV pumping, Compton cooling, photoelectric heating, atomic cooling, $H_2$ cooling, and chemical heating and cooling. Compared to the original network, we also include three-body reactions involving $H_2$ (Glover & Abel 2008; Forrey 2013), $H_2^+$ collisional dissociation by H and $H^−$ collisional detachment by He (Glover & Savin 2009). We note that these reactions are included for completeness. In this study, we consider relatively high metallicity gas in the initial conditions, setting $Z = 0.5Z_⊙$. In these conditions, $H_2$ is expected to mainly form on dust, modelled here assuming a fixed dust-to-gas ratio $D$ linearly scaling with $Z$, i.e. $D/D_⊙ = Z/Z_⊙$, where $D_⊙ = 0.00934$. The $H_2$ formation rate on dust is computed following Jura (1975) as

$$ R_f(H_2) = 3.5 \times 10^{-17} n_H n_{H,\text{tot}} Z/Z_\odot C_\rho \text{ cm}^{-3} \text{s}^{-1}, \quad (3) $$

where $n_{H,\text{tot}} = n_H + n_{H^+} + 2n_{H_2} + 2n_{H_2^+} + n_{H^-}$ is the total hydrogen nuclei number density, and $C_\rho = (\rho_g^2)/(\rho_g)^2$ is the clumping factor, accounting for the enhanced $H_2$ formation in the high-density regions unresolved in the simulation. According to Gnedin et al. (2009), the value of $C_\rho$ should be in the range 3–10, although their simulations show that even a value of 100 could not be enough to recover the relation between SFR and molecular gas surface densities. Across the paper, we use two different clumping factor models, a constant $C_\rho = 1$, and the case in which the unresolved gas above $\rho_{SF, thr}$ follows the same log-normal PDF used for SF, giving

$$ C_{\rho} = \exp(\sigma_{s}^{2}) = 1 + b^{2}M^{2}. \quad (4) $$

We also include $H^−$ photodetachment, though its effect should be irrelevant at high metallicity. This is not valid for very low metallicity and will be investigated in a future work.

The dissociation of $H_2$ occurs via two main mechanisms, the excitation to the vibrational continuum of an excited electronic state and the Solomon process (see Bovino et al. 2016), when dissociating and ionizing radiation are present.

To model the ionizing radiation, we consider two different contributions: a uniform extragalactic background, following the model by Haardt & Madau (2012), and the local radiation from young massive stars. We consider eight energy bins, ranging from 0.75 to 1000 eV, to cover the characteristic energies for $H^−$ photodetachment, $H_2$ and $H_2^+$ dissociation, and the different hydrogen and helium ionization states, following a similar approach to Katz et al. (2017). We report the bins in Table 1, with a description of the main processes occurring in every bin.

However, in reality, the extragalactic background is unable to reach the innermost parts of a galaxy, because of the shielding by the intervening gas. In order to account for this effect, we consider that the UV background is attenuated by a factor $\exp(-\tau)$, where $\tau = \sum_i \sigma_{i,\text{bin}} N_i$ is the optical depth for every energy bin, $\sigma_{i,\text{bin}}$ is the cross-section in the bin for the $i$-th species (pre-computed by KROME), and $N_i$ is the $i$-th species column density. We approximate $N_i$ assuming each gas cell as surrounded by gas with similar properties, which gives $N_i \sim n_i \lambda$, where $n_i$ is the $i$-th species number density and $\lambda$ is the absorption length, which we define equal to the Jeans length $\lambda_J$ of the gas cell. This choice is motivated by the results in Safranek-Shrader et al. (2017), where the authors explicitly show

<footnote n="1">¹ We use an approximate value for $\alpha_{vir}$, which can be obtained under the assumption of a spherical homogeneous cloud with diameter $L$, i.e. $\alpha_{vir} = 5\sigma_s^2 L/(6GM_{cloud})$, where $\sigma_v = L||\nabla \otimes v||$ is the velocity dispersion in the cloud and $M_{cloud}$ is the cloud mass. In the simulations, $L$ is computed as the grid-equivalent cell size, i.e. $L = [4\pi/(3N_{ngb})]^{1/3}h \approx 0.5h$ with $h$ the particle smoothing length and $N_{ngb} = 32$ the number of neighbours used in the kernel estimate, and $M_{cloud}$ is the particle mass.</footnote>

<footnote n="2">² We note that there are a few typos in the reaction rates of Bovino et al. (2016), and we refer to Capelo et al. (2017) for details. We also note that we removed the temperature limit at 100 K in reaction 11 of Bovino et al. (2016).</footnote>
<table><caption>**Table 1.** Energy bins for photoheating processes included in the simulations. We use eight bins in the range 0.75–1000 eV to track the different ionization states of the followed species. Columns 1–2 show the minimum and maximum bin energies, respectively, and column 3 shows the main process occurring in each bin. The second bin is actually split in three sub-bins to improve accuracy.</caption><thead><tr><th>E<sub>min</sub> (eV)</th><th>E<sub>max</sub> (eV)</th><th>Process</th></tr></thead><tbody><tr><td>0.75</td><td>1.70</td><td>H<sup>-</sup> photodetachment</td></tr><tr><td>1.70</td><td>11.20</td><td>H<sub>2</sub><sup>+</sup> dissociation</td></tr><tr><td>11.20</td><td>13.60</td><td>H<sub>2</sub> dissociation (Solomon)</td></tr><tr><td>13.60</td><td>15.20</td><td>H ionization</td></tr><tr><td>15.20</td><td>24.59</td><td>H + H<sub>2</sub> ionization</td></tr><tr><td>24.59</td><td>54.42</td><td>H + He + H<sub>2</sub> ionization</td></tr><tr><td>54.42</td><td>100.6</td><td>H + He + He<sup>+</sup> + H<sub>2</sub> ionization</td></tr><tr><td>100.6</td><td>1000.0</td><td>H + He + He<sup>+</sup> + H<sub>2</sub> ionization</td></tr></tbody></table>

that the Jeans length is a better approximation compared to the gradient-based shielding lengths.<fnref n="3" /> In addition, we also consider local radiation from young massive stars, assuming that every stellar particle samples a full stellar population following a Chabrier (2003) initial mass function (IMF) and using the updated stellar population synthesis models by Bruzual & Charlot (2003, hereafter BC03; see Section 2.4 for the model details). At high densities, gas self-shields from radiation, thus we need to include self-shielding by H₂ and dust to properly track the abundances in the high-density regions without overestimating the effective radiative flux. Both terms are implemented as a sub-grid recipe following Richings, Schaye & Oppenheimer (2014).

Metal cooling rates, instead, are provided by look-up tables pre-computed with the photoionization code CLOUDY (Ferland et al. 2013) by Shen, Wadsley & Stinson (2010), and Shen et al. (2013), as a function of density, temperature and redshift, including the extragalactic UV background by Haardt & Madau (2012). Metal cooling rates are provided assuming solar abundances and then linearly rescaled with the total metallicity, which is followed as a passive scalar.

Cooling is not allowed below 10 K, to avoid inconsistencies in the cooling rates, which are not well defined below this temperature.

## 2.4 Supernova and wind feedback

Since we cannot resolve single stars in our simulations, we assume that our stellar particles correspond to an entire stellar population, following a Chabrier (2003) IMF. According to stellar evolution, we consider three different processes: type II supernovae (SNe), type Ia SNe and low-velocity winds by stars in the range 1–8 M☉.

(i) After ~4 Myr, the most massive stars start to explode as type II SNe. For stars with mass between 8 and 40 M☉, we release $E_{SN} = 10^{51}$ erg per SN via thermal injection on to the gas cells enclosed in the stellar particle kernel sphere,<fnref n="4" /> together with mass and metals. Due to the limited resolution, the energy-conserving phase of the bubble expansion cannot be followed in the simulation and the additional energy released would be rapidly lost because of radiative cooling, resulting in an ineffective SN feedback. In order to avoid

gas from rapidly getting rid of this additional energy, we implement a delayed-cooling prescription, where gas cooling is inhibited for the time needed to reach the momentum-conserving phase (Stinson et al. 2006). In particular, we shut off cooling for the survival time of the blast-wave $t_{max} = 10^{6.85} E_{51}^{0.32} n_{Htot}^{0.34} \bar{P}_{04}^{-0.70}$ yr, where $E_{51}$ is the released energy in units of $10^{51}$ erg and $\bar{P}_{04} = 10^{-4} P_0 k_B^{-1}$, with $P_0$ the ambient pressure and $k_B$ the Boltzmann constant. This model has been proven to better reproduce the KS relation and the outflow mass-loading factor in isolated galaxy simulations compared to other more physically motivated models, as stated by Rosdahl et al. (2017). The only limitation is that it produces ‘unphysical’ high temperatures in a region of the density–temperature diagram where cooling is expected to be effective. We assume that type II SNe return all their mass but 1.4 M☉ and we follow metal production (via iron and oxygen yields whose production rates are thought to be metallicity independent) using the tabulated results by Woosley & Heger (2007) fitted by Kim et al. (2014). The total metal mass can then be defined as

$$M_Z = 2.09 M_{\text{oxygen}} + 1.06 M_{\text{iron}}, \quad (5)$$

where $M_{\text{oxygen}}$ and $M_{\text{iron}}$ are oxygen and iron mass, respectively, computed convolving the yield function with the IMF. For stars above 40 M☉, we assume a black hole would form via direct collapse (although we do not simulate sink formation), thus we release neither mass nor metals.

(ii) Type Ia SNe occur in evolved binary systems, when one of the stars has become a white dwarf and has accreted enough mass to exceed the Chandrasekhar mass limit (Chandrasekhar 1931). Type Ia SNe explode according to a distribution of delay times, as described by Maoz, Mannucci & Brandt (2012), scaling as $t^{-1}$ between 0.1 and 1 Gyr after a burst of SF. Type Ia SNe leave no remnant and release to the environment 1.4 M☉, $M_{\text{iron}} = 0.63 M_{\text{☉}}$ and $M_{\text{oxygen}} = 0.14 M_{\text{☉}}$. Since type Ia SNe occur a long time after the burst of SF, these events are usually located far away from the progenitor molecular cloud, and they are not clustered as type II SNe. In this case, we do not shut off cooling, as described in Stinson et al. (2006).

(iii) Because of the low mass, stars between 1 and 8 M☉ evolve on longer time-scales compared to their massive counterparts and do not explode as SNe. During their evolution, however, they release part of their mass as stellar winds. We model stellar winds by injecting only mass and metals on to the neighbouring cells of the stellar particle. Assuming the initial–final mass relation for white dwarfs by Kalirai et al. (2008), the mass-loss for these stars can be computed as $w_m = 0.394 + 0.109m M_{\text{☉}}$, where m is the mass of the star. Also in this case, we convolve $w_m$ with the assumed IMF to compute the net mass-loss as a function of stellar age. We assume that stellar winds only carry the progenitor star metallicity, hence they do not enrich the interstellar medium.

By taking into account these processes, the stellar feedback in our simulations is able to return 42 per cent of the stellar particle mass in a Hubble time, prolonging SF even in cases of no fresh gas inflows (e.g. Ciotti et al. 1991; Leitner & Kravtsov 2011; Voit & Donahue 2011).

## 2.5 Stellar radiative feedback

In addition to the previously described processes, we also include radiative feedback from the stellar population. This form of feedback is particularly important for young stars, because of the strong

<footnote n="3">³ Though better results are obtained when a temperature cap of 40 K to the Jeans length is used, the largest differences are observed in the CO distribution, which is not followed in this study.</footnote>

<footnote n="4">⁴ For the stellar particles, we define the kernel size using 64 neighbours, in order to guarantee that the volume around the star is fully covered.</footnote>
ionizing and dissociating flux produced, which affects the abundance of molecular gas. We assume that every stellar particle spawned in the simulation corresponds to an instantaneous burst of SF and we model the evolution of the stellar luminosity $L_*$ with age and metallicity using the BC03 models. We consider here two different sub-grid models accounting for this contribution, to compare them with on-the-fly RT calculations and identify the most accurate to be used in cosmological simulations, where RT is still prohibitive.

In model (a), we loop over all the active gas cells, collecting all the stars within a maximum distance $R_{max}$ from the cell, defined as the critical radius at which the surrounding medium would become optically thick, i.e. the optical depth $\tau \sim 1$. In order to compute this radius, we consider the local cell properties and integrate the gas density $\rho_g(R)$ outwards along the direction R of maximally decreasing density, obtaining

$$ \tau = \frac{\sigma_{\text{eff}}}{m_{\text{H}}} (\rho_{\text{g}} R_{\text{max}} - |\nabla \rho_{\text{g}}| R_{\text{max}}^2 / 2), \quad (6) $$

where $\sigma_{\text{eff}}$ is the minimum cross-section for the photoionization processes considered, mH is the hydrogen mass, and $\nabla \rho_g$ is the density gradient. Although $\sigma_{\text{eff}}$ is different for every process considered, we approximate it using a unique value $\sigma_{\text{eff}} = 10^{-20} \text{ cm}^2$, small enough to avoid an unphysical underestimation of $R_{max}$, and at the same time large enough to avoid looking for all the stars in the galaxy. Among the two solutions of equation (6), we take the smallest $R_{max}$. When a solution to the equation cannot be found, we assume $R_{max} = \max(5h, \rho_g/|\nabla\rho_g|)$, where h is the kernel size of the gas cell and $\nabla\rho_g$ is the gas density gradient around it.<fnref n="5" /> After collecting the stellar particles within $R_{max}$, we compute the total radiative flux in each bin which reaches the target gas cell as

$$ F = \sum_{i} \frac{L_{*,i}}{4\pi d_i^2} \exp(-\tau_i), \quad (7) $$

where $L_{*,i}$ and $d_i$ are the luminosity and the distance of the i-th stellar particle, respectively, and $\tau_i$ is the optical depth of the intervening gas between the cell and the source. Since we cannot estimate the real optical depth between the cell and the source, we approximate $\tau_i$ using the average species abundances around the stellar particle and an absorption length defined as the minimum between $d_i$ and $l_{\text{Sob}} = \rho_*/|\nabla\rho_*|$, where $\rho_*$ and $\nabla\rho_*$ are the gas density and the gas density gradient around the star, respectively. This model can be considered an improved version of that in Christensen et al. (2012), where shielding was not considered.

In model (b), instead, we benefit from the tree structure used for gravity computation, storing the total luminosity and the luminosity-weighted centre of the stellar sources in every node, as in Hopkins et al. (2014, 2017). The stellar luminosity is attenuated at the source using the Sobolev approximation, giving $L_{\text{eff},*} = L_* \exp(-\tau)$, where $\tau$ is defined as in model (a). Then, during the gravity force calculation, we compute the total flux f reaching a gas cell collecting the contribution of all the stellar particles in the simulated volume. Finally, we attenuate the resulting flux at the absorption point using, as for the extragalactic UV background, the Jeans length of the cell.<fnref n="6" />

## 2.6 Radiative transfer with GIZMO

For the RT simulations, we employ the momentum-based RT method with the M1 closure scheme (Levermore 1984). The implementation directly follows that of Rosdahl et al. (2013), suitably modified for a mesh-less structure.<fnref n="7" /> In particular, the Riemann problem is solved using the Harten, Lax and van Leer (HLL) solver (Harten, Lax & van Leer 1983), which better models beams and shadows compared to the simpler and more diffusive Global Lax Friedrichs (Lax 1954) function (obtained by setting the HLL eigenvalues to $\pm c$, where c is the speed of light in vacuum), but increases the asymmetry of the ionization front around isotropic sources. The RT calculations are performed using an operator-splitting approach in three steps: injection, transport and thermochemistry. During the injection step, the radiative energy emitted by the stellar source during the current time-step is distributed among the cells within the kernel volume of the stellar particle, as done for SN feedback. During the transport step, photon-advection is performed by solving the Riemann problem among the interacting cells, as done for the hydrodynamics, and finally, for the thermochemistry step, we couple the code with KROME, computing the radiative flux in a cell (for each energy bin) as

$$ F = \frac{c N_{\nu}}{h^3}, \quad (8) $$

where $N_\nu$ is the number of photons per bin in the cell. The energy bins (see Table 1) and stellar spectra used in the RT simulation are the same used for the sub-grid models (a) and (b).

The time-step for the RT calculations is constrained by the Courant condition

$$ \Delta t_{\text{RT}} \le C_{\text{CFL}} \frac{h}{c}, \quad (9) $$

where $C_{\text{CFL}} = 0.2$ is the Courant factor and h the particle/cell size. Compared to the common hydrodynamic case (where $v_{\text{max}} \sim 1000 \text{ km s}^{-1}$), the time-step condition for RT can be hundred-thousand times more severe, making the simulations much more computationally expensive. The simplest solution to this issue is the reduced speed of light approximation, introduced by Gnedin & Abel (2001), where the Courant condition is relaxed by replacing c with $\tilde{c} = f_c c$. The exact value of $f_c$ depends on the boundary conditions of the problem and the associated typical time-scales and can be expressed as

$$ f_c = \min(1; \sim 10t_{\text{cross}}/\tau_{\text{sim}}), \quad (10) $$

where $t_{\text{cross}} = r_S/c$, with $r_S = [3\dot{N}/(4\pi\alpha^B n_H^2)]^{1/3}$ the Strömgren (1939) sphere radius, $\dot{N}$ the ionizing photons injection rate and $\alpha^B \sim 2.6 \times 10^{-13} \text{ cm}^3 \text{ s}^{-1}$ the case-B recombination rate at $T \sim 10^4 \text{ K}$. $\tau_{\text{sim}}$ is the shortest relevant time-scale of the simulation. In our simulations, we assume that the typical radiation source is a young massive star ($\dot{N} \sim 1-5 \times 10^{48} \text{ s}^{-1}$), surrounded by high-density gas ($n_H \gtrsim 10 \text{ cm}^{-3}$), and we require $\tau_{\text{sim}} \sim 1 \text{ Myr}$ to resolve the stellar evolution of these massive stars. Putting everything together, we obtain a minimum value of $f_c \sim 3 \times 10^{-4}$. However, in our simulations, we assume $f_c = 10^{-3}$ as a good compromise between higher accuracy and larger computational cost.

densities to be slightly higher than in reality. Nevertheless, since a simple correction factor is not available, it represents a simple and reasonable approximation that avoids time-consuming calculations.

<footnote n="5">⁵ We choose the factor 5 because it represented the best compromise between limiting the computational overhead and collecting radiation from large enough distances.</footnote>

<footnote n="6">⁶ As Hartwig et al. (2015) pointed out, the Jeans length tends to overestimate the $H_2$ column densities ($N_{H_2}$) by up to an order of magnitude for $N_{H_2} \lesssim 10^{16} \text{ cm}^{-2}$ and to underestimate it for larger $N_{H_2}$. Hence, we expect our $H_2$</footnote>

<footnote n="7">⁷ A detailed description of how to implement radiative transfer methods in mesh-less codes can be found at https://physastro.pomona.edu/wp-content/uploads/2016/07/David-Khatami-Final-Thesis.pdf.</footnote>
<table><caption>**Table 2.** Description of the full simulation suite. We show in the first column the name of the run and in the second the presence or not of the extragalactic UV background. The third column reports the model used (if any) for stellar radiative feedback (see Section 2.5) and the last one the clumping factor $C_\rho$ used in the run.</caption><thead><tr><th>Run</th><th>UV background</th><th>Rad. feedback</th><th>C<sub>ρ</sub></th></tr></thead><tbody><tr><td>G_mol</td><td>No</td><td>/</td><td>1</td></tr><tr><td>G_UV</td><td>Yes</td><td>/</td><td>1</td></tr><tr><td>G_RT_1</td><td>Yes</td><td>RT</td><td>1</td></tr><tr><td>Gojsc</td><td>Yes</td><td>Model (a)</td><td>Variable</td></tr><tr><td>G_RAD_T</td><td>Yes</td><td>Model (b)</td><td>Variable</td></tr><tr><td>G_RT</td><td>Yes</td><td>RT</td><td>Variable</td></tr></tbody></table>

For an exhaustive description of the RT implementation in GIZMO and numerical tests, we refer to Hopkins et al. (in preparation).

## 2.7 Initial conditions

In this study, we follow the evolution of an isolated galaxy at $z=3$. The galaxy is a composite system of a DM halo, a gaseous and stellar disc, and a stellar bulge. The DM halo is described by a spherical Navarro-Frenk-White (NFW; Navarro, Frenk & White 1996) density profile within the virial radius, $R_{vir}$, and by an exponentially decaying NFW profile outside $R_{vir}$, as described in Springel & White (1999). The stellar and gaseous disc is described by an exponential surface density profile and by an isothermal sheet (Spitzer 1942), whereas the stellar bulge is described by a spherical Hernquist (1990) profile. The galaxy has a virial mass $M_{vir} = 2 \times 10^{11} M_{\odot}$ (Adelberger et al. 2005) and a virial radius $R_{vir} = 45$ kpc, with 2 and 0.4 per cent of its virial mass in the disc and bulge, respectively. The mass fraction of gas in the disc is 60 per cent, typical of high-redshift galaxies (Tacconi et al. 2010). The modelled system is similar in structure to the main galaxy in the gas-rich galaxy merger simulated in Capelo et al. (2015). This same galaxy model has also been used for a study on the effects of non-equilibrium metal (C, O and Si) cooling (Capelo et al. 2017). The mass resolution is 1594 $M_{\odot}$ for gas and stellar particles and $1.97 \times 10^5 M_{\odot}$ for DM particles, corresponding to a total number of particles in the initial conditions of 1.5 M in gas, 1.5 M in DM, 1 M in the stellar disc and 500k in the stellar bulge.

At the beginning of the simulations, the instantaneous cooling of the gas in the central regions of the galaxy would result in the vertical collapse of the disc and into an initial burst of SF. To prevent this spurious effect, we shut off SF for the first 10 Myr and let pre-existing stars explode as SNe, assuming their stellar ages are uniformly distributed in the range [0–2] Gyr. The SF density threshold is set to a total hydrogen density $n_{Htot} = 1 \text{ cm}^{-3}$ to avoid SF in very low density regions where the kernel size would extend to very large scales uncorrelated with typical GMC sizes, as stated in Section 2.2. The galaxy is evolved for 400 Myr, a time long enough to follow the evolution after the initial relaxation, but still sufficiently short to allow us to neglect environmental effects like accretion flows and mergers. We run a suite of six simulations where we include different sub-grid prescriptions which affect the evolution of molecular gas, to assess the relative contribution of the different processes and to assess the best sub-grid model for stellar radiative feedback which better matches RT calculations. The full simulation suite is reported in Table 2.

# 3 TESTS OF THE CHEMISTRY NETWORK

We discuss here the results of the runs G_mol and G_UV, aimed at testing the coupling between GIZMO and the chemistry package KROME.

## 3.1 Test 1: Non-equilibrium chemistry in GIZMO and molecular gas evolution

We describe here the results of G_mol, where we follow the evolution of molecular gas when only collisional dissociation/ionization and recombination processes are considered. In this case, we expect the cold gas to always be fully molecular, even at low densities. On the other hand, the only available sources for gas heating are SNe and adiabatic compression. All the analyses at a fixed time are carried out at $t = 350$ Myr and not at the final time of the simulation, to show the galaxy properties at the median snapshot used for the SF law, obtained by stacking over a period of 100 Myr ($300 \text{ Myr} \le t \le 400 \text{ Myr}$) around $t = 350$ Myr.

In Fig. 1, we report the gas properties averaged over 100 Myr around $t = 350$ Myr, using logarithmic bins 0.04 decade wide on both axes, except for the bottom-right plot, where we use linear bins of width 0.05 on the y-axis. In the top-left panel, showing the density-temperature diagram of the gas, we observe that the gas clusters in two distinct regions, one around $T \sim 10^4$ K, extending from $n_{Htot} \sim 0.01 \text{ cm}^{-3}$ up to $n_{Htot} \sim 10 \text{ cm}^{-3}$, and the other in the range $100 \text{ K} \lesssim T \lesssim 5000 \text{ K}$, for $n_{Htot} \gtrsim 1 \text{ cm}^{-3}$. The first one corresponds to the low to intermediate density gas which was heated up by SNe and has cooled down via Ly $\alpha$ and metal cooling to the critical temperature where molecular gas cooling and metal cooling start to dominate the cooling rate. The second, instead, is the region where gas efficiently cools down to the very low temperatures where SF occurs. Finally, we also observe a high-temperature tail ($T \gtrsim 3 \times 10^4$ K) for $0.01 \text{ cm}^{-3} \lesssim n_{Htot} \lesssim 100 \text{ cm}^{-3}$, corresponding to the gas heated up by SNe and for which cooling is still disabled.

The other three panels in Fig. 1, instead, report the hydrogen species abundances, i.e. the neutral hydrogen $f_H$, ionized hydrogen $f_{H^+}$, and molecular hydrogen $f_{H_2}$ mass fractions (normalized to the total hydrogen nuclei abundance) as a function of the total hydrogen density $n_{Htot}$. We clearly see that most of the gas is neutral/ionized, and only a small fraction becomes molecular. However, the apparent necessary condition to become fully molecular is the low temperature, independent of the density, as clearly visible in the bottom-right panel. We also find that a molecular mass fraction of 10 per cent can already be obtained at fairly low densities around few $\text{cm}^{-3}$, although the transition to fully molecular gas is observed at $n_{Htot} \sim 100 \text{ cm}^{-3}$.

In Fig. 2, we show the projected gas density maps of the galaxy in the face-on (top panels) and edge-on (bottom panels) views at $t = 350$ Myr, where the left-hand panels correspond to the total hydrogen gas and the right-hand ones to the $H_2$ component only. The molecular gas density is higher in the higher density spiral arms of the galaxy, but a not negligible fraction is also observed in the low-density medium, in agreement with the results in Fig. 1.

Thanks to the detailed chemistry model used, which allows us to follow the molecular gas abundance, we can also compare the relation between the SFR and the gas density in the galaxy, both in total gas ($H + H_2$) and $H_2$ only. To better compare our simulations with observations, we measure both SFR and gas densities in patches of 750 pc side length, matching the resolution in Bigiel et al. (2010). In addition, instead of using the instantaneous SFR computed from the gas properties, we compute the stellar luminosity of our stellar
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Average gas properties for the G_mol run around *t* = 350 Myr. The top-left panel shows the gas distribution in the density–temperature plane. The other three plots, instead, report the neutral (*f*<sub>H</sub>), ionized (*f*<sub>H</sub><sup>+</sup>) and molecular hydrogen (*f*<sub>H<sub>2</sub></sub>) mass fractions as a function of the total hydrogen density *n*<sub>Htot</sub>. Although most of gas above 100 cm<sup>-3</sup> is fully molecular as expected, high values of *f*<sub>H<sub>2</sub></sub> are also visible at lower densities, because of the missing dissociating radiation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Projected gas density maps for our G_mol run at *t* = 350 Myr, face-on (top panels) and edge-on (bottom panels). The left-hand panels show the total hydrogen density and the right-hand panels the H<sub>2</sub> component only. Molecular gas forms almost everywhere in the galaxy, where the gas is cold enough. The absence of the UV background allows H<sub>2</sub> to form also at low densities, unless the gas is heated up by SNe.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** SF law for our G_mol run in the last 100 Myr of the simulation, plotted as blue and red line contours. We show the relation for the total (H + H₂) gas in the left-hand panel and that for H₂ only in the right-hand panel and compare them with the observational data by Bigiel et al. (2010, green and red) and Schruba et al. (2011, purple). The orange lines in the right-hand panel correspond to 0.1, 1 and 10 Gyr depletion times (defined as Σ\* / ΣH₂). Our data well match the observational data in total gas, whereas in H₂ the slope is steeper, with the contours extending to depletion times of 0.1 Gyr also at moderately low densities of about 1 M⊙ pc⁻².</figcaption></figure>

particles in the GALEX FUV band (1250–1750 Å), using the BC03
models, and convert it into an SFR using the relation by Salim et al.
(2007), to be consistent with the approach used by Bigiel et al.
(2010).<fnref n="8" /> Since the SFR in the galaxy does not vary significantly
during the run (see left-hand panel of Fig. 10), we could stack the
data from 10 different snapshots in the last 100 Myr of the run
(300 Myr ≤ t ≤ 400 Myr), every 10 Myr, to increase the number
of data points. The results are reported in Fig. 3, where we show
the total gas (H + H₂) KS relation in the left-hand panel and the
H₂ SF law in the right-hand panel, and compare our data with the
observations by Bigiel et al. (2010) and Schruba et al. (2011). These
are plotted in filled contours of more than 1, 2, 5 and 10 points per
bins, binning in both gas and SFR density with 0.1 decade-wide
bins. The green contours correspond to the Bigiel et al. (2010) data
inside r<sub>25</sub>, the red ones to the Bigiel et al. (2010) data outside r<sub>25</sub>
and the purple ones to Schruba et al. (2011). Our data points, instead,
are plotted as line contours of more than 2 and 10 points per bin,
using 0.2 decade-wide bins in both gas and SFR density bins.

To highlight the SFR dependence on the radial distance from
the galaxy centre, we split our point distribution in two data sets,
corresponding to the patches inside (blue contours)/outside (red
contours) 4 kpc (corresponding to the radius where ~98 per cent of
SF occurs). The simulated data points match well the observational
data, except for the highest densities (corresponding to the central
2 kpc of the galaxy), where SFR is higher. On the other hand, in the
right-hand panel our points lay above the observed relation, and the
slope is also different. While the different slope is consistent with

a too large molecular gas fraction at low density, which moves the
data points to the right and steepens the relation, at intermediate to
high densities the molecular gas fraction is too low for the measured
SFR, probably because the resolution is too low and underestimates
the real H₂ formation rate in the high-density regions. As we will
show also for the other runs, this is a common behaviour, which can
be solved using a sub-grid clumping factor, as we will describe in
Section 4.3.

## 3.2 Test 2: Non-equilibrium chemistry and photochemistry under a uniform extragalactic background

We now describe the results of the G_UV run, where we also in-
clude photochemistry. In this case, the additional heating source for
gas comes from a uniform extragalactic UV background (Haardt &
Madau 2012) at z = 3. Because of the gas shielding prescription
implemented, we expect only the gas at low density to be ionized,
and almost no difference compared to the G_mol run at interme-
diate and high densities. In Fig. 4, we show the gas distribution in
the density–temperature plane and the H₂ mass fraction. We ob-
serve in the left-hand panel that the gas with n<sub>Htot</sub> ≲ 10<sup>-2</sup> cm<sup>-3</sup> is
heated up above 10<sup>4</sup> K, like in G_mol, except that here the heating
is mainly due to the UV background instead of the SNe. On the
other hand, at n<sub>Htot</sub> ≳ 10<sup>-2</sup> cm<sup>-3</sup>, we observe a larger mass of gas at
low temperature. The reason for this difference is that the external
UV background heats the gas, slowing down gas collapse and frag-
mentation, reducing the SFR and the subsequent SN heating. The
absence of the UV background in G_mol allows gas to fragment
and form stars more easily also in the galaxy outskirts. As a conse-
quence, SNe explosions in a very low density medium become very
effective at sweeping away gas, heating it and reducing the density
even more. In the right-hand panel, we clearly see that the gas be-
comes fully molecular above n<sub>Htot</sub> = 100 cm<sup>-3</sup>, but no molecular
gas is observed at low density, unlike in the G_mol run, because of
the enhanced dissociation rate due to the external UV background.

<footnote n="8">⁸ To test the possible bias in the measure using this approach, we compared the resulting SFR with that obtained by averaging the stellar mass formed over 10, 50 and 100 Myr, respectively, and with the instantaneous SFR estimated from the gas density, finding almost perfect agreement, with small differences only at low SFRs (Σ\* ≲ 5 × 10⁻³ M☉ yr⁻¹ kpc⁻²). Still, this was expected, because of the quiet SFR history of the galaxy.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Gas distribution in the density–temperature plane (left-hand panel) and H₂ mass fraction as a function of $n_{H\text{tot}}$ (right-hand panel) for the G_UV run averaged over 100 Myr around $t = 350$ Myr. Compared to G_mol, cold gas is observed down to $n_{H\text{tot}} \sim 10^{-2} \text{ cm}^{-3}$. In addition, the low-density tail of fully molecular gas observed in G_mol is not present anymore in G_UV (right-hand panel), because of the presence of dissociating radiation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Fig. 3 for the G_UV run. In the left-hand panel, the agreement with observational data is still good, with narrower contours overlapping with the Bigiel et al. (2010) data within $r_{25}$. In the right-hand panel, instead, the slope has slightly improved, but the normalization is offset by a factor of a few.</figcaption></figure>

Fig. 5 shows the SF law for the G_UV run, obtained with the same procedure used for G_mol. In the left-hand panel, we clearly see that the inclusion of the UV background decreases the spread in the relation at low to intermediate gas densities, without changing the slope of the relation. Also in this case, our data points lie on the observed relation, both for the inner and outer part of the galaxy. Again, at high densities we still do not observe a bending, similar to the G_mol run. In the right-hand panel, instead, the trend is significantly different, with a narrower and less steep distribution compared to the observed one, especially in the central 4 kpc of the galaxy.

## 4 THE EFFECT OF STELLAR RADIATIVE FEEDBACK

We describe here the results of our fiducial runs Giemsa, G_RAD_S, G_RAD_T and G_RT, where we include radiative feedback from stars and a variable clumping factor, comparing the same properties discussed in the previous section. Unlike in previous studies, where the clumping factor was constant for all the gas, here we

self-consistently compute $C_\rho$ assuming the same PDF used for SF, but only for gas above $\rho_{\text{SF, thr}}$.

Despite the ever increasing computational power, the resolution achieved in our simulations is still prohibitive for cosmological simulations. Appendix B presents analogues of the simulations described here, but at low resolution, showing that, provided that the gas particle mass does not exceed the typical mass of a GMC, results are consistent with those discussed in this section.

### 4.1 Molecular gas distribution in the galaxy

We compare here the molecular gas distribution in the galaxy at $t = 350$ Myr. We first discuss the qualitative behaviour from the projected gas density maps in Figs 6 and 7. The small differences in the radiative feedback models can significantly affect the SF history of the galaxy, producing visible differences in the galaxy maps, which however do not have a strong impact on the conclusions.

In the face-on view (top panels of Fig. 6), Giemsa exhibits a slightly smoother gas distribution compared to G_RAD_S and G_RT, with a less defined spiral pattern. In the bottom panels
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Projected gas density maps of total hydrogen in the Giemsa S, Giemsa T and G_RRT runs at t = 350 Myr, face-on (top panels) and edge-on (bottom panels). The gas distribution is smoother in Giemsa S compared to Giemsa T and G_RRT, with less defined spiral arms. Though not shown here, this difference persisted in other snapshots. In the edge-on view, instead, the slightly thicker disc in G_RRT is simply due to the feedback resulting from a burst of SF occurred around t = 300 Myr.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Projected gas density maps of the H₂ component for our Giemsa S, Giemsa T and G_RRT runs at t = 350 Myr, face-on (top panels) and edge-on (bottom panels). In Giemsa_S, H₂ can also form at slightly lower densities compared to the other two runs, due to an underestimation of the radiative flux in this regime. G_RRT, instead, is closer to G_RRT, thanks to the gravity tree-based model which accurately accounts for all stellar sources in the galaxy, although there is still more H₂ than in G_RRT.</figcaption></figure>

(edge-on view), the differences are smaller. Giemsa_S shows lower densities at higher distances from the galaxy plane, compared to the other runs, although the difference is negligible in terms of mass. G_RRT, instead, displays a slightly thicker disc amongst the three

runs, compatible with the stronger feedback following a burst of SF which occurred at t ~ 300 Myr (see Fig. 10).

For the molecular gas distribution (Fig. 7), the differences are more evident. Giemsa_S generally forms more H₂ compared to
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Gas properties for the G RAD_S (left-hand column), G RT (middle column) and G RAD_T (right-hand column) runs averaged over 100 Myr around t = 350 Myr. We report molecular hydrogen (top row) and neutral hydrogen (middle row) mass fractions as a function of nHtot. The bottom panels, instead, correspond to the gas distribution in the density–temperature plane. All the runs show similar features, with small differences in the temperature for warm gas around 10^4 K and in the abundance of molecular gas. The strongest difference is in the abundance of atomic hydrogen in G RAD_S, where the underestimation of the stellar radiation at intermediate densities, when τ = 1 from equation (6) is never reached, results in a massive tail extending down to fH ~ 10^3.</figcaption></figure>

the other runs, whereas G_RT and G RAD_T, instead, exhibit a
more peaked distribution, with H₂ formation more vigorous in the
high-density regions of the galaxy spiral arms and suppressed in the
lower density regions.

After a first qualitative comparison, we now quantitatively assess
the differences amongst the models in Fig. 8.

In the top row, we compare the H₂ mass fraction as a function of the total hydrogen density for the three runs. Giemsa's G RAD_S forms significantly more molecular gas mass at intermediate densities (n<sub>Htot</sub> ~ 1–10 cm<sup>-3</sup>), compared to G RT and G RAD_T, which show instead a very similar behaviour. Anyway, the transition to fully molecular gas is the same for all these three runs and corresponds to n<sub>Htot</sub> ~ 10 cm<sup>-3</sup>. In the middle row, the differences are more evident. While G RAD_T and G RT show very similar features except for a slightly larger amount of atomic hydrogen in G RAD_T at n<sub>Htot</sub> = 1–3 cm<sup>-3</sup>, G RAD_S exhibits a massive tail extending down to f<sub>H</sub> ~ 10<sup>-3</sup>, consistent with the more efficient conversion of H in H<sub>2</sub>. Finally, in the bottom row, we observe very small differences, with G RAD_S showing slightly larger amounts of cold gas (below T ~ 100 K) and warm gas (around T = 10<sup>4–5</sup> K). G RT and G RAD_T, instead, are in even better agreement, with only a little more warm gas in the range 1 cm<sup>-3</sup> < n<sub>Htot</sub> < 100 cm<sup>-3</sup>.

We now compare the H₂ column density with local observations. Since observations at Z ~ 0.5 Z⊙ are not available, we report data for the Milky Way disc and halo (at solar metallicity) and the Large Magellanic Cloud (LMC, 0.3 Z⊙), which should bracket our simulation results. We compute the column density of both atomic and molecular hydrogen, respectively N<sub>H</sub> and N<sub>H₂</sub>, in a cylinder of 10 kpc radius and 1 kpc height along 25 randomly distributed line of sights. We report in Fig. 9 the logarithmic H₂ column density fraction as a function of the total logarithmic column density. The results are binned in 200 log-spaced bins along both axes, where the colour coding represents the point density in each bin. The cyan stars correspond to the LMC data by Tumlinson et al. (2002), the blue squares to the Milky Way disc data by Wolfire et al. (2008), and the green diamonds to the Milky Way halo from the FUSE survey. The simulation results agree very well with observations, overlapping both LMC and Milky Way data for a large fraction of the available data range. This is consistent with our expectations, considering that the initial metallicity was Z = 0.5 Z⊙ and SNe enriched the galaxy during the 400 Myr of evolution, moving it closer to the solar metallicity data. Nevertheless, the lower metallicity data for the LMC are still in reasonable agreement with our data, lying at the edge of the distribution. In this case, the different local
<figure><img src="image_12.png" /><figcaption>**Figure 9.** H₂ column density fraction in our Giemsa (left-hand panel), G-RRT (central panel) and G-RRT (right-hand panel) runs at t = 350 Myr, compared with local observations of molecular clouds in the Milky Way disc and halo, and in the LMC. The observational data very well bracket our simulation results, overlapping with them in most of the available column density range. There are small differences among the three runs, especially at low column densities, where Giemsa shows a slightly higher H₂ column density fraction consistent with the enhanced formation of H₂ already observed in Fig. 8. G-RRT, instead, exhibits a denser concentration at high H₂ fractions, consistent with the stronger shielding at high gas densities.</figcaption></figure>

radiation model considered clearly affects the results. In particular, Giemsa exhibits a generally higher H₂ fraction at all column densities, because of the underestimation of the radiation flux at intermediate gas densities, as already discussed. Compared to G_RRT, the distribution in Giemsa also extends to larger column densities. A similar behaviour is also observed in G_RRT. For column densities between $10^{20}$ and $10^{21}$ cm⁻², instead, G_RRT shows a similar trend to G_RRT, but with a slightly sparser distribution. A critical feature of both G_RRT and G_RRT is that they seem to never reach H₂ fractions higher than 80–90 per cent, contrary to G_RRT, where the suppressed radiation boosts the formation of H₂, shifting the data upwards.

In summary, in G_RRT the assumption of a maximum search radius allowed us to accurately account for the radiation close to the stellar sources, but at the cost of neglecting part of the stellar flux for gas at intermediate or low densities, when $\tau = 1$ could not be reached. This inevitably resulted in a slight overestimation of the molecular gas fraction at lower densities compared to the other models. In G_RRT, instead, the attenuation at both the source and the absorption point enabled us to consider all the stellar sources in the galaxy, but with the drawback of overestimating the attenuation when the star–cell distance was shorter than the sum of the two absorption lengths. This weakened the effect of radiation in the very high-density regions, whereas at intermediate or low density the approximation was much more accurate. In spite of these limitations, the differences among G_RRT, G_RRT and G_RRT remained small, and only affected the results at a few per cent level in the total molecular gas and stellar mass formed. A crucial role in suppressing the differences is played by the clumping factor, which globally shifts the transition to the fully molecular regime at lower densities, where the differences in the radiative model were stronger. In conclusion, both sub-grid models considered here represent valid alternatives to expensive RT calculations: with the employed reduced speed of light of 0.001 c, the overhead of the G_RRT simulation was a factor of 2–3 compared to G_RRT and roughly a factor of 5 compared to G_RRT.

## 4.2 SF in the galaxy and the SF law

In Fig. 10, we show the SFR (left-hand panel) and the formed stellar mass (right-hand panel) for the G_RRT, G_RRT and G_RRT.

runs in the last 300 Myr of the simulations. We explicitly excluded the first 100 Myr (~ three dynamical times of the galaxy) from the comparison in order to avoid introducing a bias due to the initial relaxation phase, which was slightly different in the different runs. All the runs show a similar SFR, although the SF bursts are offset in time in the different cases. This suggests that, despite the good agreement in the chemical properties amongst the three simulations, the different models for stellar radiation significantly affect the details of the evolution, changing the actual SF history of the galaxy on time-scales of about 20–50 Myr. In particular, the SFR in G_RRT exhibits a second burst around t = 300 Myr not observed in the other runs, and G_RRT exhibits the lowest SFR amongst the three runs, which is likely related to a slightly stronger ionization close to the stellar sources compared to the other models. This difference in G_RRT is also reflected in the total stellar mass formed (right-hand panel), with G_RRT producing the lowest stellar mass compared to G_RRT and G_RRT, for which the agreement is much better. Nevertheless, the difference in the total stellar mass amongst the runs is smaller than 10 per cent. Moreover, if we also take into account the initial stellar mass of the galaxy, $M_* = 2.4 \times 10^9 M_\odot$, the difference drops to less than 1 per cent.

In Fig. 11, we show the SF law for the same three runs. In this case, instead of the line contours, we binned our data points along the gas density axis in 0.2 decade-wide bins and reported the average SFR in each bin, using the standard deviation in the bins as the y-axis error bar. In order to exclude poorly sampled bins, we removed all the bins with less than five data points. We find very good agreement with observations by Bigiel et al. (2010) for the KS relation in total gas (H + H₂) up to $\Sigma_* = 10^{-2} M_\odot \text{yr}^{-1} \text{kpc}^{-2}$, the only exception being the data points in the central 2 kpc of the galaxy which lie slightly above the orange shaded area. However, we notice that the observational data we compared with are for local galaxies, where the highest surface densities have small statistics. If we compare our results with other data sets from higher redshift galaxies (Orr et al. 2017), our data points are perfectly consistent in the entire range we sample. Both slope and normalization of the H₂ SF law are now in agreement with the observational points. However, in G_RRT, the relation is slightly shifted to slightly higher H₂ surface densities, and we also observe a slight deviation at low H₂ surface densities, due to the higher abundance of H₂.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** SF in our Giemsa, G_RRT and G_RAD_T runs in the latest 300 Myr of the simulations. We show the SFR in the left-hand panel and the total stellar mass formed in the right-hand panel. The green dotted line corresponds to G_RAD_T, the orange dashed line to G_RRT, and the blue solid line to G_RAD_S. Though the SFR shows bursts of SF of moderate intensity in all the runs, the SF history is quiet (~1 M⊙ yr⁻¹) during the entire simulations. A slightly lower SFR can be noticed in G_RAD_S, probably due to a stronger radiative flux close to the stellar sources compared to the other runs. As a consequence, the total stellar mass formed (right-hand panel) in G_RAD_S is a bit lower than in G_RRT and G_RAD_T.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** SF law for the G RAD_S, G_RRT and G_RAD_T runs in the last 100 Myr of the simulations, compared with the observational data as in Fig. 3. For the simulated data, we show here the average SFR, with its standard deviation σ reported in the vertical error bar. G RAD_S is plotted in cyan, G_RRT in blue and G RAD_T in red, respectively. The agreement in the left-hand panel is good, with all the runs perfectly lying within the Bigiel et al. (2010) contours. In the right-hand panel, instead, the differences are a bit larger, but almost everywhere within the simulations error bars. The slope of the relation is consistent with the observational data, with a small deviation in G RAD_S for very low H₂ surface densities, consistent with the higher H₂ abundance at intermediate nH₂tot in Fig. 8 (top panels). The normalization is also in good agreement with observations over more than one order of magnitude, and slightly drifting upwards only above 10 M⊙ pc⁻², corresponding to the central 2 kpc of the galaxy.</figcaption></figure>

compared to the other runs, in particular at intermediate densities
(see Fig. 8, top-left panel).

**4.3 The role of the clumping factor**

We now compare the G_RT_1 and the G_RT runs, where we use
the same model for stellar radiation, but with a different clump-
ing factor. We show, in Fig. 12, the projected H₂ density maps
for G_RT (right-hand panels), compared with those of G_RT_1
(left-hand panels), face-on and edge-on. Expectedly, much more

H₂ forms in G_RT, because of the enhanced formation rate, but
the spiral arm pattern is very similar in the two runs, suggesting
that the effect on the global dynamics is almost negligible. None
the less, most of H₂ still resides in the dense regions of the galaxy
spiral arms, where the gas is expected to be colder and turbulent.
In the edge-on view, G_RT exhibits a thicker molecular disc, where
H₂ is abundant up to ~500 pc from the disc plane, whereas in
G_RT_1 it does not extend to more than 100–200 pc. In Fig. 13,
we compare the gas phase diagrams for the same runs. In partic-
ular, the transition from the neutral to the fully molecular phase
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Projected gas density maps of the H₂ component for the G_RT_1 (left-hand panels) and G_RT (right-hand panels) runs at *t* = 350 Myr, face-on (top panels) and edge-on (bottom panels). The clumping factor in the fiducial run enhances the amount of H₂ forming in the galaxy, producing a thicker molecular disc (bottom panels), and also a more spread distribution across the entire galaxy, although the highest concentration still resides in the dense regions of the spiral arms.</figcaption></figure>

(top panels) shifts from ~5 (right-hand panel) to ~100 cm⁻³ (left-hand panel), because of the suppressed H₂ formation rate when Cₚ = 1. This is consistent with previous results (e.g. Gnedin et al. 2009; Christensen et al. 2012), where the gas becomes fully molecular above 100 cm⁻³. In our fiducial runs, instead, the transition occurs at lower densities because of the clumping factor. However, we must keep in mind that the density measured in the simulation is the mean density of the clouds for which we assume a log-normal PDF and most of the molecular gas actually forms in the high-density tail of the PDF, above 100 cm⁻³. In the bottom panels, instead, we compare the temperature distribution of the gas. G_RT exhibits a moderately larger amount of cold gas below 100 K, consistent with the boosted H₂ formation which provides a slightly more effective cooling. Tiny differences are also visible in the warm gas at intermediate densities, where the distribution is more spread in the G_RT run, but they can be compatible with stochastic fluctuations of a few thousand particles and are not relevant.

Fig. 14, instead, shows the SF law in the total (H + H₂) and molecular gas, following the same approach used in the previous section. As expected, the KS relation in total gas is almost unaffected by the clumping factor, only depending on the net amount of neutral hydrogen. On the other hand, the molecular counterpart is significantly affected, moving by roughly an order of magnitude between the two runs. This result suggests that the resolution achieved is not sufficient to properly resolve the very high densities where the gas becomes fully molecular all over the galaxy, and a subgrid clumping factor is necessary to account for it. It is remarkable how our clumping factor model compensates for the offset without any calibration, shifting the simulation results over the observed relation.

We find that the clumping factor can vary significantly in the simulation, depending on the gas properties, from 1 to up to a maximum of about 10³, but the distribution is far from being uniform (see Fig. 15). Comparing the clumping factor distribution at different times we also found that the exact moment chosen could vary the distribution by a few per cent, if we look for example at a few tens of Myr after a starburst or a more quiescent phase. In particular, roughly 38 ± 3 per cent of the gas is typically below Cₚ = 10 and only ~10 ± 2 per cent exceeds Cₚ = 100. The median of the distribution is 18 ± 3 and the average is ~35 ± 5. Although values above 100 can appear unphysical, we emphasize that they come from gas with 25 ≲ *M* ≲ 100 (assuming *b* = 0.4), extreme but still plausible. For the sake of simplicity, we do not discuss in detail the runs G RAD_S_1 and G RAD_T_1 as the same differences observed between G_RT and G_RT_1 are reflected into the other cases.

# 5 CAVEATS OF THE STUDY

Our results showed that the relation between SF and H₂ abundance is sensitive to the physical modelling in the simulations, in particular to the stellar radiative feedback. Though not explored in this study, it is possible that the relation would also be affected by different SF and SN feedback models.

(i) *The SF model:* We used here a physically motivated parameter-free prescription which does not require a priori calibrations, unlike the standard one used in galaxy simulations, where density and (sometimes) temperature cuts are employed instead. In addition, the SF efficiency is here self-consistently computed from the gas properties, and not chosen by hand to reproduce the KS
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Gas properties for the G_RT_1 (left-hand column) and G_RT (right-hand column) runs averaged over 100 Myr around *t* = 350 Myr. We show the molecular hydrogen mass fraction in the top row and the gas temperature in the bottom row, as a function of the total hydrogen density. The critical density to become fully molecular shifts from ~5 cm⁻³ (top-left panel) to ~100 cm⁻³ (top-right panel), due to the effect of *C*ρ, which boosts H₂ formation also at lower hydrogen densities. On the other hand, the temperature distribution of the gas is almost unaffected, with only a moderately larger amount of cold gas below 100 K around *n*Htot = 1 cm⁻³.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Same as Fig. 11, but comparing the G_RT_1 and G_RT runs. While the KS relation in the left-hand panel does not show any variation relative to the clumping factor, the H₂ SF law very well matches the slope of the observed relation, but the normalization varies by roughly an order of magnitude between the two runs. This suggests that our resolution is not high enough to properly model the very high density gas where most of H₂ forms, and a clumping factor is still necessary to properly match the results.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Clumping factor distribution for G_RT at *t* = 350 Myr, binning the gas particle distribution with 0.02 decade-wide bins on both axes. For completeness, we also show the corresponding values for *M*, obtained reversing equation (4). Though a large fraction of gas clusters at low *C*ρ, consistent with a pressure-dominated regime, a moderate amount of gas below *n*<sub>Htot</sub> = 100 cm<sup>-3</sup> exhibits a turbulent-dominated state, with *M* up to 50, where the H<sub>2</sub> formation rate is strongly enhanced.</figcaption></figure>

relation. However, a main limitation of the model is that it holds only for turbulent clouds, which means that the star-forming gas must satisfy $M \ge 2$. In spite of this, we stress that this condition does not have any impact on the results and, moreover, stars are expected to physically form in turbulent molecular clouds where $M \gg 1$, as found in the nearby starburst NGC253 (Leroy et al. 2015), and not in the sub-sonic ISM; hence, the employed cut has a clear physical motivation. On top of that, we also put a density cut to prevent SF in regions where the gas is expected to be warm and not star-forming, and this could appear as an arbitrary choice. However, the motivation for this choice resides in the strong dependence of the SF efficiency on the gas density, with values as low as $10^{-4}$ to $10^{-3}$ around the SF density threshold, as can be seen from Fig. A1. For instance, at our density threshold the expected SFR is of the order of $10^{-8}$ to $10^{-7}$ M$_\odot$ yr$^{-1}$, corresponding to less than one stellar particle during the entire evolution of the galaxy. By using a higher SF density threshold, we would have introduced an additional free parameter which would have simply made the SF more bursty, but with almost no difference in the stellar spatial distribution and in the resulting KS relation. At high densities, instead, the high SF efficiency could have played a major role in reproducing the SF–H$_2$ relation. Actually, it is possible that the high-density gas was turned into stars very rapidly, and the H$_2$ associated with the gas particle instantaneously removed, reducing the total amount of H$_2$ available. This effect could have probably been attenuated using the standard SF prescription, where the SF efficiency is kept constant also for high density gas, but it is not obvious whether the final result would have been significantly different. To test this idea, we evolved our target galaxy for 200 Myr using the same model of our G RAD_T run, but with the standard SF prescription, assuming $\rho_{thr} = 100$ cm$^{-3}$ and a constant SF efficiency $\varepsilon = 0.03$. Though not shown here, the SF–H$_2$ law did not differ significantly from that of our fiducial model.

(ii) *The SN feedback model:* In this study, we used the empiri-
cal delayed-cooling model, which was shown to better reproduce
the KS relation (Rosdahl et al. 2017), at the expense of the ISM
temperature. In particular, the gas heated up by SNe would remain
warm/hot even at densities of a few cm⁻³, where radiative cooling
should be important, and this could affect the abundance of the

chemical elements. However, although the H₂ abundance for the
intermediate density gas kept hot by SNe could have been underes-
timated in our simulations, it is also true that the gas should become
fully molecular around n<sub>H,tot</sub> ~ 100 cm<sup>-3</sup>, where the amount of
warm gas is subdominant compared to the cold counterpart (see for
instance Fig. 8). In conclusion, further detailed investigations about
the effect of different SN feedback models are needed to assess
the main advantages and limitations of each model and they could
help determining the model which better reproduces observations.
Nevertheless, this is beyond the scope of the present study.

Finally, other numerical aspects could also play a role in how
well we can reproduce the molecular KS relation, like the hydro-
dynamics method and the artificial pressure floor (not employed
here). Though these effects have not been addressed here, we refer
to Capelo et al. (2017), where the same galaxy is evolved using
a different hydrodynamics scheme (SPH; with GASOLINE2 Wadsley,
Keller & Quinn 2017), a pressure floor, and with different chemistry
networks to assess the role of non-equilibrium metal cooling at low
temperature.

6 CONCLUSIONS

We introduced a model to accurately follow the chemical evo-
lution of molecular hydrogen and performed a suite of simula-
tions of an isolated galaxy at z = 3 to test its validity. To self-
consistently follow H₂, we coupled the chemistry package KROME
with the mesh-less code GIZMO. We used model 1a in Bovino et al.
(2016), including all the relevant processes regulating H₂ forma-
tion and dissociation, together with a new physically motivated
SF prescription based on theoretical models of turbulent GMCs;
delayed-cooling SN feedback; gas, dust and H₂ shielding; extra-
galactic UV background and local UV stellar radiation. The latter
is modelled with sub-grid models (Giemans and Giesy) and
compared to RT calculations (GRT). We also included a cell-by-
cell estimate of a clumping factor, the enhancement in the formation
of H₂ in the very high density cores of the GMCs, unresolved in our
simulations.

We showed that the external UV background has the only ef-
fect of destroying H₂ at very low density, and the most important
contribution to H₂ dissociation comes from stellar radiation inside
the galaxy. The two sub-grid models for local UV stellar radiation
perform similarly well against the RT simulation, and have a much
higher computational efficiency.

The KS relation in the total gas matches well observations across
several orders of magnitudes, with only a small increase of the SFR
at higher densities. As for the H₂ counterpart, the slope is generically
in good agreement with the observational data, but the normalization
is recovered correctly only when including a clumping factor to
account for unresolved dense cores in GMCs.

In summary, we showed that, with a detailed chemical network
and a physically motivated SF prescription (with a consistent sub-
grid clumping factor), we have been able to naturally reproduce the
observed relation between the H₂ and SFR surface densities for a
large range of Σ<sub>H₂</sub>. In particular, we confirm that this relation can
naturally arise when an opportune model is used, and an a priori
dependence of SF on the H₂ abundance, as assumed in several
studies (e.g. Gnedin et al. 2009; Christensen et al. 2012; Hopkins
et al. 2017; Pallottini et al. 2017) is not necessary. We have also
found that the relation seems to extend down to lower H₂ densities
maintaining the same slope, a regime not explored yet with the
available observations.
ACKNOWLEDGEMENTS

We thank the anonymous referee for the useful suggestions which
helped improve the study. We thank Philip F. Hopkins for provid-
ing the RT implementation in GIZMO and for useful discussions.
We thank Sijing Shen for having provided the metal cooling ta-
ble described in Shen et al. (2010, 2013). We acknowledge support
from the European Research Council (Project No. 267117, ‘DARK’,
AL, JS; Project no. 614199, ‘BLACK’, AL, MV). SB thanks for
funding through the DFG priority program ‘The Physics of the
Interstellar Medium’ (projects BO 4113/1-2). PRC acknowledges
research funding by the Deutsche Forschungsgemeinschaft (DFG)
Sonderforschungsbereiche (SFB) 963 and support by the Tomalla
Foundation. This work was granted access to the High Performance
Computing resources of CINES under the allocation x2016046955
made by GENCI, and it has made use of the Horizon Cluster, funded
by Institut d'Astrophysique de Paris, for the analysis of the simula-
tion results.

REFERENCES

Adelberger K. L., Steidel C. C., Pettini M., Shapley A. E., Reddy N. A., Erb D. K., 2005, ApJ, 619, 697

Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

Bigiel F., Leroy A., Walter F., Blitz L., Brinks E., de Blok W. J. G., Madore B., 2010, AJ, 140, 1194

Bovino S., Grassi T., Capelo P. R., Schleicher D. R. G., Banerjee R., 2016, A&A, 590, A15

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000 (BC03)

Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., Governato F., 2015, MNRAS, 447, 2123

Capelo P. R., Bovino S., Lupi A., Schleicher D. R. G., Grassi T., 2017, MNRAS, in press, preprint (arXiv:1710.01302)

Chabrier G., 2003, PASP, 115, 763

Chandrasekhar S., 1931, ApJ, 74, 81

Christensen C., Quinn T., Governato F., Stilp A., Shen S., Wadsley J., 2012, MNRAS, 425, 3058

Ciotti L., D'Ercole A., Pellegrini S., Renzini A., 1991, ApJ, 376, 380

Clark P. C., Glover S. C. O., Klessen R. S., Bonnell I. A., 2012, MNRAS, 424, 2599

Fedderrath C., Klessen R. S., 2012, ApJ, 761, 156

Fedderrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364

Fedderrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010, A&A, 512, A81

Feldmann R., Gnedin N. Y., Kravtsov A. V., 2011, ApJ, 732, 115

Ferland G. J. et al., 2013, RMxAA, 49, 137

Fiacconi D., Feldmann R., Mayer L., 2015, MNRAS, 446, 1957

Forbes J. C., Krumholz M. R., Goldbaum N. J., Dekel A., 2016, Nature, 535, 523

Forrey R. C., 2013, ApJ, 773, L25

Glover S. C. O., Abel T., 2008, MNRAS, 388, 1627

Glover S. C. O., Savin D. W., 2009, MNRAS, 393, 911

Gnedin N. Y., Abel T., 2001, New Astron., 6, 437

Gnedin N. Y., Tassis K., Kravtsov A. V., 2009, ApJ, 697, 55

Grassi T., Bovino S., Schleicher D. R. G., Prieto J., Seifried D., Simoncini E., Gianturco F. A., 2014, MNRAS, 439, 2386

Haardt F., Madau P., 2012, ApJ, 746, 125

Harten A., Lax P. D., van Leer B., 1983, SIAM Rev., 25, 35

Hartwig T., Glover S. C. O., Klessen R. S., Latif M. A., Volonteri M., 2015, MNRAS, 452, 1233

Heiderman A., Evans N. J., II, Allen L. E., Huard T., Heyer M., 2010, ApJ, 723, 1019

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Hernquist L., 1990, ApJ, 356, 359

Hopkins P. F., 2015, MNRAS, 450, 53

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hopkins P. F. et al., 2017, preprint (arXiv:1702.06148)

Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, MNRAS, 458, 3528

Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151

Jura M., 1975, ApJ, 197, 575

Kalirai J. S., Hansen B. M. S., Kelson D. D., Reitzel D. B., Rich R. M., Richer H. B., 2008, ApJ, 676, 594

Katz H., Kimm T., Sijacki D., Haehnelt M. G., 2017, MNRAS, 468, 4831

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kim J.-h. et al., 2014, ApJS, 210, 14

Kim J.-h. et al., 2016, Astrophys. J., 833, 202

Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250

Krumholz M. R., Leroy A. K., McKee C. F., 2011, ApJ, 731, 25

Lax P. D., 1954, Commun. Pure Appl. Math., 7, 159

Leitner S. N., Kravtsov A. V., 2011, ApJ, 734, 48

Leroy A. K. et al., 2015, ApJ, 801, 25

Levermore C., 1984, J. Quant. Spectrosc. Radiat. Transfer, 31, 149

Lupi A., Volonteri M., Silk J., 2017, MNRAS, 470, 1673

Mac Low M.-M., 2016, in Jablonka P., André P., van der Tak F., eds, Proc. IAU Symp. 315, From Interstellar Clouds to Star-Forming Galaxies: Universal Processes? International Astronomical Union, Paris, France, p. 1

Maoz D., Mannucci F., Brandt T. D., 2012, MNRAS, 426, 3282

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Orr M. et al., 2017, MNRAS, preprint (arXiv:1701.01788)

Padoan P., Nordlund Å., 2011, ApJ, 730, 40 (PN)

Pallottini A., Ferrara A., Gallerani S., Vallini L., Maiolino R., Salvadori S., 2017, MNRAS, 465, 2540

Pelupessy F. I., Papadopoulos P. P., van der Werf P., 2006, ApJ, 645, 1024

Richings A. J., Schaye J., 2016, MNRAS, 458, 270

Richings A. J., Schaye J., Oppenheimer B. D., 2014, MNRAS, 442, 2780

Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, MNRAS, 466, 11

Safranek-Shrader C., Krumholz M. R., Kim C.-G., Ostriker E. C., Klein R. I., Li S., McKee C. F., Stone J. M., 2017, MNRAS, 465, 885

Salim S. et al., 2007, ApJS, 173, 267

Schmidt M., 1959, ApJ, 129, 243

Schruba A. et al., 2011, AJ, 142, 37

Schruba A. et al., 2012, AJ, 143, 138

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, ApJ, 765, 89

Somerville R. S., Popping G., Trager S. C., 2015, MNRAS, 453, 4337

Spitzer L., Jr, 1942, ApJ, 95, 329

Springel V., 2005, MNRAS, 364, 1105

Springel V., 2010, MNRAS, 401, 791

Springel V., White S. D. M., 1999, MNRAS, 307, 162

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Strömgren B., 1939, ApJ, 89, 526

Tacconi L. J. et al., 2010, Nature, 463, 781

Thompson R., Nagamine K., Jaacks J., Choi J.-H., 2014, ApJ, 780, 145

Tomassetti M., Porciani C., Romano-Díaz E., Ludlow A. D., 2015, MNRAS, 446, 3330

Tumlinson J. et al., 2002, ApJ, 566, 857

Voit G. M., Donahue M., 2011, ApJ, 738, L24

Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS, 471, 2357

Wolfire M. G., Tielens A. G. G. M., Hollenbach D., Kaufman M. J., 2008, ApJ, 680, 384

Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269
APPENDIX A: STAR FORMATION
EFFICIENCY

In this study, we have implemented a new prescription for SF based on the theoretical model by Padoan & Nordlund (2011) and calibrated against simulations and observations by Federrath & Klessen (2012). Unlike the standard model used in numerical simulations, where the SF efficiency $\epsilon$ is constant and must be tuned to get the right KS normalization, this model self-consistently gives an SF efficiency dependent on the grav-o-turbulent state of the gas. According to Federrath & Klessen (2012), the SF efficiency can vary across several orders of magnitude and, in principle, can also exceed 100 per cent when the GMC is strongly bound ($\alpha_{vir} < 1$) and very turbulent ($M \gg 10$), with compressive modes dominating over the solenoidal modes.<fnref n="9" /> To evaluate the range spanned by the SF efficiency in our simulations and the dependence on the gas properties, we plotted in Fig. A1 the SF efficiency from our G RAD T run at $t = 350$ Myr, as a function of the total hydrogen density (left-hand panel) and gas temperature (right-hand panel), binning all the quantities in 0.02 decade-wide bins. The SF efficiency around the SF density threshold $n_{Htot,thr} = 1 \text{ cm}^{-3}$ is much lower than the typically used value of 0.01 (most of the gas has $\epsilon \lesssim 10^{-3}$) and rapidly increases with the gas density up to more than 10 per cent<fnref n="10" /> above $n_{Htot} = 10^3 \text{ cm}^{-3}$. This suggests that $\rho_{SF,thr}$ is not a crucial parameter of the model, because of the rapid decay of $\epsilon$ at very low densities. A moderate value of $n_{Htot}$ as low as $1 \text{ cm}^{-3}$ is, however, a reasonable choice to avoid spending time computing the SFR in regions where stellar particles would never form. From the SF efficiency distribution it seems that we never approach 100 per cent efficiency, even when gas is fully molecular (cf. Fig. 13). However, it could simply be a selection effect due to the fact that at very high SF efficiencies the conversion from gas into stars is almost instantaneous, yielding to the disappearance of the particle from the SF efficiency plot shown here. As for the right-hand panel, instead, the dependence on the temperature is extremely steep, with all the gas particles with $\epsilon > 10^{-5}$ being colder than $T = 10^3$ K and the hotter ones showing almost negligible $\epsilon$. The model therefore provides a natural temperature threshold for SF, without the need to impose a strict criterion as usually done with the standard prescription. Moreover, the lack of a sharp vertical cut-off in the temperature–SF efficiency plot reveals that the Mach number criterion does not play a fundamental role, being most of the cold (and highly star-forming) gas in a turbulent-dominated state.

In Fig. A2, we show the SF efficiency as a function of H₂ density (left-hand panel) and the renormalized SF efficiency $\epsilon_0 = \epsilon/f_{H_2}$ (right-hand panel), using the same binning as in Fig. A1. We overplot, in both panels, the median of the SF efficiency distribution from the left-hand panel of Fig. A1, to highlight the change in the distribution when we introduce a dependence on the H₂ abundance. The SF efficiency, mildly increasing up to $n_{H_2} \sim 0.1 \text{ cm}^{-3}$, with an average value around $10^{-4}$, shows a steeper slope at higher densities, where most of the SF occurs. $\epsilon_0$ has a very large scatter below $n_{Htot} \sim 30-40 \text{ cm}^{-3}$, where stars form in the atomic hydrogen regime with a very weak dependence on the abundance of H₂, whereas a clear overlap with the median curve (in red) is visible at

<figure><img src="image_19.png" /><figcaption>**Figure A1.** SF efficiency versus total hydrogen density (left-hand panel) and gas temperature (right-hand panel) for our G RAD T run at $t = 350$ Myr. The SF efficiency rapidly increases with increasing density and is strongly suppressed in the warm gas, consistent with the fact that SF mainly occurs in dense, cold and turbulent clouds.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure A2.** SF efficiency versus H₂ density (left-hand panel) and renormalized SF efficiency versus $n_{Htot}$ (right-hand panel) for our G RAD T run at $t = 350$ Myr. The red curve is the median of the SF efficiency from the left-hand panel in Fig. A1, reported as a function of total hydrogen density in both panels. In the left-hand panel, it is rescaled by a factor of 2 to account for the number of hydrogen nuclei in the H₂ molecule and be consistent with the x-axis of the plot. Though the SF efficiency shows a strong dependence on the H₂ abundance, there is no apparent causal connection between the H₂ and SF. When renormalized to the H₂ abundance, the increase is still present, but a large scatter appears for low densities where SF proceeds in the H-dominated regime.</figcaption></figure>

higher densities. This is expected, since we are in the regime where gas becomes fully molecular ($f_{H_2} \approx 1$) and the normalization factor has a negligible effect. Even at high $n_{H_2}$, $\epsilon_0$ never reaches a 100 per cent efficiency in our simulations, but it self-regulates to always be below few tens per cent.

As a benchmark to compare the SF prescription implemented in this study with the standard one, where stars form in gas with $n_{Htot} > n_{Htot,thr}$ with a constant $\epsilon$, we compute the average SF efficiency

<footnote n="9">⁹ This is due to the definition of the SF efficiency per free-fall time, which is normalized with the total mass and the free-fall time at the mean density of the GMC.</footnote>

<footnote n="10">¹⁰ This only occurs for a few particles at very high densities, which cannot be seen in the plot.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure A3.** Average SF efficiency as a function of the density threshold $n_{H_{tot,thr}}$. The black solid line is the result from the averaging procedure applied to the bins in the left-hand panel of Fig. A1 and the orange dashed one is the best fit with the analytic function defined in equation (A1).</figcaption></figure>

<table><caption>**Table A1.** Parameter values for the best-fitting function to $\langle\epsilon\rangle$ defined in equation (A1).</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>α</td><td>0.04415</td></tr><tr><td>β</td><td>1.642</td></tr><tr><td>φ</td><td>0.2507</td></tr><tr><td>ξ</td><td>0.004509</td></tr></tbody></table>

$\langle\epsilon\rangle$ for gas above a total hydrogen density $n_{H_{tot,thr}}$ by stacking the snapshots of the last 100 Myr of the simulation, as

$$ \langle\epsilon\rangle = \frac{\sum_{i=i_0}^{N_i} \left[ \sum_{\epsilon_j=0}^{N_j} \epsilon_{i;j} M_{i;j} \right] / t_{\text{ff}}(n_i)}{\sum_{i=i_0}^{N_j} \left[ \sum_{j=0}^{N_j} M_{i;j} \right] / t_{\text{ff}}(n_i)}, \quad (A1) $$

where *i* runs on the $N_i$ density bins, *j* runs on the $N_j$ SF efficiency bins, $i_0$ corresponds to the $n_{H_{tot,thr}}$ bin, $M_{(i;j)}$ is the gas mass in the (i; j)-th bin, $\epsilon_j$ is the SF efficiency in the j-th bin and $t_{\text{ff}}(n_i)$ is the free-fall time in the i-th bin. Fig. A3 shows the resulting $\langle\epsilon\rangle$, together with the best fit, expressed by the analytic expression

$$ \langle\epsilon\rangle = \alpha \exp\{\beta \arctan[\phi \ln(\xi n_{H_{tot,thr}})]\}. \quad (A2) $$

We report the best-fitting parameters in Table A1. The fitting function used shows a perfect agreement across four orders of magnitude with the simulation data and varies between the usual SF efficiency value of 0.01 and $\gtrsim$ 0.14. While it is common practice to increase

the SF efficiency when a high SF density threshold is used, the only way to tune it is by rescaling it to match the local KS relation. Here we provide an alternative way to do it, by means of an analytic fitting function, aimed at reproducing, on average, the same SFR but using a simpler SF prescription with constant SF efficiency.

## APPENDIX B: THE LOW-RESOLUTION CASE

We report here the results of two low-resolution runs analogous to the high-resolution ones described in the main sections of the paper, to assess how sensitive our clumping factor model is to the numerical resolution. These simulations are run using model (b) for stellar radiation, which showed a slightly better agreement with our G_RT run.

The galaxy is here modelled using 100k gas particles, 500k DM particles and 100k stellar particles in the disc and 48k stellar particles in the bulge, corresponding to a mass resolution of 23910 M$_{\odot}$ for gas, 16600 and 15940 M$_{\odot}$ for stellar disc and bulge, respectively, and $5.9 \times 10^5$ M$_{\odot}$ for DM. The softening length is rescaled here to 10 pc for stars and 6 pc (minimum) for gas.

For the sake of simplicity, we do not report here the full analysis done for the high-resolution runs, but we limit the discussion to the KS relation only. We show, in Fig. B1, the results of the two low-resolution runs, with and without our clumping factor model, obtained by binning along the gas density axis as done for Fig. 11, and we compare them with those of the corresponding high-resolution runs G_RAD_T_1 and G_RAD_T. For the KS relation in the total gas (H + H$_2$), the results are in very good agreement, with small differences very well within the error bars. For the relation in the molecular gas, while both the low-resolution (blue dashed) and the high-resolution (red dotted) runs without clumping factor are offset from the observed relation, the offset is more pronounced in the low-resolution case. When we use the clumping factor model (cyan solid and green dot-dashed lines), instead, the difference almost vanishes, becoming much smaller than the measured scatter, and, moreover, both curves become consistent with a depletion time of ~1 Gyr, as found in observations.

This suggests that, although the agreement slightly improves increasing the resolution, the model is still valid in the low-resolution case. However, to ensure the validity of the gas log-normal PDF, we must be able to resolve single GMCs, at least the more common ones. If this is not the case, both the clumping factor and the SF prescription become invalid and should not be used. We therefore suggest a minimum gas mass resolution of a few $10^4$ M$_{\odot}$, according to the observations by Roman-Duval et al. (2010) of 580 molecular clouds in the University of Massachusetts-Stony Brook Survey, where only ~15 per cent is more massive than $10^5$ M$_{\odot}$, and most of the GMCs (~42 per cent of the total) are in the range $10^4-10^5$ M$_{\odot}$.
<figure><img src="image_22.png" /></figure>

**Figure B1.** SF law for the low- and high-resolution runs in the last 100 Myr of the simulations, with and without the clumping factor. The high-resolution cases correspond to the Giemsa and G_RAD_T runs described in the previous sections. We show the relation for the total gas (H + H₂) in the top panel and that for H₂ only in the bottom panel, following the same approach of Fig. 11. All runs agree on the total gas KS relation, independently of the clumping factor and the resolution. The H₂ SF law is instead significantly offset from the observational data, for the runs without clumping factor, and this is more evident in the low-resolution run. When the clumping factor is used, the agreement improves, bringing both low- and high-resolution runs within the respective error bars.

This paper has been typeset from a TeX/LA TeX file prepared by the author.