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 formation exists both within individual galaxies and among different galaxies (the molecular Kennicutt-Schmidt (KS) relation; e.g., Kennicutt 1989Kennicutt , 1998;;Wong & Blitz 2002;Bigiel et al. 2008;Baker et al. 2023); that relation 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, hereafter GK11;Gnedin & Draine 2014, hereafter GD14;Sternberg et al. 2014, hereafter S14; see Diemer et al. 2018 for a review). However, most existing models of molecular hydrogen gas fraction are calibrated in the relatively highmass, high-metallicity regime. The H I-H 2 transition models that are formulated for lower-Z gas are often more complex and have additional assumptions and tunable parameters (Krumholz et al. 2009b, hereafter KMT09b;Krumholz 2013, hereafter K13;GD14;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 (Glover & Clark 2012a;Krumholz 2012;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 f H 2 that assume no time limit to H 2 formation systematically overpredict f H 2 in low-metallicity gas (Krumholz 2012). Conversely, as shown by the Glover & Clark (2012a) models that use nonchemical equilibrium calculations of H 2 abundance to estimate star formation rate (SFR) will underpredict the SFR, 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 z 5 and the evolution of local dwarf galaxies. It is thus important to examine and calibrate the abundance of molecular gas and star formation efficiency in this low-metallicity regime. Likewise, theoretical models of the ISM 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 magnetohydrodynamic simulations of starforming 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 decorrelation 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 H I-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 change qualitatively when gas metallicity decreases to 0.1 Z e .

The paper is organized as follows: Section 2 describes the simulation used to calibrate our f H 2 models. In Section 3, we lay out our simple models for f 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 f H 2 models, the details of which are presented in the Appendix.

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 the averaging scale (or the tuning fork; Kruijssen et al. 2019).

The simulations are carried out using the Adaptive Refinement Tree (ART) N-body+hydrodynamics code (Kravtsov 1999; Kravtsov et al. 2002;Rudd et al. 2008;GK11), with self-consistent modeling of radiative transfer (RT; Gnedin 2014) and nonequilibrium abundance of molecular hydrogen coupled to the local UV radiation field (using the six species model described in the appendix in GK11). 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 n 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 f 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 the predictions of our model.

In our fiducial simulation, the gas metallicity is initialized to have a radial profile similar to the metallicity profile observed in NGC 300, after which it evolves self-consistently. The final snapshot that we use in our analysis contains cells with gas metallicity ranging from effectively 0-4.2 Z e (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 resimulations of the same galaxy, in which gas metallicity is fixed to values approximately evenly distributed in Z log 10 : Z = 0.01, 0.03, 0.1, 0.2, 0.3, 0.6, and 1.0 Z e . We conservatively assume a dust-to-gas ratio that scales linearly with metallicity (Wolfire et al. 2008). The dust-to-gas ratio likely decreases faster than linearly at low metallicity. However, this will only make the decrease of molecular fraction stronger at lower metallicities than what we find here, thereby strengthening our conclusion that it cannot be a linear tracer of star formation in the metal-poor regime. Note that, although gas metallicity is fixed in these runs, all other processes are modeled in the same way. In particular, RT is performed, and the UV field varies spatially, reflecting the distribution of sources and absorbing gas. We 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 Z and a more flocculent density distribution at higher Z. This trend is also apparent in the fiducial run with a nonuniform metallicity distribution shown in the right column, in which inner high-Z regions are similar to the Z = 1 Z e run, while the outer lower-Z regions have a much smoother gas distribution.

We compare the projected molecular hydrogen fraction,

2 , from the single-Z runs of the simulation to observations of F H 2 in the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC; Tumlinson et al. 2002;Bolatto et al. 2011;Welty et al. 2012) in Figure 2. We select the single-Z run closest to the gas metallicities of the LMC and SMC, adopting our 0.6 Z e run as an analog for the LMC, and our 0.2 Z e run as an analog for the SMC. The location of the H I-H 2 transition is in good agreement between the observations and the simulation.

Although we do not exactly match the high F H 2 observations of the SMC from Bolatto et al. (2011) in our 0.2 Z e run, there are two effects that likely contribute to this discrepancy. The novel method for measuring S H 2 used in Bolatto et al. (2011) does not fully distinguish between H 2 and cold H I, 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

Note. Each version of the simulation was run for a sufficient time, so M H2 did not evolve significantly between snapshots. a For our fiducial mixed metallicity run, the grid cells range in Z from 10 -20 to ∼4 Z e , with a median metallicity of ∼0.4 Z e in higher-density (n H 0.1 cm -3 ) gas.

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 highdensity, highest F H 2 regions are averaged to a 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 f 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 F 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 f 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 f H 2 fit is motivated by the fact that we expect f H 2 to exhibit a fairly sharp transition at a certain density or column density and saturate at values close to some maximum value of f H ,max 2 . We thus choose a sigmoid-like function:

This form allows us to account for the fact that the maximum possible molecular hydrogen fraction f H ,max 2 varies as a function of metallicity. We parameterize this dependence as

Here R 0 = 3.5 × 10 -17 cm 3 s -1 is the rate of H 2 formation on dust grains (see Wolfire et al. 2008) and Z is the metallicity in solar units. We define x as

The functional form and parameter values in these equations were chosen so that the average trend of f H 2 with gas density in the simulation and the maximum values of f H 2 are reproduced. Equation (5) accounts for the dependence of the location and shape of the H I-H 2 transition on UV field strength and metallicity. The density at which this transition occurs is set by the value of n tr , while the metallicity-dependent prefactors are responsible for the changing slope of f H 2 versus n H . The value of n tr can be measured from the simulation and fit directly.

Given the very low molecular gas fraction at low metallicities, the parameterization using n 1/2 -the density at which = f 0.5 H 2 in GK11 and GD14 at higher metallicities does not work in our lowest-Z runs as f H 2 never reaches 0.5 (see Figure 3). Instead, we define n tr as the hydrogen number density (per cubic centimeter) 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 n H for higher metallicities, while Figure 4 shows that n tr behaves like a power law with respect to both U MW and Z. This power law at each discrete metallicity for binned values of U MW on a volumetric cell-by-cell basis can be approximated by

Here U MW is the free-spacefoot_1 UV flux relative to the MW value U MW = J 1000 /J MW , where J 1000 is the interstellar UV flux at 1000 Å, J MW = 10 6 photons cm -2 s -1 sr -1 eV -1 (Draine 1978; Mathis et al. 1983), and D 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 a and b as a function of U MW using a simple least squares fit of simulation results:

34.7 2.25 0.0199 53.9 . MW 0.32 0.3 MW 0.31

The strength of the D dependence of n tr becomes weaker at higher metallicities for all U MW . To reflect this saturation of the metallicity dependence at near-solar metallicities, we add a correction term assuming that solar metallicity corresponds to a mass fraction of 0.0199:

To avoid nonphysical, negative n tr at high Z and very-low U MW , the floor of n tr can be set explicitly. We choose a minimum value of = n 0.1 tr cm -3 given that we anticipate very little cold, dense molecular hydrogen gas at n 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 n tr are shown in Figure 4, overplotted on the measured location of this transition.

Equations (1)-( 8) can be used to estimate f H 2 in the highdensity gas (see Figure 3), as this gas constitutes most of the 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 R 0 ≈ 3.5 × 10 -17 (D/0.019) cm 3 s -1 is the rate of H 2 formation on dust grains assumed here to scale linearly with D, n H is the number density of hydrogen nuclei, and I = 4.7 × 10 -11 s -1 is the unshielded photodissociation rate in the local interstellar UV field.

We note that the volumetric f H 2 fit presented in this section is applicable at densities n H  10 3 cm -3 probed in our simulations. The upper limit on densities in the simulations is determined by the star formation and feedback model used in these simulations and not by resolution. The gas would have reached higher densities at the resolution of the simulations if not for the exponential increase in the star formation efficiency per free-fall time assumed in our model and strong stellar feedback in the form of thermal and turbulent energies and momentum injection. Star formation and feedback disperse the gas before it can reach very high densities.

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 the projected molecular fraction below. To distinguish it from the volumetric one, we denote the projected fraction as

To obtain projected molecular fractions F H 2 on different spatial scales in the simulations, we use the face-on projection of the simulated galaxy with gas properties binned on grids with physical cell sizes of 10 pc-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 U MW and Z in the computational cells enclosed in a given bin.

Projected molecular fractions as a function of column density, F H 2 (N H ), are shown in Figure 5 for two representative averaging scales-30 and 300 pc. Similarly to the volumetric molecular fractions, we use the sigmoid-like functional form of the model for F H 2 :

where

and where, as before, R 0 = 3.5 × 10 -17 cm 3 s -1 (see Wolfire et al. 2008), and x in Equation ( 11) is

1.35 0.01 10 pc 3.4 0.6 , 16 0.25 0.6 0.02 and S 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, ( ) N Z U S , , 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 -F N H H 2 parameter space was chosen because ( ) F N H H 2 is increasing sharply with increasing column density, and the F 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 N tr (Figure 6), while the slope of ( )

is set by the prefactors, which have a weak dependence on the scale S (in parsecs) and metallicity, z = Z log 10 :

w y y y S 10 , , 0.27 0.01 9.25 9.64 , exp 0.5 1.5 6.84 , 21.96 0.19 log . 17 w y S tr MW , corr 2 norm 2 norm 10

As for the volumetric model, the functional forms and their parameter values are chosen so that Equation (11) reproduces the mean trend of F H 2 as a function of N H in the simulations for each metallicity, UV flux, and scale. N corr is the correction factor introduced to address two distinct limitations of our original fit to the measured median column density of the H I-H 2 transition. Due to a paucity of grid cells at larger scales (S > 100 pc), we only explicitly fit for N tr in the S = 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 F H 2 is defined, assuming that the transition occurs at 1.65 × 10 -4 as defined by N 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:

1 0.13 log 0.1 log 10 pc . 18 corr 10 10

Comparison of Fits to the Simulation

Volumetric f H 2 Fit

The validity and accuracy of the volumetric f H 2 fitting formulae (Equations ( 1)-( 8)) can be gauged by comparing the f H 2 according to the fit to the simulation f 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 nonuniform 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 f H , max 2 condition, the model f 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 f 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, M 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 f 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 nonequilibrium 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 by the model with the total simulation H 2 mass summed across each computational cell. The left panel of Figure 7 shows the ratio M M

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 accuracy of the fit is similar at most snapshots, except for a single snapshot where the ratio increased to ≈1.8, where likely nonequilibrium 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 Krumholz (2012) and assuming a clumping factor f c = 10 due to turbulence on the scale of molecular clouds (see Appendix A.7 in GK11), the H 2 formation time for gas of Z = 0.01 Z e and number density n H = 50 cm -3 should be 210 Myr, while it is only ≈2.1 Myr for Z = Z e 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 f H 2 is generally suppressed in low-Z gas (Z  0.1 Z e ), where the time for the gas to reach chemical (H I-H 2 ) equilibrium is The median column density for which the molecular fraction is between 5 × 10 -5 and 5 × 10 -4 , which marks the location of the H I-H 2 transition, as a function of U MW and S 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 S 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 U MW , Z, and S (Equation ( 17)) are shown by the magenta symbols connected by dashed lines. - f 10 H 5 2

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 nonuniform metallicity (magenta, right) across a range of scales. longer than the timescales on which molecular clouds persist without disruption. 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 be extrapolated to the zero metallicity case. In practice, however, a very small value of metallicity should return reasonable results. Though the model is accurate when applied to a run that includes negligible but nonzero metallicity grid cells, we suggest that the model can most confidently be used for Z  10-3 Z e . These metallicities correspond to local dwarf galaxies and are roughly consistent with the gas out of which Population II stars form, making this low-metallicity model relevant for high-Z galaxies, as well.

Projected F H 2 Fit

Figure 5 shows good agreement between the projected molecular fraction, F 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 and 30 pc, both of which are fairly well populated by simulation grid cells with

- f 10 H 5 2 , even at the lowest metallicities. As with the volumetric fit, the imposed maximum fraction F H ,max 2 results in a hard upper boundary on ( ) F N H H 2

; see Section 4.1. Also, the scatter of F 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.5 across the full probed range of metallicities and averaging scales.

Discussion

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 Z e (see the Introduction).

Figure 9, however, shows that the depletion time of molecular hydrogen gas is change by nearly 3 orders of magnitude from 0.01 Z e to Z e in our simulations, while the SFR is changing only by a factor of 10 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 it effectively trades one error for another by artificially overestimating f H 2 in an environment where the fraction is actually low. A much better alternative 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 H I-H 2 transition models for the Z  0.2 Z e regime (corresponding to metallicities typical of dwarf galaxies and galaxies at high redshift), we 1)) 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 the results of the test shown in Figure 7 are not sensitive to the specific output used.

Figure 9. The depletion time (τ dep = M gas /SFR) as a function of metallicity and gas species relative to the SFR as a function of metallicity. We use M gas within 5 kpc of the simulation center. While t dep,H2 varies by a factor of 260 between Z = 0.01 and 1 Z e , the SFR (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 star formation in low-metallicity galaxies is not directly tied to H 2 abundance. We also include the inferred SFR and H 2 and H I τ dep in the LMC and SMC 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.

constructed fitting formulae for the molecular fraction as a function of gas density (or column density), local UV flux, and metallicity down to Z = 0.01 Z e . 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 Z up to the solar metallicity.

Here we compare the fitting formulae presented in this paper and a number of molecular fraction models in the literature (namely, KMT09b; GK11; K13; Gnedin 2014; S14) with simulation results. The models parameterize average f H 2 and F 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 the Appendix for details of the model implementations).

Figure 10 shows comparisons of the average f H 2 and F 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 Z in Figures 11 and 12. The very narrow distribution of f H 2 as a function of density in the KMT09b model is likely due to its neglect of dependence on U 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 H I-H 2 transition is better at high Z. This seems to be driven by their very steep functional form of f H 2 , so that only a few grid cells with high f H 2 contribute to the inferred mass at lower Z, while at higher metallicity, the transition n H is underestimated, leading to the model underprediction of M H 2 .

In the case of the projected models, KMT09b performs consistently well for Z  0.1 Z e . The shape of ( )

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 Z  0.4 Z e . GK11 is the next most accurate of the existing models. The lowest metallicity runs result in an inferred = M M 0 H 2  with this model, but this is expected behavior for the projected GK11 model, which is known to not be accurate for Z  0.01 Z e regard to the GD14 model, it is worth noting that this model includes a phenomenological account for the H 2 selfshielding due to line overlap. This results in near independence of f H 2 on dust abundance (and thus metallicity) for D  0.2D MW . The simulations used here do not include any accounting for such line overlap and thus a part of the difference between the GD14 model and simulation results at Z  0.2Z e 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 Z  0.2Z e (see Table 2).

The fit presented in this paper reproduces the total molecular hydrogen mass in the simulations considerably better than previous models. Even restricting K13 to Z  0.1 Z e , where their model results in masses within a factor of ∼2 of those measured in simulation, it appears to be a coincidence given the steep relation and few high f 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:

1. We show that the ISM in the simulated galaxy changes qualitatively when gas metallicity is varied by 2 orders of magnitude. The density distribution becomes increasingly nonuniform as metallicity increases from 0.01 Z e to Z e (see Figure 1). 2. 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 to values much less than 1 at Z  0.2 Z e (see Figures 3 and 5), while the dependent of molecular fraction on density becomes less steep. 3. We show that accurate fitting functions for volumetric and projected molecular fractions can be constructed if they account for the dependence on gas density, gas metallicity, and the strength of the ionizing UV field (Figures 4 and 6). We demonstrate that the presented fits reproduce the dependence of the H I-H 2 transition on metallicity and the overall shape of the molecular fraction-density distribution than existing models (Figures 10-12) 4. We also show that the volumetric (projected) molecular fraction fit is applied to individual cells (projected patches) and reproduces the total molecular mass in the simulated galaxies to a factor of 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 H I-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 SFR becomes a nonlinear function of molecular gas density at metallicities <0.1 Z e , due to nonequilibrium effects (see Section 5 and Figure 9). As an alternative, the star formation efficiency and depletion time can be calibrated using such simulations, and we will present such calibrations in a followup work. ) cells in the simulation. For our mixed metallicity run, we also use the median Z (∼0.5 Z e ). For K13 and S14, we assume a clumping factor of 1, and for K13, we also set ρ SD (the density of stars and dark matter) to ∼0.1 M e pc -3 , which is the median for the high molecular fraction cells in the simulation with a nonzero ρ SD . We underplot the median density in bins of f H2 in black for each run of the simulation. within the range shown here. As in Figure 11, given how few K13 model cells have - f 10 H 5 2 at densities consistent with those in the Z = 0.01 Z e simulation, we plot the phase space location of each cell individually for that model. Unlike in Figure 11, where we compute ρ SD on a per grid cell basis, we instead use a global value ∼0.1 M e pc -3 based on the median for each run.

are generally smaller than 100 pc. We note that this is not a perfect approximation since this expression for Σ comp should be averaged on ∼100 pc scales. In the volumetric case, to compute Σ H from n H , we take Σ H = (n H /cm -3 )(L cm -1 )(m p g -1 ), which we then convert to units of solar mass per square parsec. Similarly, in the K13 model, we compute Σ 0 in the same way as Σ H .

We largely follow the implementation of the K13 model. In order to compute f H 2 , we first need to compute n CNM,hydro , which relies on P th . P th is a function of R mol ( ) =f f 1 H H 2 2 . We ultimately iterate over this step 10 times, following the more efficient approach in Diemer et al. (2018), where f H 2 is initialized at 0.5 and is averaged as it advances through the iteration so that for step i, »+ - f f f 0.3 0.7

2 . We also adopt a clumping factor, f c = 1. For Figure 10, we adopt the median density ρ SD = 0.1 M e pc -2 in our simulation in regions where et al. (2018, Section C.4), we find no significant difference between using different constant values of ρ SD , but that using different values on a cell-by-cell basis (as we do in Figure 11, where ρ SD is computed directly from the simulation) has a significant effect.

A.3. S14 Model

We follow the simple Diemer et al. (2018, Section C.5) formulation for the S14 model (Bialy & Sternberg 2016). To convert the number density to column density, we assume Note. We assume the same model parameters as in Figures 11 and 12, respectively. For the projected models, we show the accuracy for the representative S = 100 pc case. We denote model and metallicity combinations that produce no molecular hydrogen with an ellipsis.

N H /cm -2 = (n H /cm -3 )(L cm -1 ) for comparison with our volumetric model. In the same vein, when comparing with our projected model, we take n H /cm -3 = (N H /cm -2 )(1.54 × 1021/cm) -1 , which corresponds to the approximate height of the disk (∼500 pc).