© The Authors 2024

1 Introduction

Molecular gas plays a major role in the interstellar medium (ISM) of galaxies, providing the necessary conditions and likely serving as the primary fuel for star formation. The cosmic molecular gas density in the Universe, as inferred from blind surveys such as the VLA CO Luminosity Density at High Redshift¹ (COLDz, Riechers et al. 2019) and the ALMA Spectroscopic Survey in the HUDF² (ASPECS, Decarli et al. 2019; Walter et al. 2020), increases by roughly an order of magnitude between red-shifts z ~ 6 and z ~ 2. This is accompanied by a similar trend in the cosmic star formation rate density that reaches its peak value at z ~ 3−2, a period known as ‘cosmic noon’ (see Madau & Dickinson 2014; Förster Schreiber & Wuyts 2020, for a review).

Investigating the build-up of the molecular gas reservoir in galaxies and its cosmic evolution is therefore pivotal to our understanding of the star formation history and galaxy assembly in the Universe. Molecular hydrogen (H₂) is the dominant molecular gas component in galaxies. However, because of the lack of a permanent dipole moment and high excitation temperatures (T ≳ 500 K) for its ro-vibrational transitions, H₂ does not emit light under typical conditions of molecular clouds (T ≲ 100 K). Therefore, the presence and mass of H₂ are routinely inferred via emission from tracers such as dust, CO, and atomic fine-structure lines.

The low-J rotational transitions of CO are the most commonly used tracers of molecular gas in galaxies (e.g. Dickman et al. 1986; Solomon & Barrett 1991; Downes & Solomon 1998; Solomon & Vanden Bout 2005; Tacconi et al. 2006, 2008, 2010; Daddi et al. 2010a,b; Genzel et al. 2010; Bolatto et al. 2013). The CO-to-H₂ conversion factor α_CO captures the relation between the observed CO J = 1 → 0 luminosity of a galaxy and the underlying molecular gas mass. Alternatively, for spatially resolved observations within the Milky Way (MW) and nearby galaxies, X_CO relates the CO intensity (W_CO) to the H₂ column density $\left(N_{{\text{H}}_2}\right)$ along the line of sight. The variation of X_CO with physical conditions, such as metallicity, stellar surface density, galactocentric distance, among others, has been extensively investigated in the MW and nearby galaxies (see Bolatto et al. 2013, for a review). At higher redshifts, however, the limited number of galaxies observed in CO poses a challenge to studying the dependence of α_CO on other galaxy properties. This is further complicated by the fact that the CO J = 1 → 0 transition is not accessible at high redshifts and observers have to rely on higher-J transitions to obtain an estimate for it. This requires knowledge about the CO excitation ladder, thereby introducing another systematic uncertainty in employing CO as a molecular gas tracer. Moreover, in some galaxies (e.g. low-metallicity dwarf galaxies), a large amount of H₂ is not traced by CO emission and is referred to as CO-dark molecular gas (Wolfire et al. 2010; Madden et al. 2020).

Other tracers of molecular gas such as the [C II] fine-structure line suffer from similar systematic effects. The nature and origin ofthis line are highly debated as it can arise from multiple ISM phases and not all of the [C II] emission of a galaxy is associated with the molecular gas phase. The presence of a [C II] halo in some galaxies extending 2−3 times farther than their rest-frame UV emission (see e.g. Fujimoto et al. 2020) hints at an extended C⁺ reservoir devoid of molecular gas. Additionally, the decrease in the [C II]/far-infrared (FIR) luminosity ratio with increasing FIR luminosity, known as the ‘[C II] deficit’, further complicates the use of this line as a reliable tracer of molecular gas across galaxies.

Atomic carbon has been proposed as another reliable tracer of molecular gas. The fine structure lines of atomic carbon are expected to trace the bulk of the molecular gas in galaxies (Papadopoulos et al. 2004; Tomassetti et al. 2014; Glover & Clark 2016) and have been used both in the local Universe as well as at high redshifts (e.g. Gerin & Phillips 2000; Ikeda et al. 2002; Weiβ et al. 2003, 2005; Walter et al. 2011; Valentino et al. 2018; Henríquez-Brocal et al. 2022). However, the use of these lines requires assumptions about the relative abundance of atomic carbon and molecular hydrogen.

Cosmological simulations are a useful tool for investigating the reliability of molecular gas tracers under different ISM conditions and galaxy environments. However, simulating the molecular gas content of galaxies is challenging as it requires modelling the various physical and chemical processes happening on a wide range of (spatial and temporal) scales. On the one hand, it is necessary to simulate galaxies in realistic environments as their ISM is affected by outflows and gas accretion from the cosmic web. On the other hand, molecular-cloud chemistry is regulated by conditions on sub-parsec scales. Early attempts at modelling H₂ assumed a local chemical equilibrium between H₂ formation and destruction (Krumholz et al. 2008, 2009; McKee & Krumholz 2010) and these models have found extensive application in cosmological simulations (Kuhlen et al. 2012, 2013; Hopkins et al. 2014; Thompson et al. 2014; Lagos et al. 2015; Davé et al. 2016). Despite their wide use, these equilibrium models do not account for the dynamic nature of the ISM and the long formation timescale for H₂ (Tielens & Hollenbach 1985). For example, Pelupessy & Papadopoulos (2009), Tomassetti et al. (2015), Richings & Schaye (2016), Pallottini et al. (2017), Schabe et al. (2020), and Hu et al. (2021) provide a detailed discussion on the effect of non-equilibrium H₂ chemistry on the integrated properties of simulated galaxies. However, cosmological hydrodynamical simulations with on-the-fly computations of the non-equilibrium chemical abundances are rare (e.g. Dobbs et al. 2008; Christensen et al. 2012; Tomassetti et al. 2015; Semenov et al. 2018; Lupi et al. 2018; Lupi 2019; Schabe et al. 2020; Katz et al. 2022; Hu et al. 2023) and often restricted to individual galaxies. Some of these studies only include a non-equilibrium chemical network for HI and H₂ but not the tracers (Lupi & Bovino 2020; Lupi et al. 2020).

Despite tremendous improvement in the spatial resolution of cosmological simulations over the last decade, the state-of-the-art today is still far from resolving the clumpy ISM in molecular clouds. A clumpy ISM would allow for pockets of high-density gas where H₂ formation would be enhanced. This enhancement is missed by simulations as the density is uniform below their resolution scale. Gnedin et al. (2009) accounted for these unresolved densities by enhancing the H₂ formation rate by an effective clumping factor C assuming a log-normal density distribution for the gas. This technique was later tested by Micic et al. (2012) in their numerical study on the effect of the nature of turbulence on H₂ formation. They found that using a clumping factor systematically overpredicts the H₂ formation rate in regions with a high molecular fraction $\left(f_{{\text{H}}_2}\gtrsim 0.5\right)$. Christensen et al. (2012) adapted this method in their smooth particle hydrodynamics (SPH) simulations to model the non-equilibrium abundance of H₂ in a cosmological simulation of a dwarf galaxy. However, similar to Micic et al. (2012), they cautioned that this approach works well for dwarf galaxies where fully molecular gas is rare but would need further modifications at high molecular fractions. Moreover, the finite resolution of simulations implies that the temperature tracked in simulations is an average over the resolution element, similar to any other quantity followed explicitly. This average fails to capture the true heterogeneous temperature distribution that is essential for determining the rates of chemical reactions taking place in the ISM.

To overcome these limitations, Tomassetti et al. (2015, hereafter T15) developed a sub-grid model that accounts for the unresolved density structure in simulations by assuming a log-normal probability distribution function (PDF) of sub-grid densities. They assigned a temperature to each sub-grid density in the PDF using a temperature-density relation from high-resolution simulations of molecular clouds (Glover & Mac Low 2007b). They obtained the H₂ abundance in each resolution element by evolving their chemical network at the sub-grid level and integrating over the density PDF. They found good agreement between different gas and stellar properties of their simulated galaxy at z = 2 and observations of high-redshift galaxies.

In this paper, we improve upon the work of T15 and introduce a new sub-grid model, HYACINTH – HYdrogen And Carbon chemistry in the INTerstellar medium in Hydro simulations – for on-the-fly computation of the non-equilibrium abundances of H₂ and associated carbon tracers (CO, C, and C⁺) within cosmological simulations of galaxy formation. The paper is organised as follows: in Sect. 2, we describe the components of HYACINTH and how it can be incorporated into cosmological simulations. A comparison of our chemical network with two more complex chemical networks and a photon-dominated region code is presented in Sect. 3.1. We further compare the chemical evolution from HYACINTH against high-resolution simulations of molecular clouds – the SILCC-Zoom simulations (Seifried et al. 2017, 2020) and the Glover & Mac Low (2011) simulations in Sect. 3.2. Although HYACINTH is primarily designed to be embedded as a sub-grid model in cosmological simulations, as an immediate application, we applied it to a galaxy simulation (from T15) in postprocessing and directly compared it with observations related to the abundances of H₂, CO, C, and C⁺ in nearby and high-redshift galaxies. These are discussed in Sect. 4. In Sect. 5, we elaborate on the caveats involved in using HYACINTH as a post-processing tool and compare this with other approaches in the literature. A summary of our findings is presented in Sect. 6.

2 Methods

At the core of HYACINTH is a PDF of sub-grid densities and a simplified chemical network. The PDF accounts for the unresolved density structure in simulations, statistically incorporating its impact on the chemistry at resolved scales. The chemical network for hydrogen and carbon chemistry is designed to be efficient and scalable for integration into large-scale galaxy formation simulations. Additionally, HYACINTH can be used as a post-processing tool as well; see Sect. 4 for a sample application. In the following, we describe the technical specifications of these two components.

2.1 The sub-grid density PDF

The density structure of molecular clouds is governed by the interplay between turbulence and self-gravity. The PDF of densities is an important statistical property that describes this structure. For instance, the mass-weighted PDF gives the probability that an infinitesimal mass element dM has a density in the range [ρ,ρ + dp]. This distribution is expected to take a log-normal shape in an isothermal, turbulent medium, not significantly affected by the self-gravity of gas (see e.g. Vázquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; McKee & Ostriker 2007, for a review). Using near-infrared dust extinction mapping of nearby molecular clouds, Kainulainen et al. (2009) found that the column density PDF in quiescent clouds is very well described by a log-normal. However, in star-forming clouds, they found large deviations from a log-normal: the PDF in these clouds has a log-normal shape only at low column densities and resembles a power law at high column densities. A similar picture is supported by numerical simulations of turbulent and self-gravitating gas (e.g. Nordlund & Padoan 1999; Klessen 2000; Glover & Mac Low 2007b; Ballesteros-Paredes et al. 2011; Ward et al. 2014). On scales where the turbulence is supersonic and self-gravity is unimportant, the gas follows a log-normal distribution. In regions where gas is collapsing under its own gravity, the density structure is log-normal at low densities and develops a power-law tail at high densities. Motivated by these findings, we have modified the original PDF of T15 and use a log-normal PDF for turbulence-dominated regions and a log-normal+power-law PDF for gravity-dominated regions.

The mass-weighted log-normal PDF of sub-grid densities n_H is given by (T15) $${\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)=\frac{1}{\sqrt{2\pi }\sigma n_{\text{H}}}\exp \left[-\frac{{\left(\ln n_{\text{H}}-\mu \right)}^2}{2{\sigma }^2}\right],$$(1)

where σ and µ are parameters that decide the width and the location of the peak of the distribution. It is often convenient to introduce the clumping factor $$C=\langle n_{\text{H}}^2\rangle /{\langle n_{\text{H}}\rangle }^2$$(2)

that captures the inhomogeneity and the degree of clumpiness in the medium. For a log-normal PDF, it is related to the parameter σ above as $C={\text{e}}^{{\sigma }^2}$. We use a constant clumping factor of 10, which has previously been shown to reproduce the observed H₂ fractions in nearby galaxies (see e.g. Gnedin et al. 2009; Christensen et al. 2012, T15). We note, however, that it is possible to adopt a more sophisticated approach by varying C with the local Mach number ℳ or the 3D velocity dispersion (see e.g. Lupi et al. 2018). The parameter µ is related to the mean density 〈n_H〉 in the region as $\mu =\ln \langle n_H\rangle +\frac{{\sigma }^2}{2}$.

For a log-normal+power-law distribution of n_H, the PDF takes the form $${\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)=\{\begin{array}{ll}\frac{Q_1}{n_{\text{H}}}\exp \left[-\frac{{\left(\ln n_{\text{H}}-{\mu }_2\right)}^2}{2{\sigma }_2^2}\right],\hfill & \text{ if }{n}_{\text{H}}\le n_{\text{tr}}\hfill \\ Q_2{\left(\frac{n_{\text{H}}}{n_{\text{tr}}}\right)}^{\alpha },\hfill & \text{ if }{n}_{\text{tr}}<n_{\text{H}}\le n_{\text{cut}}\hfill \\ 0,\hfill & \text{ if }{n}_{\text{H}}>n_{\text{cut}},\hfill \end{array}$$(3)

where α < 0 is the slope of the power law and n_tr is the density at which the power-law tail begins. The parameters µ₂ and σ₂ characterise the location of the peak and the width of the log-normal part of the PDF. These are calculated, along with constants Q₁ and Q₂, for a given 〈n_H〉, α, and n_tr to match the mean density to 〈n_H〉 and ensure the continuity, differentiability, and normalisation of the PDF. Numerical simulations of self-gravitating molecular clouds (e.g. Kritsuk et al. 2011; Federrath & Klessen 2013) have shown that a power law with an index α > −1 provides a good fit to the density distribution in regions undergoing gravitational collapse. Girichidis et al. (2014) developed an analytical model for the time evolution of the slope of the power-law tail in a spherically collapsing cloud. They found that irrespective of the initial state of the cloud, after one free-fall time, the power-law tail has a universal slope of −0.54. Their value is in good agreement with the range of values observed for star-forming clouds in Kainulainen et al. (2009) and those found by Kritsuk et al. (2011) in simulations³. Therefore, we adopt α = −0.54. Furthermore, we set n_tr equal to 10 times the mean density 〈n_H〉 (Kritsuk et al. 2011) and, in order to prevent the integral of the PDF from diverging, we impose a cut-off of n_cut = 1000 〈n_H〉, above which the PDF is set to zero.

For both PDFs, we do not vary the parameters with the spatial resolution as long as it is larger than the typical scales of density fluctuations (i.e. the scales at which clumping takes place) in molecular clouds. For a detailed discussion on the variability of the clumping factor, we refer the interested readers to Micic et al. (2012) and Schäbe et al. (2020).

Figure 1 shows the distribution of sub-grid densities n_H in two sample simulation cells with different PDFs, each with 〈n_H〉 = 100 cm⁻³. The corresponding cumulative distribution functions (CDF) are plotted in the bottom panel. We can see that the log-normal+power-law distribution spans a broader range of densities as compared to the log-normal. This information is captured by the clumping factor C, which for the log-normal+power-law PDF used here is ~300 compared to C = 10 for the log-normal.

Effect of stellar feedback on the PDF. The internal structure and lifecycle of molecular clouds are affected by star formation and stellar feedback. While supernova (SN) feedback begins to act only a few million years (3–10) after star formation, pre-SN feedback in the form of stellar winds, photoionisation, and radiation pressure, particularly the latter in dense regions, can already start to act once stars form within molecular clouds.

Recent observations of molecular clouds in nearby galaxies (Hollyhead et al. 2015; Hannon et al. 2019; Kruijssen et al. 2019; Chevance et al. 2020) have found that these effects efficiently disperse the gas within molecular clouds.

A typical molecular-cloud region would cycle through episodes of star formation followed by stellar feedback. We use a log-normal+power-law PDF for gravitationally collapsing regions before star formation. At the onset of star formation in this region, we transition to a log-normal PDF to capture the combined effects of pre-SN and SN feedback. After a period of 40 Myr (Oey & Clarke 1997; also see Chevance et al. 2023 for a recent review of molecular cloud lifetimes), we switch back to the log-normal+power-law PDF. This 40 Myr timescale includes the molecular cloud dispersal time plus the time it takes for the assembly of the next cloud, but still prior to collapse. We note that the exact value of this timescale time will vary as a function of the local density and might be overestimated for some regions.

We note that our model does not explicitly capture the collapse of gas and star formation but only approximates their effect on chemistry by modifying the PDF from a log-normal to a log- normal+power-law form. Hence, these PDFs are a tool to mimic the effects of the ‘microscopic’ (i.e. unresolved) density structure on the ‘macroscopic’ (resolved) chemistry at different stages in the lifecycle of a molecular cloud.

Fig. 1 Probability distribution functions (PDFs) and the corresponding cumulative distribution functions (CDFs) used in this study. Panel a shows the mass-weighted PDFs in sample simulation cells as a function of the sub-grid density nH, where PM(nH) dnH denotes the fraction of the total cell mass present at sub-grid densities in the range [nH, nH + dnH]. The log-normal (Eq. (1)) and the log-normal+power-law (Eq. (3)) PDFs are shown in blue and red, respectively. The sample cells have a mean hydrogen density 〈nH〉 = 100 cm-3 (shown by the dotted black line). For the log-normal+power-law PDF, the transition density ntr and the cut-off density ncut are shown by the dashed and solid black lines, respectively. Panel b shows the corresponding CDFs, that is, CM (nH) denotes the fraction of the total cell mass present at sub-grid densities below nH . For the log-normal+power-law PDF shown here, the power-law tail encloses ~91% of the total cell mass.

2.2 The chemical network

We use a simplified version of the Nelson & Langer (1999, hereafter NL99) chemical network with some modifications (described in the following subsections) from the recent work of Gong et al. (2017, hereafter G17). Our simplifications reduce the number of chemical species and reactions that we follow to retain only the dominant formation and destruction channels for H₂, CO, C, and C⁺ under the physical conditions prevalent in molecular clouds (see Table A.1 in Appendix A). This simplifies the NL99 network for easier integration into cosmological simulations, allowing on-the-fly computation of chemical abundances without a significant increase in the computational overhead of these simulations.

Numerical simulations of molecular clouds (see e.g. Glover & Mac Low 2007b; Hu et al. 2021) show that at the densities relevant for H₂ formation in molecular clouds, the temperatures are mostly ≲200 K. At these temperatures, the contribution of ionised hydrogen (H⁺) to the total hydrogen mass is expected to be negligible. Therefore, we make a further simplifying assumption that all hydrogen is either atomic or molecular. Hence, in practice, we solve the system of rate equations for only three species, namely H₂, CO, and C⁺. The electron abundance follows from change conservation (i.e. f_e- = f_C+). We compute the abundances of H and C based on the conservation of hydrogen and carbon nuclei, respectively. We assume that the total gas-phase elemental abundances of carbon and oxygen are proportional to the gas metallicity Z, i.e. f_C.tot = 1.41 × 10⁻⁴ (Z/Z_⊙) and f_O,tot = 3.16 × 10⁻⁴ (Z/Z_⊙) (Savage & Sembach 1996).

2.2.1 H₂ chemistry

We consider the following formation and destruction channels for H₂:

1. H₂ formation on dust grains $$\text{H}+\text{H}+\text{ grain }\to {\text{H}}_2+\text{grain;}$$

2. radiative recombination of H⁺₃ with an electron $${\text{H}}_3^{+}+{\text{e}}^{-}\to {\text{H}}_2+\text{H;}$$

3. collisional reactions of H⁺₃ with atomic carbon and oxygen $${\text{H}}_3^{+}+\text{C}\to {\text{CH}}_{\text{x}}+{\text{H}}_2,$$ $${\text{H}}_3^{+}+\text{O}\to {\text{OH}}_{\text{x}}+{\text{H}}_2;$$

4. photodissociation of H₂ $${\text{H}}_2+\gamma \to 2\text{H};$$

5. ionisation of H₂ by cosmic rays $${\text{H}}_2+\text{CR}\to {\text{H}}_2^{+}+{\text{e}}^{-};$$

6. collisional reaction of H₂ with C⁺ $${\text{H}}_2+{\text{C}}^{+}+{\text{e}}^{-}\to {\text{CH}}_{\text{x}}+\text{H},$$ $${\text{H}}_2+{\text{C}}^{+}+{\text{e}}^{-}\to \text{C}+2\text{H}.$$

The rate of H₂ formation on dust grains depends on the dust abundance. We adopt a metallicity-dependent dust-to-gas mass ratio (DTG) based on observational measurements and theoretical predictions of the DTG as a function of the gas-phase metal- licity in galaxies at redshifts 0 < z ≲ 5 (Péroux & Howk 2020 and Popping & Péroux 2022). The DTG in our model is given as: $${\log }_{10}(\text{DTG})=1.3{\log }_{10}\left(\text{Z}/{\text{Z}}_{\odot }\right)-2.02,$$(4)

where Z/Ζ_⊙ is the gas-phase metallicity in solar units. We use Z_⊙ = 0.02 (Karakas 2010). The above non-linear dependence of the DTG on gas metallicity reflects a variable dust-to-metals (DTM) ratio (i.e. the fraction of metals locked up in dust grains). Based on a compilation of absorption-line studies of high- redshift objects, Popping & Péroux (2022) found that the DTG does not show significant evolution for 0 < z < 5. Hence, we further assume that the same relation holds at all redshifts.

For the collisional reaction between H₂ and C⁺, NL99 only considers the first outcome. However, it has been shown by Wakelam et al. (2010) that the reaction between C⁺ and H₂ gives CH_x + H only 70% of the times. In the remaining 30% of the cases, C + 2H are formed instead. This has important consequences for the relative abundances of C and CO. The CHx formed in the first outcome can react with an O atom to form CO or could photodissociate into C and H atoms, whereas, the second outcome acts as an additional formation channel for C in our network as well as in G17. The dissociation of H₂ is carried out by photons in two narrow bands of energies in the range 11.213.6 eV (λ = 912-1108 Å), called Lyman-Werner photons. We do not explicitly include three-body interactions in our network as these are inefficient at most ISM densities resolved in cosmological simulations. However, three-body reactions are the main mechanism for H₂ formation at high redshifts z ≳ 12 (see e.g. Christensen et al. 2012; Lenoble et al. 2024).

2.2.2 H₃⁺ chemistry

The ionisation of H₂ by cosmic rays produces the H₂⁺ ion that quickly reacts with an H atom to form H₃⁺. H₃⁺ can be destroyed by reactions with e⁻, C, and O. Because of its high reactivity (Oka 2006), we assume a local (i.e. at each sub-grid density) equilibrium between its formation and destruction channels (see Eqs. (A.4)–(A.5) in Appendix A). We note that NL99 also include the reaction of H₃⁺ with CO, which is excluded from our chemical network to limit the number of chemical species and reactions. Another difference with respect to NL99 is that, following G17, we consider two outcomes for the recombination of H₃⁺ with e⁻: a) H₂ + H b) 3 H. Of these, only the first one is included in NL99.

2.2.3 CO chemistry

The formation of CO proceeds via the reaction of CH_x (OH_x ) with O (C) which is formed by the reaction of H₃⁺ with C (O) (see reactions (3)–(6) in Table A.1). CHx can additionally be formed by the collisional reaction between C⁺ and H₂ . The destruction channels for CO are dissociation into C and O by UV photons in the 912–1100 Å band and cosmic rays. The rate of photodissociation drops off exponentially with the increasing column density of H₂, CO, and dust as a result of the shielding effect of these species on the impinging UV radiation deep inside molecular clouds (see Sect. 2.2.7 for the expression of the shielding functions).

2.2.4 Grain-assisted recombination of C⁺

Following Glover & Clark (2012, hereafter GC12) and G17, we include grain-assisted recombination of C⁺ in our network in addition to its radiative recombination. This is the main channel for C⁺ recombination at solar metallicity (G17), although its importance at sub-solar metallicities would depend on the relative amount of dust to gas in the ISM (the dust-to-gas ratio, see Sect. 2.2.1). Moreover, for several ions, including C⁺, the recombination rate on dust grains is often higher than the direct radiative recombination rate, especially in star-forming regions where dust, particularly in the form of polycyclic aromatic hydrocarbons (PAHs) is highly prevalent and frequently collides with ions leading to their recombinations. In a previous study, GC12 stressed that this reaction is particularly effective when the ratio of the UV field strength to the mean density is very small. As a result of including this additional destruction (formation) channel for C⁺ (C), we expect the C/C⁺ ratio predicted by our network (as well as in G17) to be significantly higher than predicted by NL99 (see Sect. 3.1).

2.2.5 Cosmic rays

Cosmic rays (CRs) with energies ≲0.1 GeV play a critical role in initiating ion-ion chemistry deep inside dense molecular- cloud regions that are well-shielded from UV radiation (see e.g. Padovani et al. 2009). These reactions become particularly relevant at high redshifts where the cosmic ray ionisation rate per H atom (CRIR, ζH) is expected to be higher than the canonical MW value of 3 × 10⁻¹⁷ s⁻¹ because of higher star formation rates (see e.g. Muller et al. 2016; Indriolo et al. 2018, for CRIR estimates at high redshifts). We include the ionisation of H₂, the ionisation of C, and the dissociation of CO by CRs in our network while NL99 only includes the first one.

Based on absorption studies of HD in H₂ -bearing damped Lyman-α systems, Kosenko et al. (2021) found that the CRIR scales quadratically with the UV field intensity relative to the Draine field⁴ (χ, Draine 1978). Therefore, as a default choice, we adopt the following relation between ζ_H and χ: $$\frac{{\zeta }_{\text{H}}}{{\zeta }_{\text{H},\text{MW}}}={\left(\frac{\chi }{{\chi }_{\text{MW}}}\right)}^2,$$(5)

but also consider alternative options in Sect. 4 and Appendix D. We use ζ_H,MW = 3 × 10⁻¹⁷s⁻¹ and χ_MW = 1.0 for MW. We also impose upper and a lower bounds on the CRIR of 3 ×10⁻¹⁴s⁻¹ and 10⁻¹⁸s⁻¹, respectively, to avoid unreasonable CRIRs. The effect of this upper limit is investigated in Sect. 4 and Appendix D. We note that the χ in Eq. (5) is measured in the FUV band (λ = 912–2070Å), while HYACINTH requires the UV flux in the Lyman-Werner band (λ = 912–1080 Å) in Habing units as an input. In the solar neighbourhood, the mean energy densities in the two bands are related as: U_LW / U_FUV ~ 1.1 (Parravano et al. 2003).

2.2.6 The temperature–density relation

As some chemical reactions in our network have temperaturedependent rate coefficients (see Table A.1), we need to associate a temperature with each sub-grid density in the PDF. For this, we use a metallicity-dependent temperature-density relation obtained from simulations of the ISM (Hu et al. 2021).

These simulations self-consistently include time-dependent H₂ chemistry and cooling, star formation, and feedback from photoionisation and supernovae. They adopt a linear relationship between the CRIR and the UV intensity (actually both quantities scale with the star-formation-rate surface density) which differs from our quadratic scaling and thus leads to a small inconsistency in our assumptions for dense, well-shielded gas, where CRs are an important heating source. In Sect. 4, we compare the chemical abundances obtained when using the ζ_Η − χ relation from Hu et al. (2021) against those obtained from Eq. (5).

2.2.7 Shielding functions

Dense, optically thick gas can shield itself against penetrating UV radiation because of the high column densities of H₂, CO, C, and dust. We account for this by modifying the reaction rates of photoionisation and photodissociation reactions by appropriate shielding functions. The (self-) shielding of H₂ is given by (Draine & Bertoldi 1996): $$\begin{array}{cc}{f}_{s,{\text{H}}_2}\left(N_{{\text{H}}_2}\right)=\hfill & \frac{0.965}{{\left(1+x/b_5\right)}^2}\hfill \\ \hfill & +\frac{0.035}{{(1+x)}^{0.5}}\exp \left[-8.5\times {10}^{-4}{(1+x)}^{0.5}\right],\hfill \end{array}$$(6)

where $x=\frac{N_{{\text{H}}_2}}{5\times {10}^{14}{\text{cm}}^{-2}},N_{{\text{H}}_2}$ is the column density of H₂ and b₅ is the velocity dispersion of gas in km s⁻¹. Following Sternberg et al. (2014) and T15, we use a constant b₅ = 2 throughout. However, as noted in G17, the H₂ fraction $f_{{\text{H}}_2}$ is insensitive to the value of b₅ for $f_{{\text{H}}_2}\gtrsim 0.1$.

The CO shielding function f_s,CO (N_CO, $f_{s,\text{CO}}\left(N_{\text{CO}},N_{{\text{H}}_2}\right)$ accounts for both CO self-shielding and the shielding of CO by H₂ and is calculated by interpolating over N_CO and $N_{{\text{H}}_2}$ from Table 5 in Visser et al. (2009).

The shielding function for C is given by (Tielens & Hollenbach 1985) $$f_{\text{s},\text{C}}\left(N_{\text{C}},N_{{\text{H}}_2}\right)=\exp \left(-{\tau }_{\text{C}}\right)f_{\text{s},\text{C}}\left({\text{H}}_2\right),$$(7)

where $$\begin{array}{cc}& {\tau }_{\text{C}}=1.6\times {10}^{-17}\frac{N_{\text{C}}}{{\text{cm}}^{-2}},\\ & f_{s,\text{C}}\left({\text{H}}_2\right)=\frac{\exp \left(-r_{{\text{H}}_2}\right)}{1+r_{{\text{H}}_2}},\\ & r_{{\text{H}}_2}=2.8\times {10}^{-22}\frac{N_{{\text{H}}_2}}{{\text{cm}}^{-2}},\end{array}$$(8)

where N_C and $N_{{\text{H}}_2}$ are the column densities of atomic carbon and H₂ , respectively.

In addition to self-shielding and H₂-shielding, all species are also shielded against UV radiation by dust grains; the relevant shielding function is given by $$f_{\text{dust}}=\exp \left(-\gamma A_V\right),$$(9)

where γ is different for each species and listed in Table A.1 in Appendix A. The visual extinction A_V is related to the total column density of hydrogen nuclei $N_{\text{H}}=N_{\text{HI}}+2N_{{\text{H}}_2}$ along the line of sight as $$A_V=\frac{N_{\text{H}}{Z}_{\text{d}}}{1.87\times {10}^{21}{\text{cm}}^{-2}}.$$(10)

In a simulation, recording the H₂ fraction for every sub-grid density within a grid cell at each timestep is computationally expensive and impractical. Therefore, we follow the approach of T15 for distributing the available H₂ to the different subgrid densities. We assume a sharp transition from atomic to molecular hydrogen at the sub-grid density $n_{\text{H}}=n_{\text{crit},{\text{H}}_2}$, that is $f_{{\text{H}}_2}(n_{\text{H}}<n_{\text{crit},{\text{H}}_2})=0$ and $f_{{\text{H}}_2}\left(n_{\text{H}}\ge n_{\text{crit},{\text{H}}_2}\right)=1$. This sharp transition has been observed in various numerical studies (Glover & Mac Low 2007a; Dobbs et al. 2008; Krumholz et al. 2008, 2009; Gnedin et al. 2009), and occurs at densities where H₂ becomes self-shielding (Dobbs et al. 2014). Similarly, we assume that carbon transitions from ionic to atomic form above n_crit,CI and becomes fully molecular at n_crit,CO . In a given region, these critical densities depend on the density of total hydrogen and the abundance of the different species involved in the transition (see Appendix C for calculation of $n_{\text{crit},{\text{H}}_2}$ for the HI – H₂ transition).

In practice, our sub-grid model requires six input parameters - the average density of hydrogen nuclei 〈n_H〉, the gas-phase metallicity Z, the UV flux in Lyman-Werner bands in Habing units G₀, the characteristic length scale Δx of a resolution element (e.g. cell size or smoothing length), the density PDF Ƥ_M, and the time Δt over which the chemical abundances are to be evolved. When embedded as a sub-grid model in a simulation, these parameters can be obtained directly from the simulation or calculated in post-processing (e.g. the UV flux). Depending on the value of Z, the temperature-density relation assigns a (sub-grid) temperature to each sub-grid density in the PDF. The chemical network then solves the rate equations for each sub-grid density and integrates over the PDF to obtain the average abundances of H₂ , CO, C, and C⁺ with respect to the total hydrogen within the region.

3 Comparison of chemical abundances

In this section, we compare different aspects of HYACINTH to previous approaches in the literature. First, we focus on the chemical network and contrast its predictions to those from the NL99 and G17 networks. Subsequently, we compare our full implementation of the sub-grid model (including the density PDF) with the output of high-resolution simulations of individual molecular clouds (Seifried et al. 2017, 2020; Glover & Mac Low 2011).

3.1 Comparison with NL99 and G17

In this section, we compare our chemical network predictions to Figs. 2b and 3b in G17. The setup involves a onedimensional semi-infinite slab with a uniform hydrogen density n_H = 1000 cm⁻³, a solar metallicity (Z = Z_⊙), a solar dust-togas ratio of 0.01, and fixed gas and dust temperatures of 20 K and 10 K, respectively. The slab is irradiated from one side by a UV field of strength 1 in Draine units (i.e. ,χ = 1). We test for two different values of the CRIR (same as in G17): 1 × 10⁻¹⁷ s⁻¹ H⁻¹ (left column of Fig. 2) and 2 × 10⁻¹⁶ s⁻¹ H⁻¹ (right column of Fig. 2). For this comparison, we assume the same elemental abundances for carbon and oxygen as in G17, that is, f_C,tot = 1.6 × 10⁻⁴ (Z/Z_⊙) and f_O,tot = 3.2 × 10⁻⁴ (Z/Z_⊙), throughout Sect. 3.1.

For a fair comparison with the G17 results, we calculate the shielding to the incident radiation field using their approximation for an isotropic radiation field (as described in section 2.3 of G17 and originally used by Wolfire et al. 2010). Briefly, it approximates an isotropic radiation field with a unidirectional field incident at an angle of 60° with the normal to the slab. Thus, for an incident radiation field of strength χ, the effective radiation field at a (perpendicular) depth L into the slab can be expressed as $${\chi }_{\text{eff}}=(\chi /2)f_{\text{shield }}(2L),$$

where ƒ_shield(L) represents the shielding function at a depth L into the slab and is different for each chemical species (see Sect. 2.2.7). For this uniform-density slab, the column density of hydrogen, N_H, at depth L can be written as N_H(L) = n_H L. The slab is divided into 1000 layers with N_H values spaced logarithmically in the range N_H = 10¹⁷ cm^–2–10²² cm^–2 or equivalently, visual extinction A_V (Eq. (10)) in the range 5.35 × 10^–5–5.35. For each layer, the NL99 and HYACINTH networks are evolved until equilibrium.

To put our results in context, we compare them to the output of more complex models. First, we consider the extended HYACINTH chemical network⁵ that includes the non-equilibrium treatment of two additional species, namely, He⁺ and HCO⁺. These species, particularly He⁺, serve as the main destruction agents for CO in dense, shielded regions.

Moreover, we present the output of the photon-dominated region (PDR) code used in G17 which tracks the abundances of 74 species accounting for 322 chemical reactions and includes a more sophisticated treatment of radiative transfer. This PDR code is derived from Tielens & Hollenbach (1985) and updated by Wolfire et al. (2010), Hollenbach et al. (2012), and Neufeld & Wolfire (2016).

Figure 2 shows the equilibrium abundances of the chemical species as a function of A_V. The ƒH₂ versus A_V from our network is in excellent agreement with that from NL99, G17 and the PDR code. In contrast, there are noticeable differences in the C⁺ → C and C → CO transitions in the different networks. For both values of the CRIR, at low A_V, the C⁺ abundance in HYACINTH (both standard and extended), G17 and the PDR code is slightly lower than NL99, accompanied by a higher C abundance for these networks with respect to NL99. In HYACINTH, this shift results primarily from the inclusion of grain-assisted recombination of C⁺ (see Sect. 2.2.4) and an additional outcome for the C⁺ + H₂ reaction (see Sect. 2.2.1). Together, these lead to almost an order of magnitude increase in the abundance of atomic carbon at A_V < 0.5 in HYACINTH compared to NL99. GC12 observed a similar trend when comparing their networks with and without this reaction. Previous studies have reported that the chemical networks employed in PDR codes tend to produce an elevated atomic carbon abundance with respect to the atomic carbon abundance measured in MW (interstellar) clouds (Sofia et al. 2004; Sheffer et al. 2008; Wolfire et al. 2008; Burgh et al. 2010). This has been long-standing problem for several chemical networks (see e.g. Burgh et al. 2010; Liszt 2011; Gong et al. 2017).

For ζ_H = 10^–17s^–1H^–1, the C → CO transition in HYACINTH (both standard and extended), G17 and the PDR code occurs at a slightly higher A_V than in NL99. Conversely, for ζ_H = 2 × 10^–16 s^–1 H^–1, this transition occurs at a slightly lower A_V in all other approaches compared to NL99. However, this does not have any practical implications for the modelling of CO chemistry in galaxy simulations as we expect the bulk of the CO mass in a molecular cloud or a galaxy to be present at high A_V (≳2).

A noticeable difference between our network and G17 is that at high A_V (≳1), all carbon in our network is in the form of CO, while G17 predicts ~3-10% of the carbon to be in atomic form. The PDR code predicts that ≲3% of the total carbon is in atomic form at A_V ≳ 1 for both values of the CRIR. In contrast, extended HYACINTH predicts an even higher abundance of atomic carbon than G17 at A_V ≳ 1 and consequently a lower CO abundance. Therefore, it is evident that the varying complexities of the different approaches result in significant differences in the atomic carbon abundance at A_V ≥ 1. However, since the bulk of the atomic carbon mass in all networks is present at intermediate A_V (0.1 ≲ A_V ≲ 1), the contribution of the atomic carbon in G17 and the PDR code at Av ≥ 1 to the total atomic carbon mass is ≲2% and is therefore, not significant.

Overall, both the standard and extended HYACINTH networks show a good agreement with NL99, G17, and the PDR code in the respective A_V range where each of the carbon species dominates. For H₂, the agreement is excellent for all A_V, showing that hydrogen chemistry is not sensitive to the exact treatment of carbon chemistry. There are, however, noticeable differences between standard and extended HYACINTH for the carbon-based species. For instance, there is a significant difference between the abundances of C and CO at A_V ≳ 0.5. Extended HYACINTH shows a smoother C → CO transition as compared to standard HYACINTH. This transition closely matches that in G17. Nevertheless, at A_V ≳ 2, the CO abundance in extended HYACINTH is only marginally different from that in the standard one – roughly 16% (20%) less for ζ_H = 10^–17 s^–1 H^–1 (2 × 10^–16 s^–1 H^–1). This shows that standard HYACINTH provides robust CO abundances with respect to extended HYACINTH, while requiring approximately 3.3 times less computational time. Therefore, in what follows, we only consider the standard HYACINTH network.

Fig. 2 Comparison of the chemical network in HYACINTH with NL99 and G17 networks. The abundances of H2, CO, C, and C+ as a function of the visual extinction AV in a semi-infinite plane-parallel slab are shown in different panels. The left column shows the abundances for a CRIR (ζH) of 10–17 s–1 H–1 while the right one has ζH = 2 × 10–16 s–1 H–1. The blue, turquoise, and red lines represent the results from HYACINTH, extended HYACINTH, and NL99 respectively. Here extended HYACINTH refers to the HYACINTH network with additional chemical reactions for He+ and HCO+ that are not part of standard HYACINTH (see text for more details). The dashed and dotted black lines show, respectively, the abundances from the chemical network and the PDR code in G17. The slab has a uniform hydrogen density nII = 1000 cm–3, solar metallicity and solar dust abundance, and is illuminated from one side by a UV field of strength χ = 1.

3.2 Comparison with molecular-cloud simulations

Now we shift our focus to assessing the performance of the chemical network in conjunction with the sub-grid density PDF and compare against high-resolution simulations of individual molecular clouds - the SILCC-Zoom simulations (Seifried et al. 2017, 2020) with solar metallicity (Z = Z_⊙, Sect. 3.2.1) and those from Glover & Mac Low (2011, hereafter GML11) with Z = 0.1 Z_⊙ (Sect. 3.2.2). We note that neither the SILCC-Zoom nor the GML11 runs account for star formation and stellar feedback.

3.2.1 SILCC-Zoom simulations

SILCC-Zoom are adaptive mesh refinement (AMR) simulations that follow the formation of two molecular clouds extracted from the parent SILCC simulations (Walch et al. 2015) with the zoom-in technique. They reach a minimum cell size of 0.06 pc in the densest regions. The original SILCC and SILCC-Zoom simulations were run with a chemical network based on Glover & Mac Low (2007a) and Glover & Mac Low (2007b) for hydrogen chemistry and Nelson & Langer (1997) for CO chemistry. However, the SILCC-Zoom simulations used in this comparison (Seifried et al. 2020) were performed with the NL99 network for CO chemistry.

Our goal is to compare the chemical evolution from HYACINTH (employed as a sub-grid model within a hydro simulation) with that from SILCC-Zoom. For this, we build a low-resolution copy of one of the SILCC-Zoom clouds (MC1-HD) by running a simulation with our modified version of the RAMSES code (Teyssier 2002) that uses HYACINTH for evolving the time-dependent chemistry. We consider a cubic box with a side length of 175 pc split in cells with a size of 25 pc (comparable to the spatial resolution we aim to achieve in our future cosmological simulations). We set the initial conditions (gas properties and chemical abundances) by coarse graining the SILCC-Zoom snapshot after a simulated time of 1 Myr, when all the mesh refinements have been performed, so as to exclude any possible numerical effects of variable spatial resolution on chemical abundances. We then follow the evolution of the system for 3 Myr, that is the full duration of the SILCC-Zoom runs. Beyond self-gravity, our simulation accounts for an external gravitational potential as described in sections 3.1 and 3.2 of Walch et al. (2015). The DTG is set to 0.01 and the elemental abundances of carbon and oxygen are set to ƒ_C,tot = 1.41 × 10^–4 and ƒ_O,tot = 3.16 × 10^–4. The cells are irradiated with an ISRF of strength G₀ = 1.7 in Habing units and a CRIR of ζ_H = 1.3 × 10^–17s^–1H^–1.

Before proceeding, two comments are in order. First, the lack of stellar feedback in the SILCC-Zoom simulations leads to the formation of very high-density clumps. Consequently, the density PDFs within the 25 pc cells vary over time and deviate significantly from the analytical forms we adopt in HYACINTH (see also Fig. 20 in Walch et al. 2015 and Fig. 3 in Buck et al. 2022). For instance, the clumping factor in these cells spans a wide range of values and is often much higher than the values of 10 and 300 that we adopt for the log-normal (hereafter LN) and log-normal+power-law (hereafter LN+PL) PDFs, respectively. Second, our simulation does not include radiative transfer of UV radiation. Consequently, the local ISRF is not attenuated due to the shielding from surrounding gas and dust as done in the SILCC-Zoom simulations. To quantify the impact of this effect, we perform a second simulation adopting a uniform effective visual extinction A_V = 0.5 in the entire simulation volume.

Figure 3 shows the time evolution of the masses of different chemical species in our original (AV = 0, solid lines) and UV-attenuated (A_V = 0.5, dashed lines) simulations compared to SILCC-Zoom (black diamonds). The shape of the density PDF determines how quickly the molecules are produced. The formation timescales for H₂ and CO are comparable to (for the LN+PL PDF) or longer than (for the LN PDF) the duration of the simulations which explains why we find that their masses are growing throughout the simulation. The LN PDF consistently underpre-dicts $M_{{\text{H}}_2}$ with respect to SILCC-Zoom by a factor of ~2 and produces little CO within the simulation time. This is consistent with the fact that SILCC-Zoom has a much more prominent high-density tail. On the other hand, the LN+PL PDF forms H₂ at a rate that is nearly twice as high as SILCC-Zoom and nicely matches the evolution of M_CO for the first 1.5 Myr, beyond which it saturates at about 50% and 75% of the final value in SILCC-Zoom in the runs with A_V = 0 and 0.5, respectively. To give an idea of the long-term evolution of the different species, we also compute their masses at chemical equilibrium for the final snapshot of the simulation (these are indicated with filled (A_V = 0) and empty (A_V = 0.5) circles on the right-hand side of the plots). These values never differ substantially from the simulation output at 3 Myr.

Overall, despite the differences in the exact treatment of various physical processes and the chemical network in SILCC-Zoom and our simulations, the predicted masses of the chemical species are in reasonable agreement, in particular when the LN+PL PDF is used. We stress that the goal here is not to achieve a perfect match since, as we mentioned, the sub-grid density PDF in MC1-HD is quite different from the analytical models used within HYACINTH, which are meant as averages over an ensemble of turbulent molecular clouds in the presence of stellar feedback. Finally, it is worth remembering that, when using HYACINTH as a sub-grid model in cosmological simulations, we switch between the LN and LN+PL PDFs to mimic the combined effects of self-gravity in a turbulent ISM and stellar feedback as explained in Sect. 2.1. The switch takes place on much longer timescales than those investigated in this section and would likely lead to somewhat different results.

Fig. 3 Time evolution of the total mass of different chemical species contained in a (175 pc)3 region encompassing the molecular cloud MC1-HD from the SILCC-Zoom simulations (black diamonds). The curves show the corresponding results from four low-resolution hydrodynamic simulations that use HYACINTH as a sub-grid model for the chemical evolution. The colour coding distinguishes runs obtained assuming either the LN PDF (blue) or the LN+PL PDF (red). The line style indicates the strength of the assumed uniform UV ISRF: G0 = 1.7 in Habing units with no attenuation (solid) or with an effective visual extinction of AV = 0.5 (dashed). The circles on the right-hand-side of the panels show the masses derived by assuming that the final snapshot of each simulation is in chemical equilibrium. Filled and empty symbols refer to the zero attenuation and AV = 0.5 runs, respectively.

Fig. 4 Comparison of the chemical abundances from HYACINTH against those from the Z = 0.1 Z⊙ runs from GML11 as a function of time. For 〈nH〉 = 100 cm−3, very little CO is formed when using the LN PDF (not visible in the bottom-right panel). The circles on the right-hand side denote the equilibrium abundances.

3.2.2 GML11 simulations

The GML11 simulations track the thermal and chemical evolution of magnetised and supersonically turbulent gas with typical conditions found in molecular clouds. The computational domain is a cube of side length L = 20 pc divided into a fixed number of cells. We present results for the runs with a grid-cell size of 0.156 pc which generate a closely LN density PDF well converged with the exception of the high-density tail. The GML11 simulations consider three different mean densities (100, 300, and 1000 cm⁻³) and use a chemical network (consisting of 218 reactions between 32 species) adapted to model hydrogen, carbon, and oxygen chemistry in molecular clouds. A detailed comparison of the performance in modelling CO chemistry of this more extended network and the simpler one given in NL99 is presented in GC12.

As the boxlength of the GML11 simulations is similar to the cell size we aim to achieve in our future cosmological runs, we model the entire computational volume of GML11 as a single domain in HYACINTH. The DTG is set to 0.01 (Z / Z_⊙) and the elemental abundances of carbon and oxygen are set to f_C,tot = 1.41 × 10⁻⁴(Z/Z_⊙) and f_O,tot = 3.16 × 10⁻⁴(Z/Z_⊙) (same as in GML11). All hydrogen is initially in atomic form while all carbon is in ionised form. The regions are irradiated with an interstellar radiation field (ISRF) of strength G₀ = 1.7 in Habing units and a cosmic ray ionisation rate of ζ_H = 10⁻¹⁷ s⁻¹ H⁻¹. We compute the chemical evolution for 5.7 Myr (the time interval covered by the GML11 runs) and, in Fig. 4, compare the H₂ and CO abundances with those from GML11⁶ (as presented in their Fig. 5). For 〈n_H〉 = 1000 cm⁻³, both LN and LN+PL give a reasonable agreement with GML11, while for the lower densities, the $f_{{\text{H}}_2}$ from LN closely matches the $f_{{\text{H}}_2}$ from GML11, but the one from LN+PL is consistently higher. A similar behaviour is observed for CO, although LN predicts roughly an order of magnitude lower ƒ_CO at 〈n_H〉 = 300 cm⁻³. In contrast, LN+PL gives a significantly higher ƒ_CO. At 〈n_H〉 = 100 cm⁻³, hardly any CO forms in GML11 and in our LN run, while the f_CO from LN+PL comprises a small fraction of the total carbon abundance (~5%).

As in Fig. 3, on the right hand side of each panel, we show the abundances that would be obtained at chemical equilibrium. For the LN+PL PDF, these are similar to the results obtained after 5.7 Myr, indicating that the molecules form with a characteristic timescale of a few million years. On the contrary, the equilibrium abundances obtained with the LN PDF are significantly higher than those obtained after 5.7 Myr, reflecting a much slower formation rate for the molecules due to the less prominent high-density tail in the PDF.

Overall, HYACINTH with the LN PDF gives H₂ abundances in excellent agreement with the GML11 simulations. The situation is more complex for CO, as the formation timescale of the molecules is very sensitive to the high-densitytail of the sub-grid PDF (as we noticed already in Fig. 3) and some fine tuning would be needed to get a good match between the different models. For 〈n_H〉 = 100 and 300 cm⁻³, our results with the LN PDF overestimate the formation time compared to GML11 while those with the LN+PL PDF underestimate it. At 〈n_H〉 = 1000 cm⁻³, CO is always produced at a slightly faster rate in HYACINTH than in the GML11 simulations. This simple test suggests that, at such densities and on such time intervals, HYACINTH overestimates the CO abundance by a factor of ~2 with respect to GML11. We note, however, that, in a less idealised set-up that accounts for star formation and stellar feedback, the high-density regions would be exposed to much more intense UV radiation than assumed in this test (see, e.g. Fig. 6) and this would surely affect the resulting CO abundance.

4 Sample application

In this section, we apply HYACINTH to a simulated galaxy to compare our model predictions with observations of H₂, CO, C, and C⁺ abundances in nearby and high-redshift galaxies. Although our primary goal is to integrate HYACINTH as a sub-grid component within simulations, here we post-process an existing simulation from T15 at ɀ ~ 2.5 and compute the equilibrium abundances. The galaxy in T15 was simulated with a modified version of the AMR code RAMSES (Teyssier 2002) including a sub-grid model for H₂ chemistry and a method to propagate the UV radiation in the Lyman-Werner (LW) bands to nearby cells. H₂ formation takes place on the surface of dust grains and H₂ is destroyed by LW photons. In contrast, our chemical network includes additional channels for H₂ formation and destruction (see Sect. 2.2.1 and Table A.1 for a complete list). Moreover, we adopt a different DTG than T15 (log₁₀ (DTG_T15) = log₁₀ (Z/Z_⊙) − 2.0; see Sect. 2.2.1 for our DTG).

An important caveat when applying HYACINTH in postprocessing involves the CRIR. In this scenario, the UV field χ from the simulation is used as an input to compute the equilibrium abundances, meaning that it is assumed that χ stays constant until equilibrium is reached. Because of the quadratic relationship between the CRIR and χ in HYACINTH (Sect. 2.2.5), cells in the simulation that have recently undergone star formation will have an extended period of high CRIR.

Thus, ζ_H ∝ χ² is not the ideal CRIR for a post-processing application⁷ .Therefore, we adopt a fixed CRIR of 3 × 10⁻¹⁷ s⁻¹ H⁻¹. Lastly, as applying the model in post-processing does not allow us to switch between the two PDFs, here we only use the log-normal (as in T15).

Figure 5 shows the column-density maps for the species obtained via post-processing as well as the total hydrogen and H₂ column densities from the simulation. Firstly, we see a remarkable agreement between the $N_{{\text{H}}_2}$ maps from the simulation and post-processing, indicating that the differences in our computation of H₂ chemistry compared to T15, namely the additional chemical reactions and the equilibrium treatment do not significantly impact the H₂ column density. This also indicates that the H₂ in the simulated galaxy is close to equilibrium. From the distribution of the three carbon species, we see that CO is the dominant carbon component in the highest $N_{{\text{H}}_2}$ (≳10²² cm⁻²) regions and is hardly present at $N_{{\text{H}}_2}\lesssim {10}^{21}{\text{cm}}^{-2}$, where N_CO drops below 10¹⁵ cm⁻². The distribution of atomic carbon (C) shows a great similarity to that of H₂ both in its extent and the location of peaks. C⁺ is found even more extensively throughout the galaxy, resembling the spread of the total hydrogen (H_tot). This is because carbon’s ionisation energy (11.6 eV) is slightly lower than that of hydrogen, allowing C⁺ to exist in all ISM phases. We further discuss this in the context of some recent observations in Sect. 4.2.4. Another noteworthy feature is that while the density peaks in CO and C coincide with those in H₂, the local fluctuations in C⁺ are milder. As H₂ assists in shielding the former two from UV radiation, their abundances are enhanced in high $N_{{\text{H}}_2}$ regions. On the other hand, C⁺ can exist both in the atomic and molecular phases of the ISM; therefore, we do not see a strong enhancement in $N_{{\text{C}}^{+}}$ with increasing $N_{{\text{H}}_2}$.

4.1 Comparison of H₂ abundance

A comparison of the post-processed H₂ fraction, $f_{{\text{H}}_2}=2\langle n_{{\text{H}}_2}\rangle /\langle n_{\text{H}}\rangle$, with that directly obtained from the simulation is presented in Fig. 6. We stress that because of the differences in T15 and our chemical network and because we solve for equilibrium, we do not expect to exactly reproduce the dynamically evolved H₂ abundance from the simulation in every grid cell, but only obtain values similar to those in the simulated galaxy. We also calculate $f_{{\text{H}}_2}$ in each cell using two analytical prescriptions, namely ‘KMT-EQ’ (Krumholz et al. 2009) and ‘KMT-UV’ (Krumholz 2013). In the KMT-EQ relation, $f_{{\text{H}}_2}$ is determined by the gas column density and metallicity, and is independent of the strength of the ISRF, G₀, by construction. The KMT-UV relation accounts for the effect of UV radiation on $f_{{\text{H}}_2}$ and is sensitive to the ratio G₀/〈n_H〉. In the left panel, we plot the median $f_{{\text{H}}_2}$ as a function of 〈n_H〉 for the different approaches. Firstly, we see that at 〈n_H〉 ≳ 100 cm⁻³, the median $f_{{\text{H}}_2}$ from post-processing with a uniform ζ_H (solid black line) agrees very well with the simulation and the KMT-EQ prediction. At lower densities, however, the three show some differences. By 〈n_H〉 = 10 cm⁻³, the KMT-EQ prediction gradually decreases to 0, while the simulation has a $f_{{\text{H}}_2}~0.2$. The post-processed $f_{{\text{H}}_2}$ shows a similar trend at low densities but is consistently higher than the simulation. In contrast with all other approaches, KMT-UV predicts a median $f_{{\text{H}}_2}=0$ at 〈n_H〉 ≲ 100 cm⁻³. This median is dominated by the cells that have a very high G₀/〈n_H〉 value and as such the KMT-UV prescription gives $f_{{\text{H}}_2}=0$ for all these cells. At 〈n_H〉 ≳ 100 cm⁻³, the $f_{{\text{H}}_2}$ increases gradually and gives similar results as the other methods.

Additionally, to demonstrate the impact of using a variable ζ_H in post-processing, we show the $f_{{\text{H}}_2}$ obtained for three different ζ_H − χ relations: (i) ζ_H ∝ χ², but with an upper limit of 3 × 10⁻¹⁴ s⁻¹ H⁻¹ on ζ_H (dashed green curve); (ii) ζ_H ∝ χ² without any upper limit (dotted green curve); (iii) the ζ_H − χ relation from Hu et al. (2021, solid blue curve). In all three cases, the post-processed $f_{{\text{H}}_2}$ exhibits a sharp drop at 〈n_H〉 ~ 70 cm⁻³ and shows a great resemblance to the KMT-UV prediction. We see that removing the upper bound of 3 × 10⁻¹⁴ s⁻¹ H⁻¹ on the CRIR in our ζ_H ∝ χ² relation (dotted green line) results in a decrease in the median $f_{{\text{H}}_2}$ by ~ 10% at 〈n_H〉 ≥ 100 cm⁻³.

In the middle and right panels of Fig. 6, we show, respectively, the dynamical (i.e. from the simulation) and the equilibrium $f_{{\text{H}}_2}$ (with a uniform ζ_H) for each grid cell colour-coded by G₀, the strength of the UV field in the LW band in the cell. We see that, at a given 〈n_H〉, the equilibrium $f_{{\text{H}}_2}$ predicted by our model is sensitive to the value of G₀ (similar to KMT-UV), with a higher G₀ resulting in a lower $f_{{\text{H}}_2}$. In contrast, the simulated $f_{{\text{H}}_2}$ does not show a clear trend with G₀, because the UV field varies throughout the formation history of H₂ in any given region. Overall, we see that the different approaches for H₂ chemistry result in varying predictions for $f_{{\text{H}}_2}$. In particular, the abundance computed dynamically in the simulation differs from the equilibrium calculations and shows a larger scatter at all densities. Similar plots when using a variable ζ_H are shown in Fig. D.2.

Fig. 5 Column density maps of different species in the simulated galaxy after post-processing. The total hydrogen and H2 from the simulation are shown in the first two panels. The remaining panels show the column density of H2, C+ , C, and CO obtained via post-processing. The H2 from post-processing (top row, third panel) is remarkably similar to the H2 from the simulation (top row, second panel). CO dominates the carbon budget in the highest $N_{{\text{H}}_2}$ regions and is essentially absent when $N_{{\text{H}}_2}\lesssim {10}^{21}$ cm−2. Atomic carbon traces the full extent of the H2 distribution, while C+ extends farther out to regions devoid of a significant amount of H2.

Fig. 6 Comparison of the H2 fraction $\left(f_{{\text{H}}_2}=2\langle n_{{\text{H}}_2}\rangle /\langle n_{\text{H}}\rangle \right)$ as a function of 〈nH〉 using different approaches. Left: the equilibrium $f_{{\text{H}}_2}$ from postprocessing the T15 galaxy with HYACINTH (black, green, and blue) compared with the $f_{{\text{H}}_2}$ from the simulation (red) and two analytical estimates − KMT-EQ (blue) and KMT-UV (turquoise). The lines show the median value in a given 〈nH〉 bin while the shaded area encloses the 16th to 84th percentiles. The solid black line denotes the median $f_{{\text{H}}_2}$ when using a uniform CRIR of ζH = 3 × 10−17 s−1 H−1 for HYACINTH. We also show the median $f_{{\text{H}}_2}$ obtained from post-processing when using the CRIR from Sect. 2.2.5 (i.e. ζH ∝ χ2) with (dashed green line) and without (dotted green line) the upper limit on the CRIR as well as (solid blue line) the ζH − χ relation from Hu et al. (2021, i.e., ζH ∝ χ, normalised to ζH,MW = 10−16 s−1 H−1). Middle: the $f_{{\text{H}}_2}$ from the T15 simulation colour-coded by the strength of the UV field in the LW bands in Habing units, G0 (from the simulation), in each grid cell. The red line is the same as in the left panel. Right: the post-processed $f_{{\text{H}}_2}$ (using a uniform ζH) for each grid cell within the simulated galaxy colour-coded by G0. The black line is the same as in the left panel.

4.2 Comparison with observations

4.2.1 The HI - H₂ transition

We compare the HI – H₂ transition in the post-processed galaxy (with a uniform ζ_H) with the observed one obtained from the absorption spectra of distant quasars and nearby stars in the MW and Large Magellanic Cloud (LMC). For this, we compute the column densities of H₂ and total hydrogen within a position-position-velocity data cube of 156pc × 156pc × 40 km s⁻¹⁸. For a total hydrogen column density $N_{\text{H}}=N_{\text{HI}}+2N_{{\text{H}}_2}$, the H₂ fraction can be defined as $f_{{\text{H}}_2}=2N_{{\text{H}}_2}/N_{\text{H}}$. The results of this comparison are shown in Fig. 7. In the left panel, we show the regions with MW-like conditions, that is 0.6 < Z/Z_⊙ < 1.4 and 0.6 < G₀/1.7 < 1.4 and in the right panel, those with LMC-like conditions, that is 0.18 < Z/Z_⊙ < 0.42 and 6 < G₀/1.7 <14. The observed data are taken from the Copernicus survey (Savage et al. 1977), the HERACLES survey of 30 nearby galaxies (Schruba et al. 2011), and the compilation of several observations of OB stars in the Galactic disc with the Far Ultraviolet Spectrographic Explorer (FUSE) by Shull et al. (2021). The LMC data are from Tumlinson et al. (2002) using a FUSE survey. We see that the HI – H₂ transition from post-processing agrees very well with the observed relation, particularly for MW-like conditions (left panel). As evident from both the post-processed and observed data, the transition is sensitive to the metallicity and the strength of the ISRF: it shifts to higher N_H values for LMC-like conditions as compared to MW-like conditions, since the former has 10 times stronger ISRF than the MW but only about a third of the metals in the MW.

4.2.2 The relationship between ${\text{N}}_{{\text{H}}_2}$ and N_CO

Figure 8 shows the relation between the ratio of CO to H₂ column densities $\left(N_{\text{CO}}/N_{{\text{H}}_2}\right)$ and the H₂ column density $N_{{\text{H}}_2}$ in the post-processed galaxy for three different CRIRs. For comparison, we include column densities obtained from UV absorption measurements along sight lines towards diffuse and (optically) dark molecular clouds in the MW from Sheffer et al. (2008). The $N_{\text{CO}}/N_{{\text{H}}_2}$ values for a compilation of dark-cloud observations by Federman et al. (1990) are also shown. In these observations, dark clouds are defined as those with visual extinction A_V ≳ 5, while those with A_V ≲ 5 are defined as diffuse. As these are MW observations, we only plot the post-processed regions with 0.6 < Z/Z_⊙ < 1.4. Moreover, Fig. 7 of Sheffer et al. (2008) shows that the ISRF values (expressed as χ = G₀/1.7) for these observations range from ~0.5 to ~10. Therefore, for a fair comparison, we only select the post-processed regions with G₀/1.7 in the range [0.35, 13], allowing for a scatter of 30%.

At H₂ column densities above 10²¹ cm⁻², our post-processed points (for all CRIRs) span the same region in the $N_{\text{CO}}/N_{{\text{H}}_2}$ versus $N_{{\text{H}}_2}$ plane as Federman et al. (1990), albeit with significantly less scatter. At lower $N_{{\text{H}}_2}~{10}^{21}$ cm⁻², the $N_{\text{CO}}/N_{{\text{H}}_2}$ ratio shows a sharp drop when using a uniform ζ_H = ζ_{H MW}, in contrast to the gradual decline seen in the observed data at $N_{{\text{H}}_2}~5\times {10}^{20}-{10}^{21}$ cm⁻². Conversely, the $N_{\text{CO}}/N_{{\text{H}}_2}$ ratio shows a gradual decline similar to the observed data when using a variable ζ_H (CRIR). Nevertheless, both our data and observations exhibit a large scatter in the $N_{\text{CO}}/N_{{\text{H}}_2}$ values at $N_{{\text{H}}_2}\in \left[{10}^{20},{10}^{21}\right]{\text{cm}}^{-2}$. Moreover our values are in a good agreement with the observed values at these H₂ column densities. Such large variations in $N_{\text{CO}}/N_{{\text{H}}_2}$ values were also reported previously by Smith et al. (2014) in their simulation of a MW-like galaxy. At $N_{{\text{H}}_2}\lesssim {10}^{20}$, our data is sparse and some of our sight lines have a factor ~2 higher N_CO compared to those reported by Crenny & Federman (2004) and Sheffer et al. (2008). Nevertheless, the majority of our data remains consistent with the observations shown here. Finally, it is worth noting that the regions within our post-processed galaxy with MW-like metallicity predominantly resemble dark clouds while only a few inhabit the region where N_CO < 10¹⁶cm⁻² and $N_{{\text{H}}_2}<{10}^{21} {\text{cm}}^{-2}$.

4.2.3 The abundance of atomic carbon

The fine structure lines of atomic carbon CI (corresponding to the ³P₂ – ³P₁ and ³P₁ − ³P₀ transitions) are often used to infer the H₂ masses of galaxies – both in the local Universe as well as at high redshifts (e.g. Gerin & Phillips 2000; Ikeda et al. 2002; Weiß et al. 2003, 2005; Walter et al. 2011; Valentino et al. 2018; Henríquez-Brocal et al. 2022). This method requires assuming an atomic carbon abundance relative to H₂. Several observations have tried to measure this abundance using other independent estimates of the H₂ mass of a galaxy such as dust emission in the infrared or rotational lines of CO. Here we compare the atomic carbon abundance relative to H₂ for our post-processed galaxy with those found in observations. Our galaxy (post-processed with a uniform ζ_H) has M_CI = 6.29 × 10⁶ M_⊙ and $M_{{\text{H}}_2}=4.70\times {10}^{10}$ M_⊙ which implies a galaxy-integrated neutral carbon abundance relative to H₂, $X_{\text{CI}}=\frac{M_{\text{CI}}}{6M_{{\text{H}}_2}}$ of 2.23 × 10⁻⁵ or equivalently log₁₀ X_CI = −4.65. The M_CI of our galaxy is similar to that obtained by Tomassetti et al. (2014) using a different approach. They use a X_CI of 3 × 10⁻⁵ from Alaghband-Zadeh et al. (2013) and obtained M_CI = 7.92 × 10⁶ M_⊙ for their simulated galaxy at z = 2.

Walter et al. (2011) carried out a survey of CI emission in z > 2 sub-millimetre galaxies (SMGs) and quasar host galaxies using the IRAM Plateau de Bure interferometer and the IRAM 30 m telescope. They obtain a mean log₁₀ X_CI of $-{4.08}_{-0.23}^{+0.16}$ for their sample of 10 galaxies. Valentino et al. (2018) carried out a survey of [C I] (³P₁ −³ P₀) in far-infrared-selected main-sequence (MS) galaxies in the COSMOS field at z ~ 1.2 with the Atacama Large Millimeter Array (ALMA). They found a mean log₁₀ X_CI of −4.7 + 0.1 and −4.8 + 0.2, respectively, for $M_{{\text{H}}_2}$ estimates based on dust and CO measurements for their sample of 12 galaxies. For SMGs at z > 2.5, they obtain a value of −4.3 + 0.2 and −4.2 + 0.1 for $M_{{\text{H}}_2,\text{ dust }}$ and $M_{{\text{H}}_2,\text{CO}}$, respectively. For a sample of 6 MS galaxies at z = 1–3 observed as part of the ASPECS, Boogaard et al. (2020) estimated X_CI = (1.9 ± 0.5) × 10⁻⁵ (i.e. ${\log }_{10}{X}_{{\text{C}}_1}=-{4.72}_{-0.13}^{+0.10}$). Using a sample of 21 lensed starburst galaxies with [C I] detection at z ~ 1.3–3.5, Harrington et al. (2021) found an X_CI of (6.82 ± 3.04) × 10⁻⁵ (i.e. ${\log }_{10}{X}_{{\text{C}}_1}=-{4.17}_{-0.26}^{+0.16}$). A slightly lower value of ~ −4.5 is often used in observations of galaxies at z ≳ 2 (e.g. Weiß et al. 2003). Thus, broadly speaking, the X_CI in our post-processed galaxy is consistent with those obtained for star-forming galaxies at z ≳ 2. In particular, our log₁₀ X_CI = −4.65 is similar to those in MS galaxies, but lower than the values reported for starbursts and SMGs.

Fig. 7 H2 fraction as a function of the total column density of hydrogen $N_{\text{H}}=N_{\text{HI}}+2N_{{\text{H}}_2}$, where $f_{{\text{H}}_2}=2N_{{\text{H}}_2}/N_{\text{H}}$. Left: the HYACΓNTH-post-processed regions with MW-like conditions are shown in black. The observed data are from Savage et al. (1977) (blue open squares), Shull et al. (2021) (orange open squares), Schruba et al. (2011) (magenta dots). For clarity, here we only show a contour enclosing all data points for Schruba et al. (2011). Right: the HYACINTH-post-processed regions with LMC-like conditions (black) compared to the observed data from Tumlinson et al. (2002) (green open squares). For HYACINTH data, we use a uniform ζH = ζH,MW and the column densities are calculated in a position-position-velocity cube of 156 pc 156 pc × 40 km s−1, where 156 pc is the spatial resolution of the AMR grid in the simulation at the redshift of post-processing and 40 km s−1 is an upper limit on the observed velocity dispersion in the various absorption measurements used in this comparison.

Fig. 8 Ratio of the CO-to-H2 column density $N_{\text{CO}}/N_{{\text{H}}_2}$ versus $N_{{\text{H}}_2}$ for the post-processed galaxy, compared with observations of (optically) diffuse (Rachford et al. 2002; Crenny & Federman 2004; Sheffer et al. 2008) and dark (Federman et al. 1990) molecular clouds in the MW. The black, green, and blue points are from the galaxy post-processed with HYACINTH but using different CRIRs (as indicated in the legend). In each case, the points are selected to have a gas-metallicity Z ∈ [0.7, 1.3] and χ ∈ [0.35, 1.3] to match the respective values in the observed data. The column densities are calculated in a position-position-velocity cube of 156pcl56pc× 5 km s−1, to match the velocity dispersion in the absorption measurements.

4.2.4 The extended C⁺ distribution

The [C II] fine-structure line at 158 µm is one of the brightest emission lines in star-forming galaxies and an important coolant of the ISM. The morphology, extent, and kinematics of this line give crucial insights into the different physical processes in galaxies. Several recent [C II] observations with ALMA (e.g. Nesvadba et al. 2016; Carniani et al. 2018; Fujimoto et al. 2019, 2020; Herrera-Camus et al. 2021; Lambert et al. 2023) have detected the presence of an extended C+ reservoir in galaxies out to z ~ 7. These observations find that the effective radius of the [C II]-emitting region is ~ 2–3 times larger than that of the rest-frame UV-emitting region. Possible explanations for this extended emission include metal-rich outflows, mergers, connection to Lyα halo, and contribution from satellite galaxies, without any clear consensus so far regarding its source (Fujimoto et al. 2019, 2020). What causes this extended [C II] emission remains an open question in the field of galaxy formation.

Figure 9 shows the surface density profiles (normalised to one at the centre) for all gas, H₂, CO, C⁺, young stars (age ≤ 10 Myr) and all stars in our post-processed galaxy (with a uniform ζ_H) observed face on. The H₂ and CO profiles show a steep decline such that both ${\Sigma }_{{\text{H}}_2}$ and σ_CO decrease by more than an order of magnitude from the centre out to 2 kpc. The stellar surface density profile shows a similar although less steep decline. In contrast, the C⁺ profile exhibits a gradual decline, maintaining a substantially high surface density at distances extending up to ~ 10 kpc, beyond which it shows a steep decline. In this regard, our (normalised) C⁺ profile is similar to the stacked one observed in Fujimoto et al. (2019) and used to infer the presence of a C⁺ ‘halo’ (which might be an inappropriate name as, in our simulated galaxy, the extended profile arises from the gaseous disc and not from a spheroidal distribution). Young stars also exhibit a relatively slower decline albeit with strong fluctuations. This reflects that UV radiation is essential for ionising carbon. Moreover, this is consistent with the observed correlation between the [C II] line luminosity and the star formation rate (De Looze et al. 2011, 2014; Herrera-Camus et al. 2015; Schaerer et al. 2020). Also in Fig. 5, we saw that the C⁺ distribution strongly mirrors the total hydrogen distribution which is dominated by atomic and ionised hydrogen at larger distances from the centre. Since we apply HYACINTH in post-processing, all C⁺ that is farther than the other components is formed in situ and not transported by any dynamical processes such as outflows.

Fig. 9 Surface density profiles of gas (black), H2 (blue), CO (magenta), C+ (red), stars (orange), and young stars (teal, with ages less than 10 Myr) within the HYACINTH-post-processed galaxy. For each component, the profiles are normalised by their central value σ0. For each radial bin, we calculate the surface density from the face-on projection from a cylinder with height of ±5 kpc. Among all components, C+ and young stars show the slowest decline as we move away from the centre of the galaxy, while H2 and CO are more concentrated (see also Fig. 5). The dashed grey line denotes the size of the galaxy as defined in T15.

5 Discussion

It is well known that H₂ formation can be greatly enhanced by the clumpiness of the ISM on scales typically below the resolution scale of present-day cosmological simulations. Our model accounts for this unresolved density structure by assuming a sub-grid density PDF that varies with the state of star formation in each gas element. This variable PDF is designed to mimic the effect of pre-supernova feedback on the density structure of molecular clouds. In this regard, our approach is alternative to Lupi et al. (2018) who assume a log-normal PDF with the clumping factor related to the local Mach number. While Lupi et al. (2018) always assume a log-normal PDF with a continuously varying clumping factor, our approach has a binary nature, where the PDF switches between two distinct shapes – log-normal and log-normal+power-law, each with a constant clumping factor.

In Sect. 3.2, we demonstrated that the chemical abundances from HYACINTH are in a reasonable agreement with high-resolution molecular-cloud simulations. This shows that HYACINTH works reasonably well for modelling molecular hydrogen and carbon chemistry. However, our goal is not to replace extensive chemical networks that would be more accurate for modelling chemistry in ISM-scale simulations, but rather the model is tailored to follow hydrogen and carbon chemistry in large cosmological simulations that have a resolution of 20–200 pc.

In Sect. 4, we have shown a post-processing application for a galaxy simulation (at z ~ 2.5) with a spatial resolution of 156 pc. The resulting equilibrium H₂ abundance are in agreement with the commonly adopted KMT-EQ relation and the H₂ column density map shows a great similarity to the one obtained directly from the simulation (accounting for the non-equilibrium H₂ abundance). This shows that our approach gives reasonable results on larger scales (~ 100–200 pc) that barely resolve individual molecular clouds.

An important caveat about using HYACINTH is that the formation of H₂ in HYACINTH relies on the presence of dust, similar to most other approaches in the literature (e.g. KMT). Therefore, it cannot be applied to gas with arbitrarily low metallicities (Z ≲ 10⁻³ Z_⊙). As a result it cannot be applied to pristine gas such as in simulations of population III stars and early galaxy formation (z ≳ 15, see e.g. Hirano et al. 2015; Lenoble et al. 2024). Accurately modelling hydrogen chemistry in such (proto) galaxies requires accounting for three-body reactions and more importantly a high spatial resolution to resolve the densities where these reactions are efficient (n_H > 10⁸ cm⁻³). One workaround for cosmological simulations that do not resolve these early objects is to impose a metallicity floor that mimics the metal enrichment by the first stars (i.e. the Population III stars; see e.g. Kuhlen et al. 2012, 2013; Tomassetti et al. 2015; Pallottini et al. 2022).

6 Summary

We present a new sub-grid model, HYACINTH, that can be embedded into cosmological simulations for calculating the non-equilibrium abundances of H₂ and its carbon-based tracers, namely CO, C, and C⁺. HYACINTH comprises a variable sub-grid density PDF to capture the unresolved density structure in simulations and a simplified chemical network for hydrogen and carbon chemistry. These simplifications were introduced to make the network highly efficient for use in large-scale simulations. Additionally, a metallicity-dependent temperature-density relation, based on high-resolution simulations of the star-forming ISM, was used to assign a temperature to each sub-grid density.

We compared HYACINTH against more sophisticated approaches in the literature for modelling hydrogen and carbon chemistry including two extensive chemical networks (NL99 and G17) and a PDR code, using a one-dimensional semi-infinite slab setup (Fig. 2). HYACINTH reproduced the $f_{{\text{H}}_2}-A_V$ relation from these methods highlighting that H₂ chemistry is insensitive to the exact treatment of carbon chemistry. Moreover, despite its simplicity and size, HYACINTH captured the C⁺ → C and C → CO transitions very well, as predicted by more complex approaches (see Sect. 3.1 for a detailed discussion).

We compared the predictions from HYACINTH with high-resolution molecular-cloud simulations – the SILCC-Zoom simulations at Z = Z_⊙ and the GML11 simulations at Z = 0.1 Z_⊙. We find reasonable agreement between the chemical abundances from HYACINTH and those from the simulations. We further find that our log-normal+power-law PDF shows better agreement with the SILCC-Zoom simulations while the log-normal PDF performs better for the GML11 runs, in alignment with the different density structures shown by the two simulations.

Finally, we applied HYACINTH to a simulated galaxy at ɀ ~ 2.5 from T15 in post-processing and compared it directly with observations. For regions in the post-processed galaxy with MW-like conditions, the H I-to-H₂ transition in the $f_{{\text{H}}_2}-N_{\text{H}}$ plane aligns very well with that from observations. The observed $f_{{\text{H}}_2}$ and N_H were obtained from measurements of the absorption spectra of quasars and nearby stars in sightlines towards MW molecular clouds. The same is also true for LMC-like regions (Fig. 7). Additionally, the values of $N_{\text{CO}}/N_{{\text{H}}_2}$ versus $N_{{\text{H}}_2}$ in MW-like regions are consistent with observations of MW molecular clouds (Fig. 8). It is worth noting that most of our post-processed regions resemble the observed (optically) dark molecular clouds in the $N_{\text{CO}}/N_{{\text{H}}_2}-N_{{\text{H}}_2}$ plane. Furthermore, the relative abundance of atomic carbon to molecular hydrogen (X_{C I}) in our post-processed galaxy consistently matches the X_{C I} values found in star-forming galaxies at redshifts ɀ ≳ 1. Based on surface density profiles of the different baryonic components within the post-processed galaxy, we find an excess of C⁺ at large distances from the galaxy centre, similar to those found in observations (e.g. Fujimoto et al. 2019).

In a forthcoming paper (paper II; Khatri et al., in prep.), we will present a suite of cosmological simulations using HYACINTH to model the non-equilibrium abundances of H₂ and its tracers in high-redshift (ɀ ≳ 2) galaxies. This will open up the possibility of addressing fundamental questions such as what the contribution is of low-mass galaxies to the global H₂ budget at high ɀ, what regulates the molecular gas fraction in galaxies, and for what physical conditions, environments, and galaxies are CO, C, and C⁺ all reliable tracers of molecular gas.

Acknowledgements

The authors thank an anonymous reviewer and Pascal Oesch for suggestions and comments that helped improving the paper. They gratefully acknowledge the Collaborative Research Center 1601 (SFB 1601 sub-projects C5, B1, B4) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 500700252. Early stages of this work were carried out within the Collaborative Research Centre 956 (SFB 956 sub-projects C4, C5, C6), funded by the DFG – 184018867. PK is a part of the International Max Planck Research School in Astronomy and Astrophysics, the Bonn Cologne Graduate School, and a guest at the Max Planck Institute for Radio Astronomy in Bonn. DS acknowledges funding from the programme “Profilbildung 2020”, an initiative of the Ministry of Culture and Science of the State of North Rhine-Westphalia. CP is grateful to SISSA, the University of Trieste, and IFPU, where part of this work was carried out, for hospitality and support.

Appendix A Reactions in the chemical network

Table A.1 gives a list of all chemical reactions in HYACINTH. For each of these, we use updated reaction rates from UMIST (McElroy et al. 2013) and KIDA (Wakelam et al. 2012) databases.

The fractional abundance f_X of a species X relative to hydrogen is defined as f_X = n_X/n_H, except for H₂ and ${\text{H}}_3^{+}$, where $f_{{\text{H}}_2}=2n_{{\text{H}}_2}/n_{\text{H}}$ and $f_{{\text{H}}_3^{+}}=3n_{{\text{H}}_3^{+}}/n_{\text{H}}$. For H₂, CO, and C⁺, the differential equations for each sub-grid density n_H can be written as⁹ $$\begin{array}{l}\frac{\text{d}{f}_{{\text{H}}_2}}{ \text{d}t}=2k_{{\text{H}}_2,\text{gr}}{f}_{\text{HI}}{n}_{\text{H}}-k_{\gamma ,{\text{H}}_2}{f}_{{\text{H}}_2}-k_{\text{c}.\text{r}.,{\text{H}}_2}{f}_{{\text{H}}_2}\hfill \\ -n_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}+\frac{2}{3}{n}_{\text{H}}{k}_{{\text{e}}^{-},{\text{H}}_3^{+}}{f}_{{\text{e}}^{-}}{f}_{{\text{H}}_3^{+}};\hfill \end{array}$$(A.1) $$\begin{array}{l}\frac{\text{d}{f}_{\text{CO}}}{\text{d}t}=\frac{1}{3}{n}_{\text{H}}{k}_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{O}}{}_{\text{I}}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{OH}}_x}+\frac{1}{3}{n}_{\text{H}}{k}_{\text{C}1,{\text{H}}_3^{+}}{f}_{\text{CI}}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{CH}}_x}\hfill \\ +\frac{1}{2}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}-k_{\gamma ,\text{CO}}{f}_{\text{CO}};\hfill \end{array}$$(A.2) $$\begin{array}{l}\frac{\text{d}{f}_{{\text{C}}^{+}}}{\text{d}t}=-n_{\text{H}}{k}_{{\text{C}}^{+},{\text{e}}^{-}}{f}_{{\text{C}}^{+}}{f}_{{\text{e}}^{-}}-n_{\text{H}}{k}_{{\text{C}}^{+},{\text{e}}^{-},\text{gr}}{f}_{{\text{C}}^{+}}{f}_{{\text{e}}^{-}}\hfill \\ -\frac{1}{2}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}+k_{\gamma ,\text{CI}}{f}_{\text{CI}}.\hfill \end{array}$$(A.3)

The reaction rate coefficient k for each reaction is listed in Table A.1. In the above equations, ${\alpha }_{{\text{CH}}_{\text{x}}}\left({\alpha }_{{\text{OH}}_{\text{x}}}\right)$ is the branching ratio (see Sect. B.1) for the formation of CH_x (OH_x) when ${\text{H}}_3^{+}$ reacts with atomic carbon (oxygen). The species O₂, H₂O, OH, OH⁺, H₂O⁺ and H₃O⁺ are collectively referred to as OH_x in the chemical network (as originally done in NL99 based on the argument that these species have similar chemical reactions and reaction rates). Similarly, CH, CH₂, CH⁺, and ${\text{CH}}_2^{+}$ are collectively referred to as CH_x. For ${\text{H}}_3^{+}$, we assume a local equilibrium (i.e. at each sub-grid density) between the formation and destruction pathways such that: $$\begin{array}{l}\frac{\text{d}{f}_{{\text{H}}_3^{+}}}{\text{d}t}=\frac{3}{2}{k}_{\text{C},\text{r},{\text{H}}_2}{f}_{{\text{H}}_2}-n_{\text{H}}{k}_{\text{CI},{\text{H}}_3^{+}}{f}_{\text{CI}}{f}_{{\text{H}}_3^{+}}\hfill \\ -n_{\text{H}}{k}_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{OI}}{f}_{{\text{H}}_3^{+}}-n_{\text{H}}{k}_{{\text{e}}^{-},{\text{H}}_3^{+}}{f}_{{\text{e}}^{-}}{f}_{{\text{H}}_3^{+}}=0,\hfill \end{array}$$(A.4)

which implies $$f_{{\text{H}}_3^{+}}=\frac{3}{2n_{\text{H}}}\frac{k_{\text{c},\text{r},{\text{H}}_2}{f}_{{\text{H}}_2}}{k_{\text{CI},{\text{H}}_3^{+}}{f}_{{\text{C\hspace{0.17em}I}}_{}}+k_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{OI}}+k_{{\text{e}}^{-},{\text{H}}_3^{+}}{f}_{{\text{e}}^{-}}}.$$(A.5)

Similarly, an equilibrium abundance of OH_x and CH_x implies $$\begin{array}{l}\frac{\text{d}{f}_{{\text{OH}}_{\text{x}}}}{\text{d}t}=\frac{1}{3}{k}_{{\text{OI,H}}_3^{+}}{f}_{\text{OI}}{f}_{{\text{H}}_3^{+}}{n}_{\text{H}}\hfill \\ -k_{\text{CI},{\text{OH}}_x}{f}_{\text{CI}}{f}_{{\text{OH}}_{\text{x}}}{n}_{\text{H}}-k_{\gamma ,{\text{OH}}_{\text{x}}}{f}_{{\text{OH}}_{\text{x}}}=0;\hfill \end{array}$$(A.6) $$\begin{array}{l}\frac{\text{d}{f}_{{\text{CH}}_{\text{x}}}}{\text{d}t}=\frac{1}{3}{k}_{{\text{CI,H}}_3^{+}}{f}_{\text{CI}}{f}_{{\text{H}}_3^{+}}{n}_{\text{H}}+\frac{1}{2}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{n}_{\text{H}}\\ -k_{\text{OI},{\text{CH}}_{\text{x}}}{f}_{\text{OI}}{f}_{{\text{CH}}_{\text{x}}}{n}_{\text{H}}-k_{\gamma ,{\text{CH}}_{\text{x}}}{f}_{{\text{CH}}_{\text{x}}}=0.\end{array}$$(A.7)

These give $$f_{{\text{OH}}_{\text{x}}}=\frac{\frac{1}{3}{k}_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{OI}}{f}_{{\text{H}}_3^{+}}}{k_{\text{CI},{\text{OH}}_{\text{x}}}{f}_{\text{CI}}+k_{\gamma ,{\text{OH}}_{\text{x}}}/n_{\text{H}}};$$(A.8) $$f_{{\text{CH}}_{\text{x}}}=\frac{\frac{1}{3}{k}_{\text{CI},{\text{H}}_3^{+}}{f}_{\text{CI}}{f}_{{\text{H}}_3^{+}}+\frac{1}{2}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}}{k_{{\text{OI}}_{,{\text{CH}}_{\text{x}}}}{f}_{\text{OI}}+k_{\gamma ,{\text{CH}}_{\text{x}}}/n_{\text{H}}}.$$(A.9)

Appendix B Cell-averaged rate equations

Obtaining the differential equation for the evolution of each species within a given region (e.g. a cell in a cosmological simulation) requires calculating a PDF-weighted integral of each term on the right-hand side in Eqs. (A.1)–(A.3). This leads to the following set of differential equations for the cell-averaged abundances: $$\begin{array}{l}\frac{\text{d}\langle f_{{\text{H}}_2}\rangle }{\text{d}t}=\langle 2k_{{\text{H}}_2,\text{gr}}{f}_{\text{HI}}{n}_{\text{H}}\rangle -\langle k_{\gamma ,{\text{H}}_2}{f}_{{\text{H}}_2}\rangle -\langle k_{\text{c}.\text{r}.{\text{H}}_2}{f}_{{\text{H}}_2}\rangle \\ -\langle n_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}\rangle +\langle \frac{2}{3}{n}_{\text{H}}{k}_{{\text{e}}^{-},{\text{H}}_3^{+}}{f}_{{\text{e}}^{-}}{f}_{{\text{H}}_3^{+}}\rangle \\ \underset{A_1}{\underbrace{{\displaystyle {\int }_0^{n_{\text{crit, },{\text{H}}_2}}2}{k}_{{\text{H}}_2,\text{gr}}{n}_{\text{H}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}-\underset{A_2}{\underbrace{{\displaystyle {\int }_{n_{\text{crit, },{\text{H}}_2}}^{\infty }{k}_{\gamma ,{\text{H}}_2}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ -\underset{A_3}{\underbrace{{\displaystyle {\int }_{n_{\text{crit, },{\text{H}}_2}}^{\infty }{k}_{\text{c}.\text{r},{\text{H}}_2}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}-\underset{A_4}{\underbrace{{\displaystyle {\int }_{n_{{\text{crit, H}}_2}}^{n_{\text{crit,CI }}}{n}_{\text{H}}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{\text{C},\text{ tot }}{\alpha }_{{\text{CH}}_{\text{x}}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ +\underset{A_5}{\underbrace{{\displaystyle {\int }_{n_{{\text{crit, H}}_2}}^{\infty }\frac{2}{3}}{n}_{\text{H}}{k}_{{\text{e}}^{-},H_3^{+}}{f}_{{\text{e}}^{-}}{f}_{{\text{H}}_3^{+}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}};\end{array}$$(B.1) $$\begin{array}{l}\frac{\text{d}\langle f_{\text{CO}}\rangle }{\text{d}t}=\langle \frac{1}{3}{n}_{\text{H}}{k}_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{OI}}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{OH}}_x}\rangle +\langle \frac{1}{3}{n}_{\text{H}}{k}_{\text{CI},{\text{H}}_3^{+}}{f}_{\text{CI}}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{CH}}_{\text{x}}}\rangle \\ +\langle \frac{1}{2}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}\rangle -\langle k_{\gamma ,\text{CO}}{f}_{\text{CO}}\rangle \\ =\underset{B_1}{\underbrace{{\displaystyle {\int }_{n_{{\text{criti, H}}_2}}^{\infty }\frac{1}{3}}{n}_{\text{H}}{k}_{\text{OI},{\text{H}}_3^{+}}{f}_{\text{O},\text{I}}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{CH}}_{\text{x}}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ +\underset{B_2}{\underbrace{{\displaystyle {\int }_{\max \left(n_{\text{crit, },{\text{H}}_2},n_{\text{crit, CI}}\right)}^{n_{\text{crit,Co }}}\frac{1}{3}}{n}_{\text{H}}{k}_{\text{CI},{\text{H}}_3^{+}}{f}_{\text{C},\text{ tot }}{f}_{{\text{H}}_3^{+}}{\alpha }_{{\text{CH}}_{\text{x}}}\mathcal{P}{\text{M}}_{}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ +\underset{B_3}{\underbrace{{\displaystyle {\int }_{n_{{\text{crit, H}}_2}}^{n_{\text{crit,Cl }}}\frac{1}{2}}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{\text{C},\text{ tot }}{\alpha }_{{\text{CH}}_{\text{x}}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ -\underset{B_4}{\underbrace{{\displaystyle {\int }_{n_{\text{crit. Co }}}^{\infty }{k}_{\gamma ,\text{Co}}}{f}_{\text{C},\text{ tot }}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}};\end{array}$$(B.2) $$\begin{array}{l}\frac{\text{d}\langle f_{{\text{C}}^{+}}\rangle }{\text{d}t}=-\langle n_{\text{H}}{k}_{{\text{C}}^{+},e^{-}}{f}_{{\text{C}}^{+}}{f}_{{\text{e}}^{-}}\rangle -\langle n_{\text{H}}{k}_{{\text{C}}^{+},{\text{e}}^{-},\text{gr}}{f}_{{\text{C}}^{+}}{f}_{{\text{e}}^{-}}\rangle \\ -\langle \frac{1}{2}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{{\text{C}}^{+}}{f}_{{\text{H}}_2}{\alpha }_{{\text{CH}}_{\text{x}}}\rangle +\langle k_{\gamma ,\text{CI}}{f}_{\text{CI}}\rangle \\ =-\underset{C_1}{\underbrace{{\displaystyle {\int }_0^{n_{\text{crit},\text{CI}}}{n}_{\text{H}}}\left(k_{{\text{C}}^{+},{\text{e}}^{-}}+k_{{\text{C}}^{+},{\text{e}}^{-},\text{gr}}\right)f_{\text{C},\text{ tot }}{f}_{{\text{e}}^{-}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ -\underset{C_2}{\underbrace{{\displaystyle {\int }_{n_{\text{crit, },{\text{H}}_2}}^{n_{\text{crit, Cl}}}\frac{1}{2}}{n}_{\text{H}}{k}_{{\text{C}}^{+},{\text{H}}_2}{f}_{\text{C},\text{ tot }}{\alpha }_{{\text{CH}}_{\text{x}}}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}\\ +\underset{C_3}{\underbrace{{\displaystyle {\int }_{n_{\text{crit, CI }}}^{n_{\text{crit,Co }}}{k}_{\gamma ,\text{CI}}}{f}_{\text{C},\text{ tot }}{\mathcal{P}}_{\text{M}}\left(n_{\text{H}}\right)\text{d}{n}_{\text{H}}}}.\end{array}$$(B.3)

In the above equations, $n_{\text{crit},{\text{H}}_2}$, n_crit,CI; and n_crit,CO are the critical densities for the HI → H₂, C⁺ → C, and CI → CO transitions, respectively in a given cell (see Sect. 2.2.7 and Appendix C). The $f_{{\text{H}}_3^{+}}$ in these equations can be obtained from Eq. (A.5). We note that the limits of integrals A₅ and B₁ are from $n_{{\text{crit,H}}_2}$ to ∞ as the formation of ${\text{H}}_3^{+}$ relies on the presence of H₂ (see Eq. A.4), but the exact expression for $f_{{\text{H}}_3^{+}}$ would be different for n_H < n_crit,CI (where $f_{{\text{C}}^{+}}=f_{\text{C,tot}}$, ƒ_CI = ƒ_CO = 0), n_crit,CI ≤ n_H < n_crit,CO (where ƒ_CI = ƒ_C,tot, $f_{{\text{C}}^{+}}=f_{\text{CO}}=0$), and n_H ≥ n_crit,CO (where ƒ_CO = ƒ_C,tot, $f_{{\text{C}}^{+}}=f_{\text{CI}}=0$). Furthermore, the integrals A₄, B₃ and C₂ will vanish for a given cell if $n_{{\text{crit,H}}_2}>n_{\text{crit,CI}}$.

In Eqs. (B.1)–(B.3), the resulting system of coupled differential equations is solved using the implicit integrator DASSL¹⁰ (Petzold 1982). DASSL uses a variable order (between 1 and 5) backward differential formula to compute the solution of the coupled ODEs after one time step.

B.1 Branching ratios

Often there are multiple outcomes for the reaction between two species. The probability for a given outcome is represented as a branching ratio for that outcome and it denotes the fraction of times that particular outcome will occur. For example, the reaction between C⁺ and H₂ leads to the formation of ${\text{CH}}_2^{+}$, which reacts with H₂ to give ${\text{CH}}_3^{+}$. ${\text{CH}}_3^{+}$ reacts with an electron to give CH or CH₂ in 70% of the cases (these two species are referred to as CH_X in reaction 1 in Table A.1). In the remaining 30% of the cases, ${\text{CH}}_3^{+}+{\text{e}}^{-}$ gives C + 2H (reaction 2). The reaction rates in Table A.1 account for these branching ratios.

Appendix C Calculation of n_crit

The hydrogen in a region can exist in atomic (H I) and molecular (H₂) forms and their mean densities are related to the mean density 〈n_H〉 of H nuclei in the cell as $$\langle n_{\text{H}}{}_{\text{I}}\rangle +2\langle n_{{\text{H}}_2}\rangle =\langle n_{\text{H}}\rangle .$$(C.1)

The mean H₂ fraction in a region is defined as $$f_{{\text{H}}_2}=2\frac{\langle n_{{\text{H}}_2}\rangle }{\langle n_{\text{H}}\rangle }.$$(C.2)

Assuming that hydrogen shows a sharp transition from fully atomic to fully molecular at $n_{\text{H}}=n_{{\text{crit,H}}_2}$, Eq. (C.1) becomes $${\displaystyle {\int }_0^{n_{{\text{crit, H}}_2}}{\mathcal{P}}_{\text{M}}}\text{d}{n}_{\text{H}}+\underset{2\langle n_{{\text{H}}_2}\rangle /\langle n_{\text{H}}\rangle }{\underbrace{{\displaystyle {\int }_{n_{\text{crit, },{\text{H}}_2}}^{\infty }{\mathcal{P}}_{\text{M}}}}}\text{d}{n}_{\text{H}}=1.$$(C.3)

For a log-normal 𝒫_M, this becomes $$\frac{1}{2}\left[1+\mathrm{Erf}\left(\frac{\ln n_{{\text{crit,H}}_2}-\mu }{\sqrt{2}\sigma }\right)\right]+2\frac{\langle n_{{\text{H}}_2}\rangle }{\langle n_{\text{H}}\rangle }=1.$$(C.4)

If 〈n_H〉 and $f_{{\text{H}}_2}$ are known quantities, then the above equation can be written as $$f_{{\text{H}}_2}=\frac{1}{2}\left[1-\mathrm{Erf}\left(\frac{\ln n_{{\text{crit,H}}_2}-\mu }{\sqrt{2}\sigma }\right)\right],$$(C.5)

where µ = 〈n_H〉 + 0.5σ². The critical density $n_{{\text{crit,H}}_2}$ can be obtained by finding the root of the above equation. This equation is equivalent to equation (10) in T15.

For a log-normal+power-law PDF, Eq. (C.3) becomes $$\left[{\displaystyle {\int }_0^{n_{\text{trans }}}{\mathcal{P}}_{\text{M}1}} \text{d}{n}_{\text{H}}+{\displaystyle {\int }_{n_{\text{trans }}}^{n_{{\text{crit, H}}_2}}{\mathcal{P}}_{\text{M}2}} \text{d}{n}_{\text{H}}\right]+2\frac{\langle n_{{\text{H}}_2}\rangle }{\langle n_{\text{H}}\rangle }=1.$$(C.6)

We note that the second term in the brackets will be zero for $n_{{\text{crit,H}}_2}<n_{\text{trans}}$. While using an iterative root-finding method such as the Newton-Raphson method, one needs to evaluate the terms in the brackets differently depending on whether a given guess value for the root $n_{{\text{crit,H}}_2}$ is smaller or larger than n_trans. Thus, in this case, it is not possible to obtain an analytical expression relating 〈n_H〉, $f_{{\text{H}}_2}$, and $n_{{\text{crit,H}}_2}$.

Appendix D Effect of the CRIR on chemistry

In order to investigate the impact on the chemistry of our assumptions regarding the scaling of the CRIR, we extend the results presented in Sect. 4 by repeating the post-processing of the simulated galaxy and considering different options. In Fig. D.1, we plot the ratio of the CO abundance within a cell without and with the upper limit on the CRIR introduced in Sect. 2.2.5. At a given ζ_H, the effect of imposing an upper bound on ζ_H is strongest at low densities and decreases with increasing density. Consequently, the mass-weighted ratio of the CO abundance in the two cases is very close to 1 for CRIR values up to 10⁻¹¹ s⁻¹ H⁻¹, which is larger than the highest observational estimate¹¹ of 3 × 10⁻¹² s⁻¹ H⁻¹ reported in Yang et al. (2023, denoted by a blue arrow in Fig. D.1).

In Fig. D.2, we show the $f_{{\text{H}}_2}$ obtained in each cell when using a variable ζ_H in post-processing. We consider three different CRIRs: (i) ζ_H ∝χ², with an upper limit of 3 × 10⁻¹⁴ s⁻¹ H⁻¹ on ζ_H (left panel); (ii) ζ_H ∝ χ² without any upper limit (middle panel); (iii) the ζ_H − χ relation from Hu et al. (2021, right panel). The median $f_{{\text{H}}_2}$ in each case is shown by a solid black line (also shown in the left panel of Fig. 6). We find slight variations among the three cases for the individual cells. Nevertheless, at low densities (〈n_H〉 ≲ 50 cm⁻³), $f_{{\text{H}}_2}$ decreases with an increase in G₀ (the UV field in LW bands in Habing units, denoted by the colour of the points).

Fig. D.1 Ratio of the ƒCO obtained when using a ζH ∝χ2 without and with an upper limit on ζH, as a function of ζH in the post-processed galaxy. Here we only show the grid cells with ζH greater than the imposed upper limit of 3 × 10−14 s−1 H−1, since the two ƒCO are identical for lower CRIR values. Each grid cell is colour-coded by the mean hydrogen density 〈nH〉 in the cell. The red line shows the mass-weighted mean of the cells in each ζH bin. We also indicate by a blue arrow the highest observed ζH of 3 × 10−12 s−1 H−1 (from Yang etal. 2023).

Fig. D.2 Same as the left panel of Fig. 6, but for a different CRIR used in post-processing (indicated in each panel). The solid black line shows the median $f_{{\text{H}}_2}$ in a given 〈nH〉 bin.

Finally, in Table D.1, we report the total mass of the different chemical species in the post-processed galaxy. Within the accuracy that can be expected from our simplified calculations, HYACINTH provides stable predictions for the masses that are not influenced much by the assumed scaling between the CRIR and the UV flux. The quadratic scaling with no upper limit gives the most discrepant results, with a 30% reduction in the CO mass compensated by an increase in the masses of neutral and ionised atomic carbon. In contrast, the H₂ mass is reduced only by ~7%. The results obtained with our default choice and the CRIR from Hu et al. (2021) are nearly identical.

References