Tracing the total molecular gas in galaxies: [CII] and the CO-dark gas

While the CO(1-0) transition is often used to deduce the total molecular hydrogen in galaxies, it is challenging to detect in low metallicity galaxies, in spite of the star formation taking place. In contrast, the [CII] 158 micron line is relatively bright, highlighting a potentially important reservoir of H2 that is not traced by CO(1-0), but residing in the C+ - emitting regions. We explore a method to quantify the total H2 mass (MH2) in galaxies and learn what parameters control the CO-dark gas reservoir. We present Cloudy grids of density, radiation field and metallicity in terms of observed quantities, such as [OI], [CI], CO(1-0), [CII], total infrared luminosity and the total MH2 and provide recipes based on these models to derive total MH2 mass estimates from observations. The models are applied to the Herschel Dwarf Galaxy Survey, extracting the total MH2 for each galaxy which is compared to the H2 determined from the observed CO(1-0) line. While the H2 traced by CO(1-0) can be negligible, the [CII] 158 micron line can trace the total H2. 70% to 100% of the total H2 mass is not traced by CO(1-0) in the dwarf galaxies, but is well-traced by [CII] 158 micron line. The CO-dark gas mass fraction correlates with the observed L[CII]/LCO(1-0) ratio. A conversion factor for [CII] luminosity to total H2 and a new CO-to-total-MH2 conversion factor, as a function of metallicity, is presented. A recipe is provided to quantify the total mass of H2 in galaxies, taking into account the CO and [CII] observations. Accounting for this CO-dark H2 gas, we find that the star forming dwarf galaxies now fall on the Schmidt-Kennicutt relation. Their star-forming efficiency is rather normal, since the reservoir from which they form stars is now more massive when introducing the [CII] measures of the total H2, compared to the little amount of H2 in the CO-emitting region.


Introduction
Our usual view of star formation posits the molecular gas reservoir as a necessary ingredient, the most abundant molecule being H 2 . This concept is borne out through testimonies of the first stages of star formation associated with and within molecular clouds and numerous observational studies showing the observed correlation of star formation with indirect tracers of H 2 reflected in the Schmidt-Kennicutt relationship (e.g. Kennicutt 1998;Kennicutt et al. 2007;Bigiel et al. 2008;Leroy et al. 2008;Genzel et al. 2012;Kennicutt & Evans 2012;Leroy et al. 2013;Kumari et al. 2020). Witnessing the emission of H 2 associated with the bulk of the molecular clouds directly, however, is not feasible. Indeed, H 2 emits weakly in molecular clouds due to the lack of permanent dipole moment and the high temperatures 1 necessary to excite even the lowest rotational transitions (the 2 lowest transitions have upper level energies, hν/k, of 510 K and 1015 K). Thus, H 2 observations can directly trace only a relatively small budget (15% to 30%) of warmer (∼100K) molecular gas in galaxies (e.g. Roussel et al. 2007;Togi & Smith 2016), but not the larger reservoir of cold ( ∼10-20 K) molecular gas normally associated with star formation.
In spite of 4 orders of magnitude lower abundance of CO compared to H 2 , most studies rely on CO rotational transitions to quantify the properties of the H 2 reservoir in galaxies. The low excitation temperature of CO (1-0) (∼ 5 K), along with its low critical density (n crit ) for collisional excitation (∼ 10 3 cm −3 ), make it a relatively strong, easily-excited millimeter (mm) emission line to trace the cold gas in star-forming galaxies. Studies within our Galaxy have long ago established a convenient recipe to convert the observed CO (1-0) emission line to H 2 gas mass, that is, the X CO conversion factor (see Bolatto et al. 2013, for historical and theoretical development). Likewise, galaxies in our local universe have also routinely relied on CO observations to quantify the H 2 reservoir (e.g. Leroy et al. 2011;Bigiel et al. 2011;Schruba et al. 2012;Cormier et al. 2014Cormier et al. , 2016Saintonge et al. 2017). CO is also commonly used to probe the H 2 and star formation activity in the high-z (z∼ 1 -2) universe (e.g. Tacconi et al. 2010;Daddi et al. 2010;Combes et al. 2013;Carilli & Walter 2013;Walter et al. 2014;Daddi et al. 2015;Genzel et al. 2015;Kamenetzky et al. 2017;Tacconi et al. 2018;Pavesi et al. 2018;Tacconi et al. 2020), although for higher-z sources, higher J CO lines may become more readily available, but less straightforward to relate to the bulk of the H 2 mass (e.g. Papadopoulos et al. 2008;Gallerani et al. 2014;Kamenetzky et al. 2018;Vallini et al. 2018;Dessauges-Zavadsky et al. 2020). Well-known caveats can hamper accurate total H 2 gas mass determination in galaxies based on observed CO emission only, including those associated with assumptions on CO excitation properties, dynamical effects related to X CO , presence of ensembles of molecular clouds along the lines of sight in galaxies and low filling factor and abundance of CO compared to H 2 , as in low metallicity environments (e.g. Maloney & Black 1988;Glover & Mac Low 2011;Shetty et al. 2011;Narayanan et al. 2012;Pineda et al. 2014;Clark & Glover 2015;Bisbas et al. 2017).
CO, however, can fail altogether to trace the full H 2 reservoir, very importantly, in low-metallicity galaxies. Lower dust abundance allows the far ultra-violet (FUV) photons to permeate deeper into molecular clouds compared to more metal-rich clouds, photodissociating CO molecules, leaving a larger C +emitting envelope surrounding a small CO core. H 2 , on the other hand, photodissociates via absorption of Lyman-Werner band photons, which, for moderate A V , can become optically thick, allowing H 2 to become self-shielded from photodissociation (e.g. Gnedin & Draine 2014), leaving a potentially significant reservoir of H 2 existing outside of the CO-emitting core, within the C + -emitting envelope (Fig. 1). The presence of this CO-dark molecular gas (e.g. Röllig et al. 2006;Wolfire et al. 2010;Glover & Clark 2012b), requires other means to trace this molecular gas reservoir that is not probed by CO. The total H 2 mass, then, is the CO-dark gas mass plus the H 2 within the CO-emitting core. Efficient star formation in this CO-dark gas has been shown theoretically to be possible (e.g. Krumholz & Gnedin 2011;Glover & Clark 2012a).
More recently, observations of the submillimetre (submm) transitions of [C i] have been used to trace the total H 2 gas mass in galaxies at high and low-z, especially due to the advent of the spectroscopic capabilities of SPIRE (Griffin et al. 2010) on the Herschel Space Observatory (Pilbratt et al. 2010) and with ALMA (Popping et al. 2017;Andreani et al. 2018;Jiao et al. 2019;Nesvadba et al. 2019;Bourne et al. 2019;Valentino et al. 2020;Dessauges-Zavadsky et al. 2020), while theory and simulations (Papadopoulos et al. 2004;Offner et al. 2014;Tomassetti et al. 2014;Glover & Clark 2016;Bothwell et al. 2017;Li et al. 2018;Heintz & Watson 2020) as well, deem [C i] to be a viable tracer of CO-dark H 2 . Likewise, ALMA has opened the window to numerous detections of [C ii]λ158µm in the high-z universe, making [C ii], which is more luminous than [C i], a popular tracer of total H 2 in galaxies (e.g. Zanella et al. 2018;Bethermin et al. 2020;Schaerer et al. 2020;Tacconi et al. 2020). Dust mass estimates via monochromatic or full spectral energy distribution (SED) studies have also been used to quantify the total H 2 gas mass in the nearby and distant universe (e.g. da Cunha et al. 2013;Groves et al. 2015) (more background on the different approaches to uncover CO-dark gas is presented in the following Section 2).
Local and global environmental effects, including star formation properties and metallicity, require consideration as well, when choosing wisely the approach to relate total gas mass to star formation. Low metallicity environments present wellknown issues in pinning down the molecular gas mass. Detecting CO in dwarf galaxies has been notoriously difficult (e.g. Leroy et al. 2009;Schruba et al. 2012;Leroy et al. 2012;Cormier et al. 2014;Hunt et al. 2015;Grossi et al. 2016;Amorín et al. 2016, and references within), leaving us with uncertainties in quantifying the total molecular gas reservoir.
Star-forming dwarf galaxies often have super star clusters which have stellar surface densities greatly in excess of normal H ii regions and OB associations. Thus, the dearth of detectable CO raises questions concerning the fueling of such active star formation. Possible explanations may include: 1) CO is not a reliable tracer of the total molecular gas reservoir; 2) there is a real under-abundance of molecular gas indicative of an excessively high star formation efficiency; 3) H 2 has been mostly consumed by star formation; 4) H 2 has been completely dissociated in the aftermath of star formation; 5) HI is the important reservoir to host star formation.
These issues compound the difficulty in understanding the processes of both local and global star formation in lowmetallicity environments. Developing reliable calibrators for the total M H 2 in low Z galaxies is a necessary step to be able to refine recipes for converting gas into stars under early universe conditions and thus further our understanding of how galaxies evolve. Understanding the detailed physics and chemistry of the structure and properties of the gas that fuels star formation can help us discriminate strategies for accurate methods to get to the bulk of the molecular gas in galaxies.
Since carbon is one of the most abundant species in galaxies, its predominantly molecular form, CO, and its neutral and ionised forms, C 0 and C + , carry important roles in the cooling of the interstellar medium (ISM) of galaxies. The locations of the transitions from ionised carbon, to neutral carbon to CO, which are within the photodissociation region (PDR)-molecular cloud structure, depend on numerous local physical conditions, including the FUV radiation field strength (G 0 ) 1 , hydrogen number density (n H ), gas abundances, metallicity (Z), dust properties, etc. Low-metallicity star-forming environments seem to harbor unusual PDR properties compared to their dustier counterparts, with notably bright [C ii]λ158µm emission compared to the often faint CO (1-0) emission, characteristics that can be attributed Comparison of a solar metallicity (metal-rich) molecular cloud and a low-metallicity cloud impacted by the UV photons of nearby star clusters. The decrease in dust shielding, in the case of the low metallicity cloud, leads to further photodissociation of the molecular gas, leaving a layer of self-shielded H 2 outside of the small CO-emitting cores. This CO-dark H 2 can, in principle, be traced by [C ii] to the lower dust abundance along with the star formation activity and clumpy ISM (e.g. Cormier et al. 2014;Chevance et al. 2016;. In this study we determine a basic strategy to quantify the total H 2 mass in galaxies and estimate the CO-dark molecular gas mass existing outside the CO-emitting core, for a range of Z, n H and G 0 . We start with the models of Cormier et al. (2015Cormier et al. ( , 2019, which were generated using the spectral synthesis code, Cloudy (Ferland et al. 2017), to self-consistently model the mid-infrared (MIR) and far-infrared (FIR) fine structure cooling lines and the total infrared luminosity (L TIR ) from Spitzer and Herschel. This approach allows us to constrain the photoionization and PDR region properties, quantifying the density and incident radiation field on the illuminated face of the PDR, G 0 . These models exploit the fact that there is a physical continuity from the PDR envelope, where the molecular gas is predominantly COdark, into the CO-emitting molecular cloud. This dark gas reservoir can be primarily traced by the [C ii]λ158µm line while the L [C II] /L CO(1−0) can constrain the depth of the molecular cloud and the L [C II] /L TIR and L [C II] /L [O I] constrain the G 0 and density at the illuminated face of the PDR (discussed in Section 5.2). With this study, we determine a conversion factor to pin down the total H 2 , which will make use of the observed [C ii], [O i], CO (1-0) and L TIR . Application of the models to the starforming low-metallicity galaxies of the Herschel Dwarf Galaxy Survey (DGS; Madden et al. 2013), using the observed [C ii] and CO (1-0) observations, allows us to extract the total H 2 mass and determine the fraction of the total H 2 reservoir that is not traced by CO observations, as a function of metallicity and other galactic parameters.
This paper is organised as follows. In Section 2 we review the studies that have focused on uncovering the CO-dark reservoir in galaxies. The environments of the dwarf galaxies in the DGS are briefly summarised in Section 3. In Section 4 we bring in the observational context and motivation by comparing the L [C II] /L CO(1−0) and L [C II] /L FIR 2 of the full sample of the DGS to other more metal-rich galaxies from the literature. In Section 5 we present Cloudy model grids and discuss the effects of G 0 , n H and metallicity. In this section we walk through the steps to use the models to understand the mass budget of the carbon bearing species in the PDRs and to obtain the total H 2 . In Section 6 we apply the models to the particular case of II Zw 40 observations, as an example. Section 7 discusses the model results applied to the DGS targets. In this section we extract the total H 2 for these galaxies and discuss the implication of the CO-dark gas reservoirs. We provide recipes to estimate the CO-dark H 2 and the total H 2 gas reservoir from [C ii]. Possible caveats and limitations are noted in Section 8 and the results are summarised in Section 9.

CO dark gas studies
Three decades ago, the FIR spectroscopic window became readily accessible from the Kuiper Airborne Observatory (KAO), and the Infrared Space Observatory (ISO). First estimates of the CO-dark H 2 reservoir obtained from [C ii]λ158µm observations found that 10 to 100 times more H 2 was 'hidden' in the C + -emitting regions in dwarf galaxies compared to that inferred from CO (1-0) alone (Poglitsch et al. 1995;Madden et al. 1997). These early estimates, for lack of access to additional complementary FIR diagnostics to carry out detailed , relied on a number of assumptions which have, until now, limited the use of the FIR lines to robustly quantify this CO-dark molecular gas in galaxies. Recent spatially-resolved studies of [CII] in Local Group galaxies using Herschel and the Stratospheric Observatory for Infrared Astronomy (SOFIA), have also confirmed a significant presence of CO-dark gas (e.g. Fahrion et al. 2017;Requena-Torres et al. 2016;Jameson et al. 2018;Lebouteiller et al. 2019;Chevance et al. 2020). Additionally, the MIR rotational transitions of H 2 from Spitzer have been invoked to deduce a significant reservoir of H 2 in dwarf galaxies (Togi & Smith 2016). Other molecules observed in absorption, such as OH and HCO + , have also been used as tracers of dark molecular gas in our Galaxy (Lucas & Liszt 1996;Allen et al. 2015;Xu et al. 2016;Li et al. 2018;Nguyen et al. 2018;Liszt et al. 2019). However, the need for high enough sensitivity and resolution limits their use in other galaxies.
As γ rays interact with hydrogen nuclei, they can trace the total interstellar gas. Comparison of γ ray observations in our Galaxy with CO and H 0 has uncovered an important reservoir of neutral dark gas that is not traced by CO or H i (Grenier et al. 2005;Ackermann et al. 2012;Hayashi et al. 2019). While γ ray emission, in principle, can be one of the most accurate methods possible to get at the total gas mass, its use to measure the total interstellar gas mass in other galaxies, is limited for now, because of the need for relatively high resolution observations.
Comparison of Herschel [C ii] velocity profiles to H i emission (Langer et al. 2010;Pineda et al. 2013;Langer et al. 2014;Tang et al. 2016) as well as observations comparing dust, CO (1-0) and H i (Planck Collaboration et al. 2011;Reach et al. 2017), confirm that this CO-dark gas reservoir can be an important component in our Galaxy, and comparable to that traced by CO (1-0) alone. Recent hydrodynamical simulations with radiative transfer and analytic theory also demonstrate the presence of this dark gas reservoir in the Galaxy, that should be traceable via [C ii] (Smith et al. 2014;Offner et al. 2014;Glover & Clark 2016;Nordon & Sternberg 2016;Gong et al. 2018;Franeck et al. 2018;Seifried et al. 2020), often associated with spiral arms in disk galaxies. While optically-thick H i can, in principle, be the source of the CO-dark gas, there is no compelling evidence that it is a strong contribution to the dark neutral gas (e.g. Murray et al. 2018). Thus we focus here on the CO-dark H 2 gas, and the potential of C + , in particular, to quantify the H 2 reservoir.
Dust mass measurements with Planck, compared to the gas mass measurements, have uncovered a reservoir of dark gas in the Galaxy (Planck Collaboration et al. 2011;Reach et al. 2015). Dust measurements have often provided easier means, from an observational strategy, to quantify the full reservoir of gas mass of galaxies in the local universe as well as distant galaxies (e.g. Magdis et al. 2011;Leroy et al. 2011;Magnelli et al. 2012;Eales et al. 2012;Bourne et al. 2013;Sandstrom et al. 2013;Groves et al. 2015;Scoville et al. 2016;Liang et al. 2018;Bertemes et al. 2018). This approach also poses its own fundamental issues (e.g. Privon et al. 2018). The derived total molecular gas mass is a difference measure of two large quantities (H i-derived atomic gas mass and dust-derived total gas mass), which can result in large uncertainties. Dust mass determination, which is commonly derived via the modelling of the SED, requires constraints that include submm observations and necessitates dust models with assumed optical properties as a function of wavelength, composition and size distributions. Then, turning total dust mass into total gas mass requires assumptions on dust-togas mass ratios which can carry large uncertainties depending on star formation history and metallicity and can vary by orders of magnitude, as shown by statistical studies of low metallicity galaxies (e.g. Rémy-Ruyer et al. 2014;Roman-Duval et al. 2014;Galliano et al. 2018).
One contributing factor to the difficulty in being conclusive on the relationship of gas, dust and star formation in lowmetallicity conditions, in particular, is the fact that galaxies with full MIR to submm dust and gas modelling have been limited to mostly metal-rich environments and relatively higher star formation properties due to telescope sensitivity limitations. This impediment has been ameliorated with the broad wavelength coverage, the spatial resolution and sensitivity of the Herschel, allowing the accessibility of the gas and dust properties of lowmetallicity galaxies. The DGS has compiled a large observational data base of 48 low-metallicity galaxies , motivated by the Herschel PACS (Poglitsch et al. 2010) and Herschel SPIRE 55 to 500 µm photometry and spectroscopy capabilities.

Dwarf Galaxy Survey -extreme properties in low metallicity environments
The DGS targeted the most important FIR diagnostic tracers, (Cormier et al. , 2015 in a wide range of low-Z galaxies -as low as ∼1/50 Z ⊙ . Additionally, the DGS collected all of the FIR and submm photometry from Herschel to investigate the dust properties of star-forming dwarf galaxies. The infrared SEDs exhibit distinct characteristics setting them apart from higher metallicity galaxies with different dust properties, generally exhibiting overall warmer dust than normal metallicity galaxies; an obvious paucity of PAH molecules and a striking non-linear drop in dust-to-gas mass ratios at lower metallicities (Z < 0.1 Z ⊙ ) (e.g. Galametz et al. 2011;Dale et al. 2012;Rémy-Ruyer et al. 2013 [C ii] usually ranks foremost amongst the PDR cooling lines in normal metallicity galaxies, followed by the [O i]λ63µm, with the L [C II]+[O I] /L FIR often used as a proxy of the heating efficiency of the photoelectric effect (e.g. Croxall et al. 2012;Lebouteiller et al. 2012Lebouteiller et al. , 2019. Cormier et al. (2015) summarised the observed FIR fine structure lines in the DGS noting that the range of L [C II]+[O I] /L FIR of low-metallicity galaxies is higher than that for galaxies in other surveys of mostly metal-rich galaxies (e.g. Brauher et al. 2008;Croxall et al. 2012;van der Laan et al. 2015;Cigan et al. 2016;Smith et al. 2017;Lapham et al. 2017;Díaz-Santos et al. 2017), indicating relatively high photoelectric heating efficiency in the neutral gas. In star-forming dwarf galaxies, however, the [O iii]λ88µm line is the brightest FIR line, not the [C ii]λ158µm line, as noted on full-galaxy scales (e.g. Cormier et al. 2015Cormier et al. , 2019 as well as resolved scales (Chevance et al. 2016;Jameson et al. 2018;Polles et al. 2019). In more quiescent metal-poor galaxies, however, the [C ii]λ158µm line is brighter than the [O iii]λ88µm line (e.g. Cigan et al. 2016). The predominance of the large scale [O iii]λ88µm line emission, which requires an ionization energy of 35 eV, demonstrates the ease at which such hard photons can traverse the ISM on full galaxy scales in star-forming dwarf galaxies (Cormier et al. 2010(Cormier et al. , 2015, highlighting the different nature and the porous structure of the ISM of low metallicity galaxies. Since the [O iii]/[C ii] has been observed to also be similarly extreme in many high-z galaxies (e.g. Laporte et al. 2019;Hashimoto et al. 2019;Tamura et al. 2019;Harikane et al. 2020), the results presented here, where we focus on low-Z star-forming galaxies, may eventually be relevant to the understanding of the total molecular gas content and structure in some high-z galaxies.   Stacey et al. (2010) to include dwarf galaxies and more high redshift galaxies. The dwarf galaxies are from Cormier et al. (2010Cormier et al. ( , 2015 for the DGS; Grossi et al. (2016) for HeVICs dwarfs and Smith & Madden (1997). High redshift galaxies include those from Hailey-Dunsheath et al. (2010); Stacey et al. (2010); Gullberg et al. (2015). The black and blue symbols are data from the original figure of Stacey et al. (1991) for Galactic star-forming regions, starburst nuclei and non-starburst nuclei, ULIRGS (Luhman et al. 2003) and normal galaxies (Malhotra et al. 2001). The dashed lines are lines of constant L [C II] /L CO(1−0) . Note the location of the low metallicity dwarf galaxies (red squares) which show extreme observed L [C II] /L CO(1−0) values. wide variety of environments, ranging from starburst galaxies and Galactic star forming regions, to less active, more quiescent "normal" galaxies, low metallicity star forming dwarf galaxies and some high-metallicity galaxies as well as some high-z galaxies. This figure, initially presented by Stacey et al. (1991) and subsequently used as a star-formation activity diagnostic, indicates that most normal and star forming galaxies are observed to have L [C II] /L CO(1−0) ∼ 1000 to 4000, with the more active starburst galaxies showing approximately a factor of 3 higher range of L [C II] /L CO(1−0) than the more quiescent galaxies (Stacey et al. 1991Hailey-Dunsheath et al. 2010). The more active, dusty, star forming environments possess widespread PDRs exposed to intense FUV and the L [C II] /L CO(1−0) depends on the strength of the UV field and the shielding of the CO molecule against photodissociation due to H 2 and dust (e.g. Stacey et al. 1991Stacey et al. , 2010Wolfire et al. 2010;Accurso et al. 2017b).

L [C II] /L CO(1−0) in Galaxies
It was already shown that L [C II] was somewhat enhanced relative to L FIR and more remarkably so, relative to L CO in a few cases of star-forming low-metallicity galaxies compared to more metal-rich galaxies (Stacey et al. 1991;Poglitsch et al. 1995;Israel et al. 1996;Madden et al. 1997;Hunter et al. 2001;Madden et al. 2011;). In recent, larger studies, Cormier et al. (2015Cormier et al. ( , 2019 have carried out detailed modelling of the DGS galaxies and compared their overall physical properties to those of more metal-rich galaxies attributing the enhanced L [C II] /L FIR to the synergy of their decreased dust abundance and active star formation. The consequence is a high ionization parameter along with low densities resulting in a thick cloud and considerable geometric dilution of the UV radiation field. The effect over full galaxy scales is a low average ambient radiation field (G 0 ) and relatively normal PDR gas densities (often of the order of 10 3 to 10 4 cm −3 ), with relatively large C + layers. This scenario is also coherent with the L [C II] /L CO(1−0) observed in the low metallcity galaxies.
As can be seen in Fig. 2, star-forming dwarf galaxies can reach L [C II] /L CO(1−0) ratios an order of magnitude higher, or more (sometimes reaching 80,000) than their more metal-rich counterparts, on galaxy-wide scales. While CO (1-0) has been difficult to detect in star-forming low metallicity galaxies, thus shifting these galaxies well-off the Schmidt-Kennicutt relation (e.g. Cormier et al. 2014), [C ii], on the other hand, has been shown to be an excellent star formation tracer over a wide range of galactic environments, including the dwarf galaxies (e.g. de Looze et al. 2011;De Looze et al. 2014;Herrera-Camus et al. 2015). Such relatively high [C ii] luminosities in star-forming dwarf galaxies, harboring little CO (1-0), may be indicative of a reservoir of CO-dark gas, which is one of the motivations for this study.

Using spectral synthesis codes to characterise the ISM physical conditions
The strategy of this study is to first generate model grids that will help us explore how the properties of the CO-dark gas evolve as a function of local galaxy properties, such as Z, n H and G 0 , and then determine how observational parameters can be used to pin down the mass of the CO-dark gas. We will show how L [C II] /L CO(1−0) can constrain A V from the models and then the L [C II] /L TIR can narrow down the values of n H and G 0 to finally quantify the total mass of H 2 and hence, the mass of H 2 that is not traced by CO, the CO-dark gas.
5 S. C. Madden et al.: Tracing the total molecular gas in galaxies: [CII] and the CO-dark gas lines; abundance of C + , C 0 , and CO relative to C; hydrogen density (n H ) as a function of column density (N H ) and A V ; fraction of hydrogen in the ionised, atomic, and molecular form. The vertical dashed lines indicate the depth of the main phase transitions: H + to H 0 , H 0 to H 2 , CO(1-0) optical depth of 1. Masses indicated here should be adjusted for individual galaxies, scaling as L TIR (galaxy)/10 9 L ⊙ , as the source luminosity of the model is 10 9 L ⊙ (see Appendix A).

Model Parameters and Variations with Cloud Depth
We begin with the grids of Cormier et al. (2015Cormier et al. ( , 2019, who use the spectral synthesis code, Cloudy version 17.00 (Ferland et al. 2017), which simultaneously computes the chemical and thermal structure of H ii regions physically adjacent to PDRs. The central source of the spherical geometry of the Cloudy model is the radiation extracted from Starburst99 (Leitherer et al. 2010) for a continuous starburst of 7 Myr and a total luminosity 10 9 L ⊙ . The grids were computed by varying the initial density at the illuminated face of the H ii region, n H , and the distance from the source to the edge of the illuminated H ii region, the inner radius, r in . The ionization parameter (U) is deduced in the model, based on the input ionizing source, r in and n H . The models are calculated for 5 metallicity bins: Z=0.05, 0.1, 0.25, 0.5 and 1.0 Z ⊙ , 3 and n H ranging from 10 to 10 4 cm −3 . The r in values range from log(r in cm) = 20.0 to 21.3, in steps of 0.3 dex, which, for the various models, covers a range of U ∼ 1 to 10 −5 and G 0 ∼ 17 to 11481 in terms of the Habing radiation field. A density profile that is roughly constant in the H ii region and increases linearly with the total hydrogen column density (N H ) beyond 10 21 cm −2 is assumed . To ensure that all models go deep enough into the molecular phase, regardless of the metallicity, the stopping criterion of the models is set to a CO column density of 10 17.8 cm −2 (A V ∼ 10 mag) . With this criterion, the optical depth of the CO(1-0) line, τ CO , is greater than 1 in all models, which, by our definition, means that all models have transitioned in the molecular core 4 . For all calculations, the A V /N H ratio is computed self-consistently from the assumed grain-size distribution, grain types, and optical properties. For more description of the Cloudy models that generated the grids analysed here, see Appendix A. Figure 3 shows the evolution of the accumulated mass, abundances and intensities of the [C ii], [C i] and CO (1-0) as a function of metallicity, density, temperature and G 0 , from the H + region into the molecular cloud. As we want to capture the properties of the H 2 zone, we have defined the location of the hydrogen ionization front and the H 2 front in our calculations. As can be seen in the bottom subpanels in Fig. 3, the ionization front is defined where H exists in the form of H + and H 0 with a ratio of 50%/50%; the transition between H 0 and H 2 is, similarly, where each species is 50% of the total hydrogen abundance. Our model results are relatively insensitive to the transition definitions. The C + zone is defined to be where 95% of the total carbon abundance is in the form of C + and 5% in the form of C 0 and likewise the C 0 zone has 95% of the total carbon abundance in the form of C 0 (Fig. 3 second and third subpanels).
To understand the effect of Z, n H and G 0 on the zone boundaries and accumulated mass, we compare 3 different cases in Fig.  3: Z=0.1 Z ⊙ (top right panel) for a fixed n H (10 3 cm −3 ) and fixed G 0 (380); 2. Density effects: Compare n H = 10 2 and 10 4 cm −3 with a fixed Z (0.5 Z ⊙ ) and fixed G 0 (380) (middle left and right panels); 3. G 0 effects: Compare G 0 varying from 80 to 1800 with a fixed Z (0.5 Z ⊙ ) and fixed n H (10 3 cm −3 ) (bottom left and right panels).
We give here an example to interpret these plots: in the case of Z=1.0 Z ⊙ (top left panel), the transition from H + to H 0 occurs at A V ∼ 0.015 mag and the transition from H 0 to H 2 is at A V ∼ 0.5 mag, while the C + -emitting zone spans the range of A V ∼ 0.05 to 2 mag. The evolution of the accumulated mass of the hydrogen species (top left panel, top subpanel) indicates that the C + -emitting zone includes ∼ 90% of the total gas mass, in the form of both H 0 and H 2 , while CO is not yet formed at those A V values. The C 0 -emitting region begins at A V ∼ 2 mag until CO forms (in a region containing exclusively H 2 ). If we then lower the metallicity to 0.1 Z ⊙ (top right panel), we see that the evolution of the C + -emitting zone and the H 0 and H 2 transitions, in terms of A V , are the same as the 1.0 Z ⊙ case; however the N H , or cloud depth, scales differently, requiring a larger cloud depth at 0.1 Z ⊙ to reach the same A V as the 1.0 Z ⊙ cloud. The middle panels show how the cumulative [C ii], [C i], and CO (1-0) intensities scale when changing only the n H in the ionised zone from 10 2 cm −3 (middle left panel) to 10 4 cm −3 (middle right panel), while keeping G 0 and Z fixed. In this case, increasing the n H shifts the transitions from H + to H 0 and from H 0 to H 2 to lower values of A V . For the higher density case, we are quickly out of the H + region and into the atomic regime. Also, as can be seen from the center right panel, C + is emitting over much of the H 2 zone and before CO is emitted for the higher density case.
Finally, the bottom panels demonstrate the effect of changing G 0 and keeping Z and n H fixed. The increase of G 0 from 80 to 1800, for example, shifts the H + to H 0 transition to higher A V values (from A V ∼ 0.0015 to 0.02 mag) as well as a higher N H value for the H 0 to H 2 transition and with C + tracing a lower mass of H 2 .
A V /N H is computed from the dust-to-gas ratio, assumed to be scaled by metallicity in the model for the range of metallicities studied here. For an H + region, as long as dust does not significantly compete with hydrogen for ionizing photons, the N H of the H + region will remain constant. Reducing the metallicity does, however, reduce the A V corresponding to the N H in the H + region by the same factor as the metallicity, and hence, the ratio A V /N H . The size (in A V ) of the C + zone is roughly independent of metallicity (Kaufman et al. 2006), since the abundance of C + scales directly with metallicity, while the UV opacity scales inversely. For Z=0.1 Z ⊙ the carbon abundance is reduced by a factor of ten, but the decrease in dust extinction means that C + is formed over a path-length ten times longer than that for 1.0 Z ⊙ .

Hydrogen and carbon phase transitions in the model grids
Another way to visualise the model grids is shown in Fig. 4 with some observables, such as H i, The total H 2 gas mass from the model is shown in column 4 for the various metallicity bins. The C + -emitting region and some of the C 0 -emitting region, depending on the Z, n and G 0 (Fig. 3), will harbour H 2 sitting outside the CO-emitting region. Note that the models are scaled to L TIR of 10 9 L ⊙ . These diagrams can be used to bracket gas masses for a given L [C II] /L TIR and Z bin. To apply these models to observations of a specific galaxy, the masses must be scaled by the proper L TIR of the galaxy / 10 9 L ⊙ . More Cloudy model details, including application of the scaling factor, are described in Appendix A. An example of the application of the model grids to the galaxy, IIZw 40, is demonstrated in Fig.  7 and Section 6). A simple conversion of observed [C ii], [C i] or CO (1-0) to mass of total H 2 or fraction of CO-dark gas is not immediately straightforward.

Metallicity effects:
As the metallicity decreases (from bottom to top panels in Fig.4), the grids of gas mass shift to the right: more of the H 2 is associated with the C + and C 0 -emitting zones, demonstrating the overall larger CO-dark gas reservoirs predicted by the models in order to reach the molecular core. The location of both the H 2 front and the formation of CO depend on metallicity, with both zones forming at higher A V with decreasing Z (Wolfire et al. 2010). The location of the H 2 front as a function of A V depends only weakly on metallicity, scaling as ∼ ln(Z −0.75 ). This is caused by the H 2 formation rate scaling linearly with Z, while the UV dissociation rate of H 2 depends on both the self-shielding of H 2 (independent of Z) and dust extinction (which increases at a given depth with increasing Z). The CO forming zone scales as ∼ ln(Z −2 ), due to the decreased abundance of oxygen and carbon, which scales directly with Z. This has the overall effect of increasing the H 2 zone. Additionally, our calculations compute density as a function of column density. Decreasing Z increases the column density needed to obtain the same A V (also in Fig  3), therefore making the density in the H 2 zone higher at lower metallicity. These combined effects increase the mass present in the H 2 zone.

Distribution of the H 0 and H 2 phases :
The consequence of metallicity on the partition of the mass between three regions; the H 0 associated with [C ii] emission, the H 2 associated with [C ii] emission, and the H 2 associated with [C i] emission, can be seen in the first three columns of Fig.  4. The mass of H 0 associated with [C ii] increases significantly with decreasing metallicity, as expected. The beginning of the H 0 zone starts at the hydrogen ionization front which, for most of our parameter space, is dominated by hydrogen opacity since the dust does not significantly compete with the gas for ionizing photons. Therefore, the column density corresponding to the hydrogen ionization front is nearly independent of Z, but the A V of the ionization front will decrease with decreasing Z due to the lower A V /N H . The A V of the H 2 front increases slightly with decreasing Z as a result of decreased dust extinction. Overall, these two processes lead to a larger H 0 zone, and therefore increased mass with decreasing Z. The size of the H 2 zone, and the beginning of the CO-emitting region, follow similar logic.

G 0 effects:
When G 0 increases (for the same density and Z), the H 2 and CO fronts are pushed out to higher column densities, as expected, due to an increased dissociation rate. The width (in A V ) of the H 2 zone therefore shows little variation with G 0 . Increasing G 0 does have a mild effect on the mass associated with the H 2 zone, with increasing G 0 decreasing the mass in the H 2 zone (Fig.4). This is because increasing G 0 moves the H 2 zone out to a larger column density (Fig.3). Since the density in our models increases with increasing N H , the physical size of the H 2 zone will shrink. The thinner physical size leads to a smaller integrated mass, when compared to a lower G 0 calculation. This effect is small, however, when compared to variations of CO-dark mass with Z. Increasing density also slightly decreases the H 2 region mass, due to slight changes in the location of the H 2 front and the beginning of the CO formation zone.

Density effects:
In addition to the size of each region, the temperature, n H , and G 0 will dictate if the H 0 , H 2 , and C + regions will emit, and therefore trace the CO-dark gas. For [C ii], the emission is controlled by the C + column density, the n crit for [C ii] emission (3x10 3 cm −3 for collisions with H 0 ), and the excitation temperature (92 K) (Kaufman et al. 1999). We see that a density below the n crit in the H 0 region allows for efficient emission of [C ii] (Fig. 3). This is reflected in the fact that, for the lower density models shown in Fig. 4, the H i mass and the H 2 mass traced by [C ii] are often comparable. For densities beyond n crit , the emission of [C ii] is nearly independent of n H , and is therefore controlled by the temperature and column density. Since the density increases with column density, for almost all models except the lowest density of 10 cm −3 , the density in the H 2 zone eventually exceeds n crit of [C ii], while for the H 0 zone the column density is lower, leading to lower densities in this region and a larger region of the parameter space with a density lower than n crit of [C ii]. This explains why Fig. 4 has some spread in the mass for different densities, while for the H 2 , traced by [C ii], the plot is nearly constant for decreasing Z, except for the lowest densities considered in our calculations. For Z = 1 Z ⊙ , a bit more spread is seen since the column density to reach the stopping criterion, and hence the density increase, is smaller. A similar effect occurs for the mass traced by [C i], which has a n crit similar to [C ii], thereby causing almost all of the lower metallicity calculations to be independent of n H . Overall, the reservoir of M H 2 where the C + emission is originating, is systematically lower than that of the C 0 -emitting zone. This effect is caused by the fact that where C 0 is emitting, most of the hydrogen is in the form of H 2 while in the C + -emitting region, a significant amount of the hydrogen is in the form of H i, not necessarily in H 2 and the thickness of the respective layers are comparable. The plots of the 5th column in Fig. 4 show the [C ii] to [O i] ratio as a function of L [C II] /L TIR for the range of n H and G 0 . Kaufman et al. (1999) show that this ratio depends strongly on G 0 for low density, especially for the higher Z case. Then for increasing density the functional form of the ratio changes, becoming more sensitive to density, for densities greater than the n crit of [C ii]. For larger densities, the [O i] emission increases with density, while the [C ii] stays roughly constant (for a constant G 0 ), leading to a lower [C ii] to [O i] ratio. This explains the trend of this plot with decreasing Z, as the increased density in the PDR, given our density law, leads to models with an initially low density reaching the n crit of [C ii], which causes the predicted emission of [O i] to increase relative to [C ii]. Figure 4 allows us to quantify the amount of gas mass accumulated until the stopping criterion of the model is reached (ie log N(CO)=17.8; refer to section 5.1) as a function of important physical parameters, given the observed FIR spectrum. For example, in the metallicity bin, Z=1.0 Z ⊙ , given an observed L [C II] /L TIR of 0.5%, the total mass of H 2 gas ranges from 1 × 10 8 to 3 × 10 9 M ⊙ , depending on the density, while for the lowest metallicity bin shown, Z=0.05 Z ⊙ , the quantity of H 2 gas ranges from 1 to 3 × 10 10 M ⊙ for this same L [C II] /L TIR value. This shift to higher mass ranges of H 2 gas, suggests that the mass of CO-8 dark H 2 gas can be an important component in low Z galaxies, as already pointed out in previous works (e.g. Wolfire et al. 2010).
In summary, from Fig. 4, the observed L [C II] /L TIR enables us to determine a range of plausible G 0 and n H for a given metallicity. If G 0 and n H can be determined from other assumptions or observations (e.g. [C ii]/[O i] for density) then a tighter constraint for the total H 2 gas can be determined (section 5.4), eliminating some of the spread in n H and G 0 . For the lowest metallicity cases, the range of H 2 , determined by G 0 and n H , is relatively narrow, and less dependent on variations in G 0 or n H . Thus, even having only observations of [C ii] and L TIR , Fig. 4 may bring a usefully narrow range of solutions for the total H 2 for the lowest metallicity cases based on the definition of CO dark gas adopted here and used by Wolfire et al. (2010). Caution must be exercised when using L CO /M H 2 from these figures, since the M H 2 within the COemitting region, is very sensitive to the depth at which the models are stopped. If there is reason to trust that A V of ∼ 10 mag is an accurate representation of the CO-emitting zone, then the L CO /M H 2 from this figure would be justifiable. It may not be applicable, for some low Z cases, as discussed in Section 5.4.
Note how the profiles shift to the right, to higher H 2 masses at lower Z. For the higher Z case, density, in particular, plays an important role in determining the M H 2 conversion factors for these species. For low Z, the density variations become less sensitive to the M H 2 , as noted in the Z=0.1 Z ⊙ case (right panel of Fig. 5). For the case of L [C I] /M H 2 and L CO /M H 2 at lower metallicity, the density contours collapse together (Fig. 5, right panel), as is the case for the L [C II] , except for the highest n H case (n H = 10 4 cm −3 ; the green dots). In principle, [C i] should be able to quantify M H 2 , independent of density. From the right panel of [C i] appears to be a useful tracer to quantify the M H 2 , as also illustrated from hydrodynamical models (e.g. Glover & Clark 2016;Offner et al. 2014), and, as we show here at least for low Z cases, with little dependence on n H . The fact that L [C I] is much fainter than L [C II] (Fig. 5) makes it more difficult to use as a reliable tracer of M H 2 . These effects can also be seen from the grids in Fig.4.
We continue in the following sections to demonstrate the use of the observed [C ii] to determine M H 2 for specific cases and will followup on [C i] as a tracer of M H 2 in a subsequent publication.

How to determine A V and its effect on line emission. ("Spaghetti plots")
We see from Fig. 2 that L [C II] /L CO(1−0) exhibits great variations from galaxy to galaxy. When comparing the observed L [C II] /L CO(1−0) for low metallicity galaxies (ranging from ∼ 3 000 to 80 000), to the modeled L [C II] /L CO(1−0) , where the stopping criterion is A V ∼ 10 mag, we see that the models do not reach such high observed L [C II] /L CO(1−0) values. This is because when the models are stopped at A V ∼ 10 mag, too much CO has already formed for the low Z cases (Fig. 3), which require the models to stop at lower A V . The higher metallicity molecular clouds, exhibiting lower L [C II] /L CO(1−0) than the dwarf galaxies, could be better described by stopping at A V ∼ 10 mag. To explore the sensitivity of the emission as a function of depth, we can use the observed L [C II] /L CO(1−0) as a powerful constraint. Figure 6 shows The L [C II] /L CO(1−0) is a good tracer of A V for all G 0 and n H , with the range of A V becoming narrower for the spread of n H as G 0 increases, as seen in the top panels of Figure 6. The L [C I] /M H 2 begins to turn over, peaking at A V of ∼ a few mag, more or less where the [C ii] formation has decreased. This is why we see a relative flattening of the L [C I] /M H 2 for growing A V .
We stop the plots in the figures beyond where CO has formed and becomes optically thick. Otherwise this M H 2 will continue to accumulate causing the L CO /M H 2 to turn over and begin to decrease beyond the point of τ CO = 1. Having CO (1-0) observations in addition to the [C ii] observations, brings the best constraint on the A V of the cloud and thus a better quantification of the M H 2 .
Note that to construct Fig. 6, the models are stopped at log N(CO) = 17.8 (a maximum A V ∼10 mag) and the line intensities are extracted at the different depths into the cloud -at different cloud A V values. In principle, the emerging line intensities could be different depending on the stopping criterion, due to optical depth effects and cloud temperature structure. To quantified this effect, we ran grids stopping at several maximum A V values (2, 5 and 10, for example) and extracted line intensities, comparing the values we present in Fig. 6. The effects on the [O i] and [C ii] line intensities are mostly negligible (variations < 20% throughout the grid). The intensities of the grid stopping at a maximum of A V of 2 mag are generally similar to or larger than the line intensities extracted at the cloud depth of A V = 2 mag of a grid stopping at maximum A V of 10 mag, primarily due to optical depth effects on the grids run to A V of 10 mag versus A V of 2 mag. The CO (1-0) can see up to a factor of 2.5 variations (primarily at A V = 2 mag), depending on the stopping criteria and the extraction of line intensities in A V construction. Our comparison verification has assured us that we can move forward with our use of Fig. 6.
Later, when the models are applied to particular observations (Section 6), the depth of the models, i.e. the maximum A V , will be adapted for the specific metallicity case, to determine the range of associated masses in the different zones. More precise determination of the stopping criterion will also require some knowledge of the range of G 0 and n H (for example, L [C II] /L [O I] , in Fig. 4). Only then can we have the necessary ingredients to obtain the total M H 2 . Thus, determination of A V then gives us the total M H 2 .The mass of CO-dark H 2 is then the difference between the total M H 2 from the models and the H 2 determined from the observed CO (1-0) using a Galactic X CO .

Quantify the total M H 2 and CO-dark gas -how to use the models
We walk through the steps to constrain the total M H 2 and the CO-dark gas, depending on the availability of observations. We emphasise that to obtain the total molecular gas, including He, an additional factor of 1.36 (Asplund et al. 2009) should be taken into account. The best-constrained case is that for which [C ii], [O i], CO (1-0) and L FIR observations exist. We also illustrate the range of solutions that can be obtained with less observational constraints.
Since metallicity plays an important role in the models, knowledge or assumptions of Z is necessary. The parameters that define the applicable model are G 0 , n H , and the maximum A V . We outline the steps to use the models: These steps are illustrated for one galaxy, as an example, in Section 6 and Fig. 7. Having found the parameter combination(s) that reproduce the relative strength of these key emission lines and the dust continuum, we can scale the model to the absolute line strength using the panels above and obtain the mass of H 2 . Probably the most useful scaling is based on using the L [C II] /M H 2 predictions (4th row of each panel in Fig. 6) because the [C ii] line is strong and it originates throughout the parts of the model where the 5 We provide the values of these plots in table form at the CDS where the reader can select model values of Z, n H , inner radius (r in ), G 0 , A V and τ CO to obtain predicted M H 2 , [C ii], CO (1-0), [C i] and [O i] luminosities L [C II] /L CO(1−0) varies strongly. Note that our methodology using [C i] as an observational constraint, does not (directly) aid in better determining the most applicable models in the range of interest, i.e. those situations where the CO-emitting zone is reached. However, the conversion from [C i] line luminosity to molecular gas mass for the applicable models is tighter than for [C ii] (see also Fig. 5). Therefore [C i] observations will be very useful for getting the best possible measures of the total H 2 gas mass once the matching models have been determined and compared to observations. Once the total M H 2 is determined, the difference between this value and the M H 2 determined from CO (1-0) and the X CO conversion factor, will quantify the CO-dark gas reservoir.

Applying the models: the particular example of II Zw 40
Here, we use one galaxy from the DGS, II Zw 40, (Z=0.5 Z ⊙ ), to demonstrate how to determine the total M H 2 directly from the observations and the models presented above and thus the subsequent CO-dark gas mass. The relevant steps are visualised in Fig. 7  It is interesting to compare the derived H 2 gas mass with that determined using L CO and a standard (Milky Way) X CO conversion factor 6 . The standard conversion factor yields a value of 7.1 × 10 6 M ⊙ , i.e. only a few percent of the actual total H 2 gas mass determined from these models. We already know that as 6 In this study for the standard X CO conversion we use a factor of 2 × 10 20 cm −2 (K km s −1 ) −1 in terms of X CO (I CO ∝ N H2 ); 3.2 M ⊙ pc −2 (K km s −1 ) −1 in terms of α CO (I CO ∝ M H 2 ). While α CO is normally 4.3 M ⊙ pc −2 (K km s −1 ) −1 , including helium, here we do not include the helium mass (factor of 1.36) when comparing to the model output of M H 2 . Z decreases, X CO increases and calibrations for X CO based on metallicity vary wildly in the literature (e.g. Bolatto et al. 2013). Thus, comparing with a Galactic X CO really gives a lower limit on what is expected. In any case, even for the moderately low Z of 0.5 Z ⊙ , the CO-dark molecular gas will be important.

Estimating M H 2 with fewer observational constraints
Following the full example for II Zw 40 above, we consider: case a) when using [C ii], CO (1-0) and L TIR (no [O i] line observed or no other density indicator) and case b) only [C ii] and L TIR are available (CO (1-0) and [O i] are not available). In each case, we compare to the M H 2 derived from our fiducial model above (Section 6.1), where a more complete set of observations ([C ii], L TIR , [O i] and CO (1-0)) was available.
-case a: Using [C ii], L TIR and CO (1-0). The observed ratio of L [C II] /L TIR for this galaxy (Fig. 7) indicates combinations of log(n H ) and log(G 0 ), from (1;2.2) through (2.5;2.8) to (4.0;3.3). For these models A V values between 3.5 and 6 mag reproduce L [C II] /L CO(1−0) within its uncertainty. The best model, i.e. the closest predictions to the observed values, has an A V of 4.5 mag and contains 1.3 × 10 8 M ⊙ of H 2 . The H 2 gas mass in the models that satisfactorily reproduces the observations, ranges from 0.6 to 5.1 × 10 8 M ⊙ . Comparing with the range we find for the fiducial model (1.3−5.1 × 10 8 M ⊙ ), we can see that the allowed range is significantly larger: the upper value is the same, but now allows a lower end of the range. In the particular case of II Zw 40 the higher density models that contain less M H 2 gas before reaching the observed L [C II] /L CO(1−0) values can not be excluded and the range is thus expanded to lower values. -case b: Using [C ii] and L TIR as the only available observational constraints; CO (1-0) has not been observed or with limited sensitivity and [O i] or another density tracer, has not been observed. This means that we cannot exclude "normal" L [C II] /L CO(1−0) ratios and high optical depth. Without a constraint on the A V , these models can only be used to infer an upper limit on the M H 2 in the PDR by integrating the model until τ CO ≃1. The same combinations of n H and G 0 as for case a yield upper limits on H 2 gas masses of 0.7 to 13 × 10 8 M ⊙ . Removing the CO-bright H 2 gas from these models by using the predicted L CO with a Galactic conversion factor yields upper limits on the dark H 2 gas mass of 0.5 to 11 × 10 8 M ⊙ . Thus the derived upper limit of 13 × 10 8 M ⊙ is significantly above the value we derive (3.2±1.9 × 10 8 M ⊙ ) using the full set of constraints. Whether such an upper limit is useful will depend on the specific science question one is trying to address.

Quantify the total M H 2 and CO-dark gas in the Dwarf Galaxy Survey
To quantify the CO-dark gas of the DGS galaxies, we determine the predicted total M H 2 for each DGS galaxy from the models (M(H 2 ) total ). We then compare this predicted total M H 2 to the mass of the CO-bright H 2 , M(H 2 ) CO , using the observed L CO and Galactic X CO . The mass of the CO-dark gas component, M(H 2 ) dark , would then be the difference between the modelpredicted M(H 2 ) total reservoir and the observed M(H 2 ) CO : We note n H values ranging from 10 0.5 − 10 3 cm −3 and G 0 values of 10 2 − 10 3 that are determined from the Cormier et al. (2019) model solutions for the DGS galaxies, for metallicities ranging from near solar to ≈ 1/50 Z ⊙ 7 . We then extract the total M H 2 for each galaxy from the corresponding model grids, applying the steps described above, to quantify the total M H 2 consistent with the model solutions and, consequently, derive the CO-dark gas reservoir for each galaxy. Various galactic properties, observed and modeled parameters, and their relationships with the total H 2 mass and CO-dark gas reservoirs, are inspected (Figs. 8, 9 and 10).
7.1. Trends of M(H 2 ) total and CO-dark gas with model parameters We notice right away (Fig. 8) that for DGS galaxies, M(H 2 ) total is always much larger than that determined using CO only, M(H 2 ) CO . The increasing M(H 2 ) total / M(H 2 ) CO signals the effect of the CO-dark gas reservoir becoming an increasingly important component of the total M H 2 , particularly in low Z environments, which has been noted in the literature (e.g. Poglitsch et al. 1995;Israel et al. 1996;Madden et al. 1997;Wolfire et al. 2010;Fahrion et al. 2017;Nordon & Sternberg 2016;Accurso et al. 2017a;Jameson et al. 2018;Lebouteiller et al. 2019). The total mass of H 2 can range from 5 to a few hundred times the H 2 determined from the CO-emitting phase. The CO-dark gas dominates the total H 2 reservoir in these galaxies. We are clearly missing the bulk of the H 2 by observing only CO. What is controlling the fraction of CO-dark gas?
The G 0 and density determined from the models do not seem to play an obvious role in driving the CO-dark gas fraction as shown in Fig. 8 panels a and b. Panel c shows the tight anticorrelation of M(H 2 ) total / M(H 2 ) CO with A V , and, as shown in panel d, A V anti-correlates with L [C II] /L CO(1−0) , underscoring the role of A V in regulating the C + -CO phase transition in the PDR (e.g. Wolfire et al. 2010;Nordon & Sternberg 2016;Jameson et al. 2018). The extreme range of observed L [C II] /L CO(1−0) in low metallicity galaxies seen in Fig. 2 is a consequence of their overall low average effective A V .

Trends of M(H 2 ) total and CO-dark gas with observables
What relationships exist between the observed quantities or measured galaxy properties and the total M H 2 and the quantity of CO-dark gas? In Fig. 9 (panel a) we see that the observed L [C II] /L CO(1−0) is an excellent tracer of the CO-dark M H 2 fraction. We fit this correlation to convert from observed L [C II] /L CO(1−0) to the M(H 2 ) total /M(H 2 ) CO : The standard deviation of this fit is 0.25 dex. It follows, from eq. 2 and eq. 3, that the ratio of the mass of CO-dark gas to CObright gas is, therefore: M(H 2 ) dark / M(H 2 ) CO = 10 −3.14 × [L [C II] /L CO(1−0) ] 1.09 -1.0 (4) We find a very tight correlation between L [C II] and total H 2 gas mass (Fig. 9, panel b) with a standard deviation of 0.14 dex. The C + luminosity alone can pin down the total M H 2 , making the [C ii]λ158µm a valuable tracer to quantify the molecular gas mass in galaxies. Our modelling results find about a factor of 3 larger mass of total H 2 associated with L [C II] , and hence a larger reservoir of COdark gas mass, than that determined empirically from Zanella et al. (2018). Zanella et al. (2018) determine the M H 2 from an assumed CO-to-M H 2 conversion factor which is lower then that found in this study. We find a trend of the CO-dark M H 2 fraction growing as the metallicity decreases ( Fig. 9 panel c).
Our study of the lowest-metallicity objects is, however, limited by the lack of robust CO detections in these extreme environments, thus limiting our knowledge of the behavior of Z with the CO-dark gas mass or total M H 2 at the lowest Z end. The COdark gas mass fraction does climb steeply as the Z decreases, even for moderately low-Z galaxies.
We see a weak trend of increasing fraction of CO-dark gas mass as the observed L [C II] /L TIR increases ( Fig. 9 panel d). For the relatively narrow range of L [C II] /L TIR there is a wide range of CO-dark gas fraction. The fraction of CO-dark gas, which covers almost 2 orders of magnitude, does not show a trend with the ∼ 1 order of magnitude range of surface density of SFR (Σ SFR ) in our galaxy sample ( Fig. 9 panel e). This is consistent with the effect of increasing L [C II] /L TIR (larger fraction of CO dark gas) being correlated with the increasing L [C II] (Fig. 9 panel d) and less so with direct effects of L TIR .

Consequences of the CO-dark gas fraction on the Schmidt-Kennicutt relation and the X CO conversion factor
What is the consequence of the presence of this reservoir of COdark H 2 on Σ SFR and the surface density of gas (Σ gas ) in galaxies, as described in the relationships of Kennicutt (1998) and Bigiel et al. (2008)? In Fig. 10 (panel a), we determine the Σ M H 2 within the CO-emitting region (black squares) and find their positions well-off the Σ M H 2 -Σ SFR relationships, as found from Cormier et al. (2014), which may be suggesting much higher Σ SFR for their M H 2 . Once we take into account the total M H 2 determined from [C ii] and the modelling, which now includes the CO-dark M H 2 as well as the CO-bright M H 2 (Fig. 10, panel a; red dots), we see the locations of the galaxies shift to the right, lying between both Σ M H 2 -Σ SFR relationships. The CO-dark gas is an important component to take into account in understanding the star formation activity in dwarf galaxies. While the data are limited, we find that taking into account the total M H 2 , the star-forming dwarf galaxies do show a similar relation as the star-forming disk galaxies. Thus they are not necessarily more efficient in forming stars. It has been shown from simulations (e.g. Glover & Clark 2012b;) that star formation can proceed without a CO prerequisite as well as without H 2 . The relationship observed between H 2 and star formation can be due to the cloud conditions providing the ability to shield themselves from the UV radiation field, thereby, allowing them to cool and to form stars. In this process of shielding, at least H 2 formation can also proceed, as well as CO in well-shielded environments. The conditions required to reach the necessary low gas temperatures set the stage for efficient C + cooling. CO may accompany star formation but does not have a causality effect. Thus it may not come as a surprise to see star-forming dwarf galaxies, harboring a dearth of CO, showing a similar relationship as that of the more metal-rich disk galaxies seen in the Schmidt-Kennicutt relationship.
With our determination of total H 2 we can now give an analytic expression to convert from observed CO to a total M H 2 conversion factor, α CO , and its relationship with Z. Here, again, the total M H 2 now includes the CO-dark gas plus the CO-bright H 2 mass. We find from Fig. 10 (panel b): with a standard deviation of 0.32 dex. We find a steeply rising α CO as the metallicity decreases. For example, at Z= 0.2 Z ⊙ , α CO is about 2 orders of magnitude greater than that for the Galaxy. A strong dependence of the CO-to-H 2 conversion factor on metallicity is not unexpected. Since the survival of molecules depends on how unsuccessful UV photons are in penetrating molecular clouds and photodissociating the molecules, extinction plays an important role in this process. Thus, the lower dust abundance in the low metallicity cases, comes into play. In Fig. 10 (panel b) we show comparison of our derived metallicity dependence of α CO with others in the literature. Schruba et al. (2012) have determined α CO from the observed SFR scaled by the observed L CO and a constant depletion time. This approach assumes that the efficiency of conversion of H 2 into stars is constant within different environments. Our determination of α CO and that of Schruba et al. (2012) 8 , are quite comparable, given the spread of the observations, low number statistics, as well as possible uncertainties in metallicity calibrations. The α CO from our study seems to climb even more steeply toward lower Z. However, the lower end of the metallicity space, where mostly only upper limits in CO (1-0) observations exist, is not pinned down robustly. It is well shifted up from the Glover & Mac Low (2011) relationship with CO-to-M H 2 conversion factor, which is determined from hydrodynamical simulations. A shallower slope is found for star-forming low metallicity galaxies from Amorín et al. (2016) who derive a metallicity-dependent α CO considering the empirical correlations of SFR, CO depletion time scale and metallicity. Other α CO -Z scaling functions, such as Arimoto et al. (1996) and Wolfire et al. (2010) fall near that of Amorín et al. (2016). By studying how the molecular gas depletion times vary with red-shift and their relation to the star-forming main sequence, Genzel et al. (2015) have determined a scaling of α CO taking into account CO and dust-based observations (Fig. 10, panel b). While this conversion factor is shallower to that we find from our study of local low-Z dwarf galaxies, they note that their study, which is based on massive star-forming galaxies, is probably not reliable for Z < 0.5 Z ⊙ . Accurso et al. (2017a) determine a comparable scaling of α CO with Z as Genzel et al. (2015) using similar surveys, but includes a second order dependence on distance from the star-forming main sequence in their α CO . They caution against using this relation below Z∼ 0.1 Z ⊙ , where it has not been calibrated. Bolatto et al. (2013) take a thorough look at numerous observations and theory and different metallicity regimes and propose a Z-dependent α CO which is similar to Genzel et al. (2015) for near-solar Z galaxies, but curves to steeper α CO for lower Z cases, while Tacconi et al. (2018) propose a compromise between Bolatto et al. (2013) and Genzel et al. (2015) (Fig.  10, panel b).
As soon as the ISM of the galaxy is more metal-poor, given the observed CO which, even for moderately metal-poor galaxies, is difficult to obtain, the conversion factor from CO-to-H 2 quickly grows. Even at 20% Z ⊙ , the CO conversion factor is already 1000 times that of our Galaxy. The true reservoir of M H 2 may have been severely underestimated so far at low Z and these new relations can quantify that.

Possible caveats and limitations
While [C ii] emission can be a convenient tool to quantify a reservoir of molecular gas that is not traced by CO (1-0), there are some caveats and limitations to this study.
-Lower metallicity bound: We emphasize that these relationships to determine total gas mass presented here have only been studied for the star-forming dwarf galaxies of the DGS with metallicities as low as ≈ 1/50 Z ⊙ . The models have been applied for the low metallicity galaxies for which CO is observed and this has limited the derived α CO and L [C II]to-M H 2 conversion factor only to metallicities as low as Z ∼ 0.2 Z ⊙ , even though [C ii] has been detected in DGS galaxies below this metallicity. -Limited range of model parameters: The models have been applied over a range of n H = 10 to 10 4 cm −3 and over a range of inner radii, corresponding to log G 0 ∼ 1 to 4. These results need to be studied over broader ranges of galactic properties to be applied with confidence beyond this study. For example, faint, low metallicity, more quiescent dwarf galaxies can also habor star formation rates lower than those found for the Schmidt-Kennicutt relation (e.g. Roychowdhury et al. 2009;Cigan et al. 2016), in contrast to the DGS sample ( Fig.  10 panel a). Some of these more quiescent dwarf galaxies are at the lowest metallicity range explored here and do not have CO detections, but may also be harboring some COdark molecular gas. Likewise, more massive, CO-rich galax- ies have yet to be tested with this model. -Galactic size scales: These conclusions have been drawn on unresolved galaxy scales. The reliability of [C ii] to trace M H 2 on resolved scales using these models has not been tested. A similar approach was carried out at 10 pc resolution in the 30Doradus region of the low metallicity LMC (Z = 0.5 Z ⊙ ), where > 75% of the H 2 was CO-dark and traced by [C ii] (Chevance et al. 2020). Considering that the [C ii] emission is expected to be more extended than the CO (1-0) (e.g. Jameson et al. 2018;Chevance et al. 2020), our simple scheme using Cloudy and representing galaxies with a single H ii region + molecular cloud in one dimension, does not account for realistic geometry and for the physical mixing of different physical components and thus, mixing G 0 and n H conditions. While this study points to the usefulness of [C ii] as a tracer of the molecular gas in low-Z galaxies, enlarging the galactic parameter space will be needed to apply these findings in a more general sense to a wide variety of galaxies, especially galaxies more massive than the dwarf galaxies in this study. Here we only focus on galaxies of the local universe. However, this is already a potential step forward in the possibility of estimating the gas mass from high-z galaxies using ALMA to access the [C ii]λ158µm, in the cases where CO (1-0) may be faint, perhaps due to lower metallicity.

Summary and Conclusions
This study is motivated by the extreme L [C II] /L CO(1−0) values observed for low metallicity galaxies -almost reaching 10 5 , on global scales, which can be up to an order of magnitude higher than dustier star-forming galaxies. The bright [C ii]λ158µm lines observed in low metallicity galaxies have been challenging to reconcile with the faint or undetectable CO (1-0) in light of their star formation activity. We have investigated the effects of metallicity, gas density and radiation field on the total molecular gas reservoir in galaxies and quantified the mass of H 2 not traced by CO (1-0) -the CO-dark molecular gas which can be traced by [C ii]. Cloudy grids traversing these physical parameters are inspected to understand the behavior of observed quantities, such as [O i], [C i], CO (1-0), [C ii] and L FIR in terms of M H 2 , as a function of A V , metallicity and n H . In principle, [C i] can be an important tracer of the CO-dark molecular gas, and this will be further investigated in a follow-up study. However, due to its higher luminosity, [C ii] is an ideal tracer of the molecular gas. We give recipes on how to use these models to go from observations to total M H 2 . We apply the models to the Herschel Dwarf Galaxy Survey, extracting the total M H 2 for each galaxy. The CO-dark M H 2 is then determined from the difference between the total M H 2 and the CO-bright H 2 , traced by the observed CO (1-0). Our findings indicate that CO (1-0) in the dwarf galaxies traces only a small fraction of the H 2 , if any, while, the total M H 2 is dominated by the CO-dark gas which can be uncovered by [C ii] observations.
We have determined a L [C II] -to-M H 2 conversion factor: M(H 2 ) total = 10 2.12 × [L [C II] ] 0.97 with a standard deviation of 0.14 dex. Following from this, we give a new CO-to-M H 2 conversion factor that takes into account the total M H 2 -both the CO-dark and CO-bright gas -given by application of the models: α CO = 10 0.58 × [Z/Z ⊙ ] −3.39 with a standard deviation of 0.32 dex.
Comparisons with the fraction of CO-dark gas in the low metallicity galaxies and their galactic properties reveal these findings: 1. The fraction of CO-dark gas is correlated with L [C II] /L CO(1−0) . There is a tight correlation between the [C ii]λ158µm and the total M H 2 over the range of low metallicity galaxies of the DGS. 2. The effective A V from our models is anticorrelated with L [C II] /L CO(1−0) and hence the CO-dark gas fraction. Hence the consequence of the effective low A V overall, in low metallicity galaxies is the extreme L [C II] /L CO(1−0) observed in star-forming dwarf galaxies. 3. The SFR, n H and G 0 do not individually control the CO-dark gas mass fraction. 4. The CO-dark gas accounts for most (> 70%) of the total H 2 , over the wide range of galaxy properties of the Herschel Dwarf Galaxy Survey. This study consists of galaxies with L [C II] /L TIR ranging between 0.1 to 0.5%, metallicity values from ∼ 0.02 to 0.6 Z ⊙ , L [C II] /L CO(1−0) values corresponding to 2 orders of magnitude of M H 2 and over 3 orders of magnitude of L FIR and stellar mass. 5. Our new CO (1-0)-to-M H 2 conversion factor as a function of metallicity, which also accounts for the CO-dark M H 2 , is steeper than most others in the literature, but resembles that of Schruba et al. (2012). 6. When taking into account this significant CO-dark H 2 reservoir, the star-forming dwarf galaxies, which were above the Schmidt-Kennicutt relation when the observed CO (1-0) only was used to deduce the H 2 reservoir, now shift to the normal relation of Σ SFR -Σ M H 2 found for disk, star-forming galaxies. 7. [C i] shows signs of being a convenient tracer of the M H 2 , particularly at low Z, where it can be independent of n H and G 0 . However, [C i] is intrinsically less luminous than [C ii]. The utility of [C i] as a tracer of total M H 2 for a variety of galactic conditions, requires further follow-up investigation.
We conclude that the low L CO /SFR and the high L [C II] /L CO(1−0) observed in low metallicity dwarf galaxies can be explained by the photodissociation of CO, and signals the presence of a prominent reservoir of CO-dark H 2 . Observations of a larger number of extremely low-Z galaxies is necessary to pin down the lowest Z end of the α CO -Z and the L [C II] -M H 2 relation as well as expanding the studies of parameter space and spatial scales in galaxies. Likewise, robust comparison of other methods to determine M H 2 and when these different methods are applicable, will be the subject of followup studies.
Acknowledgements. The authors wish to thank the anonymous referee for a number of suggestions that have improved the presentation of this study. We thank A. Poglitsch for valuable discussions that also helped to improve the quality of the paper. We acknowledge support from the DAAD/PROCOPE projects 57210883/35265PE and from the Programme National Physique et ionised region (at the start of the H ii region), r in , which, for these models is varied from log(r in cm) = 20.0 to 21.3, in steps of 0.3 dex. The initial hydrogen density (n H ) at the illuminated edge of the cloud is varied from 10 to 10 4 cm −3 . While the initial ionization parameter, U, can be an input command only in the intensity case, for our luminosity case, we chose to effectively have a range of U in the luminosity case, by varying r in as an input parameter, which in turn varies U, which is deduced within Cloudy for each model: where Q(H) is the number of hydrogen-ionizing photons emitted by the central source, and c is the speed of light. The intensity of the FUV radiation at the PDR front, G 0 , deduced in our Cloudy models, ranges from ∼ 17 to 11481 in units of the Habing radiation field (1.6 × 10 −3 ergs cm −2 s −1 ) (Habing 1968). Five metallicity bins were calculated for the models: Z=0.05, 0.1, 0.25, 0.5 and 1.0 Z ⊙ . The dust and PAH properties used in this model are described in Cormier et al. (2019). The opacity curves of the SMC are used. The abundance of PAHs is further reduced by metallicity to the power of 1.3, characterizing the prominent drop in PAH abundance at lower metallicity (e.g. Rémy-Ruyer et al. 2014. Tests performed in Cormier et al. (2019), inspecting the sensitivity of the PAH abundance, concluded that the PAH abundance is not an important factor in the outcome of the model results since the grain abundance is always larger than the alreadyreduced PAH abundance.
In Cloudy it is possible to choose a constant density throughout the cloud or to assume that the total pressure is constant. Cormier et al. (2019) have looked into different approaches of the density law, for a constant pressure case versus a constant density case and have found a compromise of these 2 extreme cases that works best in predicting the observations for the DGS sources. An intermediate case for the assumed density profile is constructed which is constant in the H ii region and increases linearly with the hydrogen column density in the neutral gas. The different cases of constant density, constant pressure and smoothly-increasing density law were tested in Cormier et al. (2019) where, for example, the constant density models predict less [O i] emission and the constant pressure models predict more [O i] emission than the smoothly-increasing density law adopted, which was found to successfully describe the observations.
A full galaxy has numerous stellar clusters and H ii regions and ensembles of PDR/molecular clouds. To go from a model single cluster-plus-cloud system to a representation of a full galaxy with numerous cluster-plus-cloud systems, we create a "unit model" with total source luminosity = 10 9 L ⊙ . The "unit model" of 10 9 L ⊙ chosen here is arbitrary but serves as a representative order of magnitude of luminosity of the observed galaxies. Thus all of the line luminosities predicted by the model correspond to a model for which L = 10 9 L ⊙ . The output masses are also determined for a model with source L = 10 9 L ⊙ . Since mass scales with luminosity, the masses determined for an individual galaxy, having L TIR (galaxy), are scaled by L TIR (galaxy) / 10 9 L ⊙ . This approximation assumes an energy balance between UV-optical and infrared and that L TIR scales similarly as the luminosity in the hydrogen-ionizing energy range .
We refer to Cormier et al. (2019) for further details and input parameter studies of a similar model used for this study.
16 S. C. Madden et al.: Tracing the total molecular gas in galaxies: [CII] and the CO-dark gas behavior (last column) in terms of metallicities of Z=0.05, 0.1, 1.0 Z ⊙ (from top to bottom) and for a range of G 0 and n H , the initial hydrogen density. The colour coding in each figure is log n H /cm −3 which increases from 1.0 (green) to 4.0 (purple). A range of G 0 is set by varying log r in /cm from 21.3 at the top right of the grids (large dots) to 20.0 at the bottom left of the grids (smaller dots). These r in values cover a range of log G 0 of 1.25 to 4.06 (Section 5). The Cloudy models are run with a source luminosity of 10 9 L ⊙ . Thus, the output masses should be scaled likewise (see Appendix A for details). The models are run to log N(CO)/cm −2 = 17.8 (A V ∼ 10 mag) for these grids. (1-0) (609 µm) and L CO -to M H 2 conversion factors for Z=1.0 Z ⊙ (left) and Z=0.1 Z ⊙ (right) for the particular model using a cloud depth of log N(CO)/cm −2 = 17.8. The colour coding refers to density values ranging from log n H /cm −3 = 1 (green) to 4 (purple). G 0 increases with increasing symbol size: 5, 20, 100, 200 and 500. To obtain the total molecular gas mass, which would include He, an additional factor of 1.36 (Asplund et al. 2009) should also be included for the total mass. Note that M H 2 in the CO-emitting region is very sensitive to the depth at which the models are stopped which corresponds to A V ∼10 mag here.