11institutetext:  Owens Valley Radio Observatory, California Institute of
Technology, MC 249-17, Pasadena, CA 91125, USA
11email<EMAIL_ADDRESS>22institutetext: Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA
91109-8099, USA
22email<EMAIL_ADDRESS>33institutetext: California Institute
of Technology, TAPIR, Mailcode 350-17, Pasadena, CA 91125, USA

# Constraining the H2 column densities in the diffuse interstellar medium
using dust extinction and H i data

R. Skalidis , 1122 P. F. Goldsmith 22 P. F. Hopkins 33 S. B. Ponnada 33

(Received; accepted)

###### Abstract

Context. Carbon monoxide (CO) is a poor tracer of H2 in the diffuse
interstellar medium (ISM), where most of the carbon is not incorporated into
CO molecules unlike the situation at higher extinctions.

Aims. We present a novel, indirect method to constrain H2 column densities
($N_{{H}_{2}}$) without employing CO observations. We show that
previously–recognized nonlinearities in the relation between the extinction,
$A_{V}({H}_{2})$, derived from dust emission and the H i column density
($N_{{H{\sc I}}}$) are due to the presence of molecular gas.

Methods. We employ archival $N_{{H}_{2}}$ data, obtained from the UV spectra
of stars, and calculate $A_{V}({H}_{2})$ towards these sight lines using 3D
extinction maps. The following relation fits the data: log $N_{{H}_{2}}$ =
$1.38742\leavevmode\nobreak\
\left(\log{A_{V}(H_{2}})\right)^{3}-0.05359\leavevmode\nobreak\
\left(\log{A_{V}(H_{2}})\right)^{2}+0.25722\leavevmode\nobreak\
\log{A_{V}(H_{2})}-20.67191$. This relation is useful for constraining
$N_{{H}_{2}}$ in the diffuse ISM as it requires only $N_{{H{\sc I}}}$ and dust
extinction data, which are both easily accessible. In $95\%$ of the cases, the
estimates produced by the fitted equation have deviations under a factor of
$3.5$. We construct a $N_{{H}_{2}}$ map of our Galaxy and compare it to the CO
integrated intensity ($W_{CO}$) distribution.

Results. We find that the average ratio ($X_{CO}$) between $N_{{H}_{2}}$ and
$W_{CO}$ is approximately equal to $2\times 10^{20}$ cm-2 (K km s-1)-1,
consistent with previous estimates. However, we find that the $X_{CO}$ factor
varies by orders of magnitude on arcminute scales between the outer and the
central portions of molecular clouds. For regions with $N_{{H}_{2}}$ $\gtrsim
10^{20}$ cm-2, we estimate that the average H2 fractional abundance,
$f_{{H}_{2}}$= 2 $N_{{H}_{2}}$ / (2$N_{{H}_{2}}$ \+ $N_{{H{\sc I}}}$), is
$0.25$. Multiple (distinct) largely atomic clouds are likely found along
high–extinction sightlines ($A_{V}$ $\geq 1$ mag), hence limiting
$f_{{H}_{2}}$ in these directions.

Conclusions. More than $50\%$ of the lines of sight with $N_{{H}_{2}}$ $\geq
10^{20}$ cm-2 are untraceable by CO with a $J$ = 1–0 sensitivity limit
$W_{CO}$ = $1$ K km s-1.

###### Key Words.:

ISM: abundances – (ISM:) dust, extinction – ISM: structure – Galaxy:
abundances – (Galaxy:) local interstellar medium – methods: data analysis

##  1 Introduction

Molecular hydrogen is the most abundant molecule in the Universe, rendering it
one of the most, if not the most, important chemical constituents. Molecular
hydrogen’s lack of a permanent electric dipole moment results in its
rotational levels being connected only by quadrupole transitions. Their
weakness, together with the large spacing of its rotational levels, results in
rotational emission from H2 being extremely weak at the temperatures typically
encountered in molecular regions and is only rarely observed in infrared
emission.

An alternative method to observe molecular hydrogen is through absorption
lines; molecular hydrogen has several Lyman–Werner absorption lines at UV
wavelengths (11.2 - 13.6 eV). The column density of H2 can be determined by
measuring the depths of the Lyman–Werner absorption lines against UV–bright
background sources (Savage et al., 1977; Savage & Sembach, 1996; Cartledge et
al., 2004; Shull et al., 2021; Mangum & Shirley, 2015).

The direct measurement of $N_{{H}_{2}}$, either through emission or
absorption, requires space observatories. The MIRI spectrometer onboard JWST
(Rieke et al., 2015; Bouchet et al., 2015; Kendrew et al., 2015; Boccaletti et
al., 2015; Wells et al., 2015; Glasse et al., 2015; Gordon et al., 2015;
Wright et al., 2023) is sensitive to H2 rotational transitions (e.g., Armus et
al., 2023), while the absorption lines can be detected in spectra obtained
with FUSE or HST (Savage et al., 1977; Bohlin et al., 1978; Shull et al.,
2000; Sofia et al., 2005; Sheffer et al., 2008; Rachford et al., 2002, 2009;
Gillmon et al., 2006; Shull et al., 2021). The necessity of space
observatories makes $N_{{H}_{2}}$ measurements expensive and thus relatively
sparse. For this reason, various indirect techniques have been developed to
estimate $N_{{H}_{2}}$. We summarize the most widely applied indirect methods
in the following sections.

###  1.1 Dust intensity at far infrared wavelengths

Emission from dust at far–infrared wavelengths ($\lambda\gtrsim 100$ $\mu$m)
traces cold molecular gas; the dust in dark and giant molecular clouds emits
thermal radiation as a gray–body. Using multi–wavelength observations of dust
emission, the dust spectral energy distribution (SED) can be obtained and
fitted analytically with the following free parameters (e.g., Hildebrand,
1983; Draine et al., 2007; Paradis et al., 2010, 2023; Kirk et al., 2010;
Martin et al., 2012; Planck Collaboration et al., 2014; Schnee et al., 2006;
Juvela et al., 2018): 1) dust temperature, and 2) dust grain emissivity,
which, assuming a dust-to-gas mass ratio, can be converted to $N_{{H}_{2}}$.
Dust SED fitting was the most common strategy applied to analyze Herschel
images to obtain $N_{{H}_{2}}$ towards several nearby interstellar clouds
(e.g., André et al., 2010; Miville-Deschênes et al., 2010; Palmeirim et al.,
2013; Könyves et al., 2010, 2015; Cox et al., 2016).

###  1.2 CO - H2 conversion factor ($X_{CO}$)

Carbon monoxide is an abundant molecule in the ISM. It is abundant in regions
where $N_{{H}_{2}}$ $\geq 10^{20}$ cm-2, as self-shielding there is sufficient
to prevent rapid photo-dissociation of CO by the background ISM radiation
field. Thus, when there is CO, there is also H2, although the opposite is not
always true (e.g., Grenier et al., 2005; Planck Collaboration et al., 2011b;
Barriault et al., 2010; Langer et al., 2010, 2014, 2015; Velusamy et al.,
2010; Pineda et al., 2013; Goldsmith, 2013; Skalidis et al., 2022; Madden et
al., 2020; Kalberla et al., 2020; Murray et al., 2018; Lebouteiller et al.,
2019; Glover & Smith, 2016; Seifried et al., 2020; Li et al., 2015).

The CO integrated intensity ($W_{CO}$) can be used as a proxy for
$N_{{H}_{2}}$ by employing the $X_{CO}$ factor, defined as

$N_{H_{2}}=X_{CO}W_{CO}.$ (1)

In our Galaxy, the average value of $X_{CO}$ is $\langle X_{CO}\rangle\approx
2\times 10^{20}$ cm-2 (K km s-1)-1 (Bolatto et al., 2013). When $X_{CO}$ is
measured, factor of two deviations are expected from the average value due to
the assumptions employed to estimate $N_{{H}_{2}}$ by various methods (Dame et
al., 2001; Bolatto et al., 2013).

The conversion of $W_{CO}$ to $N_{{H}_{2}}$ using the Galactic average
$X_{CO}$ is meaningful only in a statistical sense, and when it is applied on
Galactic scales. Within individual clouds in the interstellar medium (ISM),
$X_{CO}$ can vary by orders of magnitude (e.g., Pineda et al., 2008; Lee et
al., 2014; Pineda et al., 2013; Ripple et al., 2013) due to the sensitivity of
$X_{CO}$ to local ISM conditions, including gas density, turbulent line width,
and the interstellar radiation field (Visser et al., 2009; Glover & Clark,
2016; Gong et al., 2017).

###  1.3 H2 not traced by CO (CO-dark H2)

Diffuse regions of the ISM (defined as having $A_{V}$ $\lesssim 1$ mag, and
density $n\lesssim 100$ cm-3) may contain a significant fraction of molecular
hydrogen but relatively little CO compared to the abundance of carbon (Grenier
et al., 2005). The reason is that the formation of H2 formation takes place at
smaller extinctions than does that of CO due to the rapid onset of
self–shielding at lower values of the extinction.

In the transition phase where the H2 density has started rising but that of CO
is still small, oxygen and carbon are largely in atomic form. In this region,
then, H2 is primarily associated with ionized (C+) or neutral (C i) carbon,
and not with CO. Molecular hydrogen in regions that are undetectable in CO are
referred to as CO-dark H2. The mass of CO-dark H2 is estimated to be $\sim
30\%$ of the total (primarily H2) mass (Pineda et al., 2010; Kalberla et al.,
2020).

###  1.4 Dust extinction residuals

Dust in the ISM is well mixed with gas (Boulanger et al., 1996; Rachford et
al., 2009). Thus, the extinction and reddening produced by dust is
proportional to the total (atomic plus molecular hydrogen) gas column along
the LOS (Bohlin et al., 1978). All established methodologies (Sect.6) concur
that in regions where the gas is atomic, the dust reddening, $E(B-V)$, is
linearly correlated with $N_{{H{\sc I}}}$ (Liszt & Gerin, 2023; Lenz et al.,
2017; Nguyen et al., 2018; Shull & Panopoulou, 2023); changes in the dust or
hydrogen content cause $N_{{H{\sc I}}}$ / $E(B-V)$ to vary (Liszt, 2014b;
Shull & Panopoulou, 2023). At high Galactic latitudes, Lenz et al. (2017)
found that the linear correlation between $E(B-V)$ and $N_{{H{\sc I}}}$ holds
for $N_{{H{\sc I}}}$ $\lesssim 4\times 10^{20}$ cm-2. At higher H i column
densities, $E(B-V)$ increases non-linearly with $N_{{H{\sc I}}}$, behavior
plausibly attributed to the presence of molecular gas. This implies that the
residuals between the total dust extinction (or reddening) and the extinction
induced by dust mixed with H i gas should probe $N_{{H}_{2}}$ (Planck
Collaboration et al., 2011b; Paradis et al., 2012; Lenz et al., 2015; Kalberla
et al., 2020).

This strategy for estimating $N_{{H}_{2}}$ can only be applied if two
conditions are satisfied: 1) the $N_{{H{\sc I}}}$ / $E(B-V)$ ratio must be
constrained, and 2) a mapping function that converts the extinction residuals
to $N_{{H}_{2}}$ exists. $N_{{H{\sc I}}}$ / $E(B-V)$ has been constrained in
the past by various authors (e.g., Liszt, 2014b; Lenz et al., 2017; Nguyen et
al., 2018; Shull et al., 2021). The major challenge is to find an accurate way
to convert the extinction residuals to $N_{{H}_{2}}$. This is the focus of
this work.

###  1.5 The present work

In this paper we present a novel method to convert the dust extinction
residuals to $N_{{H}_{2}}$ and then construct a full-sky $N_{{H}_{2}}$ map of
our Galaxy without using CO observations. We use our $N_{{H}_{2}}$ estimates
to constraint some of the molecular gas properties of our Galaxy.

We calculate the residuals between the total extinction and the extinction
expected from dust mixed with H i gas towards lines of sight (LOSs) with
$N_{{H}_{2}}$ measurements; these $N_{{H}_{2}}$ measurements have been
obtained from absorption lines (Sect. 2). We find that the extinction
residuals are well correlated with $N_{{H}_{2}}$, and fit them with a
polynomial function. Then, we calculate the extinction residuals using full-
sky $N_{{H{\sc I}}}$ and $E(B-V)$ maps (Sect. 3). Finally, we convert the
extinction residuals to $N_{{H}_{2}}$ using our fitted function in order to
construct a full–sky $N_{{H}_{2}}$ map at $16\arcmin$ resolution.

We compare our constructed $N_{{H}_{2}}$ map with the CO intensity map of Dame
et al. (2001) (Sect. 4), reproducing the previously–reported global value of
$X_{CO}$. Our $N_{{H}_{2}}$ estimates are consistent with a previous work,
which combines dust extinction residuals and the $X_{CO}$ factor (Sect. 4.2).
We construct a full-sky map of the molecular fractional abundance, and show
that molecular hydrogen is less abundant than atomic hydrogen (Sect. 5). In
Sect. 5.2, we determine the proportion of CO–dark molecular hydrogen in
relation to the total amount of molecular hydrogen. We explore the $N_{H}$ /
$A_{V}$ ratio and compare against previous results (Sect. 6). In Sects. 7, and
8 we discuss our results in the context of the existing literature and present
our conclusions.

Figure 1: Extinction of dust mixed with molecular gas as a function of
$N_{{H}_{2}}$. The H2 column densities were determined from UV observations
towards stars, and the $A_{V}({H}_{2})$ were determined from the non-
linearities of the $E(B-V)$ versus $N_{{H{\sc I}}}$ relation (Eq. 3). Green
stars identify outlying points which were not considered in the fitting. At
high column densities, $A_{V}({H}_{2})$ is strongly correlated with
$N_{{H}_{2}}$, indicating that the nonlinear increase of dust reddening with
respect to $N_{{H{\sc I}}}$ is induced by the presence of molecular gas. The
orange curves were randomly sampled from the posterior distribution of the
MCMC fitting, representing the fitting uncertainties. The green curve
corresponds to the mean posterior profile that we employ to convert
$A_{V}({H}_{2})$ to $N_{{H}_{2}}$ (Eq. 6). Figure 2: Same as Fig. 1 but for
different $N_{{H{\sc I}}}$ /$E(B-V)$ ratios. The black points with errorbars
correspond to our formal measurements using the $N_{{H{\sc I}}}$ /$E(B-V)$
constraint of Lenz et al. (2017). Blue and cyan points correspond to
measurements obtained when we assume that $N_{{H{\sc I}}}$ / $E(B-V)$ $\times
10^{-21}$ is equal to 10 and 8 cm-2 mag-1 respectively. Approximately all data
points show statistically consistent results. Figure 3: Top panel: Mollweide
projection of our $N_{{H}_{2}}$ map. The dark blue color indicates regions
with $19\lesssim\log$ $N_{{H}_{2}}$ ( cm-2) $\lesssim 20$. Bottom panel:
$W_{CO}$ map from Dame et al. (2001). Our $N_{{H}_{2}}$ map was constructed by
calculating $A_{V}({H}_{2})$, using $E(B-V)$ and $N_{{H{\sc I}}}$ data, and
then converting to $N_{{H}_{2}}$ using the relation given in Eq. 7. The
$N_{{H}_{2}}$ map shows several relatively large angular scale structures. As
these are CO-dark H2 regions, they are undetectable in the $W_{CO}$ map.
(Sect. 5.2). In both maps, longitude increases to the left from the central
(vertical) line at $\ell=0\degr$, the grid spacing, shown by the white solid
lines, is $30\degr$, the angular resolution is $16\arcmin$, and NSIDE=1024.

##  2 Probing molecular gas using dust extinction

For atomic gas, Lenz et al. (2017) found that dust reddening correlates with
the column density of the atomic hydrogen as $N_{{H{\sc I}}}$ / $E(B-V)$
$=8.8\times 10^{21}$ cm-2 mag-1. This is our reference value and we show that
our results remain statistically the same when we consider variations about
this ratio (Sect. 2.2.2). When the gas becomes molecular, $E(B-V)$ increases
non-linearly with respect to $N_{{H{\sc I}}}$.

We assume a constant total-to-selective extinction $R_{V}=$ $A_{V}$ / $E(B-V)$
= 3.1 (Cardelli et al., 1989; Schlafly & Finkbeiner, 2011) and calculate the
visual extinction of dust mixed with atomic gas,

$A_{V}({H}\,{\sc I})=0.35\times 10^{-21}\leavevmode\nobreak\ N_{H_{\sc I}}.$
(2)

$R_{V}$ fluctuates throughout our Galaxy (e.g., Peek & Schiminovich, 2013;
Zhang et al., 2023). However, the uniformity of $R_{V}$ is implicit in our
method, because we employ dust reddening maps that have been constructed using
far-infrared multi-wavelength dust intensity data (Sect. 3, Appendix A). These
maps directly constrain the dust emission properties, such as dust temperature
and opacity, when some dust model is fitted to the data, and indirectly
$E(B-V)$, by assuming some extinction law and $R_{V}=3.1$ (e.g., Schlegel et
al., 1998). In addition, most of the existing constraints of $N_{{H{\sc I}}}$
/ $E(B-V)$ also rely on these dust reddening maps. Therefore, the uniformity
of $R_{V}$ inevitably propagates in our methods and accounts for some of our
uncertainties, which are calculated in Sect. 2.4.

We define the residuals between the total dust extinction and $A_{V}({H}{\sc
I})$ as

$A_{V}\left(H_{2}\right)=A_{V}({H})-A_{V}(H\,{\sc I}),$ (3)

where $A_{V}({H})$ denotes the total extinction from dust associated with both
atomic and molecular hydrogen. On the left hand side, we use H2 in the
parenthesis because we hypothesize, and verify below (Fig. 1), that these
residuals trace the visual extinction of dust mixed with molecular hydrogen.

Eq. 3 is the cornerstone of our methodology for estimating $N_{{H}_{2}}$.
$A_{V}({H}{\sc I})$ can be estimated using $N_{{H{\sc I}}}$ measurements (Eq.
2) from surveys such as HI4PI (HI4PI Collaboration et al., 2016), while
$A_{V}({H})$ can be extracted from publicly available full-sky dust extinction
maps (e.g., Schlegel et al., 1998; Planck Collaboration et al., 2016). The
next step in our strategy is to find an empirical relation for $N_{{H}_{2}}$
as a function of $A_{V}({H}_{2})$. Before proceeding we discuss some
pathologies in the definition of $A_{V}({H}_{2})$.

###  2.1 Pathologies in the definition of $A_{V}({H}_{2})$

Dust extinction is a positive definite quantity, but $A_{V}({H}_{2})$ can
become negative if $A_{V}({H}{\sc I})$ $>$ $A_{V}({H})$. This situation can be
encountered when either the signal-to-noise ratio (S/N) of the measurements is
low or if the value of $N_{{H{\sc I}}}$ / $E(B-V)$ deviates from its global
value; variations of this ratio propagate to $A_{V}({H}_{2})$, due to Eq. 3
(Sect. 2.2.2). We present a thought experiment to demonstrate how the use of a
global $N_{{H{\sc I}}}$ / $E(B-V)$ (Eq. 2) can affect the estimation of
$A_{V}({H}_{2})$.

Consider an ISM cloud with a local $N_{{H{\sc I}}}$ / $E(B-V)$ ratio smaller
than the Galactic value. We consider that the imaginary cloud has a given
$A_{V}({H})$, and gas is atomic. Thus, the total extinction of the cloud will
be produced by dust mixed with H i gas. If we used the local $N_{{H{\sc I}}}$
/ $E(B-V)$ of the cloud, we would derive that $A_{V}({H})$ $=$ $A_{V}({H}{\sc
I})$, hence $A_{V}({H}_{2})$ $=0$. However, in Eq. 3 we use the global value
from Lenz et al. (2017), which would yield for the cloud of our thought
experiment that $A_{V}({H})$ $<$ $A_{V}({H}{\sc I})$, and hence
$A_{V}({H}_{2})$ $<0$.

Thus, $A_{V}({H}_{2})$ can become negative due to the use of a global
$N_{{H{\sc I}}}$ / $E(B-V)$, which limits the detectability of our method.
However, this is something that we cannot overcome because local measurements
are limited. For this reason, for every measurement with $A_{V}({H}_{2})$
$<0$, we consider gas to be fully atomic, and it is not considered in the
analysis below.

###  2.2 Mapping $A_{V}({H}_{2})$ to $N_{{H}_{2}}$

We calculated $A_{V}({H}_{2})$ using Eq. 3 towards several lines of sight with
reliable direct $N_{{H}_{2}}$ determinations. The direct $N_{{H}_{2}}$
measurements were obtained by observing the Lyman–Werner absorption lines
towards background sources in our Galaxy. We extracted the $N_{{H}_{2}}$
measurements, with their corresponding observational uncertainties, from the
catalogues of Sheffer et al. (2008), Gillmon et al. (2006), Gudennavar et al.
(2012), and Shull et al. (2021); these catalogues also provide estimates for
$N_{{H{\sc I}}}$, which allow us to obtain $A_{V}({H}{\sc I})$ (Eq. 2). For
measurements with asymmetric observational uncertainties, we calculated their
average value.

The majority of the sources used for estimating $N_{{H}_{2}}$ are nearby – at
distances less than 1 kpc – and close to the Galactic plane. Therefore, the
column densities obtained from the absorption lines of these stars trace the
integrated gas abundances up to the distance of each object, and not the total
gas column. Therefore, we also need to calculate the dust extinction up to the
distance of each star ($d_{\star}$). Below, we explain how we obtained
fractional extinctions up to the distance of each star using the publicly
available 3D extinction map of Green et al. (2019).

For each object with $N_{{H}_{2}}$ measurement, we extracted $d_{\star}$ from
the catalogue of Bailer-Jones et al. (2021), which uses parallaxes from the
Gaia Data Release 3 (Gaia Collaboration et al., 2021). We then calculated the
dust extinction integral of each star as:

$E(B-V)_{\star}=\int_{0}^{d_{\star}}dz\leavevmode\nobreak\ E(B-V)(z),$ (4)

where $z$ is the line-of-sight distance between the observer and the star at
distance $d_{\star}$, and $E(B-V)(z)$ denotes the dust reddening at a given
distance $z$. We extracted the $E(B-V)(z)$ values from the 3D map (Bayestar19)
of Green et al. (2019) 111http://argonaut.skymaps.info/. This map provides a
sample of $E(B-V)(z)$ profiles drawn from a posterior distribution; for the
calculation of $E(B-V)_{\star}$ we used the most probable extinction profile.
Dust reddening values from the map of Green et al. (2019) are calibrated for
the SDSS bands. We converted to dust reddening by multiplying the output of
their map by 0.884, which is the standard color correction recommended by
those authors. Then, we converted $E(B-V)_{\star}$ to $A_{V}({H})$ by
multiplying by $R_{V}$.

Both $A_{V}({H})$ and $A_{V}({H}{\sc I})$ are constrained for the LOSs where
$N_{{H}_{2}}$ measurements exist from spectroscopic observations. This allows
us to calculate $A_{V}({H}_{2})$ for each LOS and compare with $N_{{H}_{2}}$.
Several objects with spectroscopically constrained $N_{{H}_{2}}$ are located
in the anti-center of our Galaxy, which are not covered by the extinction map
of Green et al. (2019). These measurements were not considered in the analysis
below.

Fig. 1 displays $A_{V}({H}_{2})$ versus $N_{{H}_{2}}$. We do not show
measurements with log $A_{V}({H}_{2})$ (mag) $\lesssim-1.2$ because they are
non-detections. We also note that some measurements, mostly with low S/N,
yielded negative $A_{V}({H}_{2})$. These were also discarded.

For $\log$ $N_{{H}_{2}}$ (cm-2) $\lesssim 19$, $A_{V}({H}_{2})$ remains
constant with respect to $N_{{H}_{2}}$. Most of the measurements there have
S/N $<$ 3 222The uncertainties are shown in logarithmic scale, while S/N
refers to a linear scale., and hence we cannot confidently determine whether
the flatness of the profile is physical or induced by noise in the data. We
observe a strong correlation between $A_{V}({H}_{2})$ and $N_{{H}_{2}}$ for
$\log$ $N_{{H}_{2}}$ $\geq$ 20\. This correlation strongly supports the
initial hypothesis that the extinction residuals, found in the total dust
extinction when the extinction induced by dust mixed with H i is subtracted,
yield the extinction of dust mixed with H2 (as given by Eq. 3).

The various quantities involved in our analysis ($E(B-V)$, $N_{{H}_{2}}$,
$N_{{H{\sc I}}}$) are measured with beams of differing size: a) $N_{{H}_{2}}$
is measured from UV absorption lines with a “pencil beam” defined by the star,
b) several $N_{{H{\sc I}}}$ measurements have been constrained by fitting the
Ly$\alpha$ profiles (e.g., as in Shull et al., 2021), which are also
characterized by pencil beams, while some of the $N_{{H{\sc I}}}$ measurements
were obtained from H i emission line data, where the resolution is $16\arcmin$
(HI4PI Collaboration et al., 2016), and c) the resolution of the 3D extinction
map of Green et al. (2019) varies from $3.4\arcmin$ to $13.7\arcmin$. Some of
the scatter of the points in Fig. 1 is a result of this effect, but we note
that the beam mismatch seems to be less important for diffuse and extended ISM
clouds (e.g., Pineda et al., 2017; Murray et al., 2018). Although we cannot
make a correction for this issue, we include it when we calculate the
confidence levels of our method (Sect. 2.4).

####  2.2.1 Observational uncertainties in $A_{V}({H}_{2})$

There are two different sources of uncertainty that affect the estimates of
$A_{V}({H}_{2})$: 1) Uncertainties coming from the 3D extinction map of Green
et al. (2019), and 2) distance uncertainties of stars.

To evaluate the dust extinction uncertainties, we sampled 300 random values
from the posterior distribution of the Green et al. (2019) map at the most
probable distance of each star ($d_{\star}$). We calculated the extinction
spread ($\sigma_{ext}$) at $d_{\star}$.

Parallax uncertainties affect the distance estimates of a star. We consider
the distance upper and lower limits denoted as $d_{\star,+}$ and $d_{\star,-}$
respectively. From the map of Green et al. (2019), we calculated the most
probable extinction at $d_{\star,+}$, and at $d_{\star,-}$, denoted as
$E(B-V)$ ($d_{\star,+}$) and $E(B-V)$ ($d_{\star,-}$),respectively. We then
calculated the differences $\sigma_{d+}$ = $E(B-V)$ ($d_{\star,+}$) - $E(B-V)$
($d_{\star}$) and $\sigma_{d-}$ = $E(B-V)$ $d_{\star}$ \- $E(B-V)$
($d_{\star,-}$). The average ($\sigma_{d}$) between $\sigma_{d+}$ and
$\sigma_{d-}$ measures the extinction uncertainties due to parallax
uncertainties. We take the total dust extinction uncertainty to be the
quadratic sum of $\sigma_{ext}$ and $\sigma_{d}$. These uncertainties are
shown as the error bars in Fig. 1 and discussed in the following.

####  2.2.2 Uncertainties from $N_{{H{\sc I}}}$ / $E(B-V)$ variations

At high Galactic latitudes $\left(|b|\geq 60\degr\right)$, Lenz et al. (2017)
found that the correlation between $E(B-V)$ and $N_{{H{\sc I}}}$ is very
tight, with intrinsic variations only $0.015$ mag. However, as shown recently
by Shull & Panopoulou (2023), the uncertainties in the $E(B-V)$ versus
$N_{{H{\sc I}}}$ relation are larger than what is inferred by Lenz et al.
(2017), with variations being dominated by systematic uncertainties in the
employed datasets.

Lenz et al. (2017) used the $E(B-V)$ map of Schlegel et al. (1998), and
$N_{{H{\sc I}}}$ estimated using only emission line data from HI4PI
Collaboration et al. (2016), assuming that the line is optically thin.
However, the $E(B-V)$ values of Planck Collaboration et al. (2014) are
systematically larger than the values of Schlegel et al. (1998) at high
Galactic latitudes (Shull & Panopoulou, 2023). Thus, if the Planck map were
employed, the ratio between $N_{{H{\sc I}}}$ and $E(B-V)$ would be biased
toward smaller values than what would have been inferred from the Schlegel et
al. (1998) map. In addition, beam effects in the H i emission data can induce
significant scatter in the obtained $N_{{H{\sc I}}}$ constraints. Several of
these uncertainties are summarized in Shull & Panopoulou (2023), who found
that $N_{{H{\sc I}}}$ / $E(B-V)$ $\times 10^{-21}$ ranges from 8 to 10 cm-2
mag-1; this range is orders of magnitude larger than what is claimed by Lenz
et al. (2017), and represents the formal uncertainty in $N_{{H{\sc I}}}$ /
$E(B-V)$. We explored how $N_{{H{\sc I}}}$ / $E(B-V)$ variations, which
propagate to our calculations through Eq. 2, affect our results.

Fig. 2 shows $A_{V}({H}_{2})$ versus $N_{{H}_{2}}$ for every measurement,
except for the outliers (green stars), shown in Fig. 1. Black points
correspond to $A_{V}({H}_{2})$ when we use the Lenz et al. (2017) relation,
$N_{{H{\sc I}}}$ / $E(B-V)$ $\times 10^{-21}$ = $8.8$ cm-2 mag-1; the
errorbars are calculated as explained in Sect. 2.2.1. The blue and cyan points
correspond to the obtained $A_{V}({H}_{2})$ when we assume that $N_{{H{\sc
I}}}$ / $E(B-V)$ $\times 10^{-21}$ is equal to 10 and 8 cm-2 mag-1
respectively.

For $\log$ $N_{{H}_{2}}$ $>20.56$, the colored (black and cyan) points are
very close to the black points. In this regime, gas is mostly molecular, and
thus the contribution of the dust mixed with atomic gas to the total dust
extinction is minor. Thus, the uncertainty of $N_{{H{\sc I}}}$ / $E(B-V)$ has
a weak effect in high-extinction regions. When $N_{{H}_{2}}$ is low,
deviations between black and colored points become more prominent, because the
relative abundance of atomic gas is higher there. However, even there, the
offset between the black and the colored points is statistically insignificant
for the vast majority of the points (Fig. 2). There are only six exceptions,
but this is only a minor fraction of the sample.

We conclude that our analysis remains robust in relation to variations of
$N_{{H{\sc I}}}$ / $E(B-V)$, which implies that uncertainties coming from the
3D extinction maps (Green et al., 2019) or Gaia distances (Sect. 2.2.1) are
more important. However, variations in $N_{{H{\sc I}}}$ / $E(B-V)$ likely
contribute to the scatter of the observed $A_{V}({H}_{2})$ – $N_{{H}_{2}}$
relation (Fig. 1). Our estimated confidence intervals include all the
aforementioned uncertainties (Sect. 2.4).

###  2.3 Fitting the data

We fitted a polynomial to the $A_{V}({H}_{2})$ – $N_{{H}_{2}}$ relation using
Bayesian analysis. We assumed uniform priors, and calculated the following
likelihood:

$\ln{\mathcal{L}}=-\sum_{i=1}^{N}\left[\ln{\sqrt{2\pi}}+\ln{\sigma}+\frac{\left(y_{i}-\tilde{y}\right)^{2}}{2\sigma^{2}}\right],$
(5)

where $\sigma$ is the quadratic sum of the observational uncertainties of
$\log$ $A_{V}({H}_{2})$, and $\log$ $N_{{H}_{2}}$, $y_{i}$ corresponds to the
measured $A_{V}({H}_{2})$, $\tilde{y}$ is the intrinsic model, which we assume
to be a polynomial, and $i$ is the measurement index.

We sampled the parameter space with the ”emcee” Markov Chain Monte Carlo
(MCMC) sampler (Foreman-Mackey et al., 2013). Measurements shown as green
stars in Fig. 1 were not included in the fit because they are variable sources
(Sect. 2.2.1). However, the fit changes only slightly even if they are
included.

The best–fit polynomial equation is

$\displaystyle\log A_{V}(H_{2})=\leavevmode\nobreak$ $\displaystyle
0.03077\leavevmode\nobreak\ \left(\log
N_{H_{2}}\right)^{3}-1.60469\leavevmode\nobreak\ \left(\log
N_{H_{2}}\right)^{2}$ (6) $\displaystyle+27.80868\leavevmode\nobreak\ \log
N_{H_{2}}-161.06073,$

with uncertainties for the fitted coefficients (from highest to smallest power
of the logarithm) being $0.00769$, $0.42823$, $7.88171$, and $49.91439$. The
spread of the fit, shown by the orange lines in Fig. 1, grows for $\log$
$N_{{H}_{2}}$ (cm-2) $<19$, making our fit unreliable in that range. This
determines the sensitivity threshold of our method.

In practice, $N_{{H}_{2}}$ is usually the unknown, while $A_{V}({H}_{2})$ can
be measured from dust extinction and $N_{{H{\sc I}}}$ observational data.
$N_{{H}_{2}}$ can be estimated from $A_{V}({H}_{2})$ using the following
analytical relation, which also fits our data well,

$\displaystyle\log\leavevmode\nobreak\ N_{H_{2}}=\leavevmode\nobreak$
$\displaystyle 1.38742\leavevmode\nobreak\ \left(\log
A_{V}({H_{2}})\right)^{3}-0.05359\leavevmode\nobreak\ \left(\log
A_{V}({H_{2}})\right)^{2}+$ (7) $\displaystyle+0.25722\leavevmode\nobreak\
\log A_{V}({H_{2}})-20.67191,$

The fitting uncertainties of the coefficients, from highest to smallest power
of the logarithm, are $0.17822$, $0.17277$, $0.05248$, and $0.01555$. The
above equation maps $A_{V}({H}_{2})$ to $N_{{H}_{2}}$ and we use it to
construct a full-sky $N_{{H}_{2}}$ map of our Galaxy (Sect. 3).

####  2.3.1 Fitting different polynomials

We experimented with the degree of the polynomial used to fit Eq. 7. For even
polynomials, we found that $N_{{H}_{2}}$ is reduced at high $A_{V}({H}_{2})$,
because the coefficient of the highest power of the fitted polynomial is
always negative, hence making the curve concave downward. When we forced the
coefficient of the highest-power term to be positive (by employing this
constraint in the prior distribution), the reduction of $N_{{H}_{2}}$
disappeared. However, we wish to keep a flat distribution of priors over the
entire domain range (uninformative priors). We thus decided to work with
polynomials whose highest-power term is odd because they did not show any
reduction of $N_{{H}_{2}}$ at high $A_{V}({H}_{2})$.

A linear polynomial would not reproduce the non-linear behaviour of the data
at $\log$ $N_{{H}_{2}}$(cm-2) $\lesssim 20$, and $\log$ $A_{V}({H}_{2})$ (mag)
$\lesssim-1$ (Fig 1). As a result, linear models would underestimate
$N_{{H}_{2}}$ in low-extinction regions. We conclude that a cubic polynomial
(as shown in Eq. 7) is the simplest function that accurately represents the
data, while higher-degree polynomials overfit the data.

####  2.3.2 Outliers in the $A_{V}({H}_{2})$ \- $N_{{H}_{2}}$ relation

In the fitting process, we excluded the four measurements shown as green stars
(Fig. 1). According to the SIMBAD database, one of the outliers (identified as
HD200120) is a Be star – this object also appears in the Be star catalogue of
the IAU (Jaschek & Egret, 1982), hence further supporting the validity of the
spectral classification of this object – while there is a T-Tauri (HD40111)
and a $\beta$ Cephei variable (HD172140) star. Both Be and T-Tauri stars can
show some variability in both photometric and spectroscopic observations: Be
stars eject material, due to their rapid rotation, while the brightness of
T-Tauri objects can vary significantly, due to their high accretion rate, even
within months. $\beta$ Cephei are B-type stars whose surface pulsates owing to
their rapid rotation. The last object that was considered an outlier (HD
164816) is a Be star according to the SIMBAD database. However, accurate
spectroscopic observations suggest that it is a binary system (Sota et al.,
2014). In either case, some degree of variability is expected, which explains
why this measurement is an outlier. We consider the $N_{{H}_{2}}$ measurements
towards these stars to be untrustworthy and did not include them in the fit.
But even if we include these measurements, the fit changes insignificantly.

###  2.4 Confidence intervals in the estimated $N_{{H}_{2}}$

Our fitted model (Eq. 7) allows us to estimate $N_{{H}_{2}}$ using H i and
$E(B-V)$ data. But, the estimates do not come without uncertainties.

The uncertainties of the fit (shown by the orange lines in Fig. 1) are much
smaller than the variance of the points about the fitted curve. These measured
variations could be induced by several factors, such as variations in
$N_{{H{\sc I}}}$ / $E(B-V)$ (Sect. 2.2.2), $R_{V}$ variations, metallicity
gradients, H i emission line optical depth effects (Sect.7), or due to the
beam mismatch of the employed data, unaccounted for in our model. Therefore,
the fitting uncertainties alone overestimate our confidence on the
$N_{{H}_{2}}$ estimates derived by Eq. 7.

To place reasonable confidence levels, we assumed that the depicted
$A_{V}({H}_{2})$ spread (Fig. 1) represents the intrinsic spread in our
Galaxy. We calculated the linear scale ratio ($r$) between the observed
$A_{V}({H}_{2})$ with the value obtained from Eq. 7 for every data point with
S/N $\geq 3$. From the cumulative distribution function of $r$ – the
distribution of $r$ is asymmetric with skewness equal to 1.38 – we obtained
the following probabilities: for $68\%$ of the measurements
$r\leavevmode\nobreak\ \epsilon\leavevmode\nobreak\ [0.75,2.48]$, for $95\%$
$r\leavevmode\nobreak\ \epsilon\leavevmode\nobreak\ [0.46,3.54]$, and for
$99\%$ $r\leavevmode\nobreak\ \epsilon\leavevmode\nobreak\ [0.30,4.78]$. This
implies that in $95\%$ of the cases, the ”true” $N_{{H}_{2}}$ should not
deviate by a factor greater than $3.5$ from the values obtained from the fit
(Eq. 7).

Figure 4: Mollweide projection of our constructed $X_{CO}$ map. Our global
value for $X_{CO}$ is $2\times 10^{20}$ K km s-1 cm-2, which is consistent
with previous estimates. $X_{CO}$ varies by orders of magnitude within a few
tens of arcminutes. In the surroundings of molecular regions, $X_{CO}$ $\sim
10^{21}$ cm-2 (K km s-1)-1, while in the inner parts of clouds it decreases to
$X_{CO}$ $\sim 5\times 10^{19}$ cm-2 (K km s-1)-1. In the Galactic center,
$X_{CO}$ increases significantly due to the high optical depths in that
region. The angular resolution is $16\arcmin$, and NSIDE=1024.  Figure 5: 2D
probability density functions relating the different parameters studied in
this work (Sect. 4.1). Figure 6: $N_{{H}_{2}}$ map constructed using the
method presented here (upper panel) and using that in Kalberla et al. (2020)
(lower panel). The resolution of both maps is $2\degr$ and NSIDE=1024.

##  3 Full-sky $N_{{H}_{2}}$ map of our Galaxy

We use the $N_{{H{\sc I}}}$ data from the HI4PI survey (HI4PI Collaboration et
al., 2016) and the dust extinction map of Schlegel et al. (1998) – the values
of this map are multiplied by 0.884 to account for systematic uncertainties
(Schlafly et al., 2010; Schlafly & Finkbeiner, 2011) – to construct a full-sky
$N_{{H}_{2}}$ map of our Galaxy. For comparison, we use the Planck extinction
map. Overall, the Planck map shows enhanced extinction at high Galactic
latitudes; this is also noted by Shull & Panopoulou (2023). The various
estimated quantities for the molecular hydrogen do not change by more than a
factor of two when the Planck map is used (Appendix A).

For the construction of the full-sky $N_{{H}_{2}}$ map, we calculated
$A_{V}({H}_{2})$ using the $N_{{H{\sc I}}}$ and the $E(B-V)$ maps (Eq. 3), and
then converted to $N_{{H}_{2}}$ (Eq. 7). The resolution of the extinction map
($6.1\arcmin$) is higher than that of the $N_{{H{\sc I}}}$ data
($16.2\arcmin$). We thus smoothed the extinction map to the resolution of the
H i data. Fig. 3 (top panel) shows a mollweide projection of our constructed
$N_{{H}_{2}}$ map for the Milky Way. We only visualize regions with log
$N_{{H}_{2}}$ $\geq 19$ because below this limit our calibration relation is
dominated by noise (Fig. 1).

We identify several molecular clouds of the Gould Belt, including the Polaris
Flare $\left(\ell,b\sim 125\degr,+30\degr\right)$, Taurus $\left(\ell,b\sim
170\degr,-15\degr\right)$, and Orion $\left(\ell,b\sim
209\degr,-19\degr\right)$. The North Polar Spur appears as a large-scale
molecular feature centered at Galactic coordinates $\ell,b\sim
15\degr,15\degr$ and extending up to $b\sim 45\degr$. We observe many small-
angular-scale diffuse molecular clouds in the Southern Galactic hemisphere,
while in the Northern hemisphere molecular structures seem to be coherent over
scales of many degrees.

We hardly see any molecular gas structures with $\log$ $N_{{H}_{2}}$ ( cm-2)
$\geq 19$ above the Galactic plane, for $|b|\geq 60\degr$. This, however,
seems to be the case only when we use the Schlegel et al. (1998) extinction
map. When we construct the $N_{{H}_{2}}$ map using extinctions from Planck
Collaboration et al. (2020), we see some molecular structures extending up to
the Northern and Southern poles, and weaker emission from both extended and
small-scale clouds at higher Galactic latitudes (Fig. 15). In both
$N_{{H}_{2}}$ maps (Figs. 3, and 15), the majority of the molecular gas seems
to lie within $|b|\lesssim 45\degr$.

In the Galactic plane, extinctions are larger than the extinction range used
for the fitting of the $N_{{H}_{2}}$ – $A_{V}({H}_{2})$ relation (Eq. 7).
Therefore, the estimated H2 column densities in the Galactic plane were
derived by extrapolating the fitted relation to larger values of $A_{V}$.

##  4 Comparison with previous results

###  4.1 Comparing $N_{{H}_{2}}$ and $W_{CO}$

Fig. 3 shows our $N_{{H}_{2}}$ map (top panel) and the CO (J=1-0) integrated
intensity, $W_{CO}$, map (bottom panel) from Dame et al. (2001). The two maps
are well correlated, but the observed structures in the $N_{{H}_{2}}$ map are
more extended than in $W_{CO}$. In this section our analysis is focused on CO-
bright H2 regions, while in Sect. 5.2 we study the properties of CO-dark H2,
which extends to higher Galactic latitudes than CO-bright H2.

Fig. 4 shows the full-sky $X_{CO}$ (Eq. 1) map obtained with our $N_{{H}_{2}}$
measurements and $W_{CO}$ from Dame et al. (2001). We calculated $X_{CO}$
toward LOSs with $W_{CO}$ $\geq 1$ K km s-1, which is the noise level in the
CO data. We find that the average value of $X_{CO}$ is approximately equal to
2 $\times$ 1020 cm-2 (K km s${}^{-1})^{-1}$, which matches exactly with the
previously reported value (Bolatto et al., 2013; Liszt & Gerin, 2016).

We observe that in many regions $X_{CO}$ varies by orders of magnitude within
a few arcminutes. The small-scale variations $X_{CO}$ are more prominent in
regions above the Galactic plane and in the Galactic plane but only for
$270\degr>\ell>90\degr$. $X_{CO}$ tends to be high ($\log$ $X_{CO}$ $\gtrsim
20.5$) in the surroundings of molecular regions, while it decreases
significantly ($\log$ $X_{CO}$ $\approx 19.5$) in the inner parts of molecular
clouds.

$X_{CO}$ is enhanced in the surroundings of clouds because CO molecules are
more effectively photodissociated there. The low CO abundance implies that
$W_{CO}$ is small, hence $X_{CO}$ becomes high. On the other hand, in the
inner parts of the molecular clouds, the total column density is greater,
reducing the photodestruction rate of molecules due to both shielding by dust
and by self-shielding, and the CO abundance increases. The enhanced CO
abundance implies that $W_{CO}$ increases and thus $X_{CO}$ decreases. In the
Galactic plane, $X_{CO}$ becomes large because $N_{{H}_{2}}$ is maximum there
and $W_{CO}$ saturates, due to the high optical depth of the $J=1-0$
transition usually employed to determine $W_{CO}$.

The observed behaviour is consistent with past observations, and numerical
simulations (Bolatto et al., 2013). Our results demonstrate that CO is a good
tracer of $N_{{H}_{2}}$ only for LOSs where CO is the dominant form of carbon,
and the CO rotational line is not saturated.

We explore quantitatively the correlation of the CO observables ($X_{CO}$,
$W_{CO}$) with $A_{V}$ and $N_{{H}_{2}}$. Fig. 5 shows the 2D probability
density functions of all possible combinations of $A_{V}$, $X_{CO}$, $W_{CO}$,
and $N_{{H}_{2}}$. Measurements with $N_{{H}_{2}}$ $>10^{21}$ cm-2, and
$A_{V}$ $>5.5$ mag have been derived by extrapolating our fitted equation,
hence no strong conclusions are made for this regime. The following paragraphs
summarize the results of the comparison.

$A_{V}$ – $X_{CO}$: The majority of points are clustered around $X_{CO}$
$\approx 10^{20}$ cm-2 (K km s${}^{-1})^{-1}$. The relation looks relatively
flat in this regime, although there is a slight anticorrelation; this is more
evident in the low-probability contours, which are concave upward. The
anticorrelation is due to the presence of CO-dark H2 gas in low-extinction
regions (Seifried et al., 2020; Borchert et al., 2022). For $\log$ $A_{V}$
$\gtrsim 1$, $X_{CO}$ increases with $A_{V}$ because the CO line saturates.
The critical extinction that marks the CO line saturation is log $A_{V}$
$\approx 0.6\Rightarrow A_{V}\approx 4$ mag, which is consistent with previous
estimates toward the Perseus molecular cloud (Pineda et al., 2008; Lee et al.,
2014). The local ISM conditions, such as density and intensity of the
background radiation field, can change the value of $A_{V}$ where CO
saturates. For this reason, there is an extinction range where a non-linear
increase of $X_{CO}$ is observed (e.g., Lombardi et al., 2006). Similar
behaviour between $X_{CO}$ and $A_{V}$ is also observed in numerical
simulations (Borchert et al., 2022).

$A_{V}$ – $W_{CO}$: Most of the points lie in the range $0.3\lesssim$ $\log$
$W_{CO}$ (K km s-1) $\lesssim 1$. $W_{CO}$ varies by orders of magnitude for
$\log$ $A_{V}$ $\lesssim 1$, while the dispersion in $W_{CO}$ decreases
significantly for $\log$ $A_{V}$ $\geq 1$. The scaling between $A_{V}$ and
$W_{CO}$ transitions at log $A_{V}$ $\approx 0.6$, and $W_{CO}$ converges to a
maximum value $W_{CO}$ $\approx 100$ K km s-1. This convergence is due to the
increase in the optical depth of the CO line, which becomes optically thick,
and is consistent with previous observations (Pineda et al., 2008; Lee et al.,
2014).

$A_{V}$ – $N_{{H}_{2}}$: $N_{{H}_{2}}$ increases non-linearly with respect to
$A_{V}$ for $\log$ $A_{V}$ $\lesssim 0$, because H i is transformed to H2; the
transition takes place for log $N_{H}$ (cm-2) in the range $\approx$ 20.1 -
20.8 (Savage et al., 1977; Gillmon et al., 2006; Bellomi et al., 2020). If we
assume that $A_{V}$ / $N_{H}$ $=5.34\times 10^{-22}$ cm2 mag (Bohlin et al.,
1978), then we obtain that $-1.2\lesssim$ $\log$ $A_{V}$ (mag) $\lesssim-0.5$,
which is consistent with the extinction range where we observe the nonlinear
increase in $N_{{H}_{2}}$. For $\log$ $A_{V}$ $\gtrsim 0$, the correlation
between $N_{{H}_{2}}$ and $A_{V}$ becomes quasi-linear. The positive
correlation between $A_{V}$ and $N_{{H}_{2}}$ reflects the fact that we see
more gas when there is more dust and that the hydrogen is almost entirely
molecular. We explore the dust-to-gas ratio in more detail in Sect. 6. In
agreement with previous studies (Liszt & Gerin, 2023), we find that only a few
sightlines have appreciable molecular gas content, $\log$ $N_{{H}_{2}}$ (
cm-2) $\gtrsim 19$ when $\log$ $A_{V}$ (mag) $\lesssim-0.5$.

$N_{{H}_{2}}$ – $W_{CO}$: The observed behavior in this figure is similar to
that of $A_{V}$ – $W_{CO}$. For $\log$ $N_{{H}_{2}}$ (cm-2) $>22$, $W_{CO}$
converges to a maximum value because the line saturates. For log $N_{{H}_{2}}$
(cm-2) $<21.5$, the two quantities are weakly correlated: $W_{CO}$ varies by
four orders of magnitude for $20\lesssim$ log $N_{{H}_{2}}$ $\lesssim 21$. The
rapid increase in $W_{CO}$ as function of $N_{{H}_{2}}$ reflects the
conversion of atomic carbon to CO and the resulting rise in fractional CO
abundance. The most commonly-found sightlines have $W_{CO}$ $\approx 10$ K km
s-1 and $N_{{H}_{2}}$ $\approx 10^{21}$ cm-2, which explains why the global
$X_{CO}$ is close to $10^{20}$ cm-2 (K km s-1)-1.

$W_{CO}$ – $X_{CO}$: For $\log$ $W_{CO}$ $\lesssim 0.6$, $X_{CO}$ is
anticorrelated with $W_{CO}$ due to Eq. 1. For $\log$ $W_{CO}$ $\approx 0.6$,
the CO line saturates and $X_{CO}$ becomes independent of $W_{CO}$, which
converges to 100 K km s-1.

Figure 7:  Difference between our $N_{{H}_{2}}$ estimates and those of
Kalberla et al. (2020).

###  4.2 Comparison to the $N_{{H}_{2}}$ map of Kalberla et al. (2020)

Kalberla et al. (2020) used the non-linear deviations in the $E(B-V)$ –
$N_{{H{\sc I}}}$ relation to estimate $N_{{H}_{2}}$. Our method is similar to
that of Kalberla et al. (2020), although with some key differences that we
discuss below.

####  4.2.1 Summary of Kalberla et al. (2020) method

H i emits at various Doppler-shifted frequencies (velocities). Each velocity
component is considered to be a distinct cloud along the LOS. Kalberla et al.
(2020) decomposed the H i emission spectra into distinct Gaussian components
333The decomposition of the emission spectrum into Gaussians has been widely
applied in the past (e.g., Miville-Deschênes et al., 2002, 2017; Heiles &
Troland, 2003a; Murray et al., 2014; Kalberla & Haud, 2015, 2018). There are
some uncertainties that can affect the decomposition of the emission spectra:
1) Contributions from non-Gaussian intensity peaks cannot be adequately
captured. 2) Decomposing a spectrum into Gaussians is an ill-defined problem
because different numbers of Gaussians can fit a spectrum equally well. 3)
Distinct clouds that emit at frequencies with small Doppler shifts can appear
as a single component in emission. Some of the aforementioned problems have
been treated with the development of sophisticated decomposition techniques
that require minimum user input (Lindner et al., 2015; Marchal et al., 2019;
Riener et al., 2019). . Each Gaussian is characterized by an amplitude and a
dispersion: the amplitude is proportional to the total H i emitting gas, while
the spread ($\sigma_{u}$) represents the internal velocity dispersion of the H
i gas, consisting of a thermal and turbulent component.

Kalberla et al. (2020) assigned an effective temperature ($T_{\rm D}$) - in
their paper this is referred to as the Doppler temperature - to each
decomposed Gaussian component, where $T_{\rm D}$ represents the total width of
the Gaussian component. They calculated the dust extinction expected from each
H i Gaussian component, by assuming a constant dust-to-gas ratio (Lenz et al.,
2017). Then, they added the $E(B-V)$ contribution from all Gaussian components
to estimate the total dust extinction induced by dust mixed with H i gas. They
compared their estimated H i \- based extinction to the total extinction,
using the map of Schlegel et al. (1998), and attributed the residuals to the
presence of molecular gas.

####  4.2.2 Similarities and differences with our method

The similarities between our method and that of Kalberla et al. (2020) are the
following: 1) both methods rely on the assumption of a universal dust-to-gas
ratio and on the hypothesis that nonlinear deviations in the $E(B-V)$ /
$N_{{H{\sc I}}}$ relation trace the molecular gas content. 2) Both methods
employ a mapping function that converts the $E(B-V)$ / $N_{{H{\sc I}}}$ non-
linearities to $N_{{H}_{2}}$.

The major difference between our method and that of Kalberla et al. (2020)
lies on how the mapping function was obtained. Kalberla et al. (2020) derived
their mapping function by minimizing the extinction nonlinearities from the
linear relation through bootstrapping. We derived our mapping function from
independent data by using $N_{{H}_{2}}$ measurements from UV spectra of
background objects (Sect. 2). In addition, Kalberla et al. (2020) fitted their
calibration function in regions where CO is undetectable. In CO-emitting
regions, they estimated $N_{{H}_{2}}$ by assuming a constant $X_{CO}$. On the
other hand, our method does not require any assumption concerning $X_{CO}$.
Finally, we note that Kalberla et al. (2020) decomposed the H i components
based on their spread ($T_{D}$), while we treated all the H i components
uniformly.

####  4.2.3 Comparison of the maps

Fig. 6 shows the $N_{{H}_{2}}$ map of Kalberla et al. (2020) together with our
map at $2\degr$ resolution. We only include regions with log $N_{{H}_{2}}$
(cm-2) $\geq{20}$, because below that threshold both methods might be
susceptible to noise artefacts.

The $N_{{H}_{2}}$ structures are more extended in the map of Kalberla et al.
(2020) than ours. This difference could be due to degeneracies in the
definition of $T_{\rm D}$, which consists of a non–thermal (turbulent) and
thermal component. The method of Kalberla et al. (2020) employs $T_{\rm D}$ as
a proxy for the gas kinetic temperature, which is correlated with the
molecular abundance. This approximation is accurate because the distribution
of sonic Mach numbers in the diffuse ISM has a well-defined peak 444The sonic
Mach number ($\mathcal{M}_{s}$) is defined as $\mathcal{M}_{s}=\delta
u_{turb}/c_{s}$, where $\delta u_{turb}$, and $c_{s}$ correspond to the
turbulent and thermal sound speed respectively. The Doppler temperature is
$T_{D}=21.86\leavevmode\nobreak\ \delta u^{2}$ (Payne et al., 1980), where
$\delta u$ represents the broadening of the H i emission line, which consists
of a thermal and turbulent component, hence $T_{D}\propto\delta
u^{2}_{turb}+\delta u^{2}_{thermal}$. In the diffuse ISM, the average sonic
Mach number is $\langle\mathcal{M}_{s}\rangle=3.1$ (Heiles & Troland, 2003b).
The above equations suggest that $\langle T_{D}\rangle\propto
3.1\leavevmode\nobreak\ \langle c^{2}_{s}\rangle+\langle\delta
u^{2}_{thermal}\rangle$. Both $c_{s}$, and $u_{thermal}$ vary with the square
root of the gas temperature, and thus $\langle T_{D}\rangle\propto\langle
T\rangle$. This implies that the average $T_{D}$ probes the average gas
temperature in the diffuse ISM. In the absence of a characteristic sonic Mach
number, the turbulent and thermal components of $T_{D}$ become inseparable.
(Heiles & Troland, 2003b). Although $T_{\rm D}$ represents the average gas
temperature in the diffuse ISM, deviations toward individual regions are
inevitably present.

We define the difference in the column densities of our map and the map of
Kalberla et al. (2020) as:

$\Delta\left(\log N_{H_{2}}\right)=\log N_{H_{2}}^{\rm here}-\log
N_{H_{2}}^{\rm Kalberla+2020}.$ (8)

Fig. 7 shows a full-sky map of $\Delta\left(\rm logN_{H_{2}}\right)$,
calculated only in regions where both maps yield $N_{{H}_{2}}$ $\geq 10^{20}$
cm-2. The majority of the pixels have $\Delta\left(\rm logN_{H_{2}}\right)<0$,
while in the Galactic plane we find that $\Delta\left(\rm
logN_{H_{2}}\right)>0$.

We focus on Galactic latitudes $|b|\geq 10\degr$, because typical extinctions
in the Galactic plane exceed the $A_{V}$ range used for the calibration of the
two methods. Fig. 8 shows the distribution of $\Delta\left(\rm
logN_{H_{2}}\right)$ of all pixels at $|b|\geq 10\degr$. The distribution
peaks at -0.5, which implies that the $N_{{H}_{2}}$ estimates in the map of
Kalberla et al. (2020) are, on average, $\sim 3$ times larger than ours. Their
estimates are consistent within the uncertainties of our method (Sect. 2.4).

The method of Kalberla et al. (2020) was calibrated in CO-dark H2 regions,
while in CO-bright regions they estimated $N_{{H}_{2}}$ by assuming that
$X_{CO}$ = $4\times 10^{20}$ cm-2 (K km s-1)-1. Although this value is
consistent with the range of local $X_{CO}$ variations, $1.7-4.1\times
10^{20}$ cm-2 (K km s-1)-1, it is a factor of two larger than the Galactic
average (Bolatto et al., 2013). If Kalberla et al. (2020) had adopted the
global value of $X_{CO}$, then our estimates would differ, on average, by less
than a factor of two. A careful visual comparison of the $\Delta\left(\rm
logN_{H_{2}}\right)$ (Fig. 7) and the $W_{CO}$ map (Fig. 3, bottom panel)
shows that the $N_{{H}_{2}}$ estimates of the two methods agree well,
$\Delta\left(\rm logN_{H_{2}}\right)\approx 0$, in regions where CO is
undetected.

Figure 8: Distribution of the difference of our $N_{{H}_{2}}$ estimates and
those of Kalberla et al. (2020) (Eq. 8) for regions with $|b|\geq 10\degr$.
Our $N_{{H}_{2}}$ estimates are, on average, three times smaller than those
inferred by Kalberla et al. (2020). This difference likely results in large
part from the $X_{CO}$ value employed in the method of Kalberla et al. (2020).
Figure 9: All-sky image of logarithm of fractional H2 abundance,
$f_{{H}_{2}}$. The mean value of log $f_{{H}_{2}}$is equal to -0.6, with
standard deviation equal to 0.23 Figure 10: Distribution of values of H2
fractional abundance. Our mean $f_{{H}_{2}}$ = 0.25 (with 1$\sigma$ range in
linear scale 0.15 – 0.43) is consistent with Liszt & Gerin (2023), and Shull
et al. (2021), who found the average $f_{{H}_{2}}$ to be 0.2 – 0.4 and 0.20,
respectively. Figure 11: Cumulative distribution function of $f_{{H}_{2}}$.
Atomic hydrogen is more abundant than molecular: $50\%$, and $80\%$ of the
LOSs have $f_{{H}_{2}}$ $\lesssim 25\%$, and $f_{{H}_{2}}$ $\lesssim 40\%$
respectively, while only $\sim 3\%$ have $f_{{H}_{2}}$ $\geq 80\%$.

##  5 Characterization of the total gas properties of our Galaxy

Characterizing the total gas properties of our Galaxy is an important but
challenging task. It is hard to constrain $N_{{H}_{2}}$, which is usually
estimated using CO observations by assuming some global (and constant)
$X_{CO}$ value. However, $X_{CO}$ varies by orders of magnitude between the
various regions (Sect. 4.1), and in addition a significant portion of the
molecular gas lies in diffuse clouds that are largely devoid of CO (CO-dark
H2).

Dust traces the total hydrogen column density, including atomic and molecular
in either CO-dark or CO-bright form. For the $N_{{H}_{2}}$ estimate, we
employed the non-linearities in the dust extinction (Sect. 2). Thus our
$N_{{H}_{2}}$ map (Fig. 3) probes the total H2 column (CO-dark and CO-bright),
which enables us to explore the various gas properties irrespective of the gas
phase (atomic or molecular, CO-dark or CO-bright H2). In the following
sections, we use our constructed $N_{{H}_{2}}$ map to explore the relative
abundances between H i and H2 gas, and we constrain the sky distribution of
CO-dark and CO-bright H2.

###  5.1 Relative abundance of atomic and molecular gas

The molecular fractional abundance ($f_{{H}_{2}}$) is defined as:

$f_{{H}_{2}}=\frac{2N_{{H}_{2}}}{N_{{H}\,{\sc I}}+2N_{{H}_{2}}}$ (9)

When $f_{{H}_{2}}$ $\rightarrow 0$, then the total gas column is mostly in
atomic form, while when $f_{{H}_{2}}$ $\rightarrow 1$ gas is $100\%$
molecular.

Some fraction of the dust is expected to be mixed with ionized hydrogen, which
is omitted in Eq. 9. The contribution of $\rm H^{+}$ in dust extinction, if
any, is only expected in low-extinction regions ($A_{V}$ $\lesssim 0.093$ mag,
Liszt, 2014a). However, this is below the extinction range where our method is
applicable, which is log $N_{{H}_{2}}$ $\geq 19$ (Fig. 1) or log $A_{V}$
$\geq-0.5$ (Fig. 2). Thus, we expect that our estimated $f_{{H}_{2}}$ remains
statistically unchanged even if some ionized hydrogen is present.

We used the $N_{{H}_{2}}$ values from the map we constructed (Fig. 3, upper
panel) and the $N_{{H{\sc I}}}$ data from the HI4PI survey to construct a
full-sky $f_{{H}_{2}}$ map (Fig. 9). We observe that close to the Galactic
center ($\ell\lesssim 60\degr$ and $\ell\gtrsim 300\degr$, $b=0\degr$), gas
tends to be fully molecular. On the other hand, $f_{{H}_{2}}$ dramatically
decreases in the outer parts of the Galactic plane
($300\degr\gtrsim\ell\gtrsim 60\degr$, $b=0\degr$), meaning that the relative
abundance of H i increases.

The $f_{{H}_{2}}$ full-sky map also shows several large-scale structures above
the Galactic plane with enhanced molecular abundances ($\log$ $f_{{H}_{2}}$
$\gtrsim 0.5$ ): 1) the Taurus molecular cloud at $\ell,b\sim
170\degr-15\degr$, 2) the Polaris Flare at $\ell,b\sim 125\degr,30\degr$, 3)
the North Polar Spur at $\ell,b\sim 27\degr,10\degr$, and 4) Orion at
$\ell,b\sim 170\degr,-15\degr$. Below we argue that the proximity of these
molecular structures imply that the scale height of $N_{{H}_{2}}$ of our
Galaxy should be less than 300 pc.

Zucker et al. (2021) used the 3D extinction map of Leike et al. (2020) to
derive the distances of Taurus and Orion, which are close to $150$ and $400$
pc respectively. The distance estimates to the Polaris Flare are more
uncertain, ranging from 100 pc up to 400 pc (Heithausen & Thaddeus, 1990;
Zagury et al., 1999; Brunt, 2003; Schlafly et al., 2014; Panopoulou et al.,
2016). The North Polar Spur is a gigantic radio loop and there is a debate
regarding its distance. Some Galactic models suggest that this structure is a
few kpc away from the Sun. However, recently Panopoulou et al. (2021) compared
sub-mm and radio polarization data and found that the maximum distance of the
North Polar Spur is approximately equal to 400 pc from the Sun. Altogether,
all the aforementioned large-scale H2 structures above the Galactic plane are
relatively close to the sun.

To be conservative, we assume that the maximum distance of the aforementioned
molecular clouds is 500 pc. Assuming that $b=45\degr$, which is the maximum
Galactic Latitude where molecular clouds are observed in Fig. 3, we obtain a
vertical distance from the Galactic plane equal to $350$ pc. This distance
barely exceeds the outer edge of the Local Bubble (Pelgrims et al., 2020).
Therefore, most of the H2 gas of our Galaxy lies within a vertical distance of
350 pc of the Galactic plane. This is consistent with previous constraints of
the H2 scale height of our Galaxy (Heyer & Dame, 2015; Marasco et al., 2017).

Fig. 10 shows the distribution of the logarithm of $f_{{H}_{2}}$. The
distribution is nearly symmetric with an average equal to
$\langle$$f_{{H}_{2}}$$\rangle$ $\approx 25\%$, which is consistent with
previous constraints derived with UV spectra of sparsely located point sources
(Shull et al., 2021). From the cumulative distribution function of
$f_{{H}_{2}}$, shown in Fig. 11, we calculate that a large fraction ($80\%$)
of our Galaxy with $N_{{H}_{2}}$ $\geq 10^{20}$ cm-2 has $f_{{H}_{2}}$
$\lesssim 40\%$. This indicates that either the majority of the molecular gas
of our Galaxy lies in diffuse molecular clouds or that atomic clouds are more
abundant than molecular clouds. Below we argue that the latter is more
probable.

We display $f_{{H}_{2}}$ versus $A_{V}$ in Fig. 12. For $-0.5\lesssim$ $\log$
$A_{V}$ $\lesssim-0.2$, $\log$ $f_{{H}_{2}}$ increases from -2.0 ($1\%$ ) to
-0.5 ($30\%$), because the atomic gas transitions to molecular in this range
of extinction, given an ambient radiation field of $\sim 1$ Habing. The
observed non-linear increase of $f_{{H}_{2}}$, is similar to the predicted
behaviour of semi-analytical models of the H i $\rightarrow$ H2 transition
(Draine & Bertoldi, 1996; Krumholz et al., 2008, 2009; McKee & Krumholz, 2010;
Sternberg et al., 2014, 2023; Bialy & Sternberg, 2016). For $-0.2\lesssim$
$\log$ $A_{V}$ $\lesssim 0.7$, we observe that $\log$ $f_{{H}_{2}}$ slightly
decreases from -0.5 ($30\%$) to -0.7 ($20\%$). Although $A_{V}$ increases by a
factor of 5, $f_{{H}_{2}}$ decreases by a factor of 1.5. This suppression of
$f_{{H}_{2}}$ could be due to young stars that are illuminating the
surrounding material, hence responsible for photo-destruction of the molecular
gas (Madden et al., 2020). We consider this scenario unlikely because the
decline of $f_{{H}_{2}}$ starts at $A_{V}$ $\approx 1$ mag which is typical of
translucent (non star-forming) clouds, such as the Polaris Flare.

We conclude simply that there are more atomic than molecular clouds in most of
the LOSs of our Galaxy, even where we observe molecular gas. The concatenation
of atomic clouds increases $N_{{H{\sc I}}}$, and consequently $A_{V}$, but not
$N_{{H}_{2}}$, hence the reduction of $f_{{H}_{2}}$.

This suppression of $f_{{H}_{2}}$ is more prominent close to the Galactic
plane: the majority of LOSs with $\log$ $f_{{H}_{2}}$ $\approx-1$ and $\log$
$A_{V}$ $\approx 0.7$ are observed at $|b|\approx 5\degr$. However, even in
regions at high Galactic latitudes ($|b|>40\degr$), which cover a large
portion of the sky, we observed a similar trend. For $|b|>40\degr$, there are
several LOSs with $\log$ $f_{{H}_{2}}$ $\approx-0.5$ and $\log$ $A_{V}$
$\approx 0$.

This suppression is expected from the following approximate calculation. The
average number of atomic clouds per LOS, as indicated by H i emission line
data, is close to three for $|b|>40\degr$ (Panopoulou & Lenz, 2020). On the
other hand, the number of molecular clouds per LOS is a maximum two because:
1) molecular gas is associated with the H i cloud having an exceptionally
large column, usually observed in H i emission data at low velocities measured
with respect to the local standard of rest (LSR), and 2) molecular gas
associated with H i clouds at intermediate LSR velocities tend to have lower
column densities than clouds with lower velocities and only occasionally show
any sign of H2 (e.g. Lenz et al., 2015; Röhser et al., 2016a, b). Thus, we
expect the ratio of the number of molecular clouds to atomic clouds to be a
maximum of $\approx 0.67$, implying that the high-Galactic-latitude sky has
more atomic than molecular clouds. The majority of LOSs with $\log$
$f_{{H}_{2}}$ $\approx 0$ and $\log$ $A_{V}$ $\gtrsim 0.7$ are toward the
Galacic midplane ($\ell\lesssim 30\degr$ and $\ell\gtrsim 330\degr$,
$b=0\degr$).

Our conclusion about the decrease of $f_{{H}_{2}}$ at high $A_{V}$ is also
evident in Fig. 20 of Planck Collaboration et al. (2011a). These authors
compare $f_{{H}_{2}}$ and $N_{H}$, where the $N_{{H}_{2}}$ measurements used
for the computation of $f_{{H}_{2}}$ have been obtained from UV spectra
(Rachford et al., 2002, 2009; Gillmon et al., 2006; Wakker, 2006). This
suggests that despite the significant beam difference ($16\arcmin$ in our
study, and pencil beams for the UV data), the observed trend is robust. The
concatenation of multiple clouds along the LOS can significantly complicate
observational studies of the H i $\rightarrow$ H2 transition (Browning et al.,
2003).

Figure 12: Molecular fractional abundance as a function of dust extinction.
The nonlinear increase of $f_{{H}_{2}}$ at $-0.5\lesssim$ $\log$ $A_{V}$ (mag)
$\lesssim-0.2$ corresponds to the atomic to molecular gas transition. At
larger $A_{V}$, $f_{{H}_{2}}$ decreases by a factor of $1.5$, most probably
because even along high-extinction sightlines there exist more atomic than
molecular clouds. Figure 13: Upper panel: Full-sky CO-dark $N_{{H}_{2}}$ map,
showing sightlines with $N_{{H}_{2}}$ $\geq 10^{19}$ cm-2 but devoid of CO
($W_{CO}$ $\leq 1$ K km s-1). Bottom panel: CO-bright $N_{{H}_{2}}$ map,
corresponding to $N_{{H}_{2}}$ $\geq 10^{19}$ cm-2 and $W_{CO}$ $>1$ K km s-1.
The CO-dark H2 spans an area approximately equal to $17\%$ of the total sky.
This is a factor of two larger than that of the CO-bright H2, and extends up
to high Galactic Latitude. On the other hand, CO-bright H2 is mostly
concentrated towards the Galactic plane. Moving away from the Galactic plane,
we observe that H2 transitions from CO-bright to CO-dark form at $|b|\approx
5\deg$. Although we still observe a few CO-bright H2 at higher latitude above
the Galactic plane, most of the molecular gas is in CO-dark form there.

###  5.2 The properties of CO-dark and CO-bright H2

There is a significant fraction of molecular gas that is not traced by CO.
This CO-dark H2 is observed in diffuse ISM clouds where the H i $\rightarrow$
H2 transition is ongoing. In regions where the H i $\rightarrow$ H2 transition
is well advanced, carbon monoxide formation can also occur. Between the
extinctions that characterize the predominance of H2 and CO, carbon
transitions from C+ to C i. For typical ISM conditions, the onset of the H2
formation takes place at $A_{V}$ $\lesssim 0.5$ mag, the C+ $\rightarrow$ C i
transition at $1\lesssim$ $A_{V}$ $\lesssim 3$ mag, while CO becomes the
dominant carrier of carbon at $A_{V}$ $\geq\ 3$ mag. During the transition
from C+ to C i, H2, whose electronic transitions due to photon absorption
(Lyman-Werner lines) start at 11.18 eV, has already started forming but carbon
is still in atomic form. Sufficient shielding is required to reduce the energy
of incident photons bellow the CO photo-dissociation threshold (11.09 eV),
hence allowing the buildup of the abundance of CO. Until this energy threshold
is reached, the H2 absorption lines are detectable but without any
corresponding CO emission line. In this case, molecular hydrogen is primarily
mixed with atomic carbon.

The integrated intensity of the [C ii] line (at 158 $\mu$m) is proportional to
the gas density squared (Goldsmith et al., 2018) at the densities
characteristic of diffuse and translucent clouds. Towards several clouds with
untraceable CO, the observed [C ii] intensity, even assuming all carbon is in
the form of C+, is larger than that expected solely from C-H i collisions
(Langer et al., 2010). The enhancement of the [C ii] line indicates the
presence of H2 gas.

The above discussion makes it clear that our ability to trace CO-dark H2
relies on the observational uncertainties in the CO line. Even in diffuse
regions, there should be some CO emission, although very weak. Strictly
speaking CO-dark H2 is located in regions where the atomic carbon is more
abundant than CO, or in terms of observables, when the [C ii] or [C i] lines
are stronger than CO. In this work, we adopt a detectability threshold of the
CO line in order to define CO-dark H2. The detectability limit is determined
by the noise in the CO data of Dame et al. (2001), which is 1 K km s-1.

We adopt the definition that molecular gas is in CO-dark form when the
following conditions are satisfied: 1) $\log$ $N_{{H}_{2}}$ $\geq 19$ cm-2,
and 2) $W_{CO}$ $\leq 1$ K km s-1. On the other hand, we define molecular gas
that is traceable by CO, CO-bright H2, as: 1) $\log$ $N_{{H}_{2}}$ $\geq 19$
cm-2, and 2) $W_{CO}$ $>1$ K km s-1.

Fig. 13 shows our inferred CO-dark (upper panel), and CO-bright (lower panel)
$N_{{H}_{2}}$ maps. We observe that CO-dark H2 extends to high Galactic
latitudes ($|b|\approx 60\degr$), and has generally lower column densities
than CO-bright H2, because CO-dark H2 is characteristic of regions having
lower densities than CO-bright H2. CO-dark H2 increases towards lower Galactic
Latitudes, but in the Galactic plane almost all of the molecular gas is in CO-
bright form due to the enhanced densities, and column densities there. CO-
bright H2 lies close to the Galactic plane and only a few CO-emitting clouds
can be seen at higher latitudes ($|b|\approx 30\degr$).

We estimate that CO-dark H2 covers $\sim 17\%$ of the total sky area, while
CO-bright H2 $\sim 9\%$ of the sky. These values, however, should be
considered with some caution because the $W_{CO}$ map of Dame et al. (2001) is
considered to be complete only for $|b|\leq 32\degr$, as indicated by the
comparison between the CO intensities with far-infrared, and H i
data555https://lambda.gsfc.nasa.gov/product/foreground/fg_wco_info.html. Thus,
the $W_{CO}$ map may miss some CO emission at high Galactic latitudes
($|b|>32\degr$), although we do not expect the emission at such high $b$ to be
significant. Our argument is supported by the CO survey of Dame & Thaddeus
(2022), which is complete for the entire Northern sky, and shows that very few
LOSs have detectable CO emission – assuming the sensitivity of the Dame et al.
(2001) survey – at $|b|>40\degr$. Therefore, our inferred fractional sky
coverage of the CO-dark, and CO-bright H2 components should be considered as
upper and lower limits respectively. Finally, we find that $66\%$ of the
sightlines with $\log$ $N_{{H}_{2}}$ (cm-2) $\geq 19$ are in CO-dark form,
while the remaining $34\%$ of sightlines are in CO-bright form. If in the
definition of CO-dark H2 we increase the threshold to $\log$ $N_{{H}_{2}}$
(cm-2) $\geq 20$, we obtain that $57\%$ of the LOSs are in CO-dark and $43\%$
in CO-bright form.

Figure 14: Total hydrogen column density versus dust reddening. The colored
curves correspond to previously constrained dust-to-gas ratios. The blue, red,
green, and magenta curves correspond to $N_{H}$ / $E(B-V)$ = 5.1 (Kalberla et
al., 2020), 5.8 (Bohlin et al., 1978), 7.0 (Planck Collaboration et al.,
2014), and 8.8 (Lenz et al., 2017) $\times 10^{21}$ cm-2 mag-1 respectively.

##  6 Dust-to-gas ratio

$N_{H}$ / $E(B-V)$ is a key ratio for ISM studies. Existing constraints on the
Galactic average $N_{H}$ / $E(B-V)$ vary significantly (Liszt & Gerin, 2023).
The constraint of Bohlin et al. (1978), $N_{H}$ / $E(B-V)$ $=5.8\times
10^{21}$ cm-2 mag-1, has been the gold standard for several decades. But,
several recent works (Lenz et al., 2017; Planck Collaboration et al., 2014;
Shull et al., 2021; Liszt & Gerin, 2023; Liszt, 2014a; Gillmon et al., 2006;
Rachford et al., 2009) have found larger values, although Kalberla et al.
(2020) quote a smaller value for the Galactic average of $N_{H}$ / $E(B-V)$.
However, when Kalberla et al. (2020) excluded regions close to the Galactic
plane, they found a larger value than Bohlin et al. (1978).

Variations in the $N_{H}$ / $E(B-V)$ ratio can be induced by different
factors: 1) environmental ISM conditions (Paradis et al., 2012), 2) dynamical
effects (e.g., imperfect coupling between dust and gas motions, Squire &
Hopkins, 2018; Hopkins & Squire, 2018; Hopkins et al., 2020, 2022), 3) the
V–band could be less than optimal for measuring dust extinction (Jura, 1980;
Butler & Salim, 2021).

Atomic clouds at high Galactic Latitudes, where $N_{H}$ $\approx$ $N_{{H{\sc
I}}}$, have larger dust-to-gas ratio, $N_{H}$ / $E(B-V)$ $=8-10\times 10^{21}$
cm-2 mag-1 (Liszt, 2014b; Lenz et al., 2017), than the Galactic plane, where
$N_{H}$ / $E(B-V)$ lies in the range $5.8-6.1\times 10^{21}$ cm-2 mag-1
(Bohlin et al., 1978; Rachford et al., 2002, 2009; Shull et al., 2021).
Uncertainties in the employed techniques (UV or radio spectroscopy, far-
infrared dust emission) could also induce variations in the constrained
$N_{H}$ / $E(B-V)$ (Shull & Panopoulou, 2023).

The only direct way to constrain $N_{H}$, hence the dust-to-gas ratio, is by
the use of spectra of background objects (e.g., Bohlin et al., 1978; Sheffer
et al., 2008; Gillmon et al., 2006; Rachford et al., 2009; Shull et al.,
2021). However, obtaining spectra is time consuming and only a limited number
of LOSs can be targeted, usually toward relatively bright stars and quasars.
The sky is sparsely sampled, which implies that the $N_{H}$ / $E(B-V)$
estimates could be biased as spectroscopic surveys usually focus on the
Galactic plane. Indirect methods for constraining the dust-to-gas ratio
require some assumptions and approximations, but they have the advantage of
uniformly sampling the Galaxy (Paradis et al., 2012; Planck Collaboration et
al., 2011b, 2014; Kalberla et al., 2020).

In Fig. 14 we present our $N_{H}$ / $E(B-V)$ constraint and compare it to
previous work. The figure shows the total hydrogen column density, calculated
as $N_{H}$ = $N_{{H{\sc I}}}$ \+ 2$N_{{H}_{2}}$, including data covering the
full celestial sphere; the data employed are the $N_{{H{\sc I}}}$ full-sky map
from HI4PI Collaboration et al. (2016), and our constructed $N_{{H}_{2}}$ map
(Fig. 3). The colored dashed lines correspond to previous constraints.

All the depicted lines are consistent, within variations, with our data.
However, the green and red lines, which correspond to $7.0\times 10^{21}$ cm-2
mag-1 (Planck Collaboration et al., 2014), and $5.8\times 10^{21}$ cm-2 mag-1
(Bohlin et al., 1978) respectively, agree more closely with the mean behaviour
of our data.

We can approximate, to first order, the $N_{H}$ versus $E(B-V)$ relation with
a line. But, a closer inspection to Fig. 14 suggests that some non-linearities
arise. For $E(B-V)$ $\leq 0.1$ mag, the average behaviour of our data agrees
well with the green line describing the results of Planck Collaboration et al.
(2014), while for $E(B-V)$ $>0.1$ mag the red profile of Bohlin et al. (1978)
seems to be a better fit. This indicates that at larger $E(B-V)$ the dust-to-
gas ratio decreases, which could happen due to dust coagulation.

##  7 Discussion

###  7.1 A historical perspective on the extinction residuals

Boulanger et al. (1996) found that dust intensity correlates with the H i
integrated intensity ($W_{{H}\,_{\sc I}}$). The correlation is linear for
$W_{{H}\,_{\sc I}}$ $\leq\ 300$ K km s-1 and non-linear for $W_{{H}\,_{\sc
I}}$ $>300$ K km s-1. Assuming that the H i line is optically thin the H i
column density is given by $N_{{H{\sc I}}}$ = 1.823 $\times 10^{18}$
$W_{{H}\,_{\sc I}}$ (Draine, 2011), yielding a transition at $N_{{H{\sc I}}}$
$\approx 5\times 10^{20}$ cm-2. Boulanger et al. (1996) argued that the non-
linear excess between dust intensity and $W_{{H}\,_{\sc I}}$ is due to the
presence of molecular gas that is only probed by dust intensity which traces
both the atomic and molecular gas phase, while the 21 cm emission line probes
only the atomic phase. Their argument was based on the results of Savage et
al. (1977), who found, using the H2 absorption lines, that the molecular gas
fraction rapidly increases when $N_{{H{\sc I}}}$ $\approx 5\times 10^{20}$
cm-2. An alternative explanation of the non-linear relationship between dust
intensity and $W_{{H}\,_{\sc I}}$ for $W_{{H}\,_{\sc I}}$ $>300$ K km s-1 is
that at some point the emission line becomes optical thick and $W_{{H}\,_{\sc
I}}$ misses a significant faction of the H i due to self-absorption, and
reduction of the observed intensity (Burton et al., 1992).

The same debate is found in more recent papers. Fukui et al. (2014), and Fukui
et al. (2015) explored the correlation between $W_{{H}\,_{\sc I}}$ and the
dust opacity ($\tau_{353}$), but they also added one more free parameter in
the analysis: the dust temperature ($T_{d}$). The dust temperature is
anticorrelated with the total gas column. Warm dust traces the outer regions
of ISM clouds which are primarily atomic, while cold dust traces the inner,
primarily molecular, regions of the clouds. Fukui et al. (2015) showed that
the correlation between $W_{{H}\,_{\sc I}}$ and $\tau_{353}$ is linear when
$T_{d}$ $>20$ K. When dust is cold ($T_{d}$ $<20$ K), $\tau_{353}$ increases
non-linearly with respect to $W_{{H}\,_{\sc I}}$. Fukui et al. (2015) argued
that this nonlinear excess of $\tau_{353}$ with respect to $W_{{H}\,_{\sc I}}$
is consistent with the presence of optically thick H i gas, and although some
H2 gas might be present, it is less abundant than H i. At sufficiently high
column densities, the increase in $\tau_{353}$ would be dominated by the
formation of the molecular gas, which is traceable by the CO emission line.
However, Fukui et al. (2015) suggest that the excess of $\tau_{353}$ at low
$W_{{H}\,_{\sc I}}$ is not due to CO-dark H2, but due to optically thick H i.

Murray et al. (2018) addressed the question of whether the H i optical depth
effects can explain the excess of $\tau_{353}$ when $W_{{H}\,_{\sc I}}$ is
high. They compiled a large catalogue of $\tau_{{H}\,{\sc I}}$ absorption
data, which provide the only direct way to estimate $\tau_{{H}\,{\sc I}}$.
Their sample consisted of archival data towards 151 LOSs that are covered
within the GALFA survey footprint. Following the analysis of Fukui et al.
(2015), Murray et al. (2018) found that as $\tau_{{H}\,{\sc I}}$ increases,
the H i gas that is missed by the 21 cm emission line is less than a factor of
1.3; similar results have been found by Liszt (2019). This means that the
increase in $\tau_{{H}\,{\sc I}}$ cannot adequately explain the nonlinear
relation between $W_{{H}\,_{\sc I}}$ and $\tau_{353}$.

We found that nonlinear deviations in the extinction with respect to
$N_{{H{\sc I}}}$ correlate well with $N_{{H}_{2}}$. Our findings (Fig. 1)
support the claims of Boulanger et al. (1996), and Murray et al. (2018) – that
extinction residuals are induced by molecular hydrogen. However, given the
spread in our obtained $A_{V}$ versus $N_{{H}_{2}}$ relation (Fig. 1), some,
but not dominant, contribution of optically thick H i gas is also plausible.

###  7.2 Variations of $X_{CO}$

The $X_{CO}$ factor is an important observable for estimating the molecular
gas content of molecular clouds and entire galaxies (Lee et al., 2014; Sun et
al., 2018, 2020). $X_{CO}$ is expected to vary as a function of the local ISM
properties, such as metallicity, density, and intensity of the UV background
field (Gong et al., 2017, 2018). Observations of local ISM clouds suggest that
the $X_{CO}$ factor is relatively constant, with variations not exceeding a
factor of two (Bolatto et al., 2013; Liszt & Gerin, 2016).

The constancy of $X_{CO}$ indicates that ISM clouds in our Galaxy have
generally similar properties, such as temperature, and surface densities. This
uniformity has important consequences regarding the nature of ISM clouds. For
example, it explains the relation between non–thermal velocity linewidths and
cloud size – known as Larson’s relation (Larson, 1981) – as the outcome of
clouds with magnetically–supported boundaries (Mouschovias & Psaltis, 1995). A
plethora of theoretical works address the question of why ISM clouds in the
Milky Way are so similar (e.g., Narayanan & Hopkins, 2013; Clark & Glover,
2015).

We have found that CO-emitting regions are most probably observed when
$20.5\lesssim$ $\log$ $N_{{H}_{2}}$ (cm-2) $\lesssim 21$, and $3\lesssim$
$W_{CO}$ $\lesssim 8$ K km s-1 (Fig. 5). Our results, thus suggest, in
agreement with previous studies, that our Galaxy favors clouds with relatively
uniform CO and surface density properties. This is also supported by our
inferred $f_{{H}_{2}}$ distribution (Fig. 10) which is nearly symmetric and
narrow; the 1 $\sigma$ of the distribution is 0.15 - 0.43. However, in
individual molecular clouds (on angular scales of a few arcminutes), $X_{CO}$
can vary by an order of magnitude from the outer to the inner parts, as atomic
gas transitions to fully molecular form.

###  7.3 Molecular hydrogen column densities as a function of the
Galactocentric distance

The Galactic plane contains most of the molecular gas of our Galaxy. Several
authors employed data from the CO survey of Dame et al. (2001) to model the H2
surface as a function of distance from the Galactic center, but omitted any
contribution from CO-dark H2 (Bronfman et al., 1988; Wouterloot et al., 1990;
Nakanishi & Sofue, 2006, 2016; Marasco et al., 2017; Miville-Deschênes et al.,
2017).

A common conclusion from these works is that the H2 surface density decreases
nonlinearly at large Galactocentric distances, although the different studies
differ on the exact cutoff distance. It would be interesting to see how the
effect of the CO-dark H2 could affect the relative abundance between the H i
and H2 surface densities as a function of Galactocentric distance. The
presence of CO-dark H2 could yield more extended H2 surface density profiles,
as shown by Pineda et al. (2013) who estimated that the H2 gas in CO-dark form
corresponds to $\sim 20\%$ of the total molecular abundance at 4 kpc and $\sim
80\%$ at 10 kpc.

Our $f_{{H}_{2}}$ map (Fig. 9) shows that in the Galactic plane gas is mostly
molecular for $\ell\lesssim 60\degr$ and $\ell\gtrsim 300\degr$, while it
becomes fully atomic for $90\degr\lesssim\ell\lesssim 270\degr$. Despite the
presence of CO-dark H2, the overall contribution of molecular gas at large
Galactic latitudes is minimal; the abundance of CO-dark H2 becomes maximum at
$150\degr\lesssim\ell\lesssim 210\degr$, but there $f_{{H}_{2}}$ is only
$10\%$. Thus, the CO-dark H2 abundance seems to increase with Galactocentric
distance, but this only affects the exact distance where H i starts becoming
more abundant than H2.

###  7.4 CO-dark H2 estimates in the Solar neighborhood

Existing constraints on the relative abundance of CO-dark H2 in relation to
total gas imply that a significant amount of molecular gas exists in CO-dark
form (Planck Collaboration et al., 2011b; Paradis et al., 2012; Kalberla et
al., 2020). As stated by Li et al. (2018), the derived abundances of CO-dark
H2 depend on the detectability threshold of the CO line. This was demonstrated
by Donate & Magnani (2017) who performed a high-sensitivity CO survey toward a
region with existing, lower sensitivity CO data. The (new) high-sensitivity
survey showed that the relative abundance of CO-dark with respect to the total
H2 gas was approximately two times smaller than what had been inferred by the
(previous) low-sensitivity CO survey. That said, the definition of CO-dark H2
is vague, making the comparison of the various results in the literature
challenging.

Our definition of CO-dark and CO-bright H2 components differs from several
past works (e.g., Paradis et al., 2012; Kalberla et al., 2020). In these
works, the decomposition is based on the models of photo-dissociation regions
(PDRs, Wolfire et al., 2010; Velusamy et al., 2010). In PDR models, molecular
gas phases are organized into layers with CO-dark corresponding to a diffuse,
and hot layer that surrounds the dense, and cold CO-bright H2 layer.

We define CO-dark H2, molecular hydrogen in LOSs devoid of CO. Using this
definition, CO-dark and CO-bright H2 never co-exist (Fig. 13). However, our
definition allows us to emphasize that a significant fraction of the molecular
sky ($\sim 60\%$) is untraceable in CO maps, when a detectability threshold
equal to 1 K km s-1 is employed (Sect. 5.2). Our estimates should be treated
as upper limits, given the incompleteness of the Dame et al. (2001) CO survey
at high Galactic latitudes. However, we do not expect significant variations
from these estimates because CO emission is minimal in that portion of the
Milky Way.

##  8 Conclusions

Molecular hydrogen is the most abundant molecule in the ISM, but it is hard to
directly observe it and infer its column density. Several indirect methods
have been used to probe $N_{{H}_{2}}$, usually by observing the CO emission
lines and employing a CO – H2 conversion factor ($X_{CO}$, Eq. 1). However, CO
misses the diffuse molecular hydrogen that lies in translucent clouds, where
densities are not sufficient to allow rapid formation of CO and column
densities are insufficient to provide effective shielding against
photdissociation.The result is CO-dark H2. We present a new indirect method
for estimating $N_{{H}_{2}}$ using $N_{{H{\sc I}}}$ and dust extinction data.

We show that the nonlinear increments of dust extinction with respect to
$N_{{H{\sc I}}}$, $A_{V}({H}_{2})$, observed for $N_{{H{\sc I}}}$ $\geq
3\times 10^{20}$ cm-2 (Lenz et al., 2017), are due to the presence of
molecular hydrogen. We measure the implied extinction due to molecular
hydrogen, $A_{V}({H}_{2})$, using publicly available extinction maps (Green et
al., 2019), toward LOSs with direct $N_{{H}_{2}}$ measurements that have been
obtained with spectroscopic data. We show that $A_{V}({H}_{2})$ correlates
with $N_{{H}_{2}}$ (Fig. 1). We obtain a best-fit model for the
$A_{V}({H}_{2})$ – $N_{{H}_{2}}$ relation (Eq. 7). The fitted relation is the
basis of our methodology for estimating $N_{{H}_{2}}$.

We employed the following assumptions, which are motivated by the literature:
1) $N_{{H{\sc I}}}$ / $E(B-V)$ $=8.8\times 10^{21}$ cm-2 mag-1 (Lenz et al.,
2017), 2) $R_{V}=3.1$ (Cardelli et al., 1989; Schlafly et al., 2014). We used
the fitted relation (Eq. 7) to construct a full-sky $N_{{H}_{2}}$ map at
$16\arcmin$ resolution (top panel in Fig. 3); our $N_{{H}_{2}}$ map does not
employ CO observations, and hence traces both the CO-dark (diffuse), and CO-
bright H2. The accuracy of the obtained $N_{{H}_{2}}$ estimates is better than
a factor of two, three, and five with a probability $68\%$, $95\%$, and $98\%$
respectively.

We compared our inferred $N_{{H}_{2}}$ estimates with those of Kalberla et al.
(2020) and found good agreement (Sect. 4.2); there is some discrepancy which
can be explained by the assumptions of the employed methods. We also compared
our $N_{{H}_{2}}$ map with the CO total intensity obtained from the composite
survey of Dame et al. (2001) (Sect. 4.1). We constructed a map of $X_{CO}$,
which is the ratio between $N_{{H}_{2}}$ and $W_{CO}$ (Fig. 4). We estimate
that the Galactic value of $X_{CO}$ is approximately equal to $2\times
10^{20}$ cm-2 (K km s-1)-1, which is consistent with previous constraints
(Bolatto et al., 2013; Lada & Dame, 2020). However, we found that $X_{CO}$ can
vary by orders of magnitude on angular scales of a few arcminutes.

In the diffuse (peripheral) parts of ISM clouds $X_{CO}$ $>10^{21}$ cm-2 (K km
s-1)-1, because a large portion of the molecular hydrogen there is not traced
by CO, hence the large conversion factor. On the other hand, in the inner
(denser) parts of the cloud, where CO is well mixed with molecular hydrogen,
$X_{CO}$ decreases to $\sim 5\times 10^{19}$ cm-2 (K km s-1)-1. Toward LOSs
with $A_{V}$ $\gtrsim 10$ mag – mostly encountered in the Galactic plane – the
CO line saturates and $X_{CO}$ $>10^{21}$ cm-2 (K km s-1)-1 (Fig. 5).

We used our $N_{{H}_{2}}$ map and archival $N_{{H{\sc I}}}$ data to construct
a map (Fig. 9) of the molecular fractional abundance ($f_{{H}_{2}}$), which is
the ratio between $N_{{H}_{2}}$ and the total gas column (Eq. 9). For our
Galaxy, we find an average $f_{{H}_{2}}$ equal to $25\%$ (Fig. 11), which is
consistent with previous estimates (Bellomi et al., 2020; Shull et al., 2021).
In the Galactic plane ($b=0\degr$), for $90\degr\lesssim\ell\lesssim
270\degr$, $f_{{H}_{2}}$ decreases abruptly from unity to 0.1. This agrees
with previous constraints of the molecular gas distribution of our Galaxy
(Miville-Deschênes et al., 2017; Marasco et al., 2017). We find that
$f_{{H}_{2}}$ approaches unity, implying that gas is fully molecular, only
toward $3\%$ of the LOSs with $N_{{H}_{2}}$ $\geq 10^{20}$ cm-2. A large
fraction of the sky ($\sim 66\%$) with $N_{{H}_{2}}$ $\geq 10^{19}$ cm-2 is
undetected in CO maps, when the sensitivity in the CO data is $W_{CO}$ = 1 K
km s-1 (Sect. 5.2).

We explored the correlation between $N_{H}$ and $E(B-V)$ (Sect. 6). Our data
favor a $N_{H}$ / $E(B-V)$ ratio close to $7\times 10^{21}$ cm-2 mag-1 for
$E(B-V)$ $\lesssim 0.1$ mag, consistent with the constraint of Planck
Collaboration et al. (2014). For $E(B-V)$ $\gtrsim 0.1$ mag, our data fit
better with the constraint of Bohlin et al. (1978), $N_{H}$ / $E(B-V)$
$\approx 5.8\times 10^{21}$ cm-2 mag-1. All the constructed maps will be made
publicly available.

###### Acknowledgements.

We are grateful to the reviewer, J. M. Shull, for a very constructive report
which significantly improved the manuscript. We would like to thank M. Heyer,
P. Kalberla, K. Tassis, S. Stahler, and G. V. Panopoulou for useful comments
on the draft and M. A Miville-Deschenês, C. Murray, S. Bialy, A. Bracco, and
B. Hensley for fruitful discussions. We also thank P. Kallemi for the graphic
design in part of the figures. This work was supported by NSF grant
AST-2109127. This work was performed in part at the Jet Propulsion Laboratory,
California Institute of Technology, under contract with the National
Aeronautics and Space Administration (80NM0018D0004). The authors acknowledge
Interstellar Institute’s program “II6” and the Paris-Saclay University’s
Institut Pascal for hosting discussions that nourished the development of the
ideas behind this work.

## References

  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  * Armus et al. (2023) Armus, L., Lai, T., U, V., et al. 2023, ApJ, 942, L37
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Barriault et al. (2010) Barriault, L., Joncas, G., Lockman, F. J., & Martin, P. G. 2010, MNRAS, 407, 2645
  * Bellomi et al. (2020) Bellomi, E., Godard, B., Hennebelle, P., et al. 2020, A&A, 643, A36
  * Bialy & Sternberg (2016) Bialy, S. & Sternberg, A. 2016, ApJ, 822, 83
  * Boccaletti et al. (2015) Boccaletti, A., Lagage, P. O., Baudoz, P., et al. 2015, PASP, 127, 633
  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, Annual Review of Astronomy and Astrophysics, 51, 207
  * Borchert et al. (2022) Borchert, E. M. A., Walch, S., Seifried, D., et al. 2022, MNRAS, 510, 753
  * Bouchet et al. (2015) Bouchet, P., García-Marín, M., Lagage, P. O., et al. 2015, PASP, 127, 612
  * Boulanger et al. (1996) Boulanger, F., Abergel, A., Bernard, J. P., et al. 1996, A&A, 312, 256
  * Bronfman et al. (1988) Bronfman, L., Cohen, R. S., Alvarez, H., May, J., & Thaddeus, P. 1988, ApJ, 324, 248
  * Browning et al. (2003) Browning, M. K., Tumlinson, J., & Shull, J. M. 2003, ApJ, 582, 810
  * Brunt (2003) Brunt, C. M. 2003, ApJ, 584, 293
  * Burton et al. (1992) Burton, W. B., Elmegreen, B. G., Genzel, R., Pfenniger, D., & Bartholdi, P. 1992, in Astronomical Society of the Pacific Conference Series, Vol. 31, Astronomical Society of the Pacific Conference Series, 327
  * Butler & Salim (2021) Butler, R. E. & Salim, S. 2021, ApJ, 911, 40
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Cartledge et al. (2004) Cartledge, S. I. B., Lauroesch, J. T., Meyer, D. M., & Sofia, U. J. 2004, ApJ, 613, 1037
  * Clark & Glover (2015) Clark, P. C. & Glover, S. C. O. 2015, MNRAS, 452, 2057
  * Cox et al. (2016) Cox, N. L. J., Arzoumanian, D., André, P., et al. 2016, A&A, 590, A110
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * Dame & Thaddeus (2022) Dame, T. M. & Thaddeus, P. 2022, ApJS, 262, 5
  * Donate & Magnani (2017) Donate, E. & Magnani, L. 2017, MNRAS, 472, 3169
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Draine & Bertoldi (1996) Draine, B. T. & Bertoldi, F. 1996, ApJ, 468, 269
  * Draine et al. (2007) Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fukui et al. (2014) Fukui, Y., Okamoto, R., Kaji, R., et al. 2014, Astrophysical Journal, 796
  * Fukui et al. (2015) Fukui, Y., Torii, K., Onishi, T., et al. 2015, ApJ, 798, 6
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gillmon et al. (2006) Gillmon, K., Shull, J. M., Tumlinson, J., & Danforth, C. 2006, ApJ, 636, 891
  * Glasse et al. (2015) Glasse, A., Rieke, G. H., Bauwens, E., et al. 2015, PASP, 127, 686
  * Glover & Clark (2016) Glover, S. C. O. & Clark, P. C. 2016, Monthly Notices of the Royal Astronomical Society, 456, 3596
  * Glover & Smith (2016) Glover, S. C. O. & Smith, R. J. 2016, MNRAS, 462, 3011
  * Goldsmith (2013) Goldsmith, P. F. 2013, ApJ, 774, 134
  * Goldsmith et al. (2018) Goldsmith, P. F., Pineda, J. L., Neufeld, D. A., et al. 2018, ApJ, 856, 96
  * Gong et al. (2018) Gong, M., Ostriker, E. C., & Kim, C.-G. 2018, ApJ, 858, 16
  * Gong et al. (2017) Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38
  * Gordon et al. (2015) Gordon, K. D., Chen, C. H., Anderson, R. E., et al. 2015, PASP, 127, 696
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Grenier et al. (2005) Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005, Science, 307, 1292
  * Gudennavar et al. (2012) Gudennavar, S. B., Bubbly, S. G., Preethi, K., & Murthy, J. 2012, ApJS, 199, 8
  * Heiles & Troland (2003a) Heiles, C. & Troland, T. H. 2003a, ApJS, 145, 329
  * Heiles & Troland (2003b) Heiles, C. & Troland, T. H. 2003b, ApJ, 586, 1067
  * Heithausen & Thaddeus (1990) Heithausen, A. & Thaddeus, P. 1990, ApJ, 353, L49
  * Heyer & Dame (2015) Heyer, M. & Dame, T. M. 2015, ARA&A, 53, 583
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Hopkins et al. (2022) Hopkins, P. F., Rosen, A. L., Squire, J., et al. 2022, MNRAS, 517, 1491
  * Hopkins & Squire (2018) Hopkins, P. F. & Squire, J. 2018, MNRAS, 480, 2813
  * Hopkins et al. (2020) Hopkins, P. F., Squire, J., & Seligman, D. 2020, MNRAS, 496, 2123
  * Jaschek & Egret (1982) Jaschek, M. & Egret, D. 1982, in Be Stars, ed. M. Jaschek & H. G. Groth, Vol. 98, 261
  * Jura (1980) Jura, M. 1980, ApJ, 235, 63
  * Juvela et al. (2018) Juvela, M., He, J., Pattle, K., et al. 2018, A&A, 612, A71
  * Kalberla & Haud (2015) Kalberla, P. M. W. & Haud, U. 2015, A&A, 578, A78
  * Kalberla & Haud (2018) Kalberla, P. M. W. & Haud, U. 2018, A&A, 619, A58
  * Kalberla et al. (2020) Kalberla, P. M. W., Kerp, J., & Haud, U. 2020, A&A, 639, A26
  * Kendrew et al. (2015) Kendrew, S., Scheithauer, S., Bouchet, P., et al. 2015, PASP, 127, 623
  * Kirk et al. (2010) Kirk, J. M., Polehampton, E., Anderson, L. D., et al. 2010, A&A, 518, L82
  * Könyves et al. (2015) Könyves, V., André, P., Men’shchikov, A., et al. 2015, A&A, 584, A91
  * Könyves et al. (2010) Könyves, V., André, P., Men’shchikov, A., et al. 2010, A&A, 518, L106
  * Krumholz et al. (2008) Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2008, ApJ, 689, 865
  * Krumholz et al. (2009) Krumholz, M. R., McKee, C. F., & Tumlinson, J. 2009, ApJ, 693, 216
  * Lada & Dame (2020) Lada, C. J. & Dame, T. M. 2020, ApJ, 898, 3
  * Langer et al. (2015) Langer, W. D., Goldsmith, P. F., Pineda, J. L., et al. 2015, A&A, 576, A1
  * Langer et al. (2010) Langer, W. D., Velusamy, T., Pineda, J. L., et al. 2010, A&A, 521, L17
  * Langer et al. (2014) Langer, W. D., Velusamy, T., Pineda, J. L., Willacy, K., & Goldsmith, P. F. 2014, A&A, 561, A122
  * Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
  * Lebouteiller et al. (2019) Lebouteiller, V., Cormier, D., Madden, S. C., et al. 2019, A&A, 632, A106
  * Lee et al. (2014) Lee, M.-Y., Stanimirović, S., Wolfire, M. G., et al. 2014, ApJ, 784, 80
  * Leike et al. (2020) Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138
  * Lenz et al. (2017) Lenz, D., Hensley, B. S., & Doré, O. 2017, ApJ, 846, 38
  * Lenz et al. (2015) Lenz, D., Kerp, J., Flöer, L., et al. 2015, A&A, 573, A83
  * Li et al. (2018) Li, D., Tang, N., Nguyen, H., et al. 2018, ApJS, 235, 1
  * Li et al. (2015) Li, D., Xu, D., Heiles, C., Pan, Z., & Tang, N. 2015, Publication of Korean Astronomical Society, 30, 75
  * Lindner et al. (2015) Lindner, R. R., Vera-Ciro, C., Murray, C. E., et al. 2015, The Astronomical Journal, 149, 138
  * Liszt (2014a) Liszt, H. 2014a, ApJ, 783, 17
  * Liszt (2014b) Liszt, H. 2014b, ApJ, 780, 10
  * Liszt (2019) Liszt, H. 2019, ApJ, 881, 29
  * Liszt & Gerin (2023) Liszt, H. & Gerin, M. 2023, ApJ, 943, 172
  * Liszt & Gerin (2016) Liszt, H. S. & Gerin, M. 2016, A&A, 585, A80
  * Lombardi et al. (2006) Lombardi, M., Alves, J., & Lada, C. J. 2006, A&A, 454, 781
  * Madden et al. (2020) Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141
  * Mangum & Shirley (2015) Mangum, J. G. & Shirley, Y. L. 2015, PASP, 127, 266
  * Marasco et al. (2017) Marasco, A., Fraternali, F., van der Hulst, J. M., & Oosterloo, T. 2017, A&A, 607, A106
  * Marchal et al. (2019) Marchal, A., Miville-Deschênes, M.-A., Orieux, F., et al. 2019, A&A, 626, A101
  * Martin et al. (2012) Martin, P. G., Roy, A., Bontemps, S., et al. 2012, ApJ, 751, 28
  * McKee & Krumholz (2010) McKee, C. F. & Krumholz, M. R. 2010, ApJ, 709, 308
  * Miville-Deschênes et al. (2002) Miville-Deschênes, M. A., Boulanger, F., Joncas, G., & Falgarone, E. 2002, A&A, 381, 209
  * Miville-Deschênes et al. (2010) Miville-Deschênes, M. A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104
  * Miville-Deschênes et al. (2017) Miville-Deschênes, M.-A., Murray, N., & Lee, E. J. 2017, ApJ, 834, 57
  * Mouschovias & Psaltis (1995) Mouschovias, T. C. & Psaltis, D. 1995, ApJ, 444, L105
  * Murray et al. (2014) Murray, C. E., Lindner, R. R., Stanimirović, S., et al. 2014, ApJ, 781, L41
  * Murray et al. (2018) Murray, C. E., Peek, J. E. G., Lee, M.-Y., & Stanimirović, S. 2018, ApJ, 862, 131
  * Nakanishi & Sofue (2006) Nakanishi, H. & Sofue, Y. 2006, PASJ, 58, 847
  * Nakanishi & Sofue (2016) Nakanishi, H. & Sofue, Y. 2016, PASJ, 68, 5
  * Narayanan & Hopkins (2013) Narayanan, D. & Hopkins, P. F. 2013, MNRAS, 433, 1223
  * Nguyen et al. (2018) Nguyen, H., Dawson, J. R., Miville-Deschênes, M. A., et al. 2018, ApJ, 862, 49
  * Palmeirim et al. (2013) Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38
  * Panopoulou et al. (2021) Panopoulou, G. V., Dickinson, C., Readhead, A. C. S., Pearson, T. J., & Peel, M. W. 2021, ApJ, 922, 210
  * Panopoulou & Lenz (2020) Panopoulou, G. V. & Lenz, D. 2020, ApJ, 902, 120
  * Panopoulou et al. (2016) Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517
  * Paradis et al. (2012) Paradis, D., Dobashi, K., Shimoikura, T., et al. 2012, A&A, 543, A103
  * Paradis et al. (2023) Paradis, D., Mény, C., Demyk, K., Noriega-Crespo, A., & Ristorcelli, I. 2023, A&A, 674, A141
  * Paradis et al. (2010) Paradis, D., Veneziani, M., Noriega-Crespo, A., et al. 2010, A&A, 520, L8
  * Payne et al. (1980) Payne, H. E., Salpeter, E. E., & Terzian, Y. 1980, ApJ, 240, 499
  * Peek & Schiminovich (2013) Peek, J. E. G. & Schiminovich, D. 2013, ApJ, 771, 68
  * Pelgrims et al. (2020) Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17
  * Pineda et al. (2008) Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481
  * Pineda et al. (2010) Pineda, J. L., Goldsmith, P. F., Chapman, N., et al. 2010, ApJ, 721, 686
  * Pineda et al. (2017) Pineda, J. L., Langer, W. D., Goldsmith, P. F., et al. 2017, ApJ, 839, 107
  * Pineda et al. (2013) Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103
  * Planck Collaboration et al. (2014) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11
  * Planck Collaboration et al. (2011a) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011a, A&A, 536, A24
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A138
  * Planck Collaboration et al. (2011b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011b, A&A, 536, A19
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A12
  * Rachford et al. (2009) Rachford, B. L., Snow, T. P., Destree, J. D., et al. 2009, ApJS, 180, 125
  * Rachford et al. (2002) Rachford, B. L., Snow, T. P., Tumlinson, J., et al. 2002, ApJ, 577, 221
  * Rieke et al. (2015) Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584
  * Riener et al. (2019) Riener, M., Kainulainen, J., Henshaw, J. D., et al. 2019, A&A, 628, A78
  * Ripple et al. (2013) Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296
  * Röhser et al. (2016a) Röhser, T., Kerp, J., Ben Bekhti, N., & Winkel, B. 2016a, A&A, 592, A142
  * Röhser et al. (2016b) Röhser, T., Kerp, J., Lenz, D., & Winkel, B. 2016b, A&A, 596, A94
  * Savage et al. (1977) Savage, B. D., Bohlin, R. C., Drake, J. F., & Budich, W. 1977, ApJ, 216, 291
  * Savage & Sembach (1996) Savage, B. D. & Sembach, K. R. 1996, ARA&A, 34, 279
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schlafly et al. (2010) Schlafly, E. F., Finkbeiner, D. P., Schlegel, D. J., et al. 2010, ApJ, 725, 1175
  * Schlafly et al. (2014) Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Schnee et al. (2006) Schnee, S., Bethell, T., & Goodman, A. 2006, ApJ, 640, L47
  * Seifried et al. (2020) Seifried, D., Haid, S., Walch, S., Borchert, E. M. A., & Bisbas, T. G. 2020, MNRAS, 492, 1465
  * Sheffer et al. (2008) Sheffer, Y., Rogers, M., Federman, S. R., et al. 2008, ApJ, 687, 1075
  * Shull et al. (2021) Shull, J. M., Danforth, C. W., & Anderson, K. L. 2021, ApJ, 911, 55
  * Shull & Panopoulou (2023) Shull, J. M. & Panopoulou, G. V. 2023, arXiv e-prints, arXiv:2310.12205
  * Shull et al. (2000) Shull, J. M., Tumlinson, J., Jenkins, E. B., et al. 2000, ApJ, 538, L73
  * Skalidis et al. (2022) Skalidis, R., Tassis, K., Panopoulou, G. V., et al. 2022, A&A, 665, A77
  * Sofia et al. (2005) Sofia, U. J., Wolff, M. J., Rachford, B., et al. 2005, ApJ, 625, 167
  * Sota et al. (2014) Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10
  * Squire & Hopkins (2018) Squire, J. & Hopkins, P. F. 2018, ApJ, 856, L15
  * Sternberg et al. (2023) Sternberg, A., Bialy, S., & Gurman, A. 2023, arXiv e-prints, arXiv:2308.13889
  * Sternberg et al. (2014) Sternberg, A., Petit, F. L., Roueff, E., & Bourlot, J. L. 2014, The Astrophysical Journal, 790, 10
  * Sun et al. (2020) Sun, J., Leroy, A. K., Schinnerer, E., et al. 2020, ApJ, 901, L8
  * Sun et al. (2018) Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172
  * Velusamy et al. (2010) Velusamy, T., Langer, W. D., Pineda, J. L., et al. 2010, A&A, 521, L18
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323
  * Wakker (2006) Wakker, B. P. 2006, ApJS, 163, 282
  * Wells et al. (2015) Wells, M., Pel, J. W., Glasse, A., et al. 2015, PASP, 127, 646
  * Wolfire et al. (2010) Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  * Wouterloot et al. (1990) Wouterloot, J. G. A., Brand, J., Burton, W. B., & Kwee, K. K. 1990, A&A, 230, 21
  * Wright et al. (2023) Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003
  * Zagury et al. (1999) Zagury, F., Boulanger, F., & Banchet, V. 1999, A&A, 352, 645
  * Zhang et al. (2023) Zhang, R., Yuan, H., & Chen, B. 2023, ApJS, 269, 6
  * Zucker et al. (2021) Zucker, C., Goodman, A., Alves, J., et al. 2021, ApJ, 919, 35

##  Appendix A Full-sky $N_{{H}_{2}}$ properties of our Galaxy using the
Planck extinction map

Figure 15: Full-sky $N_{{H}_{2}}$ map of our Galaxy constructed with the
extinction map of Planck Collaboration et al. (2020). Resolution and grid
properties are the same as in Fig. 3.

We constructed a full-sky $N_{{H}_{2}}$ map of our Galaxy using the extinction
data of Planck Collaboration et al. (2020) (Fig. 15). We compared against the
$N_{{H}_{2}}$ map that we constructed with the extinction data of Schlegel et
al. (1998) (Fig. 3).

Fig. 15 shows more extended low-$N_{{H}_{2}}$ regions at high Galactic
latitudes than Fig. 3. This happens because $A_{V}$ in Planck Collaboration et
al. (2020) tends to be, overall, larger than in Schlegel et al. (1998). As we
show below, the relatively enhanced $A_{V}$ values of Planck Collaboration et
al. (2020) may affect some of the estimated molecular gas properties of our
Galaxy.

We use the $N_{{H}_{2}}$ data shown in Fig. 15 to calculate the Galactic
average $X_{CO}$, and $f_{{H}_{2}}$. We obtain that
$\langle$$X_{CO}$$\rangle\approx 4\times 10^{20}$ cm-2 (K km s-1)-1, and
$\langle$$f_{{H}_{2}}$$\rangle\approx 30\%$, which are $\sim 2$ larger than
what we inferred in Sects. 4.1, and 5, but both consistent with the range of
values in the literature. In addition, the estimated sky fraction of CO-dark
with respect to the total H2 becomes $\sim 79\%$, while that of CO-bright
$\sim 21\%$ for LOSs with $N_{{H}_{2}}$ $\geq 10^{20}$ cm-2. In this case the
relative abundance of CO-dark H2 is almost two times larger than what we get
from the Schlegel et al. (1998) map (Sect. 5.2). We conclude that differences
in the extinction maps can induce at most factor of two deviations in the
estimated molecular hydrogen properties.