Modeling Molecular Hydrogen in Low-metallicity Galaxies

We use a suite of hydrodynamics simulations of the interstellar medium (ISM) within a galactic disk, which includes radiative transfer, a nonequilibrium model of molecular hydrogen, and a realistic model for star formation and feedback, to study the structure of the ISM and H2 abundance as a function of local ISM properties. We show that the star formation rate and structure of the ISM are sensitive to the metallicity of the gas with a progressively smoother density distribution with decreasing metallicity. In addition to the well-known trend of the H I–H 2 transition shifting to higher densities with decreasing metallicity, the maximum achieved molecular fraction in the ISM drops drastically at Z ≲ 0.2 Z ⊙ as the formation time of H2 becomes much longer than a typical lifetime of dense regions of the ISM. We present accurate fitting formulae for both volumetric and projected fH2 measured on different scales as a function of gas metallicity, UV radiation field, and gas density. We show that when the formulae are applied to the patches in the simulated galaxy, the overall molecular gas mass is reproduced to better than a factor of ≲1.5 across the entire range of metallicities and scales. We also show that the presented fit is considerably more accurate than any of the previous fH2 models and fitting formulae in the low-metallicity regime. The fit can thus be used for modeling molecular gas in low-resolution simulations and semi-analytic models of galaxy formation in the dwarf and high-redshift regimes.


INTRODUCTION
The cold, dense tail of the multiphase interstellar medium (ISM) is generally home to cold atomic gas (e.g., Wolfire et al. 2003) and molecules, such as CO, HCN, H 2 , etc., which play an important role in the thermodynamics of gas in this phase (e.g., Omont 2007;Draine 2011;Galli & Palla 2013).At the same time, molecular gas is one of the very few direct observational probes of this tail (see Carilli & Walter 2013;Saintonge & Catinella 2022, for reviews).Given that stars also form in high-density regions, empirical studies of molecular gas are intricately tied to studies of how star formation occurs in galaxies (e.g., Kennicutt & Evans 2012).
Empirically, it is established that a fairly tight relation between surface densities of molecular gas and star forma-Corresponding author: Ava Polzin apolzin@uchicago.edution exists both within individual galaxies and among different galaxies (the molecular Kennicutt-Schmidt relation; e.g., Kennicutt 1989Kennicutt , 1998;;Wong & Blitz 2002;Bigiel et al. 2008;Baker et al. 2023), the relation that is now well understood theoretically (e.g., Semenov et al. 2019).Observationally, that correlation is tied to CO, but the assumption is that CO traces H 2 reasonably well, with the possible exception of extreme starbursts (Meier et al. 2001;Schruba et al. 2012;Carilli & Walter 2013;Madden 2022).
The existence of such a tight correlation was used as a basis for modeling star formation in galaxy formation simulations (e.g., Robertson & Kravtsov 2008;Gnedin et al. 2009;Jaacks et al. 2013;Christensen et al. 2014) and semi-analytic models (e.g., Popping et al. 2014) and motivated development of theoretical models of molecular gas (e.g., Leroy et al. 2008;Krumholz et al. 2009a;Gnedin & Kravtsov 2011;Gnedin & Draine 2014;Sternberg et al. 2014, see Diemer et al. 2018 for a review).However, most existing models of molecular hydrogen gas fraction are calibrated in the relatively high mass, high metallicity regime.The HI-H 2 transition models that are formulated for lower  gas are often more complex and have additional assumptions and tunable parameters (Krumholz et al. 2009b;Krumholz 2013;Gnedin & Draine 2014;Bialy & Sternberg 2016).Furthermore, most models estimate the abundance of H 2 assuming chemical equilibrium whereby the process of molecular gas formation is not limited in time.
The low-metallicity regime is different.It is generally expected that the formation time of H 2 increases with decreasing metallicity.At the same time, the lifetime of dense regions of the ISM is finite due to a combination of shearing forces and effects of stellar feedback.If the lifetimes of dense ISM regions are shorter than the characteristic H 2 formation time, the molecular fraction in low-metallicity gas may never reach high values, which means that stars in such regions must form from the largely atomic gas (Krumholz 2012;Glover & Clark 2012a;Hu et al. 2016).Indeed, physically star formation can occur in purely atomic gas because cooling and other processes driving the formation of star-forming regions are only mildly affected by the presence of molecular gas (Glover & Clark 2012b).This implies that chemical equilibrium models of  H 2 that assume no time limit to H 2 formation systematically overpredict  H 2 in low-metallicity gas (Krumholz 2012).Conversely, as shown by Glover & Clark (2012a) models that use non-chemical equilibrium calculations of H 2 abundance to estimate star formation rate will underpredict star formation rate, if star-forming regions have low molecular fractions but otherwise form stars with a regular efficiency.
Low metallicities are relevant for modeling the two regimes of galaxy evolution that are at the current frontiers of extragalactic research: the earliest stages of evolution of massive galaxies at  ≳ 5 and evolution of local dwarf galaxies.It is thus important to examine and calibrate the abundance of molecular gas and star formation efficiency in this lowmetallicity regime.Likewise, theoretical models of the interstellar medium and star formation in galaxy formation simulations may potentially be tested by contrasting their results with observational estimates of the H 2 abundance and star formation in dwarf galaxies (e.g., Bolatto et al. 2011;Jameson et al. 2016) and galaxies at high redshifts.
In this paper, we examine the abundance of molecular hydrogen -the dominant mass component of molecular gasusing a suite of realistic simulations of a dwarf galaxy's ISM across a wide range of metallicities.These simulations use a generalized star formation prescription that is not based on the local H 2 abundance (Semenov et al. 2021).Instead, the star formation prescription in the simulations is based on the results of high-resolution magneto-hydrodynamic simulations of star-forming regions (Padoan et al. 2012;Semenov et al. 2016).Most importantly, the model reproduces the abundance and spatial distribution of molecular and atomic gas in NGC 300, and the observed de-correlation between cold molecular gas and clusters of young stars as a function of scale in this galaxy (Semenov et al. 2021).This gives credence to the model as a benchmark that can be used to calibrate star formation and molecular gas abundance in the ISM as a function of its properties.
We present fitting formulae for both volumetric and projected molecular hydrogen fractions that depend on the gas density, gas metallicity, and local ionizing UV field.We show that the fits reproduce known trends in the location and shape of the HI-H 2 transition with metallicity and accurately reproduce both the form of the dependence of molecular fraction on the (volume or column) density and the total molecular gas mass in simulations.In addition, we demonstrate that the structure and behavior of the ISM changes qualitatively when gas metallicity decreases to ≲ 0.1 ⊙ .
The paper is organized as follows.In Section 2, we describe the simulation used to calibrate our  H 2 models.In Section 3, we lay out our simple models for  H 2 to be used in both a volumetric and projected case, and in Section 4, we present tests of the accuracy of the models.Finally, in Section 5, we discuss the implications for models of galaxy formation and compare them against existing  H 2 models, details of which are presented in the Appendix A.

SIMULATIONS
We conduct our analysis using a suite of simulations of an isolated disk galaxy that is initialized with structural properties similar to those observed in the dwarf galaxy NGC 300 (Semenov et al. 2021).We refer the reader to that paper for a detailed description of the simulation setup and implementations of various included physical processes.Below we summarize the key aspects of the simulations.
The fiducial simulation has been shown to reproduce details of the star formation and atomic and molecular gas distributions in NGC 300, including the observed spatial decorrelation of cold gas with recent star formation as a function of averaging scale (or the "tuning fork"; Kruijssen et al. 2019).
The simulations are carried out using the ART -body+hydrodynamics code (Kravtsov 1999;Kravtsov et al. 2002;Rudd et al. 2008;Gnedin & Kravtsov 2011) with selfconsistent modeling of radiative transfer (RT; Gnedin 2014) and non-equilibrium abundance of molecular hydrogen coupled to the local UV radiation field (using the 6-species model described in the Appendix of Gnedin & Kravtsov 2011).With the inclusion of RT and a realistic ISM structure shaped by star formation and feedback, as well as the simulation's maximum resolution of ∼ 10 pc (the average grid cell size is ∼ 22 pc when  H ≥ 0.1 cm −3 ), the simulation offers a highly realistic model for the formation/destruction of molecular hydrogen gas.Star formation in the simulation is not tied to  H 2 but follows the prescription introduced in Semenov et al. (2016).This implementation uses a dynamical model for unresolved turbulence to predict locally variable star formation efficiency instead of assuming a constant tunable value.As was shown in Semenov et al. (2017Semenov et al. ( , 2019Semenov et al. ( , 2021)), modeling star formation efficiency based on local properties of turbulence is important for reproducing the linear molecular KS relation on kiloparsec scales and the spatial decorrelation between young stars (UV sources) and molecular gas regions on sub-kiloparsec scales.It was also shown that this model can reproduce star formation and molecular gas properties both in Milky Way-sized galaxies and in a dwarf galaxy like NGC 300.Results of the analyses presented in this paper should therefore be generally applicable to a wide range of regular galaxies with similar chemical and physical properties.We note, however, that the ISM in strongly starbursting galaxies can have a considerably different density distribution and both the abundance of molecular gas and star formation may behave differently in these environments compared to predictions of our model.
In our fiducial simulation, the metallicity of gas is initialized to have the radial profile similar to the metallicity profile observed in NGC300, after which metallicity is evolved selfconsistently.The final snapshot that we use in our analysis contains cells with gas metallicity ranging from effectively 0 to 4.2  ⊙ (the former corresponds to the halo gas, while the latter gas is in the regions newly enriched by supernova ejecta).
In order to examine the role of metallicity in determining the molecular gas fraction and star-forming gas fraction, we ran a suite of seven re-simulations of the same galaxy, in which gas metallicity is fixed to values approximately evenly distributed in log 10 :  = 0.01, 0.03, 0.1, 0.2, 0.3, 0.6, 1.0  ⊙ .Note that, although gas metallicity is fixed in these runs, all other processes are modelled in the same way.In particular, radiative transfer is performed and UV field varies spatially, reflecting the distribution of sources and absorbing gas.We  mixed  a 901 1.9 × 10 7 1.1 × 10 9 1.9 × 10 9 For our fiducial mixed metallicity run, the grid cells range in  from 10 −20 to ∼ 4  ⊙ , with a median metallicity of ∼ 0.4  ⊙ in higher density ( H ≥ 0.1 cm −3 ) gas.
use the variation of the UV flux within a run to study the dependence of the molecular fraction on this flux at a given gas metallicity.Details of these runs are presented in Table 1.
Figure 1 presents face-on and edge-on views of the gas distribution in three of these runs.The gas density distribution varies substantially with metallicity, with a more homogeneous density distribution at low  and a more flocculent density distribution at higher .This trend is also apparent in the fiducial run with a non-uniform metallicity distribution shown in the right column, in which inner high- regions are similar to the  = 1  ⊙ run, while the outer lower- regions have a much smoother gas distribution.
We compare the projected molecular hydrogen fraction,  H 2 ≡  H 2 /( HI + 2 H 2 ), from the single  runs of the simulation to observations of  H 2 in the Large and Small Magellanic Clouds (Tumlinson et al. 2002;Bolatto et al. 2011;Welty et al. 2012) in Figure 2. We select the single  run closest to the gas metallicities of the LMC and SMC, adopting our 0.6  ⊙ run as an analog for the LMC and our 0.2  ⊙ run as an analog for the SMC.The location of the HI-H 2 transition is in good agreement between the observations and the simulation.
Although we do not exactly match the high  H 2 observations of the SMC from Bolatto et al. (2011) in our 0.2  ⊙ run, there are two effects that likely contribute to this discrepancy.The novel method for measuring Σ H 2 used in Bolatto et al. (2011) does not fully distinguish between H 2 and cold HI, potentially yielding a slightly higher molecular fraction that is ultimately more reflective of the fraction of cold neutral gas.Additionally, even though the highest resolution of our simulation grid cells is 10 pc, the effective resolution is several times this, which means that we are not sensitive to features below this effective resolution scale.It is then possible that the less prevalent high density, highest  H 2 regions are averaged to somewhat lower molecular fraction.

MODELING THE MOLECULAR GAS FRACTION
In this section, we present two versions of fitting formulae suitable for application in different regimes: fits to a volumetric  H 2 fitted to the cell-by-cell distribution of the molecular gas fraction in simulations, which can be used in high-resolution simulations (Section 3.1), and fits to projected molecular fraction  H 2 , fitted to projected 2D maps averaged on different spatial scales, which can be used in low-resolution simulations and semi-analytic models (Section 3.2).

Volumetric 𝑓 H 2 model
We model the molecular hydrogen fraction as a function of hydrogen gas density, metallicity, and UV field strength in individual simulation grid cells.The functional form of  H 2 fit is motivated by the fact that we expect  H 2 to exhibit a fairly sharp transition at a certain density or column density and saturate at values close to some maximum value  H 2 ,max .
We thus choose a sigmoid-like function: . (1) This form allows us to account for the fact that the maximum possible molecular hydrogen fraction  H 2 ,max varies as a function of metallicity.We parameterize this dependence as where and Here  0 = 3.5 × 10 −17 cm 3 s −1 is the rate of H 2 formation on dust grains (see Wolfire et al. 2008).
We define  as where The functional form and parameter values in these equations were chosen so that the average trend of  H 2 with gas density in the simulation and the maximum values of  H 2 are reproduced.Equation 5 accounts for the dependence of the location and shape of the HI-H 2 transition on UV field strength and metallicity.The density at which this transition occurs is set by the value of  tr , while the metallicity-dependent prefactors are responsible for the changing slope of  H 2 vs.  H .
The value of  tr can be measured from the simulation and fit directly.
Given the very low molecular gas fraction at low metallicities, the parameterization using  1/2 -the density at which  H 2 = 0.5 in Gnedin & Kravtsov (2011) and Gnedin & Draine (2014) at higher metallicities does not work in our lowest- runs as  H 2 never reaches 0.5 (see Figure 3).Instead, we define  tr as the hydrogen number density (cm −3 ) at which the molecular hydrogen fraction is 5 × 10 −4 , which characterizes the transition even in this low metallicity regime.
Figure 3 shows that the atomic-to-molecular transition occurs at lower  H for higher metallicities, while Figure 4 shows that  tr behaves like a power law with respect to both  MW and .This power law at each discrete metallicity for binned values of  MW on a volumetric cell-by-cell basis can be approximated by Here  MW is the free-space1 UV flux relative to the MW value  MW =  1000 / MW , where  1000 is the interstellar UV flux at 1000 Å  MW = 10 6 photons cm −2 s −1 ster −1 eV −1 (Draine 1978;Mathis et al. 1983), and  is the dust-to-gas ratio, which we assume to be equal to the unnormalized mass fraction of heavy elements in the gas.
We then determine fit parameters  and  as a function of  MW using a simple least squares fit of simulation results:

MW
The strength of the -dependence of  tr becomes weaker at higher metallicities for all  MW .To reflect this saturation of the metallicity dependence at near-solar metallicities, we add a correction term assuming that solar metallicity corresponds to mass fraction of 0.0199: To avoid non-physical, negative  tr at high  and very low  MW , the floor of  tr can be set explicitly.We choose a minimum value of  tr = 0.1 cm −3 given that we anticipate very little cold, dense molecular hydrogen gas at  H ≲ 0.1 cm −3 , but this can be set even lower without affecting the accuracy of the overall model.The results of this fit to  tr are shown in Figure 4 overplotted on the measured location of this transition.
Equations 1-8 can be used to estimate  H 2 in the highdensity gas (see Figure 3), as this gas constitutes most of gas mass in galaxies.The molecular fraction in the low-density unshielded gas does not require a fitting formula and can be obtained by simply equating the H 2 formation and photodissociation rates (see, e.g., Section 3.2 in Wolfire et al. 2008): where  0 ≈ 3.5 × 10 −17 (/0.019)cm 3 s −1 is the rate of H 2 formation on dust grains assumed here to scale linearly with ,  H is the number density of hydrogen nuclei, and  = 4.7 × 10 −11 s −1 is the unshielded photodissociation rate in the local interstellar UV field.

Model for projected molecular fraction
In observations, low-resolution simulations, and semianalytic models one often needs to work with the projected mass densities and we thus present fitting formulae for projected molecular fraction below.To distinguish it from the volumetric one we will denote the projected fraction as To obtain projected molecular fractions  H 2 on different spatial scales in the simulations, we use the face-on projection of the simulated galaxy with gas properties binned on . Molecular hydrogen gas fraction (  H 2 ) in computational grid cells of different hydrogen gas density,  H , for each run in our suite.The values in simulation cells are shown in blue, with the model described in Equation 1 overplotted in pink.The model explicitly captures the behavior of  H 2 at high  H owing to the metallicity-and density-dependent cap value  H 2 ,max (see Equation 2).The model was not designed to reproduce the tail of high  H 2 at lower densities for the reasons discussed in Section 4.1.The  tr model results as a function of metallicity and  MW (Equation 7) are shown by magenta symbols connected by the dashed lines.
grids with physical cell sizes 10 pc to 1 kpc.The gas surface (column) density is computed simply as the gas mass (atom number) in each bin divided by its area.The UV flux and metallicity in each bin are estimated as the gas-density weighted averages of  MW and  in the computational cells enclosed in a given bin.
Projected molecular fractions as a function of column density,  H 2 ( H ), are shown in Figure 5 for two representative averaging scales -30 pc and 300 pc.Similarly to the volumetric molecular fractions, we use the sigmoid-like functional form of the model for  H 2 : where and and where, as before,  0 = 3.5 × 10 −17 cm 3 s −1 (see Wolfire et al. 2008), and  in Equation 11is and  is the resolution of the projected map from the simulation, i.e. the scale on which surface densities and fractions are averaged.The transition column density,  tr (,  MW , ), is defined as the median column density of all cells with molecular hydrogen fractions between 5×10 −5 and 5×10 −4 .This region in the  H 2 − H parameter space was chosen because  H 2 ( H ) is increasing sharply with increasing column density and  H 2 range is sufficiently low to be used at very small metallicities.As in the volumetric model, the location of the transition is set by  tr (Figure 6), while the slope of  H 2 ( H ) is set by the pre-factors, which have a weak dependence on the scale  (in parsecs) and metallicity,  = log 10 : As for the volumentric model, the functional forms and their parameter values are chosen so that equation 11 reproduces the mean trend of  H 2 as a function of  H in the simulations for each metallicity, UV flux and scale. corr is the correction factor introduced to address two distinct limitations of our original fit to the measured median column density of the HI-H 2 transition.Due to a paucity of grid cells at larger scales ( > 100 pc), we only explicitly fit for  tr in the  = 10, 30, and 100 pc cases, which means that our parameterization does not account for the behavior on larger spatial scales.In addition, given that  H 2 is defined assuming that the transition occurs at 1.65 × 10 −4 as defined by  tr , while the median molecular fraction between 5 × 10 −5 and 5 × 10 −4 varies subtly with scale and metallicity.We capture these effects using the following functional form: The validity and accuracy of the volumetric  H 2 fitting formulae (Equations 1-8) can be guaged by comparing the  H 2 according to the fit to the simulation  H 2 for all grid cells in the simulation with molecular fractions larger than 10 −5 .We do this for all fixed metallicity runs and for the fiducial run with non-uniform metallicity, which was not used in deriving the fit.Note also that the latter includes cells with metallicities outside of the range within which we calibrated the model.
Figure 3 shows a good qualitative agreement between the molecular fractions produced by the model fits and the simulation results for each run.Because we impose a strict  H 2 ,max condition, the model  H 2 distribution has a sharp boundary at the highest fraction values at each density.In principle, one can introduce scatter around the model relations to reproduce the tail of high  H 2 cells at small densities.However, we do not think it is worthwhile for two reasons.First, as we show next, the presented fit accurately recovers molecular mass,  H 2 , for the galaxy at all metallicities.This means that the molecular mass in this tail is fairly small.Second, the bulk of the high  H 2 gas at low densities could be due to numerical diffusion of radiation and molecular gas around star-forming regions and thus may be a non-equilibrium artifact.
The accuracy of the fit in reproducing the molecular content of the ISM in the simulated galaxy can be assessed by comparing the total H 2 mass estimated with when the molecular fraction in each computational cell is estimated using the fit to the actual mass in the simulation.The left panel of Figure 7 shows the ratio  H 2 ,sim / H 2 ,mod as a function of metallicity for the volumetric model as magenta open circles and indicates that the fit is accurate in predicting H 2 mass to ≲ 25%.
To test whether results depend on the time snapshot of the simulation we estimated the ratio of the actual-to-modelpredicted molecular mass at a series of different simulation snapshots.These estimates are shown in Figure 8, which shows that generally, the fit accuracy is similar at most snapshots, except for a single snapshot where the ratio increased to ≈ 1.8 one where likely non-equilibrium processes related to a local starburst changed the ISM significantly.
The figure shows that variation of the mass ratio from snapshot to snapshot increases significantly in lower metallicity runs.This is likely related to the rapidly increasing H 2 formation time with decreasing metallicity.Indeed, using  chem = 105 (/ ⊙ ) −1 ( H /cm −3 ) −1 (10/   ) Myr from Krumholz (2012) and assuming a clumping factor   = 10 due to turbulence on the scale of molecular clouds (see Appendix A.7 of Gnedin & Kravtsov 2011), the H 2 formation time for gas of  = 0.01  ⊙ and number density  H = 50 cm −3 should be 210 Myr, while it is only ≈ 2.1 Myr for  =  ⊙ gas of the same density.Given that the typical lifetime of a molecular cloud is significantly shorter, ∼ 5 − 15 Myr (Semenov et al. 2017), it is not surprising that  H 2 is generally suppressed in low- gas ( ≲ 0.1  ⊙ ), where the time for the gas to reach chemical (HI-H 2 ) equilibrium is longer than the timescales on which molecular clouds persist without disrup- Given that we only fit explicitly for averaging scales  between 10 and 100 pc due to the number of available grid cells, we do not include the measured values for scales between 300 and 1000 pc.The model results as a function of   , , and  (Equation 17) are shown by the magenta symbols connected by dashed lines.
tion.In low-metallicity runs, H 2 abundance is much more susceptible to disruption of individual star-forming regions (which are also fewer), which leads to larger variations of the H 2 abundance.We note that the fit should not to be extrapolated to the zero metallicity case.In practice, however, a very small value of metallicity should return reasonable results.

The projected 𝐹 H 2 fit
Figure 5 shows good agreement between the projected molecular fraction,  H 2 , estimated using fit (equations 11-18) and simulation results across metallicities and across averaging scales.We show this explicitly for two representative scales of 300 pc and 30 pc, both of which are fairly wellpopulated by simulation grid cells with  H 2 ≥ 10 −5 even at the lowest metallicities.As with the volumetric fit, the imposed maximum fraction  H 2 ,max results in a hard upper boundary on  H 2 ( H ) (see Section 4.1).Also, the scatter of  H 2 in the model at a given number density is somewhat smaller than in the simulations.As we argued in the discussion of the volumetric model comparisons above, the additional scatter can be added to the model, but the amount of molecular gas associated with the tails of the distribution is fairly small.Indeed, Figure 7 shows that the total molecular gas mass estimated using the fit formulae is accurate to better than a factor of ≲ 1.5across the full probed range of metallicities and averaging scales.

Implications for galaxy formation modeling
In simulations and analytical models of galaxies, molecular gas is sometimes used as a proxy for star-forming gas, which is motivated by the observed constant depletion times of molecular gas in normal (non-starburst) galaxies of metallicities ≈ 0.1 − 1 ⊙ (see Introduction).
Figure 9, however, shows that the depletion time of molecular hydrogen gas is expected to change by nearly three orders of magnitude from 0.01 ⊙ to  ⊙ in our simulations, while the SFR is changing only by a factor of ten over the same metallicity range.Thus, according to our simulations the star formation rate at low metallicities is a nonlinear function of the molecular mass, which implies that the fraction of stars forming in atomic gas increases with decreasing metallicity.
As discussed by Krumholz (2012), the use of chemical equilibrium models to estimate H 2 abundance for star formation rate calculations can partly compensate for this trend and will result in higher SFR at low metallicity.However, this is hardly justified because this effectively trades one error for another by artificially overestimating  H 2 in the environment where the fraction is actually low.A much better alternative The ratio of H 2 mass measured in the simulations using the grid cells with  H 2 ≥ 10 −5 to the H 2 mass predicted by the model using the densities, metallicities, and UV fluxes in the individual simulation cells as a function of metallicity for both the volumetric case (magenta, left) and the scale-dependent projected measurements.Notably, the model is accurate to better than a factor of ≈ 1.5 in most cases for both the fixed metallicity runs and the fiducial run with non-uniform metallicity (magenta, right) across a range of scales. is to switch to a star formation prescription that captures the metallicity dependence of the actual star formation.

Comparisons with other models
Motivated by the paucity of reliable HI-H 2 transition models for the  ≲ 0.2  ⊙ regime (corresponding to metallicities typical of dwarf galaxies and galaxies at high redshift), we constructed fitting formulae for the molecular fraction as a function of gas density (or column density), local UV flux, and metallicity down to  = 0.01 ⊙ .This metallicity range includes the smallest metallicities observed in galaxies in the ultra-faint regime (Hidalgo 2017;Simon 2019).Simultaneously, as we showed above these fitting formulae perform well at higher  up to the solar metallicity.
Here we compare the fitting formulae presented in this paper and a number of models of molecular fraction in the literature (namely, Krumholz et al. 2009b;Gnedin & Kravtsov 2011;Krumholz 2013;Sternberg et al. 2014;Gnedin 2014) with simulation results.The models parameterize average  H 2 and  H 2 as a function of gas density, metallicity, and UV field strength (the KMT09b model only accounts for the dependence on column density and metallicity; see Appendix A for details of the model implementations).
Figure 10 shows comparisons of the average  H 2 and  H 2 to the molecular fraction as a function of volumetric and projected densities in simulations.In order to have a more direct comparison between existing models and the fits presented in this study, we also show predictions of each model when applied to the individual computational grid cells in our simulation for our single metallicity runs with the lowest and highest  in Figures 11 and 12 .The depletion time ( dep =  gas /SFR) as a function of metallicity and gas species relative to the star formation rate (SFR) as a function of metallicity.We use  gas within 5 kpc of the simulation center.While  dep,H 2 varies by a factor of 260 between  = 0.01  ⊙ and 1  ⊙ , the star formation rate (defined here by the mass of stars formed over the last 10 Myr) only varies by a factor of 11.This is indicative of the fact that the star formation in low metallicity galaxies is not directly tied to H 2 abundance.We also include the inferred SFR and H 2 and HI  dep in the Large and Small Magellanic Clouds from Jameson et al. (2016) as squares, which we correct by a factor of 1.35 for the presence of He, and the inferred SFR and  dep for NGC 300 digitized from Kruijssen et al. (2019) as triangles.
of  H 2 as a function of density in the KMT09b model is likely due to its neglect of dependence on  MW , which drives the scatter in the simulation and in other models, including the model calibrated in this study.
Given that most existing models have been calibrated in the high-metallicity regime, it is not surprising that their accuracy improves with increasing metallicity (see Table 2).Interestingly, the accuracy of the K13 model is better at low metallicities for volumetric grid cells, but the match to the shape and behavior of the HI-H 2 transition is better at high .This seems to be driven by their very steep functional form of  H 2 , so that only a few grid cells with high  H 2 contribute to the inferred mass at lower , while at higher metallicity, the transition  H is underestimated, leading to the model underprediction of  H 2 .
In the case of the projected models, KMT09b performs consistently well for  ≳ 0.1  ⊙ .The shape of  H 2 ( H ) is very close to what we observe in the simulation for these metallicities and the transition between atomic and molecular gas is consistent with the location of the transition in the simulation for  ≳ 0.4  ⊙ .GK11 is the next most accurate of the existing models.The lowest metallicity runs result in an inferred  H 2 = 0  ⊙ with this model, but this is expected behavior for the projected GK11 model, which is known to not be accurate for  ≲ 0.01  ⊙ (Gnedin & Kravtsov 2011).
With regard to the GD14 model, it is worth noting that this model includes a phenomenological account for the H 2 self-shielding due to line overlap.This results in near independence of  H 2 on dust abundance (and thus metallicity) for  ≲ 0.2 MW .The simulations used here do not include any accounting for such line overlap and thus a part of the difference of GD14 model and simulation results at  ≲ 0.2 ⊙ may be due to this difference.Otherwise, the GK11 and GD14 models are quite similar and thus the accuracy of the GK11 model should be comparable to predictions of the GD14 model without the line overlap effect.The similarity of these two models is the reason why their estimated molecular mass is similar at  ≲ 0.2 ⊙ (see Table 2).
The fit presented in this paper reproduce the total molecular hydrogen mass in the simulations considerably better than previous models.Even restricting K13 to  ≲ 0.1  ⊙ , where their model results in masses within a factor of ∼two of those measured in simulation, it appears to be a coincidence given the steep relation and few high  H 2 cells.

CONCLUSIONS
In this work we presented fitting formulae for volumetric (equations 1-8) and projected (equations 11-18) molecular gas fractions.The fits consist of a set of simple scalings calibrated to reproduce mean trends measured in simulations of a realistic dwarf galaxy similar to NGC 300 (Semenov et al. 2021).Both volumetric and projected fits parameterize the molecular fraction as a function of gas density, gas metallicity, and the strength of the local free space ionizing UV field.
Our main results and conclusions are as follows: • We show that the interstellar medium in the simulated galaxy changes qualitatively when gas metallicity is varied by two orders of magnitude.The density distribution becomes increasingly non-uniform as metallicity increases from 0.01 ⊙ to  ⊙ (see Figure 1).
• In addition to the well-known trend of the HI-H 2 transition shifting to higher densities with decreasing metallicity, the maximum achieved molecular fraction in the interstellar medium drops drastically to values much less than 1 at  ≲ 0.2  ⊙ (see Figures 3 and 5), while the dependent of molecular fraction on density becomes less steep.
• We show that accurate fitting functions for volumetric and projected molecular fractions can be constructed if they account for the dependence on gas density, • We also show that the volumetric (projected) molecular fraction fit is applied to individual cells (projected patches) reproduces the total molecular mass in the simulated galaxies to a factor ≲ 1.25 (≲ 1.5) across the entire explored range of galaxy metallicities (Figures 7  and 8).This is considerably better than the estimates using existing models of the molecular hydrogen fraction and HI-H 2 transition (see Figures 10 and 11 and Table 2).
The presented model should be useful in modeling molecular gas abundance in simulations that do not include explicit modeling of H 2 and low-resolution simulations and semianalytical models.However, we argue that star formation modeling in simulations should not be based on molecular gas fraction because our simulation results indicate that star formation rate becomes a non-linear function of molecular gas density at metallicities < 0.1 ⊙ due to non-equilibrium effects (see Section 5 and Fig. 9).As an alternative, the star formation efficiency and depletion time can be calibrated using such simulations and we will present such calibrations in the follow-up work.
The simulations used in this work were carried out on the Midway cluster maintained by the University of Chicago Research Computing Center.A.K. was supported by the National Science Foundation grants AST-1714658 and AST-1911111 and NASA ATP grant 80NSSC20K0512.V.S. is grateful for the support provided by Harvard University through the Institute for Theory and Computation Fellowship.Software: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), yt (Turk et al. 2011), WebPlotDigitizer (Rohatgi 2022) APPENDIX A. IMPLEMENTATION OF OTHER HI-H 2 MODELS Given that there are different formulations of the HI-H 2 transition models which we use in comparisons in Section 5, we offer details of how  H 2 and  H 2 were computed as a function of density and column density for each model.

Figure 1 .
Figure1.Face-on (top)  and edge-on (bottom) gas density slices of our simulated NGC 300-like galaxy in our runs with the lowest and highest  fixed metallicity as well as our fiducial simulation with variable metallicity.

Figure 2 .
Figure 2. Molecular hydrogen gas fraction in projected patches with  = 10 pc as a function of gas column density.We compare the  = 0.2  ⊙ and  = 0.6  ⊙ simulation runs to observations of  H 2 in the Small Magellanic Cloud (SMC; top) and Large Magellanic Cloud (LMC; bottom) respectively.

Figure 4 .
Figure 4. Values of  tr estimated for each metallicity and  MW bin as a function of  MW are shown for each of the seven runs with fixed metallicities (blue symbols connected by solid lines).The shaded region corresponds to the 16th and 84th percentile of  tr in each bin.The  tr model results as a function of metallicity and  MW (Equation7) are shown by magenta symbols connected by the dashed lines.

Figure 5 .
Figure 5. Molecular hydrogen gas fraction ( H 2 ) in projected patches of size  of different hydrogen gas column density,  HI +  H 2 , for each run in our simulation suite at two representative scales,  = 30 pc and  = 300 pc.The simulation values are shown in blue, while the model described in Equation 1 is overplotted in pink.

Figure 6 .
Figure6.The median column density for which the molecular fraction is between 5 × 10 −5 and 5 × 10 −4 , which marks the location of the HI-H 2 transition, as a function of  MW and  for each of our single metallicity runs.The measured values are shown by the dark blue symbols connected by lines.Given that we only fit explicitly for averaging scales  between 10 and 100 pc due to the number of available grid cells, we do not include the measured values for scales between 300 and 1000 pc.The model results as a function of   , , and  (Equation17) are shown by the magenta symbols connected by dashed lines.

Figure 8 .
Figure 8.The ratio of the H 2 mass in the simulation and the H 2 mass inferred from the volumetric model (Equation1) as a function of snapshot time.The most advanced snapshot used in calibrating our model is shown by the magenta point.We show the evolution of this ratio for the two lowest metallicity and two highest metallicity runs we use in this work.The figure indicates that results of the test shown in Figure7are not sensitive to the specific output used.
Figure 9.The depletion time ( dep =  gas /SFR) as a function of metallicity and gas species relative to the star formation rate (SFR) as a function of metallicity.We use  gas within 5 kpc of the simulation center.While  dep,H 2 varies by a factor of 260 between  = 0.01  ⊙ and 1  ⊙ , the star formation rate (defined here by the mass of stars formed over the last 10 Myr) only varies by a factor of 11.This is indicative of the fact that the star formation in low metallicity galaxies is not directly tied to H 2 abundance.We also include the inferred SFR and H 2 and HI  dep in the Large and Small Magellanic Clouds fromJameson et al. (2016) as squares, which we correct by a factor of 1.35 for the presence of He, and the inferred SFR and  dep for NGC 300 digitized fromKruijssen et al. (2019) as triangles.

Table 1 .
Basic properties of the snapshots (within 15 kpc of the simulation center), which were fit to construct our models.Each version of the simulation was run for sufficient time that  H 2 was not evolving significantly between snapshots.RunTime/Myr  H 2 / ⊙  HI / ⊙  gas / ⊙