SILCC-Zoom: H$_2$ and CO-dark gas in molecular clouds -- The impact of feedback and magnetic fields

TL;DR

This study uses high-resolution simulations to analyze the distribution of CO-dark molecular gas in molecular clouds, highlighting the limitations of projected observations and proposing a new method to estimate H$_2$ mass more accurately.

Contribution

It provides new insights into the properties and distribution of CO-dark gas and introduces a novel approach to determine H$_2$ mass using combined CO emission and visual extinction maps.

Findings

CO-dark gas fractions range from 40% to 95%.

Projected visual extinction can misrepresent local gas conditions.

Proposed method reduces H$_2$ mass uncertainty to within a factor of 1.8.

Abstract

We analyse the CO-dark molecular gas content of simulated molecular clouds from the SILCC-Zoom project. The simulations reach a resolution of 0.1 pc and include H$_2$ and CO formation, radiative stellar feedback and magnetic fields. CO-dark gas is found in regions with local visual extinctions $A_\text{V, 3D} \sim$ 0.2 - 1.5, number densities of 10 - 10$^{3}$ cm$^{-3}$ and gas temperatures of few 10 K - 100 K. CO-bright gas is found at number densities above 300 cm$^{-3}$ and temperatures below 50 K. The CO-dark gas fractions range from 40% to 95% and scale inversely with the amount of well-shielded gas ($A_\text{V, 3D}$ $\gtrsim$ 1.5), which is smaller in magnetised molecular clouds. We show that the density, chemical abundances and $A_\text{V, 3D}$ along a given line-of-sight cannot be properly determined from projected quantities. As an example, pixels with a projected visual extinction of $A_\text{V, 2D} \simeq$ 2.5 - 5 can be both, CO-bright or CO-dark, which can be attributed to the presence or absence of strong density enhancements along the line-of-sight. By producing synthetic CO(1-0) emission maps of the simulations with RADMC-3D, we show that about 15 - 65\% of the H$_2$ is in regions with intensities below the detection limit. Our clouds have $X_\text{CO}$-factors around 1.5 $\times$ 10$^{20}$ cm$^{-2}$ (K km s$^{-1}$)$^{-1}$ with a spread of up to a factor $\sim$ 4, implying a similar uncertainty in the derived total H$_2$ masses and even worse for individual pixels. Based on our results, we suggest a new approach to determine the H$_2$ mass, which relies on the availability of CO(1-0) emission and $A_\text{V, 2D}$ maps. It reduces the uncertainty of the clouds' overall H$_2$ mass to a factor of $\lesssim$ 1.8 and for individual pixels, i.e. on sub-pc scales, to a factor of $\lesssim$ 3.

Tables

Table 1. Table 1: Overview of the simulations giving the run name, the starting time of the zoom-in procedure t 0 subscript 𝑡 0 t_{0} , the time t end subscript 𝑡 end t_{\mathrm{end}} up to which the clouds are evolved, whether magnetic fields (B, no B) or radiative feedback (FB, no FB) are included, and the underlying reference.

run	$t_0$ (Myr)	$t_{\mathrm{end}}$ (Myr)	run type	Ref.
MC1-HD	11.9	$t_0$ + 4.0	no B, no FB	(1)
MC2-HD	11.9	$t_0$ + 4.0	no B, no FB	(1)
MC1-HD-FB	11.9	$t_0$ + 4.0	no B, FB	(2)
MC2-HD-FB	11.9	$t_0$ + 4.0	no B, FB	(2)
MC3-MHD	16.0	$t_0$ + 5.5	B, no FB	(3)
MC4-MHD	16.0	$t_0$ + 5.5	B, no FB	(3)

Equations

$$A_{\mathrm{V},\mathrm{i}}=(N_{\mathrm{H},\mathrm{tot},\mathrm{i}}\times 5.348\times 10^{-22}{\mathrm{cm}}^2)\mathrm{mag}.$$

$$A_{\mathrm{V},3D}=\frac{-1}{\gamma }\text{ln}\left(\frac{1}{48}\sum _{i=1}^{48}\text{exp}(-\gamma A_{\mathrm{V},\mathrm{i}})\right),$$

$$B_x=B_{x,0}\sqrt{\rho (z)/{\rho }_0},$$

$$\text{DGF}=1-\frac{\frac{M_{\mathrm{CO}}}{28m_{\mathrm{p}}}\times \frac{1}{2.8\times 10^{-4}}}{\frac{M_{{\mathrm{H}}_2}}{2m_{\mathrm{p}}}},$$

$$\mathrm{PDF}(x)=\frac{\mathrm{d}m}{\mathrm{d}x}\frac{1}{M_{\mathrm{tot}}}.$$

$$f_{{\mathrm{H}}_2}=\frac{m_{{\mathrm{H}}_2,\mathrm{cell}}}{m_{\mathrm{H},\mathrm{cell},\mathrm{tot}}}\text{and}{f}_{\mathrm{CO}}=\frac{\frac{12}{28}{m}_{\mathrm{CO},\mathrm{cell}}}{m_{\mathrm{C},\mathrm{cell},\mathrm{tot}}},$$

$${\text{DGF}}_{\mathrm{cell}}=1-\frac{f_{\mathrm{CO}}}{f_{{\mathrm{H}}_2}}.$$

$$A_{\mathrm{V},2D}=(N_{\mathrm{H},\mathrm{tot}}\times 5.348\times 10^{-22}\mathrm{c}{\mathrm{m}}^2)\mathrm{mag}$$

$${\text{DGF}}_{2D}=1-\frac{N_{\mathrm{CO}}\times \frac{1}{2.8\times 10^{-4}}}{N_{{\mathrm{H}}_2}},$$

$$\Delta f_{{\mathrm{H}}_2}(x)=\frac{M_{{\mathrm{H}}_2}^x}{M_{{\mathrm{H}}_2}^{\mathrm{CO}}+M_{{\mathrm{H}}_2}^x}=\frac{M_{{\mathrm{H}}_2}^x}{M_{{\mathrm{H}}_2}}.$$

$$N_{{\mathrm{H}}_2}=X_{\mathrm{CO}}\times I_{\mathrm{CO}}.$$

$$N_{{\mathrm{H}}_2}=10^{20.7}\text{cm−2}\times {\left(\frac{I_{\mathrm{CO}}}{\text{K km s−1}}\right)}^{1/5}$$

$$N_{{\mathrm{H}}_2}=10^{20.75}\text{cm−2}\times {\left(\frac{N_{\mathrm{H},\mathrm{tot}}}{10^{21.5}\text{cm−2}}\right)}^{1.3}=10^{20.75}\text{cm−2}\times {\left(\frac{A_{\mathrm{V},2D}}{1.7}\right)}^{1.3}$$

$$N_{{\mathrm{H}}_2}=0.5\times N_{\mathrm{H},\mathrm{tot}}=0.935\times 10^{21}\text{cm−2}\times A_{\mathrm{V},2D}$$

$$X_{\mathrm{CO}}=X_{\mathrm{CO},0}\times {\left(\frac{⟨A_{\mathrm{V},2D}⟩}{3.5}\right)}^{-3.5}.$$

$$\text{DGF}=1-\frac{⟨f_{\mathrm{CO}}⟩}{⟨f_{{\mathrm{H}}_2}⟩}.$$

Full text

††pubyear: 2019††pagerange: SILCC-Zoom: H2 and CO-dark gas in molecular clouds – The impact of feedback and magnetic fields–A.3

1 Introduction

Molecular clouds (MCs) are the densest structures of the interstellar medium (ISM) and are defined as those regions in which hydrogen exists in its molecular form, H2. Due to the low temperatures of a few 10 K and its vanishing permanent dipole moment, H2 and thus MCs are observable only indirectly e.g. by molecules which trace the presence of H2. One of the most frequently used molecules is CO (e.g. Wilson et al., 1970; Scoville & Solomon, 1975; Larson, 1981; Solomon et al., 1987; Dame et al., 2001; Bolatto et al., 2013, and many more, but see Dobbs et al. 2014 for a review).

However, it has been shown that CO requires a more efficient shielding of the interstellar radiation field (ISRF) than H2 to form (van Dishoeck & Black, 1988; Wolfire et al., 2010). Thus, CO is not a perfect tracer of H2 gas as it misses out a considerable fraction of molecular gas. In this context, the terminology of “CO-poor” or “CO-dark” gas, i.e. gas, where H2 is present but no CO, was established (Lada & Blitz, 1988; van Dishoeck, 1992; Grenier et al., 2005). Using gamma ray emission Grenier et al. (2005) conclude that more than 30% of the H2 gas is CO-dark (see also e.g. Ackermann et al., 2012; Donate & Magnani, 2017).

Since then, CO-dark gas has been subject to a number of studies. By means of dust extinction measurements, CO-dark gas fractions of up to several 10% were found (e.g. Lee et al., 2012; Planck Collaboration et al., 2015), similar to the results obtained from gamma ray emission studies. Furthermore, observations of the [C II] 158 $\mu$m line suggest the presence of CO-dark H2 gas in order to consistently explain the observed line intensity (Langer et al., 2010, 2014; Pineda et al., 2013). These findings were recently supported theoretically by Franeck et al. (2018), who show that up to 20% of the [C II] line emission stems from the molecular phase. Moreover, emission of atomic carbon has been suggested to be a good tracer of CO-dark molecular gas (Gerin & Phillips, 2000; Papadopoulos et al., 2004; Offner et al., 2014; Glover et al., 2015; Li et al., 2018b; Clark et al., 2019).

More recently, other tracers have also been used to study CO-dark gas, in particular the hydroxyl radical OH (e.g. Crutcher et al., 1993; Barriault et al., 2010; Cotten et al., 2012; Allen et al., 2015; Li et al., 2015; Li et al., 2018a; Ebisawa et al., 2019). In a sub-pc resolution observation of OH, Xu et al. (2016) find CO-dark gas fractions varying from 80% to 20% across the boundary of the Taurus molecular cloud. In addition, other molecules like HF, HCl and ArH+ have been suggested to be able to probe the atomic and molecular hydrogen content of MCs and thus the amount of CO-dark gas (Schilke et al., 1995, 2014; Sonnentrucker et al., 2010; Neufeld et al., 1997; Neufeld et al., 2005; Neufeld & Wolfire, 2016).

Taken together, these observations draw a clear picture of MCs where a significant fraction of H2 is not detectable in CO. In order to still be able to infer the amount of H2 gas from CO observations, a conversion factor from the observed CO luminosity into an H2 column density has been established, the so-called “$X_{\mathrm{CO}}$-factor”. Its canonical value in the MilkyWay is assumed to be about 2 $\times$ 1020 cm*-2* (K km s*-1*)-1 (see e.g. the review by Bolatto et al., 2013). However, there are significant cloud-to-cloud variations of $X_{\mathrm{CO}}$ reported in the literature in both observations of galactic and extra-galactic MCs (e.g. Blitz & Thaddeus, 1980; Scoville et al., 1987; Dame et al., 1993; Strong & Mattox, 1996; Melchior et al., 2000; Lombardi et al., 2006; Nieten et al., 2006; Leroy et al., 2011; Smith et al., 2012; Ripple et al., 2013). In addition, also metallicity variations affect the value of XCO (Glover & Mac Low, 2011; Shetty et al., 2011a; Bolatto et al., 2013). All these variations imply uncertainties of a factor of a few in the masses of H2 inferred from CO observations and thus, the $X_{\mathrm{CO}}$-factor might be applicable only for an ensemble of clouds rather than individual clouds (Kennicutt & Evans, 2012).

The large spread of the amount of CO-dark gas and the resulting $X_{\mathrm{CO}}$-factor was confirmed by a number of recent numerical simulations of MC formation (Glover & Mac Low, 2011; Smith et al., 2014; Duarte-Cabral et al., 2015; Glover & Clark, 2016; Richings & Schaye, 2016a, b; Szűcs et al., 2016; Gong et al., 2018; Li et al., 2018b). However, there are two stringent constraints on the accuracy of such numerical approaches: First, due to the highly turbulent structure of MCs and the associated mixing of molecules (Glover et al., 2010; Valdivia et al., 2016; Seifried et al., 2017), the chemical evolution of such clouds has to be modelled on-the-fly in the simulations in order to obtain an accurate picture of their – partly non-equilibrium – chemical state. Second, it can be shown numerically (Seifried et al., 2017) and analytically (Joshi et al., 2019) that a very high spatial resolution of $\lesssim$ 0.1 pc is required to obtain accurate and converged chemical abundances. In case a numerical simulation does not reach this resolution or the chemical (non-equilibrium) abundances are not modelled on-the-fly, inferred fractions of CO-dark gas and values of the $X_{\mathrm{CO}}$-factor have to be considered with caution. In addition, a potential complication arises from highly idealized initial conditions, which do not match the full complexity of real MCs (Rey-Raposo et al., 2015).

So far, only a few simulations match the aforementioned requirements. Moreover, the impact of stellar radiative feedback and of magnetic fields on the amount of CO-dark gas and the properties of CO emission has obtained very little attention so far. In Seifried et al. (2017); Seifried et al. (2019) and Haid et al. (2019) we present some of the first numerical simulations of MC formation which include an on-the-fly chemical network for H2 and CO, high spatial resolution ($\sim$ 0.1 pc), a larger-scale, galactic environment for realistic initial conditions, magnetic fields and stellar feedback. In the following we will use these simulations to investigate the impact of feedback and magnetic fields on CO-dark gas and the observable CO emission in detail.

The structure of the paper is as follows: First, we describe the initial conditions and numerical methods used for the MC simulations (Section 2). We then present our results and discuss the amount and distribution of CO-dark gas (Sections 3.1 and 3.2). Next, we investigate projection effects (Section 3.3) and the effect of CO-dark gas on the CO emission and the $X_{\mathrm{CO}}$-factor (Section 3.4). Finally, in Section 4, we develop a new approach to determine the H2 mass in MCs with a higher accuracy than via the $X_{\mathrm{CO}}$-factor before we conclude in Section 5.

2 Numerics and initial conditions

We present results of the SILCC-Zoom simulations of MC formation (Seifried et al., 2017). The simulations are preformed within the SILCC project (see Walch et al., 2015; Girichidis et al., 2016, for details) and make use of the zoom-in technique discussed in Seifried et al. (2017). The simulations are performed with the adaptive mesh refinement code FLASH 4.3 (Fryxell et al., 2000; Dubey et al., 2008) and use a magneto-hydrodynamics (MHD) solver which guarantees positive entropy and density (Bouchut et al., 2007; Waagan, 2009). We model the chemical evolution of the interstellar medium (ISM) using a simplified chemical network for H+, H, H2, C+, CO, e-, and O (Nelson & Langer, 1997; Glover & Mac Low, 2007; Glover et al., 2010), which also follows the thermal evolution of the gas including the most important heating and cooling processes. We do not apply any particular treatment for compressive shocks like shattering or sputtering of dust grains, the cooling in shocks, however, is captured self-consistently by the applied chemical network.

The interstellar radiation field (ISRF) is that of Draine (1978), i.e. $G_{0}$ = 1.7 in Habing units (Habing, 1968), and its shielding is calculated according to the surrounding column densities of total gas, H2, and CO via the OpticalDepth module (Wünsch et al., 2018) based on the TreeCol algorithm (Clark et al., 2012). For this purpose, we determine for each cell the visual extinction, $A_{\mathrm{V,i}}$, separately along 48 directions by converting the total gas column density $N_{\mathrm{H,tot}}$ into a visual extinction via (Draine & Bertoldi, 1996)

[TABLE]

Since the scheme is based on a Healpix tessellation (Górski & Hivon, 2011), all directions are equally weighted. The average local visual extinction in each cell is then obtained as

[TABLE]

with $\gamma$ = 2.5 (Bergin et al., 2004). With this definition the local attenuation factor of the ISRF due to dust, which we use for the dissociation reactions of H2 and CO, is then exp(-$\gamma$ $A_{\mathrm{V,3D}}$) (Glover et al., 2010). In addition, similar to Eq. 2 for $A_{\mathrm{V,3D}}$, we calculate the self-shielding of H2 and CO from the H2 and CO column densities, which further reduces the dissociation rates (Draine & Bertoldi, 1996; Lee et al., 1996, but see also Section 2 of Glover et al. 2010). This approach thus allows us to properly assess the dissociation of H2 and CO in the simulations due to incident UV radiation.

Furthermore, we solve the Poisson equation for self-gravity with a tree-based method (Wünsch et al., 2018) and include a background potential from the old stellar component in the galactic disc, modeled as an isothermal sheet with $\Sigma_{\mathrm{star}}$ = 30 M*☉* pc*-2* and a scale height of 100 pc.

Our setup represents a small section of a stratified galactic disc with solar neighborhood properties and a size of 500 pc $\times$ 500 pc $\times$ $\pm$ 5 kpc. The gas surface density is $\Sigma_{\mathrm{gas}}$ = 10 M*☉* pc*-2* and the initial vertical gas distribution has a Gaussian profile with a scale height of 30 pc and a midplane density of $\rho_{0}$ = 9 $\times$ $10^{-24}$ g cm*-3*. The gas is initially at rest. Near the disc midplane it has an initial temperature of 4500 K and consists of atomic hydrogen and C+. For the magnetised runs, we initialize a magnetic field along the $x$-direction as

[TABLE]

where we set the magnetic field in the midplane to $B_{x,0}$ = 3 $\mu$G in accordance with recent observations (e.g. Beck & Wielebinski, 2013).

Up to $t_{0}$ (see Table 1), we drive turbulence in the disc with supernovae (SNe). Half of the SNe are randomly placed in the $x$-$y$-plane following a Gaussian profile with a scale height of 50 pc in the vertical direction, the other half is placed at density peaks. The SN rate is constant at 15 SNe Myr*-1*, corresponding to the Kennicutt-Schmidt star formation rate surface density for $\Sigma_{\mathrm{gas}}$ = 10 M*☉* pc*-2* (Kennicutt, 1998) and assuming a standard initial mass function (Chabrier, 2001). For a single SN we inject 1051 erg in the form of thermal energy if the Sedov-Taylor radius is resolved with at least 4 grid cells. Otherwise, we heat the gas within the injection region to $10^{4}$ K and inject the momentum, which the swept-up shell has gained at the end of the Sedov-Taylor phase (see Gatto et al., 2015, for details).

The base grid resolution is 3.9 pc up to $t_{0}$. At $t_{0}$ we stop further SN explosions. We choose different regions in which MCs are about to form. These “zoom-in” regions have a rectangular shape with a typical linear extent of about 100 pc. We then continue the simulations for another 1.5 Myr over which we progressively increase the spatial resolution in these zoom-in regions from 3.9 pc to 0.12 pc assuring that the Jeans length is refined with 16 cells (Seifried et al., 2017, Table 2). In the surroundings we keep the lower resolution of 3.9 pc. Afterwards we continue the simulations with the highest resolution of 0.12 pc in the zoom-in regions.

We consider two purely hydrodynamical (HD) simulations without magnetic fields (runs MC1-HD and MC2-HD, see Seifried et al., 2017, and Table 1) and two simulations with magnetic fields (MC3-MHD and MC4-MHD, see Seifried et al., 2019). For these runs we turn off sink particle formation and stellar feedback inside the clouds. The MCs with and without magnetic fields emerge from different stratified galactic disc simulations. As magnetic fields delay the formation of dense molecular gas (Walch et al., 2015; Girichidis et al., 2018), we start to zoom in at a somewhat later time ($t_{0}$) for the magnetised runs, such that the cloud masses of a few 104 M*☉* are roughly comparable for all four clouds.

Furthermore, in order to investigate the effect of stellar radiative feedback, we rerun MC1-HD and MC2-HD including sink particles and radiative stellar feedback (runs MC1-HD-FB and MC2-HD-FB, see Haid et al., 2019). A comparison between these runs allows us to isolate the impact of feedback from that of different initial conditions. Feedback from the massive stars sets in at $t$ = 13.8 Myr and 13.6 Myr for run MC1-HD-FB and MC2-HD-FB, respectively.

In the two feedback runs, sink particles are used to model the formation of stars or star clusters and their subsequent radiative stellar feedback. The sinks form from Jeans-unstable gas once the gas density exceeds a value of 1.1 $\times$ 10*-20* g cm*-3* and are treated with a 4th-order Hermite predictor-corrector scheme (Dinnbier et al., in prep.). We assure that the cells hosting the sinks are always refined to the highest level of refinement. As time evolves, the sinks accrete gas and form stars. Every 120 M*⊙* of accreted mass, one massive star between 9 and 120 M*⊙* is randomly sampled from an initial mass function assuming a slope of -2.3 between 9 and 120 M*☉* (Salpeter, 1955).

Each massive star follows its individual, mass-dependent stellar evolutionary track (Ekström et al., 2012; Gatto et al., 2017; Peters et al., 2017) where we follow in detail the amount of photoionizing radiation released by each star (Haid et al., 2018, 2019). The radiative feedback is treated with a backwards ray-tracing algorithm TreeRay (Wünsch et al., 2018, Wünsch et al., in prep.), which efficiently uses the available octal-tree structure. The radiative transport equation is solved for hydrogen-ionizing EUV radiation assuming the On-the-Spot approximation with a temperature dependent case B recombination coefficient (Draine, 2011). The resulting number of hydrogen-ionizing photons and the associated heating rate are processed within the chemical network (Haid et al., 2018).

3 Results

Throughout the paper we refer to the time elapsed since $t_{0}$ as $t_{\mathrm{evol}}$ = $t$ - $t_{0}$. Hence, for the runs MC1-HD-FB and MC2-HD-FB, feedback sets in at $t_{\mathrm{evol}}$ = 1.9 and 1.7 Myr, respectively. In the following we constrain ourselves to the times $t_{\mathrm{evol}}$ = 2, 3, 4 and 5 Myr, with the latest time being considered only for the runs MC3-MHD and MC4-MHD, which evolve more slowly111Note that assuming a typical turbulent velocity of $\sim$ 5 km s*-1* (see Fig. 5 of Seifried et al., 2017) and an extent of $\lesssim$ 50 pc, we obtain characteristic turnover timescales of $\lesssim$ 10 Myr.. In Fig. 1 we show the hydrogen nuclei column density and the CO column density of the runs MC1-HD, MC1-HD-FB, and MC3-MHD at $t_{\mathrm{evol}}$ = 4 Myr along the $x$-direction. Overall, CO shows a significantly more compact distribution leading to the problem of CO-dark gas discussed in detail in the following. The typical CO column densities span a range from $\sim$ 1016 – 1019 cm*-2*, peaking around 1017 cm*-2*. This is comparable with actual observations of e.g. the Taurus molecular cloud (Goldsmith et al., 2008). A more detailed analysis of the CO column density distribution, however, will be presented in a subsequent paper (Borchert et al., in prep.). For further details on the dynamical evolution of the clouds we refer to the references given in Table 1. In the following we mainly focus on their chemical composition.

3.1 The CO-dark gas fraction in molecular clouds

As already visible by eye from Fig. 1, the distribution of CO is significantly more compact than the distribution of the total gas. For this reason, we first determine the global mass fraction of CO-dark gas (henceforth DGF) in our simulated MCs using the full 3D information. For this purpose we calculate the ratio of the total CO mass, $M_{\mathrm{CO}}$, to the total H2 mass, $M_{\mathrm{H_{2}}}$, in the zoom-in region. We correct for the fact that in our simulations the total fractional abundance of carbon with respect to hydrogen nuclei is 1.4 $\times$ 10*-4*, or 2.8 $\times$ 10*-4* with respect to H2 molecules (under the assumption that hydrogen is completely in its molecular form), i.e.

[TABLE]

where $m_{\mathrm{p}}$ is the proton mass. Note that this definition describes the fraction of intrinsically CO-dark gas, i.e. gas with no CO molecules but H2 molecules, which is only accessible in simulations. It thus differs from the common observational definition of a DGF, which is based on observational sensitivity limits for CO observations (Wolfire et al., 2010, but see also our Section 3.4.1).

At this point we would like to note that in observational literature, the terms ”CO-dark“, ”CO-poor“ and ”CO-faint“ appear to be used interchangeably (e.g. Lada & Blitz, 1988; van Dishoeck, 1992; Grenier et al., 2005; Bolatto et al., 2013; Velusamy & Langer, 2014). In this work, with (intrinsically) ”CO-dark“ gas we refer to the actual molecular content (Eq. 4) and with ”CO-faint“ gas to an observational definition (Section 3.4.1).

In the top panel of Fig. 2 we show the DGF as a function of time. We first focus on the runs without radiative feedback. For the hydrodynamical runs (MC1-HD and MC2-HD, black and blue solid lines), the DGF remains roughly constant over time with values of $\sim$ 0.4. In the presence of magnetic fields (runs MC3-MHD and MC4-MHD, red and green lines), however, the amount of CO-dark gas is initially significantly higher with values around 0.95 at $t_{\mathrm{evol}}$ = 2 Myr and then decreases over time to values of 0.6 and 0.85, respectively, as the clouds become increasingly denser and more CO forms. Interestingly, the evolution of $M_{\mathrm{H_{2}}}$ (middle panel of Fig. 2) is similar for all four runs. It increases over time with a spread among the simulations of about 104 M*☉, i.e. relative differences of $\sim$ 20%. Hence, the significantly higher fraction of CO-dark gas for the MHD clouds (about a factor of 2) cannot be attributed to changes in $M_{\mathrm{H_{2}}}$ but to a lower amount of CO in these runs (bottom panel of Fig. 2). We attribute this lower amount of CO to differences in the structure of the clouds. This becomes already apparent by eye when investigating Fig. 1. The MCs with magnetic fields appear to be more diffuse and filamentary than the MCs without magnetic fields. For a fixed $M_{\mathrm{H_{2}}}$, however, a more filamentary structure would result – on average – in lower visual extinctions and thus lower $M_{\mathrm{CO}}$ (Röllig et al., 2007; Glover et al., 2010). For completeness, in Fig. 18 in Appendix A, we also show the mean H2 and CO densities in the dense gas ($n$ $\geq$ 100 cm-3*) and the global mass fractions of both species in the zoom-in regions.

At this point, we emphasize that here and in Section 3.2 we consider the local visual extinction of the ISRF at each point in the cloud, i.e. $A_{\mathrm{V,3D}}$, which is calculated directly during the simulation via the OpticalDepth module (Eq. 2 and Wünsch et al., 2018). It thus gives an reasonable approximation of the local attenuation of the ISRF, but does not directly correspond to the visual extinction obtained in observations by averaging along the line-of-sight (LOS), which we consider in Sections 3.3 and 4.

In Fig. 3, we plot the mass-weighted probability density function (PDF) of log($A_{\mathrm{V,3D}}$) for the four different MCs without feedback at $t_{\mathrm{evol}}$ = 2, 3 and 4 Myr, where the mass-weighted PDF of a quantity $x$ is given by

[TABLE]

As speculated before, the magnetised clouds MC3-MHD and MC4-MHD are more diffuse objects with significantly smaller mass fractions at $A_{\mathrm{V,3D}}$ $>$ 1 than the clouds without magnetic fields (MC1-HD and MC2-HD). This is also supported by the mean densities of H2 shown in the top left panel of Fig. 18. Hence, as CO only starts to form at AV,3D $>$ 1, and H2 already at $A_{\mathrm{V,3D}}$ $\gtrsim$ 0.3 (e.g. Röllig et al., 2007; Glover et al., 2010; Bisbas et al., 2019), this explains the observed differences in the DGF. It also matches our previous findings that magnetic fields significantly hamper the formation of dense, well-shielded molecular gas (Walch et al., 2015; Girichidis et al., 2018). Our findings are, however, in contradiction to the result of 1D calculations of Wolfire et al. (2010), who claim that the amount of CO-dark gas is insensitive to the internal density – and thus $A_{\mathrm{V,3D}}$ – distribution, thus emphasising the need of 3D, MHD simulations. We note that the increase of mass at $A_{\mathrm{V,3D}}$ $>$ 1 is accompanied with an increase in $M_{\mathrm{CO}}$ (bottom panel of Fig. 2), which is particularly pronounced for MC3-MHD (red lines). We will investigate the dependence of the DGF on the shielding in more detail in the next section confirming the results shown so far.

For runs including feedback (MC1-HD-FB and MC2-HD-FB, dashed lines in Fig. 2), the DGF increases over time from $\sim$ 0.4 to $\sim$ 0.7 as CO is apparently destroyed more efficiently via photodissociation than H2 (compare middle and bottom panel). We attribute this to the fact that (i) radiative feedback from young, massive stars acts where the stars are born, i.e. preferentially in the densest regions of MCs, which are fully molecular, and (ii) the dissociation rate per molecule and per UV photon incorporated in our chemical network is a factor of 3.86 times higher for CO than for H2 (van Dishoeck & Black, 1988; Röllig et al., 2007). Hence, MCs actively forming massive stars appear to have a higher amount of CO-dark gas than their quiescent counterparts.

The overall high values of the DGF and its large spread agree well with recent observations showing DGFs in MCs of 30% and more (Grenier et al., 2005; Lee et al., 2012; Planck Collaboration et al., 2015). Our results are also in agreement with theoretical results (Wolfire et al., 2010; Smith et al., 2014; Gong et al., 2018; Li et al., 2018b), although the definition used by these authors does not exactly match the DGF as defined in Eq. 4 (see Section 3.4.1 for details).

3.2 The distribution of CO-dark and CO-bright gas

3.2.1 Dependence on the local $A_{\mathrm{V,3D}}$

As shown before, the differences in the DGF by a factor of $\lesssim$ 2 between runs with and without magnetic fields (Fig. 2) can be attributed to different morphologies of the clouds, resulting in less mass at $A_{\mathrm{V,3D}}$ $>$ 1 for the magnetised clouds (Fig. 3). For this reason, we next consider the distribution of H2 and CO relative to each other and with respect to the local visual extinction, $A_{\mathrm{V,3D}}$. For this purpose, we determine the mass fractions of H2 and CO

[TABLE]

where $m_{\mathrm{H,cell,tot}}$ and $m_{\mathrm{C,cell,tot}}$ are the total mass of hydrogen and carbon available in the considered cell, $m_{\mathrm{H_{2},cell}}$ and $m_{\mathrm{CO,cell}}$ the mass of all H2 and CO molecules, respectively, and the factor $\frac{12}{28}$ corrects for the mass of oxygen.

In Fig. 4 we show the results for run MC1-HD at $t_{\mathrm{evol}}$ = 3 Myr. The qualitative behaviour also holds for all other runs and times. As can be seen from the left panel, significant amounts of CO are only formed once $\sim$ 50% of the hydrogen is in molecular form. Both mass fractions become roughly comparable above $f_{x}$ $\simeq$ 0.8 $\pm$ 0.1 (where the subscript $x$ stands for H2 and CO, respectively). This value is a rather rough estimate, which also depends on the considered simulation and increases over time ranging from $\sim$ 0.5 – 0.8.

Considering the middle and right panel of Fig. 4, we find that for both H2 and CO $f_{x}$ = 0.8 is approximately reached at $A_{\mathrm{V,3D}}$ $\simeq$ 1.5 (see violet lines), around which we expect the DGF to drop to zero. At lower $A_{\mathrm{V,3D}}$ of a few 0.1, $f_{\mathrm{H_{2}}}$ is as high as 0.1 – 0.2 (middle panel), i.e. noticeably higher than $f_{\mathrm{CO}}$, which remains close to zero in this range (right panel) and starts to rise around $A_{\mathrm{V,3D}}$ $\simeq$ 1 in good agreement with detailed chemical PDR models (e.g. Röllig et al., 2007; Glover et al., 2010, see also Gong et al. 2018 for similar results in 3D-MHD simulations). At even lower values of $A_{\mathrm{V,3D}}$ ($\lesssim$ 0.1), neither H2 nor CO are present. We note that sometimes $f_{\mathrm{CO}}$ can be slightly higher than $f_{\mathrm{H_{2}}}$, which can be attributed to the short formation time of CO, once a sufficient amount of H2 is present (see Eq. 9 in Seifried et al., 2017, but also Nelson & Langer 1997; Glover et al. 2010).

Next, we investigate how the DGF of individual cells depends on $A_{\mathrm{V,3D}}$. The DGF in a cell is related to the mass fractions of H2 and CO (Eq. 6) as

[TABLE]

In Fig. 5 we show the $A_{\mathrm{V,3D}}$-DGFcell-phase diagram for the clouds MC1-HD, MC1-HD-FB and MC3-MHD at $t_{\mathrm{evol}}$ = 2, 3 and 4 Myr. We emphasize that for the remaining runs the phase diagrams are similar. The general shape of the distribution changes only moderately with time and for the different runs. Once molecular hydrogen starts to form ($A_{\mathrm{V,3D}}$ $\gtrsim$ 0.1), DGFcell quickly rises to $\sim$ 1. It remains high until $A_{\mathrm{V,3D}}$ $\simeq$ 1 and then drops to almost zero around $A_{\mathrm{V,3D}}$ = 1.5 – 2 as indicated by the black solid and dashed lines which denote the mean and median of the distribution. This is in good agreement with our findings that the mass fractions of H2 and CO become comparable at $A_{\mathrm{V,3D}}$ $\simeq$ 1.5 (Fig. 4). Furthermore, the observed drop of DGFcell around $A_{\mathrm{V,3D}}$ = 1.5 is also in good agreement with the findings of Xu et al. (2016) in the Taurus molecular cloud.

For the runs with stellar feedback (middle row of Fig. 5), somewhat more CO-dark gas appears at $A_{\mathrm{V,3D}}$ $>$ 1.5 at later stages (seen as a horizontal stripe in the phase diagram). The overall shape, however, remains almost unchanged. For MC3-MHD (bottom row), initially ($t_{\mathrm{evol}}$ $\leq$ 3 Myr), the drop of DGFcell seems to appear at slightly higher $A_{\mathrm{V,3D}}$. As at this evolutionary stage the structure of the cloud is significantly more diffuse (Fig. 3), we speculate that also the self-shielding due to H2 and CO, which prevents CO from being dissociated, is reduced. Hence, CO forms at higher $A_{\mathrm{V,3D}}$ than for the more compact hydrodynamical clouds. However, towards later stages, the phase-diagrams approach those from the runs without magnetic fields. Similar results are also found for MC4-MHD (not shown here).

Integrating the plots in Fig. 5 along the $x$-axis gives the amount of gas at a given DGFcell. In Fig. 6 we show the PDF of DGFcell weighted by the CO and H2 mass, respectively. For the sake of readability, we only show the lines for the runs MC1-HD, MC3-MHD, and MC1-HD-FB, but note that the remaining runs are qualitatively and quantitatively very similar. As can be seen, the vast majority (note the logarithmic scaling of the $y$-axis) of CO sits in cells with a low DGFcell. For H2, however, significant amounts of gas can be found close to DGFcell = 1 and DGFcell = 0, which is in good agreement with the partly high global DGF (Fig. 2).

To summarize, our results indicate that – despite significant variations in the absolute amount – CO-dark gas is mainly present in gas with visual extinctions 0.2 – 0.3 $<$ $A_{\mathrm{V,3D}}$ $<$ 1 – 1.5, independent of the presence or absence of magnetic fields or stellar feedback. Above $A_{\mathrm{V,3D}}$ $\simeq$ 1.5, the gas is mainly CO-bright, which happens once about 50 – 80% of both hydrogen and carbon are in molecular form.

3.2.2 Dependence on density and temperature

Next, in Fig. 7 we investigate the distribution of CO-dark and CO-bright gas in the density ($n_{\mathrm{H,tot}}$) – temperature (T) phase space for a run without feedback (MC1-HD). The findings discussed in the following are also representative for the remaining runs with feedback and magnetic fields (see Fig. 19 in the Appendix). We define the CO-dark gas (black contours) as all H2 gas in cells with DGFcell $>$ 0.5 (Eq. 7) and the CO-bright gas (green contours) as all H2 gas in cells with DGFcell $\leq$ 0.5.

The bulge of CO-bright gas sits at $n_{\mathrm{H,tot}}$ $\gtrsim$ 300 cm*-3* and temperatures below $\sim$ 50 K, although in particular for the runs with feedback (top panel of Fig. 19) some CO-bright gas can be found at temperatures up to a few 100 K. The bulge of CO-dark gas, however, occurs at densities of 10 cm*-3* $\lesssim$ $n_{\mathrm{H,tot}}$ $\lesssim$ 103 cm*-3* and temperatures of a few $10$ K $\lesssim$ $T$ $\lesssim$ a few 100 K. Our results are thus in rough agreement with the findings of Glover & Smith (2016), who find CO-dark gas to reside at temperatures above $\sim$ 30 K.

There is, however, a substantial overlap of CO-dark and -bright gas in the $n_{\mathrm{H,tot}}$-$T$-plane (see also Section 3.3). Moreover, radiative feedback (top panel of Fig. 19) even further extends the region in the $n_{\mathrm{H,tot}}$-$T$-parameter space, in which CO-dark gas is found. We speculate that this broad distribution of CO-dark gas might complicate its identification in actual observations. This is supported by recent theoretical studies showing the necessity of observing various lines like [OI], [CI] and [CII] to capture the entire CO-dark gas component of MCs (Glover & Smith, 2016; Franeck et al., 2018; Li et al., 2018b; Clark et al., 2019). In this context, we note that other tracers like ArH+, HF, or HCl have also been suggested to differentiate between the atomic and molecular (hydrogen) phase in MCs (Schilke et al., 1995, 2014; Neufeld et al., 1997; Neufeld et al., 2005; Neufeld & Wolfire, 2016), which might allow for a more accurate estimate of the H2 content of MCs.

3.3 The DGF in 2D maps

So far we have based our analysis on the local value of the visual extinction, $A_{\mathrm{V,3D}}$. This value, however, is not accessible in observations, where only a LOS-integrated column density and the corresponding (2-dimensional) visual extinction are accessible, e.g. via dust extinction measurements (e.g. Lombardi & Alves, 2001). Hence, in order to allow for a better comparison with actual observations, we integrate along the $x$-, $y$-, and $z$-direction of each zoom-in region to obtain maps of the projected total hydrogen column density, $N_{\mathrm{H,tot}}$, and convert it to a visual extinction, $A_{\mathrm{V,2D}}$, via

[TABLE]

(Bohlin et al., 1978; Draine & Bertoldi, 1996). In addition, we calculate the projected DGF for each pixel in the map similar to Eq. 4:

[TABLE]

where $N_{\mathrm{H_{2}}}$ and $N_{\mathrm{CO}}$ are the surface densities of H2 and CO in each pixel. We use a pixel size of 0.12 pc identical to the maximum resolution of the simulations. In Fig. 8 we show the resulting H2-mass-weighted $A_{\mathrm{V,2D}}$-DGF2D-phase diagram obtained from the 2D maps integrated along the $y$-direction for the runs MC1-HD, MC1-HD-FB, and MC3-MHD at $t_{\mathrm{evol}}$ = 3 Myr. Similar results are obtained for the remaining runs, times and directions.

In contrast to the results for $A_{\mathrm{V,3D}}$ (Fig. 5) and observational results of Xu et al. (2016), DGF2D drops to zero at somewhat higher values of $A_{\mathrm{V,2D}}$ $\simeq$ 2 – 4. In addition, we find a significant amount of CO-dark gas up to $A_{\mathrm{V,2D}}$ $\simeq$ 5, where one would not expect it (see Fig. 5, but also Röllig et al., 2007; Glover et al., 2010). Furthermore, there appears to be a broad distribution of DGF2D in the range of $A_{\mathrm{V,2D}}$ $\simeq$ 2.5 – 5, with both CO-dark and CO-bright gas being present. This demonstrates that, under certain circumstances, $A_{\mathrm{V,2D}}$ – which is an average quantity – can give only little insight about the actual conditions along the entire LOS. This is, however, not surprising given the overlap of CO-dark and -bright gas in the $\rho$ - $T$ - plane found in Fig. 7.

The appearance of CO-dark gas at $A_{\mathrm{V,2D}}$ $\simeq$ 5 could be understood when considering the average density along the LOS of such a pixel. Assuming a typical length of the LOS of $\sim$ 50 pc, with Eq. 8 we obtain a total hydrogen column density of $\sim$ 1022 cm*-2* and a volume density of about 60 cm*-3*. Using the relation between the density and $A_{\mathrm{V,3D}}$ found in our simulations (see Fig. 11 in Seifried et al., 2017), such a density corresponds to a typical $A_{\mathrm{V,3D}}$ $\leq$ 1. Hence, under the assumption that the gas is uniformly distributed along the LOS, we expect CO to not have formed yet, and thus DGF2D to be close to 1 at $A_{\mathrm{V,2D}}\simeq 5$.

However, the assumption of a uniform density distribution along the LOS is clearly an oversimplification. Hence, in order to fully understand the reason for the broad distribution of DGF2D around $A_{\mathrm{V,2D}}$ $\simeq$ 2.5 – 5, we investigate the distribution of the local visual extinction $A_{\mathrm{V,3D}}$ and the density along the LOS of all pixels with

2.5 $<$ $A_{\mathrm{V,2D}}$ < 5 and DGF2D > 0.9, i.e. CO-dark gas (magenta box in the left panel of Fig. 8) and 2. 2.

2.5 $<$ $A_{\mathrm{V,2D}}$ < 5 and DGF2D < 0.1, i.e. CO-bright gas (green box).

We show the mass-weighted $A_{\mathrm{V,3D}}$-PDF and $n_{\mathrm{H,tot}}$-PDF of both subsets (i) and (ii) for MC1-HD at $t_{\mathrm{evol}}$ = 3 Myr in Fig. 9; for the other runs and times we obtain qualitatively and quantitatively similar results. We find no differences in the $A_{\mathrm{V,3D}}$-PDF for the three directions considered (left panel). There is, however, a clear difference in the distributions for CO-dark and CO-bright gas. CO-dark gas (solid lines) shows a much more narrow $A_{\mathrm{V,3D}}$-distribution which peaks at $A_{\mathrm{V,3D}}$ $\sim$ 0.7 and reaches a maximum of $A_{\mathrm{V,3D}}$ $\simeq$ 1.5. For CO-bright gas (dashed lines), however, the distribution is much more wide-spread and shifted to higher $A_{\mathrm{V,3D}}$ with the peak occurring around $\sim$ 1.5.

The density PDF (right panel) shows a corresponding behaviour, which is not surprising given the tight relation between $A_{\mathrm{V,3D}}$ and the density found in our simulations (Fig. 11 in Seifried et al., 2017). For CO-dark gas, the PDF peaks around a value of $n_{\mathrm{H,tot}}$ $\sim$ 100 cm*-3*, which is comparable to the average density obtained under the assumption of a uniform gas distribution along the LOS (see above). For CO-bright gas, however, the PDF peaks at 20 – 30 times higher densities. Altogether, we can thus attribute the broad distribution of DGF2D found in Fig. 8 to pixels with different density distributions along the LOS: For CO-dark gas at $A_{\mathrm{V,2D}}$ $\simeq$ 2.5 – 5, we have a rather uniform density – and thus $A_{\mathrm{V,3D}}$ – distribution with $A_{\mathrm{V,3D}}$ $\leq$ 1 and hence very little CO, whereas H2 is present already. For the CO-bright gas, however, we have regions with strong density contrasts and locally well-shielded gas ($A_{\mathrm{V,3D}}$ $\geq$ 1.5), causing CO to form. This is in excellent agreement with complementary theoretical (Levrier et al., 2012) and observational works (Busch et al., 2019), which both find that the abundance of CO is increased by local density enhancements along the LOS.

To summarise, our results indicate that the visual extinction inferred from a LOS-averaging process ($A_{\mathrm{V,2D}}$) as naturally done in observations, is a partly misleading quantity to assess the CO content along the LOS. It should therefore be considered with caution and be complemented with actual CO observations. Furthermore, for a given $A_{\mathrm{V,2D}}$, the actual $A_{\mathrm{V,3D}}$-distribution can be relatively broad and show significant qualitative differences for different pixels in agreement with findings of Clark & Glover (2014, their Fig. 12). We emphasize that recent observations of M17 and Monoceros R2 with the SOFIA telescope also indicate significant emission in [CII] at $A_{\mathrm{V,2D}}$ $\sim$ 8 – i.e. implying a high DGF – in good agreement with our findings (Guevara et al. in prep.)

3.4 CO observations

Next, we investigate to which extent the actual H2 mass, $M_{\mathrm{H_{2}}}$, can be obtained from 12CO(1-0) line observations222In the following we drop the superscript “12“.. For this purpose, we use the freely available radiative transfer code RADMC-3D (Dullemond, 2012) to produce synthetic CO(1-0) line emission maps of our MCs along the $x$-, $y$-, and $z$-direction at the same resolution as the simulation data, i.e. 0.12 pc. We use the Large Velocity Gradient method to calculate the level population and the resulting intensity of the CO(1-0) line transitions. The molecular data, e.g. the Einstein coefficients, are taken from the Leiden Atomic and Molecular database (Schöier et al., 2005). The line emission maps cover a velocity range of $\pm$ 20 km s*-1*, which guarantees that all emission is captured properly. The channel width is 200 m s*-1*, which results in 201 channels. We show the CO-spectra of the various runs in Fig. 20 in the Appendix. Further discussion of the spectra, however, is beyond the scope of this paper and will be postponed to a subsequent publication (Nürnberger et al., in prep.).

3.4.1 The sensitivity of CO observations

As a first step we define the fraction of H2 which is in regions with CO emission below the observational sensitivity limit, following the definition of Wolfire et al. (2010) in the notation given by Smith et al. (2014, Eq. 4):

[TABLE]

Here, $M_{\mathrm{H_{2}}}^{x}$ is the mass of all H2 gas in pixels which have a CO(1-0) intensity of $I_{\mathrm{CO}}\leq x$ and $M_{\mathrm{H_{2}}}^{\mathrm{CO}}$ all H2 gas in pixels with $I_{\mathrm{CO}}>x$ such that $M_{\mathrm{H_{2}}}^{\mathrm{CO}}+M_{\mathrm{H_{2}}}^{x}=M_{\mathrm{H_{2}}}$. I.e., assuming an observational sensitivity limit of $I_{\mathrm{CO}}=x$, the fraction $\Delta f_{\mathrm{H_{2}}}(x)$ of H2 cannot be traced via CO in the observation regardless of how accurately $M_{\mathrm{H_{2}}}^{\mathrm{CO}}$ can be determined. This CO-faint H2 gas333We note that $\Delta f_{\mathrm{H_{2}}}$ has also been denoted as a dark-gas fraction (Wolfire et al., 2010; Levrier et al., 2012; Smith et al., 2014; Gong et al., 2018; Li et al., 2018b). However, as stated before, while the DGF defined here via Eq. 4 is a quantity which requires knowledge about the actual CO abundance, which is only accessible via simulation data and does not depend on any threshold, $\Delta f_{\mathrm{H_{2}}}$ relies on the observable CO intensity and thus depends on the sensitivity limit of the observation. thus amplifies the problem of intrinsically CO-dark gas, i.e. gas which has no CO molecules but H2 molecules (Section 3.1). We again note that the latter problem exists for the H2 gas inside and outside the observable region, i.e. for $M_{\mathrm{H_{2}}}^{\mathrm{CO}}$ and $M_{\mathrm{H_{2}}}^{x}$.

As shown in Fig. 10, $\Delta f_{\mathrm{H_{2}}}$ does only weakly depend on the considered direction. Furthermore, the qualitative shape of the curves seems to be similar for all six MCs: Below a threshold value of $I_{\mathrm{CO}}$ $\simeq$ a few 10 K km s*-1*, the value of $\Delta f_{\mathrm{H_{2}}}$ decreases steadily with decreasing threshold $I_{\mathrm{CO}}$, which is in agreement with recent observational results of Donate & Magnani (2017). At $I_{\mathrm{CO}}$ $\simeq$ a few 10 K km s*-1*, all curves rise quickly until they reach unity at a few 100 K km s*-1* (see also Smith et al., 2014; Li et al., 2018b).

However, there are significant differences in the absolute values of $\Delta f_{\mathrm{H_{2}}}$ for the different MCs. Assuming a typical sensitivity limit of 0.1 K km s*-1* of recent CO(1-0) observations of nearby MCs (e.g. Nieten et al., 2006; Pineda et al., 2010; Smith et al., 2012; Ripple et al., 2013; Leroy et al., 2016, indicated by the grey vertical lines in Fig. 10), we obtain $\Delta f_{\mathrm{H_{2}}}(x)$ $\simeq$ 15 – 65%. For the runs without feedback, $\Delta f_{\mathrm{H_{2}}}$ drops over time as more and more CO forms and the intensity in the individual pixels increases. Contrary to that, for the runs with feedback, $\Delta f_{\mathrm{H_{2}}}$ increases over time as CO gets destroyed by the radiation released from the forming stars. Furthermore, when comparing Fig. 2 and 10, we find a positive correlation between the DGF and $\Delta f_{\mathrm{H_{2}}}$.

As the curves in Fig. 10 are rather shallow below $I$ = 10 K km s*-1*, $\Delta f_{\mathrm{H_{2}}}$ is not very sensitive on the chosen intensity threshold. Increasing the sensitivity limit to 1 K km s*-1* typical for larger-scale CO surveys (e.g. Dame et al., 2001, and references therein) increases $\Delta f_{\mathrm{H_{2}}}$ only to values of 20 – 75%. Conversely, for a (hypothetical) 10 times lower CO(1-0) sensitivity limit of 0.01 K km s*-1*, the fraction of CO-faint H2 gas would be only marginally reduced to $\sim$ 10 – 50%. We note that similar values of $\Delta f_{\mathrm{H_{2}}}$ are found by various other authors (Wolfire et al., 2010; Levrier et al., 2012; Smith et al., 2014; Gong et al., 2018; Li et al., 2018b), although, as stated before, the authors denote it as DGFs.

The results show that in some cases there is a significant amount of (diffuse) H2 gas in low-$I_{\mathrm{CO}}$ regions, which can be problematic in actual observations: Even when neglecting complications of intrinsically CO-dark gas in regions with $I_{\mathrm{CO}}$ $>$ 0.1 K km s*-1*, i.e. assuming that the amount of H2 in those regions ($M_{\mathrm{H_{2}}}^{\mathrm{CO}}$) can be determined accurately from CO, the observations would miss a significant amount of H2 in the clouds due to the sensitivity limit.

3.4.2 The $X_{\mathrm{CO}}$-factor

In order to obtain $M_{\mathrm{H_{2}}}$ of MCs from CO(1-0) observations, typically a fixed conversion factor, the so-called $X_{\mathrm{CO}}$-factor is used, such that

[TABLE]

The canonical value of $X_{\mathrm{CO}}$ in the MilkyWay is assumed to be about 2 $\times$ 1020 cm*-2* (K km s*-1*)-1 (Dame et al., 1993; Strong & Mattox, 1996, but see also the review by Bolatto et al. 2013). However, as the DGF is varying strongly among the MCs (Section 3.1), this raises the question to what extent the $X_{\mathrm{CO}}$-factor is affected as well.

For this reason, we first investigate the relation between the H2 column density, $N_{\mathrm{H_{2}}}$, and the simulated CO(1-0) intensity for three representative clouds MC1-HD, MC1-HD-FB, and MC3-MHD along the $x$-direction at $t_{\mathrm{evol}}$ = 3 Myr (Fig. 11). Overall, we find only little variation in the general functional shape when including either magnetic fields or stellar feedback and when considering different directions or times (with the latter two not shown here), although for the runs with feedback the distribution becomes somewhat broader. For all runs we find a strong increase in $I_{\mathrm{CO}}$ above $N_{\mathrm{H_{2}}}$ $\sim$ a few $\times$ 1019 cm*-2* (see also Federman et al., 1980; Liszt & Lucas, 1998; Sheffer et al., 2008, for observational examples). Above $N_{\mathrm{H_{2}}}$ $\simeq$ 1021 cm*-2*, the maximum CO intensity saturates around $I_{\mathrm{CO}}$ $\sim$ 100 K km s*-1*, where also most of the H2 mass sits. This saturation is also seen in other theoretical and observational works (e.g. Pineda et al., 2008; Ripple et al., 2013; Smith et al., 2014; Gong et al., 2018) and can be attributed to the fact that CO(1-0) becomes optically thick already around an $A_{\mathrm{V,2D}}$ of $\sim$ 1 (Seifried et al., 2017).

For comparison, we also show the relation obtained when using a fixed $X_{\mathrm{CO}}$-factor to convert $I_{\mathrm{CO}}$ to $N_{\mathrm{H_{2}}}$ (blue lines). Overall, there is a scatter of up to several orders of magnitude around this relation, which is in agreement with our previous results (Seifried et al., 2017, but see also e.g. Smith et al. 2014; Gong et al. 2018). Furthermore, the distribution of $N_{\mathrm{H_{2}}}$ and $I_{\mathrm{CO}}$ shows an (almost) linear relation – required for XCO to be applicable – only for a small range of column densities around $N_{\mathrm{H_{2}}}$ $\simeq$ 1021 cm*-2*. For runs with stellar feedback this range is somewhat more extended to higher column densities, as here a significant amount of CO gets destroyed (Section 3.1) and thus CO(1-0) becomes optically thick, i.e. the $N_{\mathrm{H_{2}}}$-$I_{\mathrm{CO}}$ relation becomes flat, only at higher $N_{\mathrm{H_{2}}}$. Overall, however, our findings agree with previous results that on sub-pc scales, i.e. for individual pixels, the $X_{\mathrm{CO}}$-factor is not applicable (e.g. Glover & Mac Low, 2011; Shetty et al., 2011a, b; Bolatto et al., 2013; Finn et al., 2019).

In Fig. 12 we show the $X_{\mathrm{CO}}$-factor obtained by integrating both the column density and the CO(1-0) intensity over the observable regions, i.e. where $I_{\mathrm{CO}}>$ 0.1 K km s*-1*. There are clear differences between the different MCs recognisable, which directly translates into uncertainties in the inferred H2 cloud mass. Most prominently, as already found in Glover & Clark (2016) and Seifried et al. (2017), the $X_{\mathrm{CO}}$-factor of the hydrodynamical runs without feedback (left panel) increases over time with typical values around 0.5 – 2 $\times$ 1020 cm*-2* (K km s*-1*)-1.

In contrast to that, for the runs including magnetic fields (right panel), $X_{\mathrm{CO}}$ partly decreases over time with typical values from 0.8 – 4 $\times$ 1020 cm*-2* (K km s*-1*)-1 (see also Richings & Schaye, 2016a, b). However, towards later stages the values appear to converge around 1 – 2 $\times$ 1020 cm*-2* (K km s*-1*)-1. We find that more diffuse clouds (here the MHD clouds, see Fig. 3) tend to have somewhat higher $X_{\mathrm{CO}}$-factors than more compact clouds (here the HD clouds). This is in good agreement with observational findings for the Perseus, Taurus, and Orion molecular cloud (Pineda et al., 2008; Pineda et al., 2010; Ackermann et al., 2012; Lee et al., 2014) showing lower $X_{\mathrm{CO}}$-factors for denser and more compact sub-regions (see also Glover & Mac Low, 2011; Szűcs et al., 2016; Seifried et al., 2017, for similar numerical results).

The values of $X_{\mathrm{CO}}$ of the runs with feedback (middle panel of Fig. 12) are typically somewhat lower ($\lesssim$ 1 $\times$ 1020 cm*-2* (K km s*-1*)-1) than those of the runs without feedback. We attribute this to the fact that in these runs CO(1-0) becomes optically thick later (middle panel of Fig. 11) and thus the amount of CO intensity for a given amount of H2 is higher than that for a run without feedback.

Overall, our average value for $X_{\mathrm{CO}}$ is around 1.5 $\times$ 1020 cm*-2* (K km s*-1*)-1 in good agreement with other theoretical works (e.g. Glover & Mac Low, 2011; Smith et al., 2014; Duarte-Cabral et al., 2015; Glover & Clark, 2016; Richings & Schaye, 2016a, b; Szűcs et al., 2016; Gong et al., 2018; Li et al., 2018b), although these works tend to have spatial resolutions coarser than the required limit of 0.1 pc (Seifried et al., 2017; Joshi et al., 2019) or simplified descriptions of the chemical evolution. However, our results also show that the actual $X_{\mathrm{CO}}$-factor can vary by up to a factor of $\sim$ 4 in either direction for different MCs. This in turn implies the same uncertainty of a factor of $\sim$ 4 for the inferred H2 cloud masses. As e.g. the virial parameter scales with $M_{\mathrm{H_{2}}}^{2}$, this can result in uncertainties of one order of magnitude for inferred quantities. We note that the partly significant cloud-to-cloud variations of $X_{\mathrm{CO}}$ reported here are in good agreement with variations reported over decades in observations of galactic and extra-galactic MCs (e.g. Blitz & Thaddeus, 1980; Scoville et al., 1987; Dame et al., 1993; Strong & Mattox, 1996; Melchior et al., 2000; Lombardi et al., 2006; Nieten et al., 2006; Leroy et al., 2011; Smith et al., 2012; Ripple et al., 2013, but see also the review by Bolatto et al. 2013) and the theoretical works noted before. This indicates that, besides being not applicable on sub-pc scales, the $X_{\mathrm{CO}}$-factor might have its strength when being applied for an ensemble of MCs rather than individual MCs (e.g. Kennicutt & Evans, 2012).

4 Towards a new approach to determine the H2 content of molecular clouds

In the literature several ways are commonly used to estimate the (molecular) mass of MCs. As described before, the often used $X_{\mathrm{CO}}$-factor is subject to significant cloud-to-cloud variations, which in turn impose a factor-of-a-few uncertainty for the H2 mass in any MC (neglecting the CO-faint H2 gas discussed in Section 3.4.1). Furthermore, at low column densities ($N_{\mathrm{H_{2}}}$ $\lesssim$ 1021 cm*-2*), the strong correlation between CO and H2 (see Fig. 11) has been used to estimate the H2 content of MCs (e.g. Federman et al., 1980; Liszt & Lucas, 1998; Sheffer et al., 2008). This approach, however, breaks down at high column densities, where CO becomes optically thick. In this latter regime, the usage of extinction maps in combination with Eq. 8 can provide the hydrogen mass of MCs. However, this mass describes the total hydrogen mass (atomic hydrogen and H2) and should thus not be mixed up with the actual molecular H2 mass, as the gas might not necessarily be fully molecular (Fig. 13). Hence, it is questionable to which extent $A_{\mathrm{V,2D}}$ can be used to assess the molecular content/chemical state of an MC; a question, which we will investigate further below.

In the following we describe a new approach, which tries to reduce the shortcomings of the aforementioned approaches and is based on three column density regimes (or alternatively $A_{\mathrm{V,2D}}$ regimes, Eq. 8). We define these regimes in Fig. 13, which shows the relation between $N_{\mathrm{H_{2}}}$ and $N_{\mathrm{H,tot}}$ for MC1 at $t$ = 3 Myr. Above $N_{\mathrm{H,tot}}$ $\simeq$ 1021.5 cm*-2*, $N_{\mathrm{H_{2}}}$ shows a good correlation with $N_{\mathrm{H,tot}}$, whereas below there is a scatter of half an order of magnitude and more. For this reason, for regime (1) we rely on CO(1-0) observations and for the regimes (2) and (3) we rely on extinction measurements. In the following, we first describe in detail the principles of the approach, before we provide the actual numbers and interpret the results.

1. (1)

$N_{\mathrm{H,tot}}<10^{21.5}$ cm*-2*:

As shown in Fig. 11, for values between $I_{\mathrm{CO}}$ $\simeq$ 10*-4* – 10 K km s*-1*, there is a good correlation between $N_{\mathrm{H_{2}}}$ and $I_{\mathrm{CO}}$, which allows us to express $N_{\mathrm{H_{2}}}$ as a function of $I_{\mathrm{CO}}$. For this purpose, in the left panel of Fig. 14 we show the mean H2 column density, $\left\langle N_{\mathrm{H_{2}}}\right\rangle$, as a function of $I_{\mathrm{CO}}$ for all runs at all times and all three LOS considered. For a given CO intensity, the values of $\left\langle N_{\mathrm{H_{2}}}\right\rangle$ vary by about 0.1 dex, i.e. 25% in either direction. Given this good qualitative agreement among the different runs, we approximate $\left\langle N_{\mathrm{H_{2}}}\right\rangle$ by a powerlaw (see Eq. 12 below) where we focus on matching the curves in the range 10*-2* K km s*-1* $<$ $I_{\mathrm{CO}}$ $<$ 10 K km s*-1*. 2. (2)

$10^{21.5}$ cm*-2* $\leq N_{\mathrm{H,tot}}<10^{23}$ cm*-2*:

As the correlation between $I_{\mathrm{CO}}$ and $N_{\mathrm{H_{2}}}$ breaks down at $N_{\mathrm{H_{2}}}\gtrsim 10^{21}$ cm*-2* ($I_{\mathrm{CO}}$ $\gtrsim$ 10 K km s*-1*, Fig. 11), we next consider the relation between $N_{\mathrm{H_{2}}}$ and $N_{\mathrm{H,tot}}$, which is accessible through dust extinction measurements. This is shown exemplarily in Fig. 13 for MC1-HD along the $x$-direction at $t_{\mathrm{evol}}$ = 3 Myr. We find that for $N_{\mathrm{H,tot}}$ $\geq$ 1021.5 cm*-2* there is a good correlation between both quantities. For this reason, in the right panel of Fig. 14 we plot $\left\langle N_{\mathrm{H_{2}}}\right\rangle$ as a function of $N_{\mathrm{H,tot}}$ (and $A_{\mathrm{V,2D}}$) for all runs, times and three LOS considered. In the range $10^{21.5}$ cm*-2* $\leq$ $N_{\mathrm{H,tot}}$ $<$ 1023 cm*-2*, $\left\langle N_{\mathrm{H_{2}}}\right\rangle$ shows only moderate deviations between the individual runs of $\pm$ 0.2 dex, i.e. a factor of 1.6 in either direction and an almost linear trend in double-logarithmic representation. For this reason, we apply a powerlaw in this range (see Eq. 13 below), which smoothly connects to the range where $N_{\mathrm{H_{2}}}$ and $I_{\mathrm{CO}}$ are still correlated (regime (1)). 3. (3)

$N_{\mathrm{H,tot}}\geq 10^{23}$ cm*-2*:

For $N_{\mathrm{H,tot}}\geq 10^{23}$ cm*-2* ($A_{\mathrm{V,2D}}\geq 53.5$), the gas becomes (almost) fully molecular within typical deviations of $<$ 20% (black dashed curve in the right panel of Fig. 14). We thus approximate $N_{\mathrm{H_{2}}}$ by 0.5 $\times$ $N_{\mathrm{H,tot}}$ (see Eq. 14 below).

4.1 The final approach

Taking all results together, we obtain a new way to calculate $N_{\mathrm{H_{2}}}$ which is based solely on CO(1-0) and visual extinction observations:

1. (1)

For \frac{N_{\mathrm{H,tot}}}{\textrm{1 cm{}^{-2}}}<10^{21.5}\textrm{ or }A_{\mathrm{V,2D}}<1.7:

[TABLE] 2. (2)

For 10^{21.5}\leq\frac{N_{\mathrm{H,tot}}}{\textrm{1 cm{}^{-2}}}<10^{23}\textrm{ or }1.7\leq A_{\mathrm{V,2D}}<53.5:

[TABLE] 3. (3)

For \frac{N_{\mathrm{H,tot}}}{\textrm{1 cm{}^{-2}}}\geq 10^{23}\textrm{ or }A_{\mathrm{V,2D}}\geq 53.5:

[TABLE]

Our approach is applicable for a wide range column densities. In the intermediate- to high-$A_{\mathrm{V,2D}}$ range, it relies on the availability of extinction measurements in order to minimize the uncertainties introduced by the variations in $X_{\mathrm{CO}}$. At $10^{21.5}$ cm*-2* $\leq$ $N_{\mathrm{H,tot}}$ $<$ 1023 cm*-2*, $N_{\mathrm{H_{2}}}$ we introduce a novel super-linear relation between $A_{\mathrm{V,2D}}$ and $N_{\mathrm{H_{2}}}$ as the gas becomes increasingly molecular. Our approach also covers the low-$A_{\mathrm{V,2D}}$ range, which is accessible via high-sensitivity, long-exposure CO observations. Here, the weak dependence of $N_{\mathrm{H_{2}}}$ on $I_{\mathrm{CO}}$ in Eq. 12 – or inversely the strong dependence of $I_{\mathrm{CO}}$ on $N_{\mathrm{H_{2}}}$ – can be understood as a consequence of CO being underabundant (with respect to H2) at low column densities and then rapidly builds up to “catch up” with H2 towards higher column densities. This is also indicated in the left panel of Fig. 4 (see also Gong et al., 2018), where we indeed find a steep rise of $f_{\mathrm{CO}}$ above $f_{\mathrm{H_{2}}}$ $\simeq$ 0.5. Finally, we note that using somewhat different fitting ranges in Eqs. 12-14 results in changes for $N_{\mathrm{H_{2}}}$ of a few percent only.

4.2 Comparison of the new approach and previous approaches

In order to test the applicability of the suggested approach, we apply it to our simulations to determine the total H2 mass in the clouds and compare it to the classical approach of a constant $X_{\mathrm{CO}}$-factor. In addition, we compare it to the approach of Glover & Mac Low (2011), which suggest that the $X_{\mathrm{CO}}$-factor varies with the visual extinction averaged over the entire cloud, $\left\langle A_{\mathrm{V,2D}}\right\rangle$, as

[TABLE]

Motivated by the findings of Fig. 12, we use a value of 1.5 $\times$ 1020 cm*-2* (K km s*-1*)-1 for $X_{\mathrm{CO,0}}$ as well as for the approach with a constant $X_{\mathrm{CO}}$-factor. We constrain ourselves to the observable regions444We only focus on a particular threshold value as the total mass of H2 depends only weakly on it (Fig. 10). where $I_{\mathrm{CO}}$ $>$ 0.1 K km s*-1*.

In Fig. 15 we compare the inferred H2 masses to the actual H2 masses for the various MCs, directions and approaches. Overall, we find that our new approach (left panel) matches the actual H2 mass best: the estimated H2 masses range from 80% to 180% of the actual masses, i.e. deviations of a factor of 1.8 at most. This is about 2 times better than with the classical approach via a constant $X_{\mathrm{CO}}$-factor, which shows deviations of a factor of 3 – 4 in either direction (middle panel). The approach of Glover & Mac Low (2011) results in even larger deviations of up to one order of magnitude, in particular for the runs with feedback (right panel). Similar deviations, in particular the rather poor match for the approach of Glover & Mac Low (2011), were also found by Szűcs et al. (2016). We speculate that this is due to the fact that Glover & Mac Low (2011) use turbulent box simulations, which are difficult to compare to our simulations and might also be problematic when considering the convergence of the chemical abundances (Joshi et al., 2019). In addition, the authors only approximate the CO(1-0) intensity by a simplified radiative transfer scheme, which could also contribute to the observed differences.

We emphasize that our new approach works equally well for all three different situations considered in this paper, i.e. clouds with and without magnetic fields and stellar feedback. There is, however, a slight tendency for the runs without feedback to overestimate the H2 mass at early times and to underestimate it at late times (vice versa for the runs with feedback), which we investigate further below.

Next, we compare the (total) hydrogen masses of the various MCs obtained via Eq. 8 to their H2 masses, only taking into account pixels above a threshold of $A_{\mathrm{V,2D}}$ = 0.1, 1, and 10, respectively (Fig. 16). As expected, the ratio approaches unity with increasing threshold since the gas becomes increasingly more molecular (see Fig. 13). However, for typical observational sensitivity thresholds around $A_{\mathrm{V,2D}}$ = 0.1 – 1, the total mass is up to a factor of $\sim$ 7 higher than the H2 mass and should therefore not be used to estimate the molecular gas content. Even when considering only the densest parts of the clouds ($A_{\mathrm{V,2D}}$ $>$ 10), $M_{\mathrm{H_{2}}}$ would be overestimated by a few 10% up to a factor of 2. This is a consequence of the fact that the gas becomes (on average) fully molecular only above an $A_{\mathrm{V,2D}}$ of a few times 10 (our regime (3), see Fig. 13), which our approach accounts for by introducing the non-linear relation between N${}_{\mathrm{H_{2}}}$ and $A_{\mathrm{V,2D}}$ in regime (2) (Eq. 13).

4.3 The accuracy of the approach

There are two main sources of uncertainty of the suggested approach, first the simulations themselves (Section 4.4), and second the combination of different runs and times to obtain the final fits. In order to explore the effect of the latter, we calculate the ratio of the inferred and actual H2 column density in each pixel of the projected 2D-maps (again considering only pixels with $I_{\mathrm{CO}}>$ 0.1 K km s*-1*) as a function of $N_{\mathrm{H,tot}}$ ($A_{\mathrm{V,2D}}$) and plot the resulting H2-mass-weighted phase diagram of MC1-HD projected along the $z$-direction at $t_{\mathrm{evol}}$ = 3 Myr (Fig. 17). For each $N_{\mathrm{H,tot}}$-bin, we also calculate the mean of the inferred and actual H2 column densities, $\left\langle N_{\mathrm{H_{2},inferred}}\right\rangle$ and $\left\langle N_{\mathrm{H_{2},real}}\right\rangle$, respectively, and plot their ratio.

We find that our approach matches the actual H2 column density even on an individual pixel basis relatively well, which also holds for the other MCs, times and directions. In regime (1) and (2), the typical deviations are of the order of at most $\pm$ 0.5 dex, i.e. lower than a factor of $\sim$ 3. Towards the highest $N_{\mathrm{H,tot}}$ (regime (3)), where most of the H2 mass resides, the ratio is close to 1. These deviations are in rough agreement with those of the underlying correlations (Figs. 11 and 13) and the deviations of their mean values from the applied fit (Fig. 14).

Hence, overall the uncertainties in our approach appear to be dominated by variations in the intermediate column density range, i.e. regime (2) of the approach (Eq. 13), where a significant amount of the gas resides (Fig. 17). The typical deviations of the fitted from the real $N_{\mathrm{H_{2}}}$ in this column density regime (see also right panel of Fig. 14) also roughly match the overall accuracy of our approach of a factor of $\lesssim$ 1.8 for the total H2 mass. Also the time trends seen in the left panel of Fig. 15 are mainly caused by (time) variations in this regime. However, as we aim at providing a method which is valid to estimate H2 for various evolutionary stages, this uncertainty cannot be further reduced. We note, however, that the match for the $X_{\mathrm{CO}}$-factor (not shown, but see Fig. 11), would be significantly worse with typical deviations of a factor of 10 (or even more) in either direction.

4.4 Caveats

The fits to determine the H2 mass (Eqs. 12–14) are strictly seen only valid for an ISRF and a cosmic ray ionisation rate (CRIR) corresponding to solar neighborhood conditions, i.e. $G_{0}$ = 1.7 and CRIR = 1.3 $\times$ 10*-17* s*-1*. However, Clark & Glover (2015) show that the $X_{\mathrm{CO}}$-factor, and thus $I_{\mathrm{CO}}$, for MCs with masses around 104 M*☉* is barely affected when varying the ISRF and CRIR over two orders of magnitude. Also for more massive MCs ($\sim$ 105 M*☉*), they find only a relatively weak dependence on these quantities. In the densest and thus UV-shielded regions, where cosmic rays are expected to have a larger impact on the CO-H2 ratio (Bisbas et al., 2015; Bisbas et al., 2017), our approach relies on the extinction measurement, thus being less sensitive to the CRIR. Furthermore, Glover & Mac Low (2011) and Shetty et al. (2011a) show that also the influence of a moderately varying metallicity on $X_{\mathrm{CO}}$ is rather limited (see also Szűcs et al., 2016).

Uncertainties in the calculation of the UV shielding by surrounding gas can also affect the calculation of molecular abundances. The assumed converstion factor of 5.348$\times$10*-22* cm2 mag (Eq. 8) shows cloud-to-cloud variations of several 10% (see Fig. 2 of Bohlin et al., 1978), which would affect the dissociation of H2 and CO (i.e. regime (1)) but also directly the conversion of observed extinction into $N_{\mathrm{H,tot}}$ (regimes (2) and (3)). However, as shown by Glover & Mac Low (2011) and Szűcs et al. (2016), the H2 content depends only weakly on the extinction and varies only by a few 10% even when varying the ISRF by a factor of 10. Furthermore, the error introduced by the discretisation within the TreeCol/OpticalDepth algorithm into 48 directions (Eq. 2) is typically only of the order of 1 – 10% (Clark et al., 2012). Taken together, we thus speculate that the formulae given in Eqs. 12–14 might still be applicable even under somewhat different environmental conditions than in the solar neighborhood.

We also note that the used chemical network (Nelson & Langer, 1997; Glover & Mac Low, 2007; Glover et al., 2010) is simplified to allow for an efficient application in 3D, MHD simulations. However, comparison calculations with a more extended network based on Nelson & Langer (1999) and Glover et al. (2010, see the appendix of for details) show a reasonable agreement with the more simple network applied here. Furthermore, our results are also in agreement with the work of Levrier et al. (2012), which use a significantly more extended chemical network. The authors chemically post-process MHD simulations of Hennebelle et al. (2008) with the Meudon Code (Le Petit et al., 2006), which constrains them to a plane-parallel geometry and equilibrium chemistry. Despite this quite different approach, their results agree well with ours concerning the importance of density enhancements along the LOS (Section 3.3) and the observational deficit $\Delta f_{\mathrm{H_{2}}}$ (Section 3.4.1), which makes us confident about our chosen network.

Finally, we note that so far for the radiative stellar feedback we have only considered a single energy band (all photons with energies above 13.6 eV). We are currently working on including additional energy bands, which would allow for an even more detailed description of the dissociation processes of H2 and CO. Furthermore, for future simulations we plan to achieve an even higher spatial resolution in order to assure that in particular the CO content in our simulations is fully resolved (Joshi et al., 2019).

5 Conclusions

We present high-resolution (0.1 pc) simulations of molecular cloud formation including a live chemical network for H2 and CO as well as the necessary shielding processes. The simulations are part of the SILCC-Zoom project (Seifried et al., 2017) and include the galactic environment of the clouds, radiative stellar feedback and magnetic fields. We investigate six different simulations, 4 hydrodynamical runs, out of which 2 include star formation and ionisation feedback of young massive stars, and 2 magneto-hydrodynamical runs. In the simulations we can differentiate between the local visual extinction in each point, $A_{\mathrm{V,3D}}$, obtained directly from the 3D simulation data and the LOS-integrated visual extinction, $A_{\mathrm{V,2D}}$, as accessible in actual observations. In the following we list our main findings.

• •

The fraction of intrinsically CO-dark H2 gas (DGF) varies from 40% to 95%, with higher values for magnetised MCs. We show that differences in the DGF can be attributed to the structure of the clouds: clouds with a high amount of CO-dark gas have less well-shielded gas. The DGF, however, does not correlate with the total H2 mass.

• •

CO-bright gas is typically found at hydrogen nuclei densities above 300 cm*-3*, temperatures below 50 K and local visual extinctions, $A_{\mathrm{V,3D}}$ $\gtrsim$ 1.5, where 50 – 80% or more of the total hydrogen and carbon atoms are in the form of H2 and CO.

• •

CO-dark gas extends into the more diffuse (10 cm*-3* $\lesssim$ $n_{\mathrm{H,tot}}$ $\lesssim$ 103 cm*-3*) and moderately cool gas (a few $10$ K $\lesssim$ T $\lesssim$ a few 100 K). We speculate that this makes it difficult to probe the entire CO-dark gas with a single tracer. The typical $A_{\mathrm{V,3D}}$ of CO-dark gas ranges from 0.2 – 0.3 to about 1 – 1.5 independent of the presence or absence of either stellar feedback or magnetic fields.

• •

We demonstrate that with the LOS-integrated $A_{\mathrm{V,2D}}$, the conditions along the LOS cannot be determined properly. The actual distribution of the local visual extinction ($A_{\mathrm{V,3D}}$) along the LOS is broad and not in any way unique.

• •

Related to that, we show that up to $A_{\mathrm{V,2D}}$ $\simeq$ 5, pixels can be CO-bright and CO-dark, i.e. the DGF is not well constrained by the observable visual extinction. This can be attributed to different density – and thus $A_{\mathrm{V,3D}}$ – distributions along the LOS: Pixels with a high DGF have a rather uniform density distribution with $A_{\mathrm{V,3D}}$ $\lesssim$ 1 where CO is not formed. For CO-bright pixels, however, regions with strong density enhancements and locally well-shielded gas ($A_{\mathrm{V,3D}}$ $\gtrsim$ 1.5) are present along the LOS.

In addition, we produced synthetic CO(1-0) observations of our simulated molecular clouds using RADMC-3D.

• •

We show that about 15 – 65% of the H2 is in regions with CO(1-0) emission below an observational detection limit of 0.1 K km s*-1*, which amplifies the problem of intrinsically CO-dark gas in regions with detectable emission. This fraction increases only slightly to 20 – 75% when a detection limit of 1 K km s*-1* is used.

• •

We find a mean $X_{\mathrm{CO}}$-factor of $\sim$ 1.5 $\times$ 1020 cm*-2* (K km s*-1*)-1 in our simulations with significant variations of a factor up $\sim$ 4 in good agreement with other observational and theoretical works. Hence, using $X_{\mathrm{CO}}$ can result in significant errors in the estimated H2 masses of individual clouds.

In order to overcome the long-standing problem to determine the H2 content of MCs and to avoid the problem of CO-dark gas, we suggest a new approach to determine the H2 content of MCs under solar neighborhood conditions. The approach relies on observations of the CO(1-0) line transition and the visual extinction. The formulae given in Eqs. 12–14 present an approximation to the data obtained from all simulations considered here, which cover a variety of cloud conditions including and excluding both radiative stellar feedback and magnetic fields. Furthermore, the approach is applicable for a wide range of visual extinctions: from the low-extinction range ($A_{\mathrm{V,2D}}$ $\lesssim$ 1) covered by high-sensitivity, long-exposure CO observations, to the intermediate- and high-extinction range, where we introduce a novel a non-linear relation between $A_{\mathrm{V,2D}}$ and $N_{\mathrm{H_{2}}}$.

The total H2 cloud masses obtained with our new approach match the actual masses within a factor of at most 1.8 independent of whether feedback or magnetic fields are included or not. In contrast to that, the classical approach via a fixed $X_{\mathrm{CO}}$-factor results in deviations of up to a factor of 4. Moreover, our approach also allows us to calculate the H2 column density for individual pixels, i.e. on sub-pc scales, which is not possible with the $X_{\mathrm{CO}}$-factor. Here, we find typical deviations from the real H2 column density by less than a factor of 3, while the standard $X_{\mathrm{CO}}$-factor results in deviations by an order of magnitude and more.

Acknowledgements

The authors like to thank the anonymous referee for the comments which helped to significantly improve the paper. DS and SW acknowledge the support of the Bonn-Cologne Graduate School, which is funded through the German Excellence Initiative. DS, SH, SW and TGB also acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center SFB 956 “Conditions and Impact of Star Formation” (subprojects C5 and C6). SW and TGB acknowledge support via the ERC starting grant No. 679852 "RADFEEDBACK". The FLASH code used in this work was partly developed by the Flash Center for Computational Science at the University of Chicago. The authors acknowledge the Leibniz-Rechenzentrum Garching for providing computing time on SuperMUC via the project “pr94du” as well as the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu).

Appendix A Supplemental figures

A.1 Density and mass fractions

In the top row of Fig. 18 we show the mean H2 and CO densities in the dense gas, i.e only taking into account gas which has a density above 3.84 $\times$ 10*-22* g cm*-3* corresponding to a particle density of $n$ = 100 cm*-3* for $\mu$ = 2.3. The MHD clouds show somewhat lower H2 densities than the HD clouds, confirming the results of Section 3.2.1 that the MHD clouds are somewhat more diffuse. Similar holds true for the CO densities, which, due to the assumed fractional abundance of carbon atoms of 1.4 $\times$ 10*-4*, are roughly a factor of 10*-4* lower.

In the bottom row we show the mass fractions of H2 and CO in the entire zoom-in regions. The mass fractions relate to the DGF (Eq. 4) as

[TABLE]

As already indicated by the values of DGF $>$ 0 (Fig. 2) the mass fractions of CO are smaller than that of H2. Furthermore, due to the more diffuse structure of the MHD clouds, in general both mass fractions are lower than that of the clouds without magnetic fields.

A.2 Density-temperature phase diagrams

In Fig. 19 we show the H2-mass-weighted $n_{\mathrm{H,tot}}$-$T$-phase diagram of the total, CO-dark and CO-bright gas in the runs MC1-HD-FB and MC3-MHD.

A.3 CO spectra

In Fig. 20 we show the CO spectra of the entire zoom-in regions of the six different simulations for three different LOS at $t_{\mathrm{evol}}$ = 3 Myr. The spectra differ significantly between the different simulations but also for individual cloud when considering different LOS. The spectra show features which can be attributed to multiple velocity components along the LOS but also to optical depth effects towards the line center. As a detailed analysis of the spectra is beyond the scope of the paper, we defer it to a subsequent publication (Nürnberger et al. in prep.).