# Dust, CO and [C i]: Cross-calibration of molecular gas mass tracers in
metal-rich galaxies across cosmic time

L. Dunne,${}^{\\!1}$ S. J. Maddox,${}^{\\!1}$ P. P.
Papadopoulos,${}^{\\!3,2,1}$ R. J. Ivison4 and H. L. Gomez1
1School of Physics & Astronomy, Cardiff University, Queens Buildings, The
Parade, Cardiff CF24 3AA
2Dept of Physics, Section of Astrophysics, Astronomy and Mechanics, Aristotle
University of Thessaloniki, GR-54124, Greece
3Research Center for Astronomy, Academy of Athens, Soranou Efesiou 4, GR-11527
Athens, Greece
4European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748
Garching, Germany<EMAIL_ADDRESS>
( Submitted to MNRAS Main Journal, 2022 Xxxxxx; Manuscript ID: MN-22-XXXX-MJ )

###### Abstract

We present a self-consistent cross-calibration of the three main molecular gas
mass tracers in galaxies, namely the 12CO(1–0), [C
i](${}^{3}P_{1}$–${}^{3}P_{0}$) lines, and the submm dust continuum emission,
using a sample of 407 galaxies, ranging from local disks to submillimetre-
selected galaxies (SMGs) up to $z\approx 6$. A Bayesian statistical method is
used to produce galaxy-scale universal calibrations of these molecular gas
indicators, that hold over 3–4 orders of magnitude in infrared luminosity,
$L_{\rm IR}$. Regarding the dust continuum, we use a mass-weighted dust
temperature, $T_{\rm mw}$, determined using new empirical relations between
temperature and luminosity. We find the average L/$M_{\rm mol}$ gas mass
conversion factors (including He) to be $\mbox{$\alpha_{850}$}=6.9\times
10^{12}\,\rm W\,Hz^{-1}\,M_{\odot}^{-1}$, $\alpha_{\rm CO}$=4.0 M⊙ (k km s-1
pc2)-1 and $\alpha_{\rm C\,I}$=17.0 M⊙ (k km s-1 pc2)-1, based on the
assumption that the mean dust properties of the sample ($\kappa_{\rm H}$= gas-
to-dust ratio/dust emissivity) will be similar to those of local metal rich
galaxies and the Milky Way. The tracer with the least intrinsic scatter is [C
i](1–0), while CO(1–0) has the highest. The conversion factors show a weak but
significant correlation with $L_{\rm IR}$ which is not apparent when $T_{\rm
mw}$ is held constant. Assuming dust properties typical of metal-rich
galaxies, we infer a neutral carbon abundance $\mbox{$X_{\rm
C\,I}$}=\rm[C^{0}/\mbox{H${}_{2}$}]=1.6\times 10^{-5}$, similar to that in the
Milky Way. We find no evidence for bi-modality of $\alpha_{\rm CO}$ between
main-sequence (MS) galaxies and those with extreme star-formation intensity,
i.e. ultraluminous infrared galaxies (ULIRGs) and SMGs. The means of the three
conversion factors are found to be similar between MS galaxies and
ULIRGs/SMGs, to within 10–20 per cent. The overarching conclusion of our work
is that, for metal-rich galaxies, near-universal average values for
$\alpha_{\rm CO}$, $X_{\rm C\,I}$ and $\kappa_{\rm H}$ are adequate for global
molecular gas estimates within the expected uncertainties. The 1$\sigma$
scatter in our optimised values for $\alpha_{\rm CO}$, $X_{\rm C\,I}$ and
$\kappa_{\rm H}$ are 0.14, 0.11 and 0.15 dex respectively.

###### keywords:

ISM: dust, extinction; Galaxies: high redshift; Submillimetre: galaxies, ISM;
Radio lines: galaxies, ISM

††pubyear: 2022††pagerange: Dust, CO and [C i]: Cross-calibration of molecular
gas mass tracers in metal-rich galaxies across cosmic time–I

##  1 Introduction

The cosmic star-formation rate (SFR) density has declined by more than an
order of magnitude during the past $\approx 8$ Gyr of cosmic history (Lilly et
al., 1996; Madau et al., 1996; Madau & Dickinson, 2014). The driver of star
formation is the molecular gas supply in galaxies, and indeed the SFR–stellar
mass (SFR–$M_{\star}$) relationship known as the galaxy main sequence (MS) is
purely a by-product of the relationship between SFR and molecular gas (e.g.
Baker et al., 2022), for unperturbed galaxies with significant gas reserves. A
major observational goal is to produce a combined census of the molecular gas
– the ‘potential for future star formation’ – and the stellar content – the
‘record of past star formation’ – over this period (e.g. Keres et al., 2003;
Dunne et al., 2003; Dunne et al., 2011; Zwaan et al., 2004; Zafar et al.,
2013; Walter et al., 2014; Decarli et al., 2016; Saintonge et al., 2017;
Driver et al., 2018; Rhee et al., 2018; Decarli et al., 2019; Riechers et al.,
2019).

The molecular gas fraction of a galaxy is a crucial component in models of
galaxy formation (e.g. Obreschkow et al., 2009; Popping et al., 2014; Lagos et
al., 2015; Chen et al., 2018) and thus measurements of $\rm H_{2}$ and stellar
mass over large representative galaxy samples are key requirements for
understanding how galaxies have transformed from clouds of gas residing in
dark matter haloes into the regular agglomerations of stars we see in the
local Universe. While it is clear that CO(1–0)-luminous gas is the phase
linked with star formation (e.g. Wong & Blitz, 2002), observations of
molecules with higher critical densities (e.g. HCN) revealed that it is the
dense H2 gas phase ($n>10^{4}$ cm-3) that correlates most tightly and linearly
with tracers of star-formation (Gao & Solomon, 2004).

Atomic hydrogen (H i), on the other hand, constitutes a longer-term gas
reservoir for star formation, where under certain conditions of pressure, far-
UV radiation field, density and metallicity, a phase transition H i
$\rightarrow$ H2 takes place, catalysed by dust grains (e.g. Elmegreen, 1993;
Papadopoulos et al., 2002; Blitz & Rosolowsky, 2006): a picture supported by
numerous observations (e.g. Honma et al., 1995; Leroy et al., 2008; Bigiel et
al., 2008; Schruba et al., 2011). This transition occurs in the inner H i
distribution of galaxies, in the cold neutral medium (CNM: $n\sim 50$–100
cm-3, $T_{\rm kin}\sim 100$–200 k), meanwhile pure H i gas often extends many
optical radii beyond the luminous stellar disk (e.g. Péroux & Howk, 2020),
where it can be found concomitant with cold dust (e.g. Thomas et al., 2002).

Unlike H i and its hyperfine line emission at 21 cm, the H2 molecule in its
S(0):$J=2-0$ transition at 28$\mu$m (the least excitation-demanding $\rm
H_{2}$ line) is essentially invisible at temperatures typical of giant
molecular clouds (10–20 k). This is because its $\Delta E/k_{\rm
B}$$\sim$$510\,{\sc k}$, limits its excitation and detetection only to shocked
regions of molecular clouds, where gas temperatures can rise past $\sim 1000$
k, for small ($\sim 1$–2 per cent) gas mass fractions. Even then, to observe
this H2 line at 28 $\mu$m requires space-borne telescopes.

For these reasons the rotational transitions of CO (the next most abundant
molecule with $\rm[CO/H_{2}]\sim 10^{-4}$) are commonly used to trace $\rm
H_{2}$ gas, with the lowest transition (12CO $J=1$–0) being the most
established tracer. Its $E_{10}/k_{\rm B}\sim 5.5$ k ensures a well-populated
upper level even in the coldest gas, while its low critical density, $n_{\rm
cr}\sim 400$ cm-3, ensures its excitation even at low densities111Because the
CO(1–0) line is typically optically thick, with $\tau_{10}\sim 5$–10 (Bryant &
Scoville, 1996; Papadopoulos et al., 2012a, e.g.: their Eqn. 11), the
effective critical density is lower still: $n_{\rm
cr}(\beta_{10})=\beta_{10}n_{\rm crit}\sim 40$–80 cm-3, where
$\beta_{10}=(1-e^{-\tau_{10}})/\tau_{10}$ is the line escape probability..

The CO(1–0) line has significant optical depths in the typically macro-
turbulent $\rm H_{2}$ gas, though these arise locally within the velocity-
coherent gas cells allowing the CO emission to trace gas mass throughout
molecular clouds (e.g. Dickman et al., 1986). The conversion factor,
$\alpha_{\rm CO}$, in the relation $\mbox{$M_{\rm H2}$}=\mbox{$\alpha_{\rm
CO}$}\mbox{$L^{\prime}_{\rm CO}$}$ cannot be determined using standard
optically thin line formation physics due to the high line optical depths.
This created the need for a $\alpha_{\rm CO}$ calibration as soon as the
ubiquity of CO line emission in H2 clouds was established. Observational and
theoretical investigation of $\alpha_{\rm CO}$ suggests it is sensitive to
metallicity, molecular gas surface density and kinematic state in galaxies
(e.g. Pelupessy & Papadopoulos, 2009; Narayanan et al., 2011; Papadopoulos et
al., 2012a; Bolatto et al., 2013).

Three distinct problems are now recognised regarding the use of CO as a global
tracer of H2 mass in galaxies:

  1. 1.

The $\alpha_{\rm CO}$ factor is sensitive – in a highly non-linear fashion –
to the ISM metallicity and ambient far-UV radiation fields (e.g. Israel, 1997;
Pak et al., 1998; Bolatto et al., 2013).

  2. 2.

Non-self-gravitating molecular clouds – and/or very different average ISM
states in terms of average temperature and gas density range from those found
in spiral galaxies where $\alpha_{\rm CO}$ was first calibrated – can yield
systematically different $\alpha_{\rm CO}$ factors. For example, $\alpha_{\rm
CO}$$\sim 1/5-1/4\times$ Galactic was initially reported for a sample of four
ULIRGs by Downes & Solomon (1998).

  3. 3.

Elevated cosmic ray (CR) energy densities can destroy CO below a certain gas
density threshold, leaving behind more C-rich gas. This density threshold
depends on the CR energy density in a highly non-linear fashion, as explored
by Bisbas et al. (2015), who found that regions of CO suppression may occur
even in moderately enhanced CR conditions if the gas density is low, while the
very high CR energy densities expected in ultraluminous infrared galaxies
(ULIRGs) may be partly compensated by higher gas densities in such starbursts.
Modelling [C i/CO] ratios as a function of CR, turbulence, gas density and
metallicity is an active area of theoretical research (e.g. Bisbas et al.,
2015; Bisbas et al., 2017; Bisbas et al., 2021; Glover & Clark, 2016; Clark et
al., 2019a; Papadopoulos et al., 2018; Gong et al., 2020).

In the distant Universe, additional problems arise. High-redshift galaxies are
often observed solely in high-$J$ CO lines ($J=3$–2 and higher), due to the
observational challenge of observing the two low-$J$ CO lines222Prior to the
commissioning of its bands 1 and 2, low-$J$ lines from high-redshift galaxies
are inaccessible to the Atacama Large Millimetre Array (ALMA). The Jansky Very
Large Array (JVLA), the Australia Telescope Compact Array (ATCA) and the
Greenbank Telescope (GBT), have in some cases been able to access the faint
low-$J$ ($J_{\rm u}\leq 2$) CO lines, but it requires huge amounts of
observing time in the best available weather.. Using the high-$J$ lines means
that global CO($J+1,J$)/(1–0) ratios must be assumed before an $\alpha_{\rm
CO}$ factor can be used; given the wide range of CO spectral-line energy
distributions (SLEDs) found for LIRGs for $J=3$–2 and higher (Papadopoulos et
al., 2012b; Greve et al., 2014; Kamenetzky et al., 2016), these assumptions
come with large uncertainties. Finally, at the highest redshifts ($\ga$4),
low-$J$ CO lines (and dust emission) can be severely suppressed for cold gas
(and dust) reservoirs due to their low contrast against the ambient, rest-
frame cosmic microwave background (da Cunha et al., 2013; Zhang et al., 2016).

In principle, radiative transfer models of well-sampled CO (and 13CO) SLEDs
can yield $\alpha_{\rm CO}$ values appropriate for a particular galaxy (or
even galaxy class) (e.g. Papadopoulos et al., 2012b; Papadopoulos et al.,
2014; Harrington et al., 2021). Nevertheless, the size of the CO line datasets
per galaxy required to do this make it impractical (in terms of telescope
time) to obtain $M$(H2) for large galaxy samples. Amassing a large sample
typically means only one or two lines can be gathered per galaxy, and thus a
calibration of $\alpha_{\rm CO}$ and its uncertainties remains very valuable.
The only practical way to achieve this is to cross-calibrate against the other
galaxy-scale $\rm H_{2}$ mass tracers.

Large-area far-infrared (FIR) and submillimetre (submm) surveys (e.g. Armus et
al., 2009; Eales et al., 2010; Vieira et al., 2010; Kennicutt et al., 2011;
Oliver et al., 2012; Hodge et al., 2013) ushered in a new era in which submm
continuum emission from dust has been used widely as an alternative tracer of
$M_{\rm H2}$, although it has been clear that submm-derived dust masses
($\propto\mbox{$L_{850}$}$) and CO-derived molecular gas masses
($\propto\mbox{$L^{\prime}_{\rm CO}$}$) are tightly correlated ever since the
first statistical submm survey of 100 local FIR-bright galaxies (SLUGS – Dunne
et al., 2000). The first suggestions to use dust as an alternative to CO at
high redshift (e.g. Santini et al., 2010; Magdis et al., 2012; Scoville et
al., 2014) were followed quickly by work demonstrating its potential (e.g
Scoville et al., 2016; Hughes et al., 2017; Orellana et al., 2017).

An advantage of using submm continuum emission from dust as an $\rm H_{2}$ gas
tracer is that it becomes easier to measure at high redshift, because of the
negative $K$-correction (e.g. Blain & Longair, 1993), while recent
technological advances made it possible to image areas large enough to be free
of cosmic variance, leading to the FIR/submm detection of many thousands of
galaxies by the Herschel Space Observatory, for example. The use of dust as a
gas mass proxy requires an estimate of metallicity, since the dust-to-gas
ratio, $\delta_{\rm GDR}$, is roughly proportional to metallicity (e.g. Muñoz-
Mateos et al., 2009; Magdis et al., 2012; Sandstrom et al., 2013; Draine et
al., 2014). The appropriate $\delta_{\rm GDR}$ can then be applied (e.g.
Valentino et al., 2018). Whilst this requirement is often raised as a problem
regarding the use of dust as a gas mass tracer, its dependence on metallicity
is in fact weaker than that of CO333Moreover, since $\rm H_{2}$ cannot be
traced (in bulk) by any of its own lines, regardless of which other tracer (X)
is used (dust emission, CO, or 13CO, or C i line emission), it will always be
necessary to assume a $\rm[X/H_{2}]$ abundance in order to proceed to a final
$\rm H_{2}$ gas mass estimate..

For galaxies selected at FIR/submm/mm wavelengths, it is safe to assume that
the metallicity will be high, such that $\delta_{\rm GDR}$ will be broadly
similar to those found for local metal-rich spirals and the Milky Way (Dunne &
Eales, 2001; Draine, 2009; Magdis et al., 2012; Sandstrom et al., 2013;
Rowlands et al., 2014; Yang et al., 2017; Berta et al., 2021). A detailed
discussion of the advantages and disadvantages of using dust as a tracer of
gas can be found in Genzel et al. (2015) and Scoville et al.
(2017)444Continuum dust emission does not yield information on kinematics,
unlike spectral lines..

A third method of tracing molecular gas – the use of atomic carbon lines – has
come to the fore since ALMA became operational. Its promise was recognised by
Papadopoulos et al. (2004) and its first application as a tracer for molecular
gas mass in galaxies gave good results (Weiß et al., 2003; Papadopoulos &
Greve, 2004), implying that: a) the [C i](${}^{3}P_{1}$–${}^{3}P_{0}$) lines
are optically thin for the bulk of $\rm H_{2}$ gas (Pérez-Beaupuits et al.,
2015) and b) atomic carbon is present throughout CO-rich molecular cloud
volumes.

The latter contradicts the earlier simple plane-parallel PDR model where
atomic carbon (and its line emission) occupied only a thin layer, sandwiched
between $\rm C^{+}$ in the outer and CO in the inner regions of FUV-
illuminated molecular clouds (Tielens & Hollenbach, 1985). However
observations have repeatedly shown excellent concomitance of C i line emission
with CO line emission, by area and by velocity, and C i shows a tighter
correlation with 13CO than with 12CO. C i is now thought to arise from same
volume as the CO, with similar excitation conditions (e.g. Plume et al., 1999;
Ikeda et al., 2002; Schneider et al., 2003; Beuther et al., 2014; Pérez-
Beaupuits et al., 2015). Moreover, it may be that C i lines can also trace CO-
dark molecular gas, should such phase exist in galaxies in significant
amounts, e.g. due to CR-induced dissociation of CO to C (and O) (Bisbas et
al., 2015).

Despite being much fainter than the $\rm C^{+}$ line at 158 $\mu$m (the prime
ISM cooling line), atomic carbon lines do hold certain advantages, namely: a)
they solely trace $\rm H_{2}$ gas, whereas the $\rm C^{+}$ line also traces
the H i and H ii gas reservoirs, which can be significant, especially in
metal-poor systems (e.g. Madden et al., 1997; Liszt, 2011; Papadopoulos &
Geach, 2012; Pérez-Beaupuits et al., 2015; Clark et al., 2019a); b) the C i
lines can remain excited for cold gas (e.g. for [C i](1–0): $E_{10}/k_{\rm
B}\sim 24$ k) unlike the $\rm C^{+}$ line, where the $\Delta E/k_{\rm B}\sim
92$ k will keep it very faint for cold gas; c) the frequencies of the two C i
lines, at 492 and 809 GHz, remain accessible for galaxies over a much larger
redshift range (and thus cosmic volume) than the $\rm C^{+}$ line. In the
latter case, its rest-frame frequency, $\rm\nu(C^{+})\sim 1.9$ THz, means the
$\rm C^{+}$ line is observable by ALMA’s most sensitive receivers only at
$z\ga 4$.

Nevertheless, the high rest-frame frequencies of the C i line made early
observations (and thus any calibration efforts) in the local Universe very
difficult. Initially there had been relatively little observational work
outside of the Milky Way, largely confined to extreme systems such as quasars
and starburst nuclei (e.g. White et al., 1994; Weiß et al., 2005; Walter et
al., 2011). These studies advocated a higher carbon abundance for these
extreme systems, $\mbox{$X_{\rm C\,I}$}=\rm[C^{0}/H_{2}]=5$–$12\times
10^{-5}$, compared to the $\mbox{$X_{\rm C\,I}$}=1$–$2.5\times 10^{-5}$ seen
in the Milky Way (Frerking et al., 1989).

More recently, Herschel observed many local galaxies in C i, although the [C
i](1–0) line was at the edge of the observable range for the Herschel Fourier
Transform Spectrometer (FTS), such that the sensitivity was somewhat
compromised. As a result, most of the detected galaxies were either ULIRGs,
starbursts or low-metallicity dwarfs (Kamenetzky et al., 2014; Rosenberg et
al., 2015; Lu et al., 2017; Jiao et al., 2017). A small sample of normal disk
galaxies was mapped in C i (Jiao et al., 2019, hereafter J19 – see also
Crocker et al., 2019). J19 studied the spatial distribution of
$L^{\prime}_{\rm C\,I}$ and $L^{\prime}_{\rm CO}$ at a $\sim 1$-kpc scale in
15 local galaxies. They concluded that C i is a good tracer of molecular gas,
in the sense that it correlates well with CO and the ratio $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ is distributed smoothly across galaxies.
Comparing against CO(1–0) maps and the independent estimates of $\alpha_{\rm
CO}$ from Sandstrom et al. (2013), these resolved studies suggested $X_{\rm
C\,I}$=1.3–2.5$\times 10^{-5}$, similar to the range in the Galaxy, and that
found by the absorber study of Heintz & Watson (2020).

Table 1:  Samples used for our comparisons.

Sample Selection $z$ $N_{\rm CO}$ $N_{\mbox{C\,{\sc i}}}$ $N_{\rm sub}$ Notes
SF mode References name $\lambda_{\rm obs}$ ($\mu$m) (see below) high-$z$ SMG
850–2000 2–6 89 42 114 Corrected for lensing Both $a$ Local SF 0 35 19 35 C i
from FTS MS $b$ (U)LIRGs 60 0 85 19 114 C i from FTS Both $c$ $z=1$ 850 1 11
18 9 CO(2–1) MS $d$ $z=0.35$ 250 0.35 12 12 12 MS $e$ $0.04<z<0.3$ 160 0-0.3
48 0 54 VALES MS $f$

With ALMA now in routine operations, studies of C i have expanded to a broader
range of galaxies, with a greater variety of average ISM conditions, over a
wider range of redshift. These include SMGs, which lie mainly at $z>1$ (e.g.
Alaghband-Zadeh et al., 2013; Bothwell et al., 2017; Popping et al., 2017;
Oteo et al., 2018; Nesvadba et al., 2019; Dannerbauer et al., 2019; Gómez-
Guijarro et al., 2019), and main-sequence (MS) galaxies at $z=0.35$–1.2
(Valentino et al., 2018; Bourne et al., 2019; Valentino et al., 2020; Dunne et
al., 2021). C i has even been detected in the intracluster medium of the
Spiderweb galaxy cluster at $z=2.16$, as well as in several of its individual
galaxies (Emonts et al., 2018). Routine use of C i as a tracer of molecular
gas is currently limited by the lack of calibration studies to explore and
determine the values and behaviour of the parameters involved, i.e. $X_{\rm
C\,I}$ and $\alpha_{\rm C\,I}$. $\mbox{$X_{\rm C\,I}$}=3\times 10^{-5}$ has
been adopted by almost all recent studies, taken from Weiß et al. (2003),
determined from a comparison of analyses of CO and C i in the centre of M 82,
which has unusually high $\rm[C^{0}/CO]\sim 0.5$, whereas attempts to estimate
$X_{\rm C\,I}$ in other ways – e.g. from absorption studies of Gamma-ray
bursts and quasar absorbers (Heintz & Watson, 2020) – have found lower values,
consistent with the range seen in the Milky Way.

This paper presents the first dedicated cross-calibration study of the dust,
12CO(1–0) and [C i](${}^{3}P_{1}$–${}^{3}P_{0}$) emission in a sample of 407
galaxies from the literature, including MS galaxies and SMGs, such that we can
compare their properties and tracer-($\rm H_{2}$ mass) conversion factors. We
include the 250-$\mu$m-selected galaxies at $z=0.35$ observed with ALMA in all
three tracers by Dunne et al. (2021) where our method was first briefly
presented.

In §2 we describe the samples used in this analysis, the observables, and the
derived quantities. In §3 we describe the Bayesian approach for producing
optimised, self-consistent tracer-($\rm H_{2}$-mass) conversion parameters
between multiple tracers simultaneously. We then examine correlations of the
observables to look for trends in §4. In §5 we investigate the trends we have
found in the conversion factors and provide refined calibration recipes.
Finally, in §6 we discuss the results and highlight the open questions.
Throughout, we use a cosmology with $\Omega_{\rm
m}=0.27,\Omega_{\Lambda}=0.73$ and $H_{0}=71$ km s-1 Mpc-1.

##  2 Deriving observational quantities

###  2.1 Sample

Figure 1: $T_{\rm mw}$ measurements from SED fitting as a function of redshift
(top), $L_{\rm IR}$ (middle) and SED colour $L_{\rm IR}$/$L_{850}$(bottom).
The best-fit line and 2$\sigma$ error are shown as a blue line and shaded
region, the gray-dashed horizontal line shows the value of $T_{\rm mw}$=25 k
used by Sco16. There is a significant correlation with all three observables,
such that the average $T_{\rm mw}$ increases with redshift ($r=0.64$), $L_{\rm
IR}$($r=0.74$) and SED colour ($r=0.80$). Fit parameters are listed in Table
LABEL:fitsT.

The samples used in our study are those available in the literature – up-to-
date as of early 2022 – which have at least two of the three tracers: submm
dust continuum emission at $\lambda_{\rm rest}>500$$\mu$m, 12CO(1–0) or (2–1),
and [C i](${}^{3}P_{1}$–${}^{3}P_{0}$). Summarising: 326 galaxies have both CO
and submm continuum detections; 140 have both C i and submm dust continuum
detections; 109 have both C i and CO detections; 101 have all three tracers.
The sample covers the redshift range $0<z<6$, and includes galaxies lying
within 1 dex of the MS as well as extreme starbursts such as local ULIRGs and
most high-$z$ submm-selected galaxies. Full details and references are listed
in Table 1. Lensed galaxies are included only where there is an estimate of
the magnification, $\mu$, and all luminosities have been corrected by the
magnification factor. Our sample includes the galaxies from one of the most
comprehensive studies of dust as a tracer of molecular gas across cosmic time
- Scoville et al. (2016), henceforth Sco16555Although lensed galaxies were
included in their work, the luminosities were not de-magnified.. The Scoville
et al. sample has been updated as described in Appendix A.

In order to test for the effect of SF intensity or ‘SF-mode’ on any later
results, we divide the sample into two groups, referred to hereafter as ‘MS
galaxies’ and ‘SMGs’, the names of the groups are not meant to be accurate
definitions but rather a reference to familiar categories. For this
heterogeneous data-set, defining a simple criterion for two groups is not
possible, and even if it were, a fuzzy boundary would still remain due to
measurement errors and the inability to capture the complexity in a single
parameter. The extreme starburst ‘SMG’ group contains the high-redshift
submillimeter selected galaxies which were discovered in the pre-ALMA era and
as such are extreme star forming systems (else they could not have been
detected), plus the local ULIRGs and some LIRGs which have evidence for very
intense and obscured regions (e.g. NGC 4418, IC 860) where conditions are
likely to be extreme (Díaz-Santos et al., 2017; Falstad et al., 2021). The ‘MS
galaxy’ group contains the lower luminosity local disk galaxies plus the LIRGS
which are not extreme, the intermediate redshift sources selected at 250$\mu$m
from the Herschel-ATLAS – $z=0.35$ galaxies from Dunne et al. (2021) and the
$z<0.3$ VALES galaxies (Hughes et al., 2017), the $z\sim 1$ galaxies
(Valentino et al., 2018; Bourne et al., 2019) and the ASPECs sources denoted
as ‘MS’ in that survey (Boogaard et al., 2020). (Full references are provided
in Table 1.)

There are two situations where corrections to luminosities may be required:

#####  H i-dominated galaxies at low L${}_{\textbf{IR}}$.

For galaxies with a large fraction of H i within their optical disk, their
dust tracing H i rather than H2 makes a significant contribution to the submm
continuum emission. Since our intention is to provide a calibration for H2
rather than total gas, we apply a correction to $L_{850}$ for galaxies with
$\mbox{$f_{\rm H\,I}$}=\mbox{H\,{\sc i}}/\mbox{H${}_{2}$}>1$, as described in
Appendix B.1. Galaxies corrected in this way are shown as cyan diamonds in the
plots.

##### Local galaxies mapped in C i by the Herschel FTS.

The local galaxies mapped using the Herschel FTS by J19 present some complex
issues. Some do not have C i and dust continuum measurements in matched
apertures, and those same galaxies are often only detected in C i in the inner
few kpc of the galaxy, where the ratios of $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ may also be biased – for example, by a lower
$\alpha_{\rm CO}$ in galaxy centres. We discuss the issues in more detail in
§4.4 and Appendix B.2. Galaxies requiring a significant correction ($>0.1$
dex) to $L^{\prime}_{\rm C\,I}$ are labelled as C icor; they are shown in the
plots as pink diamonds, but not included in the analysis unless specified.

###  2.2 Observables

We will compare three tracers of molecular gas, where the observables (the
luminosities $L_{850}$, $L^{\prime}_{\rm CO}$ and $L^{\prime}_{\rm C\,I}$) are
empirically related to the molecular gas mass as

$\mbox{$M_{\rm
mol}$}=\mbox{$L_{850}$}/\mbox{$\alpha_{850}$}=\mbox{$\alpha_{\rm
CO}$}\mbox{$L^{\prime}_{\rm CO}$}=\mbox{$\alpha_{\rm
C\,I}$}\mbox{$L^{\prime}_{\rm C\,I}$}$ (1)

The goal of this analysis is to determine self-consistent conversion factors
$\alpha_{850}$, $\alpha_{\rm CO}$ and $\alpha_{\rm C\,I}$ and study the
physical properties they depend on, e.g. C abundance, gas-to-dust ratio
($\delta_{\rm GDR}$), dust emissivity. Our definition of the ‘observables’ is
intended to be independent of as many assumptions as possible. For CO and C i,
we use $L^{\prime}$ as defined by Solomon & Vanden Bout (2005):

$L^{\prime}=\frac{3.25\times 10^{7}}{\nu_{\rm rest}^{2}}\left(\frac{D_{\rm
L}^{2}}{1+z}\right)\left[\frac{\int_{\Delta V}Sdv}{\rm
Jy\,km\,s^{-1}}\right]\,\,\,\,{\rm K\,km\,s}^{-1}\,{\rm pc}^{2},$ (2)

where $\int_{\Delta V}Sdv$ is the velocity-integrated line flux density,
$D_{\rm L}$ is the luminosity distance (Mpc), and $\nu_{\rm rest}$ is the rest
frequency666Where we use $\nu_{\rm rest}$, Solomon & Vanden Bout (2005) use
$\nu_{\rm obs}$, hence the different exponent for $(1+z)$ cf. their equation
(3). of the transition in GHz.

Most of the galaxies in our compilation have been observed in the
12CO($J=1$–0) transition. However, some observations at high redshift target
the 12CO($J=2$–1) line. We convert $L^{\prime}_{21}$ to $L^{\prime}_{10}$
using the line luminosity ratio $R_{21}=0.8$; if instead we were to set
$R_{21}$ to unity, this would not affect any of our conclusions. We do not use
$J\geq 3$ CO lines because the uncertainties in the global excitation
corrections become too large for a useful calibration study.

We use only the C i ${}^{3}P_{1}$–${}^{3}P_{0}$ line, as it is the least
sensitive to the average excitation conditions, and correlates better with the
low-$J$ CO emission (Jiao et al., 2017; Jiao et al., 2019; Crocker et al.,
2019). Moreover, there is now evidence of strongly sub-thermal excitation for
C i(2–1) (Harrington et al., 2021; Papadopoulos et al., 2022), making it
difficult to use this line as an $\rm H_{2}$ mass tracer since its excitation
is extremely uncertain.

For the dust continuum emission, we use $L_{850}$, the luminosity at rest-
frame 850 $\mu$m:

$\mbox{$L_{850}$}=4\pi S_{\rm{\nu(obs)}}\times K\left(\frac{D_{\rm
L}^{2}}{1+z}\right)\,\,\,\,\,{\rm W\,Hz}^{-1},$ (3)

where $D_{\rm L}$ is the luminosity distance, $S_{\rm{\nu(obs)}}$ is the
observed flux density and $K$ is the $K$-correction to rest-frame 850$\mu$m,
defined as

$K=\left(\rm{\frac{353\,GHz}{\nu_{\rm
rest}}}\right)^{3+\beta}\,\left(\frac{e^{\rm{h\nu_{\rm rest}/k\mbox{$T_{\rm
d}$}}}-1}{e^{16.956/\mbox{$T_{\rm d}$}}-1}\right).$ (4)

Here, $\nu_{\rm rest}=\nu_{\rm obs}(1+z)$, $T_{\rm d}$ is the luminosity-
weighted dust temperature, from an isothermal fit to the spectral energy
distribution (SED) with the dust emissivity, $\beta$, allowed to vary between
1.8–2.0.

Sco16 assumed a $T_{\rm d}$=25 k and $\beta=1.8$, respectively, to extrapolate
(or $K$ correct) their observed submm luminosities to rest-frame 850 $\mu$m.
We make full use of the available data to refine this procedure as follows: 1)
With sufficient data, we fit the SED ourselves with $\beta=1.8$ and estimate
the rest-frame 850 $\mu$m luminosity directly from the SED fit. 2) Failing
that, we use the reported $T_{\rm d}$ in the literature to make the
extrapolation from the longest wavelength measurement available. 3) For SMGs
with insufficient data points to have had their SED fitted, we adopt their
observed average, $T_{\rm d}$=38 k (da Cunha et al., 2015). The bulk of the
high-$z$ samples now have observations between 2-3 mm with ALMA, as such the
extrapolation to rest-frame 850$\mu$m is small, even at the highest redshifts.
The shortest rest-frame wavelengths we deal with are $\lambda_{r}\sim
250$$\mu$m for sources at $z\sim 2-3$ observed at 850$\mu$m, which require
K-corrections in the range 50–140. However, the important consideration is the
potential uncertainty in that K-correction, not its absolute value. We tried
alternatively using the Sco16 method of assuming $T_{\rm d}$=25 K to
extrapolate to rest-frame $L_{850}$ and found a maximum difference of a factor
1.6, with the average being 1.15 times. The true uncertainty due to the SED
sampling and K-correction will be smaller than this, as we know from our work
in Section 2.3.1 and Figure 4 that the dust temperatures in SMG are much
higher than 25 K. We thus do not consider the extrapolation to rest-frame
850$\mu$m to be a significant source of uncertainty or bias in this analysis.

###  2.3 Physical dependencies of gas mass tracers

####  2.3.1 Dust–$\rm{H_{2}}$ calibration

Large dust grains ($a\sim 0.1$$\mu$m) in thermal equilibrium with their
incident radiation field emit as a modified black body (MBB), where the
emission is related to the mass of hydrogen as:

$\mbox{$M_{\rm H2}$}=\frac{L_{\nu}}{4\pi B(\nu,\mbox{$T_{\rm
mw}$})}\,\,\,\mbox{$\kappa_{\rm H}$}(\nu).$ (5)

The two physical quantities needed to calibrate dust continuum emission as a
tracer of gas are therefore $T_{\rm mw}$ and $\kappa_{\rm H}$. Expressing Eqn
5 in astronomical units for $\lambda=850$ $\mu$m, we can write:

$\frac{\mbox{$M_{\rm H2}$}}{[\mbox{$\rm M_{\odot}$}]}=6.14\times
10^{-14}\frac{\mbox{$\kappa_{\rm
H}$}}{[\rm{kg\,m^{-2}}]}\frac{\mbox{$L_{850}$}}{[\rm{W\,Hz^{-1}}]}\left(\frac{24.5}{\mbox{$T_{\rm
mw}$}}\right)^{-1.4},$ (6)

where we have simplified the exponential term in the Planck function as
$\sim(24.5/\mbox{$T_{\rm mw}$})^{1.4}$ for $17<\mbox{$T_{\rm mw}$}<30$ k.

The mass-weighted dust temperature, $T_{\rm mw}$, is often lower than the
luminosity-weighted dust temperature, $T_{\rm d}$, as derived from an
isothermal MBB fit to the dust SED, because warm dust outshines cold dust per
unit mass. There is an excellent discussion of this in the Appendix to Sco16
which we will not repeat here (see also Dunne & Eales, 2001).

To determine $T_{\rm mw}$, we require a multi-component MBB fit to a well-
sampled dust SED (e.g. Dunne & Eales, 2001), or an SED fit using a model that
allows a range of radiation field strengths, leading to a range of dust
temperatures (e.g. Draine et al., 2007). These methods give broadly consistent
results. As a rule of thumb, the range of $T_{\rm mw}$ in local star-forming
galaxies is 15–25 k (Dunne & Eales, 2001; Draine et al., 2007; Hunt et al.,
2015; Dale & Helou, 2002; da Cunha et al., 2008; Bendo et al., 2015; Clark et
al., 2015), increasing to 25–30 k in luminous starbursts at higher redshifts
(Rowlands et al., 2014; da Cunha et al., 2015).

As the dust ($\alpha_{850}$) factor is only weakly dependent on the assumed
temperature at rest-frame 850 $\mu$m, Sco16 and others assumed a constant
$T_{\rm mw}$=25 k. If instead the true $T_{\rm mw}$ were to be 15 [30] k, the
dust and gas mass would be under-[over-]estimated by a factor $\sim 2\times$
[$1.3\times$], which is overshadowed by the other uncertainties. On the other
hand, failure to account for any systematic trend of $T_{\rm mw}$ with another
physical parameter can introduce or mask correlations of the conversion
factors with that physical parameter.

We therefore explore the validity of assuming constant $T_{\rm mw}$ by
collating measurements of $T_{\rm mw}$ from the literature (Dunne & Eales,
2001; Hunt et al., 2019), and additionally making our own fits where possible.
In Fig. 1 we show that there are indeed strong correlations of $T_{\rm mw}$
with the observables, namely $z$, $L_{\rm IR}$ and the SED colour, $L_{\rm
IR}$/$L_{850}$. There is a clear difference between samples with low and high
SFRs, with $\langle\mbox{$T_{\rm mw}$}(\rm MS)\rangle=23.0\pm 0.4$ k while
$\langle\mbox{$T_{\rm mw}$}(\rm SMG)\rangle=30.1\pm 0.7$ k. We fit empirical
relations for the correlations in Fig. 1 (see Table LABEL:fitsT). Where there
is no direct estimate of $T_{\rm mw}$ from an SED fit, which is the case for
two-thirds of galaxies, we use these empirical relations777We restrict the
predicted $T_{\rm mw}$ such that $\mbox{$T_{\rm mw}$}\leq\mbox{$T_{\rm d}$}$.
to derive $T_{\rm mw}$ for use in our subsequent analysis. Appendix H.2
compares our approach, where we use individual estimates of $T_{\rm mw}$, to
the adoption of a constant $T_{\rm mw}$=25 k, where we will discuss those
findings in §5.

In their work on strongly lensed SMGs at high redshift, which includes dust
continuum emission as a constraint in a large-velocity gradient (LVG) model,
Harrington et al. (2021) find that $T_{\rm d}$$\sim$ $T_{\rm mw}$ for most
SMGs, with both measures of temperature higher than the $T_{\rm mw}$=25 k
commonly used in the literature. To ensure consistency with our other
estimates of $T_{\rm d}$, we fitted the Harrington et al. (2021) photometry
with three simple models: 1) an isothermal optically thin MBB; 2) an MBB with
variable optical depth and – where there were enough data – 3) a two-component
MBB. In agreement with Harrington et al., we find that a single dust
temperature adequately describes the SED of these galaxies, in contrast to
lower redshift (U)LIRGs and normal galaxies which are better fit with multiple
dust components ($T_{\rm d}$ $>$ $T_{\rm mw}$) and/or fits with higher FIR
optical depths. The temperatures from the Harrington et al. turbulence model
correlate best with our isothermal $T_{\rm d}$ meansurements for these
galaxies, and the temperatures returned when allowing variable optical depth,
$\mbox{$T_{\rm d}$}(\tau)$, are always significantly higher than those from
the Harrington et al. model. We therefore do not use optically thick fits to
yield $T_{\rm mw}$ for our high-redshift SMGs. We instead use two-component
SED fits to the lensed Planck sources and the handful of SMGs with sufficient
data for the empirical relations shown in Fig. 5.

The other key physical parameter in the dust–H2 conversion is $\kappa_{\rm
H}$, which is a combination of the dust mass absorption coefficient
($\kappa_{850}$) and $\delta_{\rm GDR}$, such that888Literature studies
generally present the dust-H2 conversion in terms of $\delta_{\rm GDR}$ for a
fixed emissivity, $\kappa_{850}$. Given the mounting evidence that
$\kappa_{850}$ varies within our own (Remy et al., 2017; Ysard et al., 2015;
Ysard et al., 2018; Köhler et al., 2015) and other galaxies (e.g. Clark et
al., 2019b, but see also Priestley & Whitworth 2020), we prefer to work with
$\kappa_{\rm H}$ to avoid projecting all the variation in the $\rm H_{2}$-dust
conversion factor onto $\delta_{\rm GDR}$. $\kappa_{\rm H}$= $\delta_{\rm
GDR}$/$\kappa_{850}$. Briefly, $\kappa_{850}$ is sensitive to the grain
composition and structure (amorphous; crystalline; coagulated; mantled), while
$\delta_{\rm GDR}$ is roughly proportional to metallicity and, for galaxies
with metallicity within a factor 2 of $Z_{\odot}$, as expected for those in
our samples, can be taken to be roughly constant, at $\delta_{\rm
GDR}$$=100$–150 (Sodroski et al., 1997; Dunne et al., 2000; Dunne & Eales,
2001; Draine et al., 2007; Muñoz-Mateos et al., 2009; Leroy et al., 2011;
Sandstrom et al., 2013; Planck Collaboration et al., 2011; Jones et al., 2017;
De Vis et al., 2021).

Fortunately, observational measures of $\kappa_{\rm H}$ are available, both
for the Milky Way and for external galaxies, with values of $\kappa_{\rm H}$
$\sim 2000$ kg m-2 in the Milky Way’s diffuse interstellar medium (ISM) and
$\kappa_{\rm H}$ $\sim 800$ kg m-2 in dense clouds. Appendix C discusses in
more detail how it is measured, and Table LABEL:kappalitT provides a
comprehensive set of observational and theoretical values for $\kappa_{\rm H}$
from the literature.

It is impossible to disentangle the effect of changing dust properties
($\kappa_{850}$) from changes in $\delta_{\rm GDR}$ in observational
determinations of $\kappa_{\rm H}$. While the decrease in $\kappa_{\rm H}$
towards denser sightlines in the Milky Way is thought to be due to the dust
grains coagulating in denser environments – a process expected to increase
their emissivity (e.g. Köhler et al., 2015) – there may also be some decrease
in $\delta_{\rm GDR}$ if the gas is accreted into dust mantles or ices (i.e.
grain growth). Both effects are to be expected (e.g. Jones et al., 2017;
Jones, 2018) and both act to decrease $\kappa_{\rm H}$. Counter to that, the
higher estimates of $\kappa_{\rm H}$ in the diffuse atomic phase (lowest
$N_{\rm H}$ sightlines at high latitudes) in the Milky Way may be due in part
to a lower dust emissivity for grains without ice mantles, where only the
refractory cores remain, subjected to harsher ultraviolet (UV) irradiation.
Additionally, there is likely a metallicity gradient at high latitudes,
leading to a higher $\delta_{\rm GDR}$, further increasing $\kappa_{\rm H}$.
There is thus a qualitative expectation that denser regions with higher
metallicity will have higher dust emissivity, $\kappa_{850}$, and lower
$\delta_{\rm GDR}$, producing a lower $\kappa_{\rm H}$. More diffuse regions
with lower metallicity will move in the opposite direction. In §5, we find
that we can constrain the range of $\kappa_{\rm H}$, at least, and therefore
the combination of $\delta_{\rm GDR}$/$\kappa_{850}$.

####  2.3.2 $\rm C\,I$–$\rm{H_{2}}$ calibration

Here, we introduce the two physical parameters the $X_{\rm
C\,I}$=$\rm[C^{0}/H_{2}]$ average abundance ratio and the average excitation
factor $Q_{10}$$=\rm N_{1}/N_{tot}$, pertinent to the use of C i as a tracer
of H2. The relationship between $M_{\rm H2}$ and the ‘observable’ – [C i](1–0)
line emission – is (in astronomical units):

$\mbox{$M_{\rm H2}$}({\rm M}_{\odot})=\frac{0.0127}{X_{\rm
C\,I}\,\,Q_{10}}\left(\frac{D_{\rm
L}^{2}}{1+z}\right)\,\,\left[\frac{\int_{\Delta V}S_{\rm[CI](1-0)}dv}{\rm
Jy\,km\,s^{-1}}\right]$ (7)

with $D_{\rm L}$ in Mpc and $\int_{\Delta V}S_{\rm[CI](1-0)}\Delta v$ in Jy km
s-1. Expressed in units of line luminosity, this becomes:

$\mbox{$M_{\rm H2}$}(\mbox{$\rm M_{\odot}$})=\frac{9.51\times
10^{-5}}{\mbox{$X_{\rm C\,I}$}\,Q_{10}}\,\mbox{$L^{\prime}_{\rm C\,I}$}.$ (8)

The excitation term, $Q_{10}$, describes the relative fraction of carbon atoms
in the $J=1$ state. Under general non-LTE conditions it is a function of both
gas density, $n$, and $T_{\rm k}$ and is derived analytically in the Appendix
to Papadopoulos, Thi & Viti (2004). A recent study of the [C i](2–1)/(1–0)
line ratio (Papadopoulos, Dunne & Maddox, 2022) finds that the C i lines are
both sub-thermally excited in the ISM of galaxies, with the [C i](2–1)
especially so (see also Harrington et al., 2021). Thus the LTE expressions for
$Q_{10}$ should not be used, nor will the C i line ratio produce an estimate
of $T_{\rm k}$ (both methods having been widely used in the literature to
date). Details for $Q_{10}$ are in the Appendix D, but in summary we find:

  1. 1.

The [C i](1–0) excitation term, $Q_{10}$, is a non-trivial function of density
and temperature, but for the range $\mbox{$T_{\rm k}$}\geq 20$ K and log
$n\geq 2.5$ – which is where the bulk of H2 in star forming galaxies is
thought to reside – $\langle\mbox{$Q_{10}$}\rangle=0.48\pm 0.08$ where the 99
per cent confidence range is quoted (see Papadopoulos et al. 2022 and Figure
13 for details).

  2. 2.

Due to a slight super-thermal behaviour, higher density, higher $T_{\rm k}$
conditions can produce similar or even lower $Q_{10}$ than lower density,
lower $T_{\rm k}$ conditions. This breaks any intuitive link between $Q_{10}$
and the ISM conditions, i.e. we do not necessarily expect a higher $Q_{10}$ in
SMGs compared to MS galaxies (see Fig. 13).

  3. 3.

As the [C i](2–1) line is even more strongly sub-thermally excited, its
$Q_{21}=N_{2}/N_{\rm tot}$ factor varies strongly999The $Q_{21}(n,T_{k})$ that
enters the estimates of molecular gas mass when the [C i](2–1) line is used
can vary almost by a factor of $\sim 5$, depending on $(n,T_{k})$.. This is
the main reason why our current study is restricted to the [C i](1–0) line.

As the [C i](${}^{3}P_{1}$–${}^{3}P_{0}$) line is optically thin for most
conditions expected in spiral disks (Weiß et al., 2005; Pérez-Beaupuits et
al., 2015; Harrington et al., 2021), the relationship between $L^{\prime}_{\rm
C\,I}$ and $M_{\rm H2}$ is proportional to $X_{\rm C\,I}$ – the abundance of
carbon atoms relative to H2. This dependence on abundance is as expected for
any method that employs tracers of $\rm H_{2}$ gas mass other than the $\rm
H_{2}$ lines themselves101010Even for optically thick tracers of $\rm H_{2}$
gas, such as CO(1–0) line emission, a $\rm[CO/H_{2}]$ abundance still enters
the method via the CO–H2 cloud volume-filling factor, $f_{\rm CO}$, albeit not
in a sensitive fashion unless a combination of strong FUV radiation and/or low
metallicities selectively dissociate CO in the outer cloud layers while
leaving the largely self-shielding $\rm H_{2}$ intact (then $f_{\rm CO}$ can
be $\ll 1$, see Pak et al. 1998 for details)..

With the excitation factor $Q_{10}$ varying no more than 16 per cent over the
typical range of H2 conditions in galaxies ($\mbox{$T_{\rm k}$}\geq 20$ K, log
$n\geq 2.5$), the major source of uncertainty in C i-based molecular gas mass
estimates (and thus the major source of scatter in the $\alpha_{\rm C\,I}$
conversion factor) is the neutral carbon abundance, $X_{\rm C\,I}$. The
relatively recent introduction of the [C i](1–0) line as a gas tracer means
that $X_{\rm C\,I}$ has not been widely explored – constraining it and
investigating any potential trends is a key outcome of our cross-calibration
work.

In the Milky Way, $X_{\rm C\,I}$ is found to vary only modestly, from
0.8–$2.2\times 10^{-5}$ (e.g. Zmuidzinas et al., 1988; Frerking et al., 1989;
Tauber et al., 1995; Ikeda et al., 2002), while a much higher value
($\mbox{$X_{\rm C\,I}$}=5\times 10^{-5}$) has been inferred for the nearby
starburst nucleus of M 82 (Schilke et al., 1993; White et al., 1994; Stutzki
et al., 1997)111111The measurement is in fact the $\rm[C^{0}/CO]$ abundance,
and a value for [CO/H2] has then to be assumed to infer $X_{\rm C\,I}$..
Thanks to ALMA, very high localised ratios of $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ (translating to high $X_{\rm C\,I}$=5-7$\times
10^{-5}$) have also been measured in extreme regions, such as the Circum-
Nuclear Disk (CND) of NGC 7469 which is believed to host an X-ray Dominated
Region (XDR) (Izumi et al., 2020) and the outflow region in NGC 6240 (Cicone
et al., 2018). More modestly elevated $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm
CO}$ ratios tend to be found in the central nuclear regions of starburst
galaxies (Jiao et al., 2019; Salak et al., 2019; Saito et al., 2020). However,
when averaged over larger kpc scale regions – the ratios become consistent
with the average global ratios measured for this sample (see Figure 5).
Independent measurements of $\mbox{$X_{\rm C\,I}$}=1.6^{+1.3}_{-0.7}\times
10^{-5}$ (for solar metallicity) were made by Heintz & Watson (2020) using UV
absorption measures for a range of absorber systems across cosmic time. Cosmic
rays (and X-rays) are expected to dissociate CO in favour of atomic carbon,
increasing $[{\rm C}^{0}/{\rm CO}]$, a hypothesis supported by both
simulations and observations (e.g. Bisbas et al., 2015; Clark et al., 2019a;
Israel, 2020; Izumi et al., 2020).

####  2.3.3 CO–$\rm{H_{2}}$ calibration

The 12CO(1–0) line is optically thick in most (but not all see Aalto et al.,
1995, e.g.) ISM conditions expected in galaxies. Unlike dust continuum
emission where optical depths build up over large columns of dust, the entire
CO line optical depth builds up within very small gas ‘cells’ ($<$0.1 pc) due
to the very turbulent nature of the velocity fields, and the small thermal
line widths (Tauber et al., 1991; Falgarone et al., 1998). This localised
nature of CO line optical depths and the macro-turbulent CO line formation
mechanism allows a great simplification of the radiative transfer models of
such lines, i.e. the use of the so-called Large Velocity Gradient (LVG)
approximation. However, it also complicates the relationship between the CO
line luminosity and the underlying $\rm H_{2}$ gas mass, making the
corresponding conversion factor, $\alpha_{\rm CO}$, dependent on the thermal
state of the gas, its average density, as well as its dynamic state.

Following Papadopoulos et al. (2012b) the $\alpha_{\rm CO}$ factor in an LVG
setting is given by:

$\mbox{$\alpha_{\rm CO}$}=2.65\frac{\sqrt{n_{\rm H2}}}{T_{\rm b}}\,K_{\rm
vir}^{-1}\,\,\,\,[\mbox{M${}_{\odot}$\,({\sc
k}\,km\,s${}^{-1}$\,pc${}^{2}$)${}^{-1}$}]$ (9)

where $n_{\rm H2}$ and $T_{\rm b}$ are the average density (in cm-3) and the
CO(1–0) brightness temperature121212Here the cloud CO-H2 volume filling factor
is set $f_{\rm CO}=1$. for the molecular cloud ensemble while $K_{\rm vir}$
describes the average dynamic state of the gas (self-gravitating clouds
$K_{\rm vir}\sim 1$, unbound clouds $K_{\rm vir}>1$). In principle, multi-
phase LVG models of CO (and 13CO) SLEDs can be used to constrain $\alpha_{\rm
CO}$, but in practice this demands large line datasets per galaxy (e.g.
Papadopoulos et al., 2014; Harrington et al., 2021), making it impractical for
use in large galaxy samples. This is why in our study $\alpha_{\rm CO}$
remains an empirical conversion factor to be (cross)-calibrated.

####  2.3.4 Conversion factors and physical parameters

The two optically thin tracers – thermal dust continuum emission and C i –
have a simple relation between the empirical ‘mass-to-light’ conversion
parameter ($\alpha_{\rm X}$) and the physical conditions in the ISM (e.g.
abundance, emissivity, temperature). We can write the empirical factors (Eqn.
1) in terms of these physical parameters as follows131313Hereafter we omit the
units for $\alpha_{\rm CO}$, $\alpha_{\rm C\,I}$ and $\alpha_{850}$.:

$\mbox{$\alpha_{850}$}=\frac{1.628\times 10^{16}}{1.36\,\mbox{$\kappa_{\rm
H}$}}\left(\frac{24.5}{\mbox{$T_{\rm mw}$}}\right)^{-1.4}\,\,\,{\rm
W\,Hz}^{-1}{\rm M}_{\rm mol}^{-1}\,\,,$ (10)

where the factor 1.36 corrects to total molecular mass, including He.

$\mbox{$\alpha_{\rm C\,I}$}=16.8\,\left[\frac{\mbox{$X_{\rm C\,I}$}}{1.6\times
10^{-5}}\right]^{-1}\left[\frac{Q_{10}}{0.48}\right]^{-1}\,\,\,\rm{\mbox{$\rm
M_{\odot}$}\,(K\,km\,s^{-1}\,pc^{2})^{-1}}$ (11)

Eqn. 11 also includes the factor 1.36 for He.

##  3 Deriving self-consistent calibration of conversion factors

We next describe how we combine the measurements of multiple gas tracers in
the most efficient way, in order to determine their cross-calibrations. Our
goal is to find the empirical conversion factors ($\alpha_{850}$, $\alpha_{\rm
C\,I}$, $\alpha_{\rm CO}$) or physical parameters ($\kappa_{\rm H}$, $X_{\rm
C\,I}$), which produce a consistent estimate for $M_{\rm H2}$ in a given
galaxy.

Our dataset provides nested samples, each with a different set of available
gas tracers. The daX sample has all three tracers available – dust continuum,
CO and C i, and contains 101 galaxies ($N_{\rm daX}$ = 101). The names and
statistics for the other samples are as follows: ad – CO and dust, $N_{\rm
ad}$ = 326; Xd – C i and dust, $N_{\rm Xd}$ = 140; Xa – C i and CO, $N_{\rm
Xa}$ = 109. The properties of these samples are in Table 2.

The best constraints at log10 $L_{\rm IR}$ $>11$141414Hereafter we will refer
to $\log_{10}$ as simply log. come from the daX sample because it has three
independent tracers of gas mass, but it lacks coverage of luminosities below
$l_{\rm IR}=11$. The ad sample is the largest and spans the widest range in
$L_{\rm IR}$, reflecting the longer time for which CO observations have been
possible for nearby galaxies.

We begin with the daX sample, to illustrate the method of optimisation for the
estimates of all three conversion factors simultaneously.151515This method was
first presented in brief in Dunne et al. (2021), where it was applied to the
sample of $z=0.35$ galaxies. There are four unknowns namely: $m={\rm
log}(\mbox{$M_{\rm H2}$})$, $X={\rm log}(\mbox{$X_{\rm C\,I}$})$, $\kappa={\rm
log}(\mbox{$\kappa_{\rm H}$})$, and $\alpha={\rm log}(\mbox{$\alpha_{\rm
CO}$})$, and three observables: $L^{\prime}_{\rm CO}$, $L^{\prime}_{\rm C\,I}$
and $L_{850}$.

With an independent measure of the true $M_{\rm H2}$, the observables would
provide direct estimates of the three conversion factors; however, the value
of $M_{\rm H2}$ is not known a priori, so we must use a probabilistic argument
based on the fact the observations do provide constraints on the relative
values of the conversion factors for each galaxy. There is thus a set of self-
consistent conversion factors which link the observables to the true $M_{\rm
H2}$, with an unknown common constant factor.

The Bayesian approach we use is described in detail in Appendix F and requires
an estimate of the intrinsic scatter for the logarithms of each of the
factors: $s_{\rm X}$, $s_{\kappa}$ and $s_{\alpha}$. The observable
luminosities relate to these factors as follows, where the coefficients of
proportionality are listed in Table 3:

$\displaystyle\frac{\mbox{$L^{\prime}_{\rm CO}$}}{\mbox{$L^{\prime}_{\rm
C\,I}$}}\propto\mbox{$\alpha_{\rm CO}$}\mbox{$X_{\rm C\,I}$},$
$\displaystyle\frac{\mbox{$L_{850}$}}{\mbox{$L^{\prime}_{\rm
C\,I}$}}\propto\mbox{$\kappa_{\rm H}$}\mbox{$X_{\rm C\,I}$},\,$
$\displaystyle\frac{\mbox{$L^{\prime}_{\rm
CO}$}}{\mbox{$L_{850}$}}\propto\frac{\mbox{$\alpha_{\rm
CO}$}}{\mbox{$\kappa_{\rm H}$}}.$ (12)

We begin by measuring the intrinsic scatter between the three pairs of
observables using an orthogonal distance regression (ODR) fitting method,
which includes the intrinsic scatter, $\lambda$, as a third parameter in the
analysis161616We need to multiply $\lambda$ from the ODR fitting routine by
$\sqrt{2}$ because we need to know the intrinsic scatter of $X-Y$ in our
dataset in order to determine the intrinsic scatter of each conversion factor
in turn. (see Appendix I for full details). The three pair variances derived
from the data are then used to estimate the intrinsic variance of the three
individual conversion factors (the derivation can be found in Appendix E). The
values of the intrinsic scatter for the parameters are given in Table 3, with
$X_{\rm C\,I}$ having the smallest scatter between galaxies. This finding is
purely empirical, requiring no assumptions about the values or trends of the
conversion factors, and as such is very interesting.

Table 2: Samples used in cross-calibration analysis.

Sample Tracers present $N$ median log $L_{\rm IR}$ daX Dust, CO and C i 101
(90) 11.65 (11.77) Xa CO and C i 109 (97) 11.66 (11.88) Xd Dust and C i 140
(128) 11.88 (12.06) ad CO and dust 326 (240) 11.54 (12.07)

Table 3:  Summary of the parameters required to reproduce this analysis.

Quantity Set C i CO Dust physical $X_{\rm C\,I}$ $\alpha_{\rm CO}$
$\kappa_{\rm H}$ empirical $\alpha_{\rm C\,I}$ $\alpha_{\rm CO}$
$\alpha_{850}$ $s_{X,\alpha,\kappa}$ log $L_{\rm IR}>11$ 0.082 0.1646 0.1339
BL 0.1125 0.1436 0.1294

Figure 2:  Tracer luminosity ratios as a function of lensing, redshift and
SFR. $p_{\rm KS}$ is the probability of no difference from a two-sided KS
test. Left column: $L_{850}$/$L^{\prime}_{\rm CO}$ – there are no indications
of any bias due to lensing or redshift. The third panel shows the effect of
high $f_{\rm H\,I}$ in MS galaxies (we have not applied the H i correction
factor to show the difference in the raw ratio). The green-shaded region
consists of MS galaxies with log $\mbox{$L^{\prime}_{\rm CO}$}>8.9$ and
$\mbox{$f_{\rm H\,I}$}<1$. The black line is for all galaxies with
$log\mbox{$L^{\prime}_{\rm CO}$}<8.9$, including the lo-VALES galaxies
suspected of having $\mbox{$f_{\rm H\,I}$}>1$. Galaxies known to have
$\mbox{$f_{\rm H\,I}$}>1$ are shown as the dotted line to illustrate the
similarity. The MS sub-group in the lower panel excludes the $\mbox{$f_{\rm
H\,I}$}>1$ and lo-VALES galaxies. Middle column: $L_{850}$/$L^{\prime}_{\rm
C\,I}$ – there are no indications of bias due to lensing or redshift. The
third panel shows the C icor galaxies from J19 compared to the rest of the MS
sub-group; even with the aperture correction applied to the C i flux (Appendix
4.4), they have a significantly different distribution of observed ratios. In
the lower panel, we exclude C icor galaxies and find no difference between
galaxies with high and low SFRs. Right column: $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ – the highest redshift SMGs ($z>3$) have higher
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ at marginal significance
($p=0.012$), but there are only 11 galaxies at $z>3$. The third panel shows
once more the difference between the C icor and other galaxies in both high
and low luminosity bins. Even though the C i and CO are measured in the same
apertures, the C icor galaxies have much lower $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$. The bottom panel shows a low signficance
difference between high and low SFR galaxies ($p=0.056$), becoming significant
when including the C icor galaxies. Means and KS results are given in Table
LABEL:hypT.

As we do not have any independent measure of the gas mass with which to
normalise our cross-calibration (four unknowns but only three measurements),
we must make an assumption about the sample average of one of the physical or
empirical conversion factors. However, with that assumption made transparent,
the individual values can always be scaled to whichever normalisation a reader
wishes to adopt. The relative values, however, are always the optimal
solution.

We choose to use the dust parameter $\kappa_{\rm H}$= $\delta_{\rm
GDR}$/$\kappa_{850}$ for this normalisation, because there are no trends of
$L_{850}$/$L^{\prime}_{\rm CO}$ or $L_{850}$/$L^{\prime}_{\rm C\,I}$ with
$L_{\rm IR}$ (see Figs 2, 14) and $\kappa_{\rm H}$ also has the best
observational constraints.

For the Xa sample, where there are no dust continuum measurements, we
normalise to $\mbox{$X_{\rm C\,I}$}^{\rm N}=1.6\times 10^{-5}$, which is the
mid-range of the values suggested by the independent study of absorption lines
by Heintz & Watson (2020).

All the galaxies in our sample are metal rich ($0.5<Z/Z_{\odot}<2$) and so we
assume that $\kappa_{\rm H}$ ($\delta_{\rm GDR}$) should be similar to that in
the Milky Way and other local disks. Throughout the rest of this work, we will
use as our reference point the mid-range of extragalactic determinations,
$\kappa_{\rm H}$= 1500–2200 kg m-2, which are consistent with measurements of
the diffuse ISM in the Milky Way. Our chosen normalisation value, then, is
$\kappa_{\rm H}$N = 1884 kg m-2, which is a good match to current theoretical
dust models (THEMIS: Jones et al., 2017, and the updated Draine et al., 2007
modified by Hensley & Draine, 2021).

For the standard Milky Way value of $\delta_{\rm GDR}$ (= 135, Jones et al.,
2017; Magdis et al., 2012), $\kappa_{\rm H}$ = 1884 kg m-2 implies that
$\kappa_{850}$ = 0.071 m2 kg-1, similar to that used in many extragalactic
studies (Dunne et al., 2000; James et al., 2002; da Cunha et al., 2013). For a
$\delta_{\rm GDR}$ fixed to 135, the range of $\kappa_{\rm H}$ in
extragalactic studies implies a range in $\kappa_{850}$ of 0.06–0.09 m2 kg-1.
Table LABEL:kappalitT lists $\kappa_{\rm H}$ values from extragalactic and
Galactic observations, as well as from theoretical dust models.

Note that the choice of normalisation does not affect any of the trends, nor
the ratio of the conversion parameters in the pairings; it merely sets the
average value of the reference calibration parameter, to which the others are
relative.

The sample mean expectation values for the other two parameters,
$\langle\mbox{$\alpha_{\rm CO}$}\rangle$ and $\langle\mbox{$X_{\rm
C\,I}$}\rangle$, are next derived from our assumed value of
$\langle\mbox{$\kappa_{\rm H}$}\rangle$, together with the mean ratios of the
observables listed in Table 3. The effective standard deviation is also
calculated – the intrinsic scatter of each parameter added in quadrature to
the measurement error for that gas tracer. For example, for CO:

$\sigma_{\rm eff}=\sqrt{s_{\alpha}^{2}+\sigma_{\rm CO}^{2}},$

where $s_{\alpha}$ is the intrinsic scatter in log($\alpha_{\rm CO}$), and
$\sigma_{\rm CO}$ is the measurement error on log($L^{\prime}_{\rm CO}$).

We can now estimate the probability of finding a particular set of conversion
factors for any given galaxy. We use $a_{i},\,i=1,2,3$ to denote the
logarithms of the three conversion factors171717For ease of representation,
$a_{850}=-\log(\mbox{$\alpha_{850}$})$., and write the mean expectation values
and effective standard deviations as $\langle a_{i}\rangle$, and
$\sigma_{i,\rm eff}$ respectively. Assuming that these follow Gaussian
distributions, the probability of finding the factors, $a_{i}$, for any galaxy
is:

$P\propto\displaystyle\prod_{i=1}^{N}\exp\left(-\frac{(a_{i}-\langle
a_{i}\rangle)^{2}}{2\sigma_{\rm i,eff}^{2}}\right)\\\
=\exp\left(-\displaystyle\sum_{i=1}^{N}\frac{(a_{i}-\langle
a_{i}\rangle)^{2}}{2\sigma_{\rm i,eff}^{2}}\right).$ (13)

Thus, the ratios of observable luminosities for any given galaxy can be used
to determine the ratios of conversion factors (Eqn. 12), and the common
scaling factor that maximises the probability in Equation 13 is the best
estimate of $M_{\rm H2}$. The derivation in Appendix F shows that this reduces
analytically to a simple inverse variance weighted mean, such that:

$\log M_{\rm H_{2}}^{\rm opt}=\frac{\sum^{N}_{i=1}(m_{i}\times
w_{i})}{\sum^{N}_{i=1}w_{i}},$ (14)

where $w_{i}=1/\sigma_{i,\rm eff}^{2}$, and $m_{i}$ is the log mass estimate
for each tracer.

$m_{i}=l_{i}+\langle a_{i}\rangle,$

where $l_{i}$ is the measured observable (log luminosity) and $\langle
a_{i}\rangle$ is the sample mean expectation value for the conversion factor.
Once the optimal mass is determined this way, we can then estimate the
corresponding optimal conversion factor on a per-galaxy basis, as:

$a_{i}=m^{\rm opt}-l_{i}.$ (15)

The error on the optimal mass is simply the error on the inverse variance
weighted mean:

$\sigma_{m}^{\rm opt}=\left(\sum^{N}_{i=1}w_{i}\right)^{-1/2},$ (16)

and the error on each of the conversion factors, accounting for co-variance
is:

$\sigma_{ai}=\sqrt{\sigma_{\rm
m^{opt}}^{2}+\sigma_{li}^{2}\left(1-\frac{2w_{i}}{\sum^{N}_{j=1}w_{j}}\right)},$
(17)

where $\sigma_{li}$ is the logarithmic measurement error on the observable
quantity, e.g. $L^{\prime}_{\rm CO}$, $L^{\prime}_{\rm C\,I}$, $L_{850}$.

By design, each tracer for a given galaxy, together with its optimised
conversion factors, will produce the same gas mass, such that $\mbox{$M_{\rm
H2}$}^{\rm CO}=\mbox{$M_{\rm H2}$}^{\rm C\,I}=\mbox{$M_{\rm H2}$}^{\rm dust}$.

##  4 Trends in the luminosity ratios

As our cross-calibration process relies on measurements of the luminosity
ratios, it is first instructive to look at the trends in these observables to
better understand any subsequent trends in the derived conversion factors.

Histograms of the tracer ratios are shown in Fig. 2, split by factors such as
lensing, redshift, SFR, and other notable quantities. The correlations of the
three tracer luminosities are shown in Fig. 3, where the various samples are
colour coded and labelled and in each panel the blue line and shaded region
represent the best fit and $2\sigma$ error interval. Fitting was performed
using our own Orthogonal Distance Regression (ODR) method, which includes $x$
and $y$ errors, intrinsic scatter as a third parameter, and covariance in
errors where required. The method is described in detail in Appendix I. Fit
parameters, slope $m$, intercept $c$, and scatter ln $\lambda$ , are listed in
Table LABEL:fitsT and statistics for the various subsets from Fig. 2 are given
in Table LABEL:hypT. It is instructive to look at these two plots together for
the same luminosity pairs.

Figure 3: Relationships between the observed luminosities. The blue line and
shaded regions in each panel show the best fit and the $2\sigma$ uncertainty
(see Table LABEL:fitsT for parameters). Top left: $L_{850}$–$L^{\prime}_{\rm
CO}$. Cyan diamonds indicate those local galaxies which have $f_{\rm H\,I}$
$>1$ and have corrected values of $L_{850}$. The best fit is shown, including
the H i-corrected set, but excluding the lo-VALES galaxies (peach circles with
log $L^{\prime}_{\rm CO}$ $<8.9$; see Appendix B.1). Choosing only galaxies
with $f_{\rm H\,I}$ $<1$ does not change the results; the fit is consistent
with a linear slope with a high degree of confidence. Top right:
$L_{850}$–$L^{\prime}_{\rm C\,I}$. The resolved C icor galaxies from J19 that
have an aperture correction (Appendix A) applied to $L^{\prime}_{\rm C\,I}$ in
order to compare to $L_{850}$ are shown as pink diamonds. The dotted blue line
shows the fit when these galaxies are included. The solid blue line is the fit
that excludes these galaxies, while the dashed blue line is the fit to the log
$L_{\rm IR}>11$ galaxies only. Bottom left: $L^{\prime}_{\rm
C\,I}$–$L^{\prime}_{\rm CO}$. The pink diamonds are the C icor galaxies, where
now the C i and CO are from matched regions, however, they are still offset to
lower $L^{\prime}_{\rm C\,I}$ for a given $L^{\prime}_{\rm CO}$. The fit
including the pink diamonds is steeper than linear, while excluding them
(solid line) returns a linear slope. Bottom right:
$\sigma^{\prime}_{850}$–$L^{\prime}_{\rm C\,I}$. We compare the 850-$\mu$m
emissivity ($\mbox{$\sigma^{\prime}_{850}$}=\mbox{$M_{\rm
dust}$}\times\mbox{$\kappa_{850}$}$) to $L^{\prime}_{\rm C\,I}$. This allows
us to show dust data in the resolved galaxies from Jiao et al. (2021) (pink
stars) measured in the same region as [C i](1–0). The dotted blue line shows
the fit including the pink stars, the solid line shows the fit excluding the C
icor galaxies and the dashed line shows the fit to log $L_{\rm IR}>11$
galaxies only. The local resolved (C icor) galaxies from J19/J21 appear to
have less C i emission for a given dust or CO luminosity, an effect which
persists even when comparing matched regions (bottom row).

###  4.1 $L_{850}$ vs. $L^{\prime}_{\rm CO}$

The left-hand column of Fig. 2 and Fig. 3(a) show the quantities $L_{850}$ vs.
$L^{\prime}_{\rm CO}$. There are no significant differences in the
distribution of $L_{850}$/$L^{\prime}_{\rm CO}$ with lensing, redshift or SF-
mode, but the observed uncorrected ratio is significantly higher for galaxies
with log $L^{\prime}_{\rm CO}$ $<8.9$ (Fig. 2; left green histogram). These
log $L^{\prime}_{\rm CO}$ $<8.9$ galaxies tend to have optical disks dominated
by atomic hydrogen ($f_{\rm H\,I}$ $>1$) and the likely increased contribution
to $L_{850}$ from dust associated with H i rather than H2 results in the
offset to higher $L_{850}$/$L^{\prime}_{\rm CO}$ ratios. For local galaxies
with $f_{\rm H\,I}$ $>1$ we apply a correction (Appendix B.1) which appears to
remove this offset (cyan diamonds in Fig. 3(a)). Galaxies with log
$L^{\prime}_{\rm CO}$ $<8.9$ from the VALES sample at $0.04<z<0.3$ (Villanueva
et al., 2017; Hughes et al., 2017) have noticeably higher
$L_{850}$/$L^{\prime}_{\rm CO}$ ratios than VALES galaxies with log
$L^{\prime}_{\rm CO}$ $>8.9$ (peach circles in Fig. 3(a) and Table
LABEL:hypT). There are no published H i measurements for VALES, but we suspect
that the lo-VALES sample with log $L^{\prime}_{\rm CO}$ $<8.9$ is likely to be
H i-rich, based on the similarity between the log $L^{\prime}_{\rm CO}$ $<8.9$
and $f_{\rm H\,I}$ $>1$ categories in the third panel. We therefore exclude
lo-VALES from our averages, as we suspect they are in need of a correction for
H i but we have no means to apply one. We also recommend that a low
$L^{\prime}_{\rm CO}$ requires careful consideration of H i-related dust. In
Fig. 3(a), the tracers show a linear dependence, regardless of exactly which
galaxies are included (Table LABEL:fitsT). The lo-VALES galaxies are excluded
from all the fits.

###  4.2 $L_{850}$ vs. $L^{\prime}_{\rm C\,I}$

The central column of Fig. 2 and Fig. 3(b) show $L_{850}$ vs. $L^{\prime}_{\rm
C\,I}$, the pair with the least scatter ($\ln\lambda$ in Table LABEL:fitsT).
There are no significant differences in the distributions of
$L_{850}$/$L^{\prime}_{\rm C\,I}$ as a function of lensing, redshift or SF-
mode. The $L_{850}$–$L^{\prime}_{\rm C\,I}$ sample has 12 galaxies with large
angular sizes from J19 that have only part of their disk detected in C i,
denoted C icor. To compare $L^{\prime}_{\rm C\,I}$ and $L_{850}$ for these
galaxies, we need to apply an aperture correction to their C i fluxes, thereby
assuming that their $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ ratios
remain roughly constant across the disk (see Appendix 4.4). Some of the
correction factors are very large (up to 0.7 dex); even after correction,
their $L_{850}$/$L^{\prime}_{\rm C\,I}$ ratios are significantly offset from
the rest of the MS sample (green histogram). Fig. 3(b) shows $L_{850}$ vs.
$L^{\prime}_{\rm C\,I}$, with the C icor galaxies as pink diamonds. The fit to
all galaxies including the C icor subset is shown as a dotted blue line, and
has a sub-linear slope $m$ ($3\sigma$), although the difference is very small
($m=0.937\pm 0.020$). Excluding the C icor galaxies from the fit leaves a
linear relationship, shown by the blue solid line. At first glance, this and
the green histogram in Fig. 2 suggest that our CO-based corrections are
insufficient; however, Fig. 3(d) shows that the problem is not simply the
assumption used to correct $L^{\prime}_{\rm C\,I}$ to match the global
$L_{850}$. This plot shows the dust emissivity cross-section,
$\sigma^{\prime}_{850}$, equivalent to $L_{850}$ but with the temperature
sensitivity removed181818$\sigma^{\prime}_{850}$ is derived from the data
provided in Jiao et al. (2021) by multiplying the dust mass in the C i
aperture by the $\kappa_{850}$ used in their method, $\kappa_{850}$ = 0.034 m2
kg-1 (Draine, 2003; Draine et al., 2007).. Importantly, this quantity is
measured in the same aperture as $L^{\prime}_{\rm C\,I}$. The C icor galaxies
are shown as pink stars, and they – along with many other low-luminosity
galaxies in J21 – still appear to have less C i emission for a given amount of
dust.

Table 4: Parameters of robust ODR fits between variables using MCMC, co-
variant errors and including intrinsic scatter, $\ln\lambda$.

Log $x$ Log $y$ Group $m$ $c$ ln $\lambda$ $r_{\rm s}$ $p$ $N$
$L^{\prime}_{\rm CO}$ $L^{\prime}_{\rm C\,I}$ BL 1.023 (0.029) $-0.93$ (0.29)
$-2.09$ (0.09) 0.94 109 $L^{\prime}_{\rm CO}$ $L^{\prime}_{\rm C\,I}$ BL+C
icor 1.078 (0.026) $-1.49$ (0.25) $-2.03$ (0.09) 0.96 121 $L^{\prime}_{\rm
CO}$ $L^{\prime}_{\rm C\,I}$ log $L_{\rm IR}>11$ 0.950 (0.035) $-0.18$ (0.37)
$-2.07$ (0.10) 0.92 97 $L^{\prime}_{\rm C\,I}$ $L_{850}$ BL 0.976 (0.024)
14.33 (0.23) $-2.11$ (0.09) 0.95 140 $L^{\prime}_{\rm C\,I}$ $L_{850}$ BL+C
icor 0.937 (0.020) 14.71 (0.19) $-2.08$ (0.08) 0.96 152 $L^{\prime}_{\rm
C\,I}$ $L_{850}$ log $L_{\rm IR}>11$ 1.024 (0.030) 13.86 (0.30) $-2.17$ (0.09)
0.94 128 $L^{\prime}_{\rm C\,I}$ $\sigma^{\prime}_{850}$ BL 1.024 (0.025)
$-2.28$ (0.24) $-2.05$ (0.08) 0.95 140 $L^{\prime}_{\rm C\,I}$
$\sigma^{\prime}_{850}$ BL+C icor 0.997 (0.021) $-2.01$ (0.20) $-2.07$ (0.08)
0.96 152 $L^{\prime}_{\rm CO}$ $L_{850}$ BL 1.003 (0.015) 13.42 (0.15) $-2.01$
(0.05) 0.96 326 $L^{\prime}_{\rm CO}$ $L_{850}$ $f_{\rm H\,I}$<1 1.002 (0.017)
13.43 (0.16) $-1.99$ (0.05) 0.96 310 $L^{\prime}_{\rm CO}$ $L_{850}$ log
$L_{\rm IR}>11$ 0.983 (0.026) 13.63 (0.26) $-1.91$ (0.06) 0.93 226 $L_{\rm
IR}$ $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ BL+C icor 0.071 (0.02)
$-1.57$ (0.23) $-1.70$ (0.09) 0.26 0.005 121 $L_{\rm IR}$ $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ BL 0.034 (0.02) $-1.11$ (0.25) $-1.70$ (0.09)
0.09 0.34 109 $T_{\rm d}$ $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ BL+C
icor 1.23 (0.23) $-2.60$ (0.35) $-$2.20 (0.12) 0.31 115 $T_{\rm d}$
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ BL 0.83 (0.28) $-2.0$ (0.4)
$-$2.00 (0.15) 0.15 0.12 103 $L_{\rm IR}$/$L_{850}$ $T_{\rm mw}$ 0.216 (0.010)
3.90 (0.12) $-3.76$ (0.15) 0.82 152 $L_{\rm IR}$ $T_{\rm mw}$ 0.070 (0.004)
0.60 (0.05) $-3.10$ (0.10) 0.80 152 $L_{\rm IR}$ $\alpha_{850}$ BL 0.045
(0.007) 12.271 (0.082) $-3.36$ (0.30) 0.46 230 $L_{\rm IR}$ $\alpha_{\rm CO}$
BL 0.59 (0.09) $-0.91$ (1.10) $-1.00$ (0.50) 0.46 230 $L_{\rm IR}$
$\alpha_{\rm C\,I}$ BL $-0.052$ (0.010) 1.896 (0.124) $-3.90$ (0.20) $-0.48$
82 $L_{\rm IR}$ $X_{\rm C\,I}$ BL $0.028$ (0.011) $-5.136$ (0.133) $-3.70$
(0.23) 0.29 0.008 82

Table 5: Two-sample KS-test result and Z-test statistic for the following
parameter pairs shown in Figs 2 and 4.

Quantity $A$ $N_{\rm A}$ $B$ $N_{\rm B}$ $\bar{A}$ $\bar{B}$ $Z(\sigma)$
$P_{\rm KS}$ $L_{850}$/$L^{\prime}_{\rm CO}$† log $L^{\prime}_{\rm CO}$<8.9 50
$L^{\prime}_{\rm CO}$>8.9 (MS) 138 $13.600\pm 0.030$ $13.420\pm 0.016$ 5.3
3e-5 $L_{850}$/$L^{\prime}_{\rm CO}$† $f_{\rm H\,I}$<1 (MS) 168 $f_{\rm
H\,I}$>1 (MS) 24 $13.456\pm 0.014$ $13.685\pm 0.042$ 5.2 3e-6
$L_{850}$/$L^{\prime}_{\rm CO}$† lo-VALES 7 $f_{\rm H\,I}$>1 17 $13.741\pm
0.061$ $13.663\pm 0.053$ 0.57 $L_{850}$/$L^{\prime}_{\rm CO}$ $f_{\rm H\,I}$<1
(MS) 168 $f_{\rm H\,I}$>1 (MS∗) 17 $13.456\pm 0.014$ $13.448\pm 0.039$ 1.0
$L_{850}$/$L^{\prime}_{\rm C\,I}$ MS 61 C icor 12 $14.108\pm 0.027$ $14.36\pm
0.036$ 5.6 6e-4 $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ $z<3$ (SMGs) 37
$z>3$ 11 $-0.686\pm 0.031$ $-0.512\pm 0.044$ 3.2 0.012 $L^{\prime}_{\rm
CO}$/$L^{\prime}_{\rm C\,I}$ MS 55 SMGs 54 $-0.743\pm 0.028$ $-0.651\pm 0.028$
2.3 0.056 $L^{\prime}_{\rm CO}$/$L^{\prime}_{\rm C\,I}$ MS+C icor 66 SMGs 54
$-0.786\pm 0.026$ $-0.651\pm 0.028$ 3.5 0.006 $T_{\rm d}$ MS 174 SMGs 160
$31.1\pm 0.4$ $38.3\pm 0.7$ 8.8 2e-12 $T_{\rm mw}$ MS 82 SMGs 52 $23.0\pm 0.4$
$30.1\pm 0.7$ 8.8 7e-13

###  4.3 $L^{\prime}_{\rm C\,I}$ vs. $L^{\prime}_{\rm CO}$

The right-hand column of Fig. 2 and Fig. 3(c) show $L^{\prime}_{\rm C\,I}$ vs.
$L^{\prime}_{\rm CO}$. There are no significant differences in
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ for strongly lensed vs. unlensed
sources, but the highest redshift, $z>3$, galaxies have higher
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ ratios at marginal significance
($p=0.012$). There are only 11 galaxies at $z>3$ and a larger sample is needed
to determine if this is a genuinely significant trend. Fig. 3(c) shows the C
icor galaxies as pink diamonds, where the C i and CO fluxes are measured in
the same apertures by J19. The solid line shows the fit to all galaxies, which
is non-linear at 3$\sigma$ significance ($m=1.078\pm 0.026$). The slope
becomes linear once the C icor galaxies are removed. The green histogram in
Fig. 2 (right) shows more clearly why we see this: the C icor galaxies have
significantly lower $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ ratios
compared to other galaxies at similar or higher luminosity. The bottom
histogram shows a marginal difference between galaxies with different SFRs
($p=0.056$) when excluding the C icor galaxies, which becomes significant when
they are included ($p=0.006$).

###  4.4 Resolved C i fluxes from Herschel FTS mapping

The local resolved galaxies observed with Herschel FTS (Jiao et al., 2019) lie
off the global trends seen in Fig. 3. There are possible physical explanations
why lower luminosity, and more quiescently star-forming galaxies might have
lower C i/CO line ratios (for example, different ISM environments in terms of
their position in the CR energy density vs average molecular gas density
diagram: see Figure 1 in Bisbas et al., 2015). Low ratios of $L^{\prime}_{\rm
C\,I}$/$L^{\prime}_{\rm CO}$ have been found in other studies, most
intriguingly in the case of the interacting LIRG NGC 6052 using ALMA
(Michiyama et al., 2020), and some high-$z$ strongly lensed sources
(Harrington et al., 2021). Such ratios tend to be unusual in higher luminosity
samples, however, whereas the resolved FTS sample has a very low average for
the C i line ratios with both CO and dust.

We recommend caution in the interpretation of the data for these resolved FTS
sources because another team subsequently presented the same data but drew
different conclusions (Crocker et al., 2019). We can therefore only note that
the C i fluxes from Herschel FTS mapping datasets are not always consistent
when analysed by different teams.191919Crocker et al. did not provide
integrated fluxes, nor a method to determine them from their published
measurements; hence, we cannot use their work directly in our analysis. Q.
Jiao has provided us with the maps used in J19, enabling us to check the
measurements independently and extend our analysis, but we have had no
responses to our requests for integrated fluxes or the details of the method
used from the Crocker team. As these resolved FTS measurements are essentially
the only source of C i data at log $L_{\rm IR}<11$, and carry a lot of weight
in $L_{\rm IR}$ and SFR correlations, we chose not to include the C icor
galaxies in the statistical analysis. If, instead, we take the J19
measurements at face value – they signpost a fundamental physical change in C
i properties, a finding which clearly warrants further study with ground-based
facilities. We discuss possible physical mechanisms for changes in the
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$and $L^{\prime}_{\rm
C\,I}$/$L_{850}$ ratios in Appendix B.2.

###  4.5 Trends with global indicators of star-formation.

Figure 4: Histograms of (top) luminosity-weighted dust temperature ($T_{\rm
d}$) from MBB fits (bottom) mass-weighted dust temperature ($T_{\rm mw}$) from
fits allowing for multiple dust temperatures. Means and KS test results are
given in Table LABEL:hypT.

Figure 5: $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ luminosity ratio
(equivalent to $\alpha_{\rm CO}$/$\alpha_{\rm C\,I}$) as a function of
luminosity-weighted (peak SED) dust temperature (left) and log $L_{\rm IR}$
(right). The yellow shaded regions in the right panel represent the 1$\sigma$
range of values found by Crocker et al. (2019) for resolved galaxies mapped
with the Herschel FTS. The vertical red and blue bars show the 1$\sigma$ range
of values for QSOs/SMGs and the Milky Way, respectively (Walter et al., 2011;
Frerking et al., 1989). The best fit to all galaxies and 2$\sigma$ confidence
region are shown as the blue line and shaded area (parameters are given in
Table LABEL:fitsT). The $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ ratio
has a significant correlation with both $T_{\rm d}$ and log $L_{\rm IR}$ when
all galaxies are included, but when excluding the C icor galaxies, the
correlation all but disappears (blue dotted line).

Finally, we check to see if any of the tracer ratios are sensitive to SFR
indicators. In our data-set the dust observables $L_{\rm IR}$ and $T_{\rm d}$
are indicators of the intensity and magnitude of star formation in galaxies
(e.g. Kennicutt, 1998; Foyle et al., 2012; Liu et al., 2021). As expected, the
distribution of $T_{\rm d}$ is very different for the MS galaxies and SMGs
(Fig. 4 and Table LABEL:hypT), reflecting the increase in the intensity of
star formation in the SMGs. The only tracer ratio sensitive to these SF
indicators is $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$, which in Figure 5
is seen to increase with $L_{\rm IR}$ and $T_{\rm d}$ when all galaxies are
considered. Such a trend was not reported for smaller samples over a more
limited range of luminosity (e.g. Jiao et al., 2017; Jiao et al., 2019),
presumably because of limited statistics. However, if the C icor galaxies
(pink diamonds) are excluded, the correlation all but disappears (blue dotted
line).

Naively, we might expect $L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ to rise
with increasing SFR intensity, due to the expected destruction of CO by cosmic
rays (CR) in high-SFR environments (e.g. Bisbas et al., 2015). For a given
range of H2 densities in the typically hierarchical molecular clouds, any
increased CR-induced ionisation rate, $\zeta_{\rm CR}$ (due to a rising
average CR energy density, $U_{\rm CR}$) will destroy CO in the lower-density,
more extended areas, while leaving CO still tracing H2 in the more compact,
denser regions (see also Figure 1 of Papadopoulos et al., 2018 for a visual
effect of this). Intriguingly, the gas density, $n$(H2), and CR ionization
rate, $\zeta_{\rm CR}$, will compete against each other in ULIRG/SMG
environments, with the higher <n> expected in their highly turbulent ISM
tending to keep the ordinary CO/C i chemistry in place, even when exposed to
the higher $\zeta_{\rm CR}$ values.202020We here assume that CR energy density
$U_{\rm CR}\propto\rm\rho_{SFR}$ and CR ionisation rate $\mbox{$\zeta_{\rm
CR}$}\propto U_{\rm CR}$. Guessing which one will win this highly non-linear
competition (see Fig 1, 8 in Bisbas et al., 2015) is dangerous in the absence
of CO and C i line data. These effects have been probed with a variety of
simulations (e.g. Bisbas et al., 2015, 2021; Clark et al., 2019a; Gong et al.,
2020) and while showing similar trends, they are not easily parameterisable in
terms of $n$(H2) and $\zeta_{\rm CR}$; one reason why such cross-calibration
efforts of the available gas mass tracers are so important.212121On an
individual galaxy basis one could assemble well-sampled CO, 13CO and C i line
SLEDs and overcome these problems with detailed analysis (e.g. Papadopoulos et
al., 2014). However, even in the ALMA era this remains very expensive in terms
of telescope time making it prohibitive for large samples of galaxies.

The other two tracer ratios show no trends with either $L_{\rm IR}$, $T_{\rm
d}$ (our proxies for SFR) – we present the relevant plots in Appendix 14 for
completeness.

##  5 Results

In this section we present the results of the optimisation method, firstly for
the daX sample, for which we have all three gas tracers, and later for the
other three samples, for which we have pairs of tracers. We investigate trends
of the conversion factors with $L_{\rm IR}$ and SFR. Mean values for the
conversion factors are listed in Table LABEL:caloptT.

Fig. 6 shows the results for the daX sample (Figs 7–9 present the same results
for each of the samples in turn). The top row of each plot shows the
distribution of the relevant physical parameter for C i and dust, and the
conversion factor for CO: $X_{\rm C\,I}$, $\kappa_{\rm H}$ and $\alpha_{\rm
CO}$. The lower left panels show the same quantities for the individual
galaxies as a function of $L_{\rm IR}$; each panel indicates a reference
measure to give context. For $X_{\rm C\,I}$, the horizontal lines indicate the
measured extremes found in the local Universe: Orion A/B clouds in the Milky
Way (Ikeda et al., 2002) and the starburst centre of M 82 (White et al.,
1994), while the grey shaded region shows the range of values inferred from
observations of GRB hosts and QSO absorbers for solar metallicity by Heintz &
Watson (2020). They use a method which does not rely on emission measures of
dust, CO or C i and so can be considered independent. For $\alpha_{\rm CO}$,
the horizontal lines indicate the typical $\alpha_{\rm CO}$ for the Milky Way
(Bolatto et al., 2013) and that commonly adopted222222In this panel,
$\alpha_{\rm CO}$ does not include the factor 1.36 for He. for ULIRGs and SMGs
(Downes & Solomon, 1998). For $\kappa_{\rm H}$, we show a shaded band
indicating the range derived for local galaxies (see Table LABEL:kappalitT),
along with lines showing the value for the most diffuse and dense sight lines
in the Milky Way (Remy et al., 2017). The right lower panels show the running
log-means as a function of $L_{\rm IR}$ to make it easier to see any trends,
and additionally includes232323As elsewhere, the empirical parameters
$\alpha_{\rm C\,I}$ and $\alpha_{850}$ include the factor 1.36 for He. the
empirical parameters, $\alpha_{\rm C\,I}$ and $\alpha_{850}$. The solid shaded
bins are the means for the grey points, which are those used to determine the
calibration; the yellow points are C icor and the semi-transparent pentagon is
the mean of those – see §4.4 for more details.

Fig. 6 shows that for galaxies with log $L_{\rm IR}>11$ there are only weak
trends of the conversion factors with $L_{\rm IR}$. While the normalisation
($\mbox{$\kappa_{\rm H}$}^{\rm N}=1884$ kg m-2) was chosen to produce average
dust properties consistent with the Milky Way and other nearby spirals, the CO
and C i conversion factors derived from the luminosity ratios also lie within
the ranges expected from independent studies. The averages at log $L_{\rm
IR}<11$ are based on only a small number of points (12) and more C i studies
are required to probe quiescent local galaxies.

Figure 6: Top row: Distributions of the three conversion factors, $X_{\rm
C\,I}$, $\alpha_{\rm CO}$ and $\kappa_{\rm H}$ for the daX sample. The overall
normalisation has been set as $\mbox{$\kappa_{\rm H}$}^{\rm N}=1884$ kg m-2,
which produces the $\mbox{$\delta_{\rm GDR}$}=135$ of the Milky Way for
$\mbox{$\kappa_{850}$}=0.071$ m2 kg-1(Jones, 2018). The intrinsic pairwise
scatter and conversion factors were estimated excluding C icor and lo-VALES
galaxies. Lower left: The optimised conversion factors as a function of
$L_{\rm IR}$. The grey shaded band represents $X_{\rm C\,I}$ at solar
metallicity using the relationship found in GRB and QSO absorber hosts by
Heintz & Watson (2020), a measure which is independent of assumptions about
$\alpha_{\rm CO}$ or $\delta_{\rm GDR}$. Lower right: The running means of the
conversion factors, with error bars of $\sigma/\sqrt{N_{\rm bin}}$. The solid
shaded bins are the means for the grey points which are those used to
determine the calibration, the yellow points are C icor and the semi-
transparent pentagon is the mean of those – see §4.4 for more details. Table
6: Mean optimised conversion factors for our various samples.

Sample Selection $N$ $X_{\rm C\,I}$ $\sigma_{\bar{X}}$ $\alpha_{\rm CO}$
$\sigma_{\bar{\alpha}}$ $\kappa_{\rm H}$ $\sigma_{\bar{\kappa}}$ $\delta_{\rm
GDR}$ /$\times 10^{-5}$ M⊙ (k km s-1 pc2)-1 m2 kg-1 daX log $L_{\rm IR}>11$ 90
$1.59^{+0.45}_{-0.38}$ 0.04 $2.66^{+0.96}_{-0.70}$ 0.10 $1990^{+738}_{-607}$
86 141 daX log $L_{\rm IR}<11$ 12 $1.18^{+0.60}_{-0.29}$ 0.13
$2.44^{+0.56}_{-0.51}$ 0.17 $1571^{+732}_{-525}$ 163 112 Xd log $L_{\rm
IR}>11$ 128 $1.59^{+0.47}_{-0.35}$ 0.04 $1946^{+654}_{-464}$ 53 138 Xd log
$L_{\rm IR}<11$ 12 $1.24^{+0.57}_{-0.27}$ 0.13 $1503^{+604}_{-369}$ 145 107 ad
log $L_{\rm IR}>11$ 240 $3.08^{+1.32}_{-0.81}$ 0.07 $1936^{+658}_{-504}$ 45
137 ad log $L_{\rm IR}<11$ 88 $3.52^{+0.95}_{-0.84}$ 0.10 $1718^{+502}_{-339}$
44 122 Xa log $L_{\rm IR}>11$ 97 $1.61^{+0.39}_{-0.31}$ 0.04
$2.57^{+0.71}_{-0.62}$ 0.08 Xa log $L_{\rm IR}<11$$+$C icor 24
$1.30^{+0.2}_{-0.23}$ 0.05 $1.88^{+0.41}_{-0.34}$ 0.10 Xa log $L_{\rm IR}<11$
12 $1.37^{+0.34}_{-0.31}$ 0.09 $2.11^{+0.18}_{-0.57}$ 0.15

###  5.1 A calibration for the gas masses

We next give a prescription for estimating gas mass, tailored to how many
tracers are available and – where appropriate – the type of galaxy being
investigated.

####  5.1.1 Dual-band

While the information content is greatest for the daX sample, which has three
tracer pairs to optimise, the method presented in §3 still improves the cross-
calibration for samples which have two tracer measurements, i.e. one pair. The
results for the Xd, Xa and ad samples are shown in Figs 7–9 and behave
similarly to the daX sample, as one would hope given that the daX galaxies are
a subset of the others. The pink diamonds in the lower-left panels in Figs 7
and 8 denote the C icor galaxies. The cyan diamonds in the lower-left panel of
Fig. 9 are galaxies with $f_{\rm H\,I}$>1 which have been corrected for the
contribution of dust mixed with the H i gas, as described in Appendix A. The
open peach circles are the lo-VALES galaxies, which we suspect to have $f_{\rm
H\,I}$>1 (see §4.1) but which we cannot correct. We do not include these in
any averages or histograms.

The method previously used in the literature (e.g. Alaghband-Zadeh et al.,
2013; Scoville et al., 2016; Orellana et al., 2017; Hughes et al., 2017;
Valentino et al., 2018) has been to assume one tracer in a pair (e.g.
$L^{\prime}_{\rm CO}$) has a known conversion ($\alpha_{\rm CO}$), then to fix
that factor for all galaxies in order to estimate the second (i.e. the one of
interest). We show this simple method alongside our optimal method as grey
lines and dashed grey error bars in the relevant panels of Figs 7–9. The
scatter in the conversion factors for the optimised estimates are governed by
the intrinsic scatter we inferred in our analysis of the data in Appendix E,
free of assumptions. In contrast, the simple method proscribes that there is
no scatter in the known conversion factor, and therefore all of the intrinsic
scatter in the luminosity ratio is attributed to the second conversion factor
of interest. The optimised method presented here does not assume an ad hoc
preference for any particular conversion factor: as it is based on empirical
variance analysis, it uses more of the available information to improve the
accuracy of the estimated conversion factors. The histograms in Figs (7–9)
show that the scatter in the factor of interest is larger when using the
simple method, and the trends in the running medians are also more
exaggerated.

Comparing the parameter estimates using three tracers to those using two
tracers for the same galaxies allows us to test the accuracy of these two-
tracer estimates. The details are in Appendix H, but in summary there is a
reasonable correlation between the three-tracer and two-tracer estimates,
without bias (Fig. 15) and an average scatter of 0.06–0.08 dex.

Figure 7: Top: Distributions of $X_{\rm C\,I}$ and $\kappa_{\rm H}$ for the Xd
sample, where the normalisation used is $\mbox{$\kappa_{\rm H}$}^{\rm N}=1884$
kg m-2. Left: Optimised conversion factors as a function of $L_{\rm IR}$ for
the galaxies in the Xd sample. The grey shaded bands represent the range of
$X_{\rm C\,I}$ ($Z=Z_{\odot}$) from a study of GRB and QSO absorber hosts by
Heintz & Watson (2020), a measure which is independent of assumptions about
$\alpha_{\rm CO}$ or $\delta_{\rm GDR}$, and also the range of $\kappa_{\rm
H}$ in nearby galaxies. Right: Running means of the conversion factors with
error bars of $\sigma/\sqrt{N_{\rm bin}}$. The solid coloured pentagons are
the means for the grey points, which are those used to determine the
conversion; the yellow points are the C icor galaxies and the semi-transparent
pentagon is the mean of those (see §4.4 and Appendix B.2). The grey dashed
error bars in the $X_{\rm C\,I}$ running mean are for the so-called simple
method, where $\mbox{$\kappa_{\rm H}$}=1884$ kg m-2 is fixed for all galaxies.

Figure 8: Top: Distributions of $X_{\rm C\,I}$ and $\alpha_{\rm CO}$ for the
Xa sample, where the normalisation has been set to $\mbox{$X_{\rm C\,I}$}^{\rm
N}=1.6\times 10^{-5}$ – the middle of the range for $Z_{\odot}$ found in GRB
hosts and QSO absorbers by Heintz & Watson (2020). Lower left: Optimised
conversion factors as a function of $L_{\rm IR}$ for the galaxies in the Xa
sample. Right: Running means of the conversion factors with error bars of
$\sigma/\sqrt{N_{\rm bin}}$. The solid coloured pentagons are the means for
the grey points which are those used to determine the conversion; the yellow
points represent the C icor galaxies and the semi-transparent pentagon is the
mean of those – see §4.4 and Appendix B.2. The grey dashed error bars in the
$\alpha_{\rm CO}$ running mean are for the so-called simple method, where
$\mbox{$X_{\rm C\,I}$}=1.6\times 10^{-5}$ is fixed for all galaxies.

Figure 9: Top: Distributions of $\alpha_{\rm CO}$ and $\delta_{\rm GDR}$ for
the ad sample, where the normalisation has been set to $\mbox{$\kappa_{\rm
H}$}^{\rm N}=1884$ kg m-2. Lower left: Optimised conversion factors as a
function of $L_{\rm IR}$ for the galaxies in the ad sample. The cyan diamonds
indicate galaxies with $f_{\rm H\,I}$>1 which have been corrected for the
contribution of dust mixed with H i, following the procedure outlined in
Appendix B.1. The open peach circles are the lo-VALES galaxies which we
suspect have $f_{\rm H\,I}$>1 – these are not included in the analysis. Right:
Running means of the conversion factors with error bars of
$\sigma/\sqrt{N_{\rm bin}}$. The coloured points reflect the optimised results
and the grey dashed error bars for the $\alpha_{\rm CO}$ running mean are for
the so-called simple method, where $\mbox{$\kappa_{\rm H}$}=1884$ kg m-2 is
fixed for all galaxies.

Thus, when multiple H2 tracers are available, we recommend the procedure
outlined in the example below.

##### Example:

Take the example of a galaxy with observations of both dust and CO. We take
the mean value for the appropriate pair combination from Table 3: $\alpha_{\rm
CO}$/$\kappa_{\rm H}$=0.00133. For our adopted sample mean normalisation of
$\kappa^{\rm N}=1884$ kg m-2, we now infer the sample mean expectation value
of $\langle\alpha\rangle=0.00133\times 1884=2.5$ (excluding He). The sample
means, $\kappa^{\rm N}$ (assumed) and $\langle\alpha\rangle$ (derived), are
next used to estimate an initial gas mass for our galaxy in each of the two
tracers, $L_{850}$ and $L^{\prime}_{\rm CO}$.

$M_{\kappa}=\kappa^{\rm N}\mbox{$L_{850}$}/4\pi B(\nu_{850},\mbox{$T_{\rm
mw}$})\\\ M_{\alpha}=\langle\alpha\rangle\mbox{$L^{\prime}_{\rm CO}$}\\\ $
(18)

Next, we calculate the effective standard deviation by adding the
observational error on the tracer luminosities in quadrature to the intrinsic
scatter for $\alpha$ and $\kappa$.

$\sigma^{\kappa}_{\rm eff}=\sqrt{s_{\kappa}^{2}+\sigma_{850}^{2}}\\\
\sigma^{\alpha}_{\rm eff}=\sqrt{s_{\alpha}^{2}+\sigma_{\rm CO}^{2}}\\\ $ (19)

where $\sigma_{850}$, $\sigma_{\rm CO}$ are the errors on
$\log(\mbox{$L_{850}$})$ and $\log(\mbox{$L^{\prime}_{\rm CO}$})$, and
$s_{\kappa}$ and $s_{\alpha}$ are the intrinsic scatter on
$\log(\mbox{$\kappa_{\rm H}$})$ and $\log(\mbox{$\alpha_{\rm CO}$})$ listed in
Table 3. The optimal H2 mass estimate is then calculated thus:

$M^{\rm
opt}=\frac{M_{\kappa}/\sigma^{2}_{\rm\kappa,eff}+M_{\alpha}/\sigma^{2}_{\rm\alpha,eff}}{1/\sigma^{2}_{\rm\kappa,eff}+1/\sigma^{2}_{\rm\alpha,eff}}$
(20)

We now work back to find the optimal conversion parameters for this galaxy:

$\mbox{$\kappa_{\rm H}$}^{\rm opt}=\kappa^{\rm N}M^{\rm opt}/M_{\kappa}\\\
\mbox{$\alpha_{\rm CO}$}^{\rm opt}=\langle\alpha\rangle M^{\rm
opt}/M_{\alpha}.\\\ $ (21)

Figure 10: Correlations of the empirical calibration factors with their tracer
luminosity. The largest sample was used in each case: ad (top and bottom); Xd
(top centre). For $\alpha_{850}$ (top), we include the relationship from
Scoville et al. (2016), exactly as quoted in that paper. Parameters for the
median posterior fit (dark blue line), accounting for co-variance in the
errors, are quoted in the legend, and we include 100 random fits from sampling
the posterior as pale blue lines to show the scatter. Further details are
given in Table LABEL:calfitsT. All tracers have a significant correlation with
their tracer luminosity.

[]figs/python/cal_empirical_daX(0.39,0.39)

[]figs/python/cal_empirical_Xa(0.33,0.33)
[]figs/python/cal_empirical_Xd_kh_fullMS(0.33,0.33)
[]figs/python/cal_empirical_aod_kh(0.33,0.33)

Figure 11:  The empirical conversion factors, $\alpha_{\rm CO}$, $\alpha_{\rm
C\,I}$ and $\alpha_{850}$, separately for MS galaxies and SMGs. The blue line
histograms represent the MS galaxies. The dark blue dashed lines show the
$f_{\rm H\,I}$>1 and C icor galaxies. We keep these separate because the
corrections are uncertain and the C icor galaxies are significantly different
in their distributions to the rest of the MS galaxies, as discussed in §4.4.
The SMG histograms are shown in red. The result of two-sample KS tests are
shown on each panel for each of the comparisons (MS galaxies vs. SMGs). The
main comparison between MS galaxies and SMGs is shown first; below, in
parentheses, are the results for comparison between the robust MS galaxies and
those requiring the corrections described above. daX: Only $\alpha_{\rm C\,I}$
shows a significant ($3.2\sigma$) difference between the two groups,
decreasing to $1.9\sigma$ if constant $T_{\rm mw}$=25 k is used. Xa: There is
a marginal difference in $\alpha_{\rm C\,I}$ (2.5–$2.9\sigma$) but no
difference in $\alpha_{\rm CO}$. There is a significant difference between the
C icor and other MS galaxies (dotted blue line and parentheses), which cannot
be attributed to uncertain aperture corrections as CO and C i are measured in
the same aperture. Xd: For $\alpha_{\rm C\,I}$ (and $\kappa_{\rm H}$ – see
Table 8) there is a marginal difference (2.6–$3.0\sigma$) between the MS
galaxies and SMGs, where these differences vanish when using fixed $T_{\rm
mw}$=25 k. There is no difference between MS galaxies and SMGs for
$\alpha_{850}$. ad: There is a small ($\sim 10$ per cent) difference in the
mean $\alpha_{\rm CO}$ between the MS galaxies and SMGs, significant at the
$3.5\sigma$ level. This drops to $2.5\sigma$ when comparing only log $L_{\rm
IR}>11$ MS galaxies, and drops further when using constant $T_{\rm mw}$=25 k.
We conclude that this is a luminosity-driven effect, not a property of the
star-formation mode per se. The difference in $\alpha_{850}$ is highly
significant (as is the diference in $\kappa_{\rm H}$, to a lesser extent – see
Table 8). Fig. 10 shows that this difference is the result of a steady change
in $\alpha_{850}$ with $L_{\rm IR}$ (or $L_{850}$), rather than being a bi-
modal split between MS galaxies and SMGs.

####  5.1.2 Single band: empirical conversion factors

Table 7: Empirical conversion factors recommended for use when only a single
gas tracer observation is in hand. These include a factor 1.36 to account for
He.

Sample $\alpha_{850}$($\times 10^{12}$) $\alpha_{\rm CO}$ $\alpha_{\rm C\,I}$
W Hz-1 $\rm M_{\odot}$-1 M⊙ (k km s-1 pc2)-1 log $L_{\rm IR}<11$ $5.8\pm 0.1$
$4.7\pm 0.1$ .. log $L_{\rm IR}>11$ $6.9\pm 0.1$ $4.0\pm 0.1$ $17.0\pm 0.3$ MS
$6.2\pm 0.1$ $4.4\pm 0.1$ $19.1\pm 0.6$ SMGs $7.3\pm 0.1$ $3.8\pm 0.1$
$16.2\pm 0.4$

Table 8: Log means and statistical tests for conversion factors for MS
galaxies and SMGs.

C i Sample $\mbox{$X_{\rm C\,I}$}{\rm(MS)}$ $\mbox{$X_{\rm C\,I}$}{\rm(SMG)}$
$\mbox{$\alpha_{\rm C\,I}$}{\rm(MS)}$ $\mbox{$\alpha_{\rm C\,I}$}{\rm(SMG)}$
$Z(\sigma)$ $P_{\rm KS}$ $P_{\rm KS}$ (25 k) daX $1.41\pm 0.07$ $1.70\pm 0.06$
$19.1\pm 0.9$ $15.8\pm 0.5$ 3.1 0.016 0.074 Xd $1.44\pm 0.05$ $1.65\pm 0.05$
$18.7\pm 0.8$ $16.4\pm 0.5$ 2.6 0.013 0.48 Xa $1.44\pm 0.04$ $1.64\pm 0.05$
$18.7\pm 0.6$ $16.4\pm 0.5$ 3.0 0.028 CO $\mbox{$\alpha_{\rm CO}$}{\rm(MS)}$
$\mbox{$\alpha_{\rm CO}$}{\rm(SMG)}$ daX $2.56\pm 0.14$ $2.61\pm 0.12$ 0.3
0.35 0.017 Xa $2.45\pm 0.12$ $2.70\pm 0.10$ 1.6 0.05 ad $3.32\pm 0.07$
$2.87\pm 0.09$ $3.8$ $0.0008$ 0.04 Dust $\mbox{$\kappa_{\rm H}$}{\rm(MS)}$
$\mbox{$\kappa_{\rm H}$}{\rm(SMG)}$ $\mbox{$\alpha_{850}$}{\rm(MS)}$
$\mbox{$\alpha_{850}$}{\rm(SMG)}$ daX $1758\pm 90$ $2039\pm 113$ $6.8\pm 0.3$
$7.7\pm 0.3$ 2.2 0.04 (0.25) 0.76 Xd $1757\pm 66$ $2025\pm 58$ $6.9\pm 0.2$
$7.4\pm 0.2$ 3.0 0.015 (0.41) 0.58 ad $1722\pm 36$ $1981\pm 65$ $6.0\pm 0.1$
$7.3\pm 0.2$ $3.6$ $0.0017$ 0.11 $(6.5)$ $(4\times 10^{-11})$

We compare MS galaxies and SMGs for each subset to look for differences in the
parameters. The MS group excludes C icor and lo-VALES galaxies, but it does
include the $f_{\rm H\,I}$>1 galaxies, after applying the correction from
Appendix B.1. The numbers in each subset are: daX: MS=46, SMG=55; Xd: MS=60,
SMG=79; Xa: MS=54, SMG=55; ad: MS=184, SMG=144. Bold indicates parameters
which are significantly different between the MS galaxies and SMGs in the
Z-test and KS tests. The C i parameters, $X_{\rm C\,I}$ and $\alpha_{\rm
C\,I}$, are simply linked due to our adoption of constant
$\mbox{$Q_{10}$}=0.48$, meaning that the distributions have the same KS
results. The dust parameters, $\kappa_{\rm H}$ and $\alpha_{850}$, are related
to each other as a function of $T_{\rm mw}$ and so they can behave
differently, e.g. $\kappa_{\rm H}$ can be indistinguishable between samples
but the $\alpha_{850}$ can be significantly different. Thus in the dust
section, there are two $P_{\rm KS}$ values: those for $\kappa_{\rm H}$ and
then, in parentheses, those for $\alpha_{850}$. The final column, $P_{\rm
KS}$(25 k), is the KS result when $T_{\rm mw}$ is fixed to 25 k, which makes
the distributions of $\kappa_{\rm H}$ and $\alpha_{850}$ identical.

Table 9: Fits to the conversion parameters and their tracer luminosities. All
quantities include He.

$y$ $x$ $m$ $c$ $r_{\rm s}$ $N$ Sample log $\alpha_{\rm CO}$ log
$L^{\prime}_{\rm CO}$ $-0.062$ (0.009) 1.25 (0.09) $-0.33$ 335 ad log
$\alpha_{850}$ log $L_{850}$ $0.052$ (0.008) 11.60 (0.18) 0.37 335 ad log
$\alpha_{\rm C\,I}$ log $L^{\prime}_{\rm C\,I}$ $-0.052$ (0.013) 1.73 (0.12)
$-0.3$ (0.0003) 140 Xd

The three empirical conversion factors, $\alpha_{850}$, $\alpha_{\rm C\,I}$
and $\alpha_{\rm CO}$, directly relate the observable tracer luminosity to a
gas mass, according to Eqn. 1. If only one tracer ($L^{\prime}_{\rm CO}$,
$L^{\prime}_{\rm C\,I}$ or $L_{850}$) is available, the empirical conversion
factor we have estimated in Table LABEL:prescT is the best choice. We adopt a
convention that the empirical parameters are referenced to $M_{\rm mol}$,
which includes a factor 1.36 for He.

Table 10: Empirical calibration factors derived from this study. $M_{\rm mol}$
columns include a factor of 1.36 for He.

$M_{\rm H2}$ $M_{\rm mol}$ Sample $N$ $\alpha_{850}$$\,(\times 10^{12})$
$\alpha_{\rm CO}$ $\alpha_{\rm C\,I}$ $\alpha_{850}$$\,(\times 10^{12})$
$\alpha_{\rm CO}$ $\alpha_{\rm C\,I}$ $\rm{W\,Hz^{-1}\mbox{$\rm
M_{\odot}$}^{-1}}$ M⊙ (k km s-1 pc2)-1 $\rm{W\,Hz^{-1}\mbox{$\rm
M_{\odot}$}^{-1}}$ M⊙ (k km s-1 pc2)-1 daX 101 $9.9\pm 0.2$ $2.6\pm 0.1$
$12.6\pm 0.4$ $7.3\pm 0.2$ $3.5\pm 0.1$ $17.2\pm 0.5$ ad (log $L_{\rm IR}>11$)
240 $9.1\pm 0.1$ $3.1\pm 0.1$ $6.7\pm 0.1$ $4.2\pm 0.1$ ad 326 $8.8\pm 0.1$
$3.2\pm 0.1$ $6.5\pm 0.1$ $4.3\pm 0.1$ Xd 140 $9.8\pm 0.3$ $12.7\pm 0.4$
$7.2\pm 0.2$ $17.3\pm 0.5$ Xa 109 $2.5\pm 0.1$ $12.5\pm 0.3$ $3.4\pm 0.1$
$17.0\pm 0.4$

Figs 6–9 show the empirical conversion factors as a function of $L_{\rm IR}$;
Fig. 10 shows the empirical conversion factors as a function of the tracer
luminosity, and Fig. 11 shows their distribution when the sample is split into
MS galaxies and SMGs (see Table 8 for details). All empirical factors show
significant but shallow correlations with the tracer luminosity. We have
carefully accounted for the co-variance between the $x$ and $y$ parameters
when fitting, so the correlations we find are not caused by the
involvement242424The inclusion of the co-variance matrix in the fit reduces
the slope (closer to zero) by 0.05. of $L_{850}$, $L^{\prime}_{\rm C\,I}$ and
$L^{\prime}_{\rm CO}$ in the derivation of $\alpha_{850}$, $\alpha_{\rm C\,I}$
and $\alpha_{\rm CO}$. Table LABEL:calfitsT lists the fit parameters. The
intrinsic scatter in all of these relationships is very small once the
measurement errors are accounted for. Correlations are also seen between
$L_{\rm IR}$ and $\alpha_{850}$, $\alpha_{\rm C\,I}$ and $\alpha_{\rm
CO}$(Fig. 7–9), albeit with more scatter.

A correlation between $\alpha_{850}$ and $L_{850}$ was also noted by Sco16
(shown as the red dashed line on our plot) but ours is somewhat shallower
($m=0.052$ vs $m=0.07$ from Sco16), although the difference is
unlikely252525Sco16 do not quote an error on their fit, so it is difficult to
be certain, but the error on Sco16 would likely be larger than ours, which
means they are consistent to within 2$\sigma$. to be significant. These
shallow but significant relationships with their tracer luminosity could be
further applied to give more accurate calibration (see Table LABEL:calfitsT).
The trends of $\alpha_{\rm C\,I}$ and $\alpha_{\rm CO}$ with their respective
tracers are explored for the first time in large numbers here.

###  5.2 Discussion of empirical factors relative to the literature.

####  5.2.1 Submillimetre dust empirical calibration, $\alpha_{850}$

The final calibration factors from this work are provided in Table LABEL:empT.
A compilation of $\alpha_{850}$ values from our optimal method262626There is
no significant difference if we fix $\alpha_{\rm CO}$ to 4.3. and those from
the literature, referenced to a common $\mbox{$\alpha_{\rm CO}$}=4.3$: the
Galactic value including He from Bolatto et al., 2013) is presented in Table
LABEL:asublitT. Literature values cover the range, $\alpha_{850}$ =
3.6–$10.1\times 10^{12}$, comfortably within the range of estimates here:
$\mbox{$\alpha_{850}$}^{\rm opt}=$ 6.5–$7.2\times 10^{12}$. The lowest value,
$\mbox{$\alpha_{850}$}=3.6^{+3.6}_{-1.9}\times 10^{12}$, comes from the local
sample of Orellana et al. (2017), who include H i as well as $M_{\rm mol}$
(from $L^{\prime}_{\rm CO}$). Their $\alpha_{850}$ refers to the total gas
mass – sensibly, since their lower luminosity sample is more H i-dominated
than the others we compare to – meaning that a lower value for $\alpha_{850}$
is required. The highest value, $\mbox{$\alpha_{850}$}=(10.1\pm 0.3)\times
10^{12}$, is from Sco16.272727The original value of
$\mbox{$\alpha_{850}$}=6.7\times 10^{12}$ quoted by Sco16 assumed that
$\mbox{$\alpha_{\rm CO}$}=6.5$ to calibrate the gas mass from CO. Re-
normalising the Sco16 result to the same $\mbox{$\alpha_{\rm CO}$}=4.3$ as our
literature comparison increases the Sco16 value to
$\mbox{$\alpha_{850}$}({\mbox{$\alpha_{\rm CO}$}=4.3})=10.1\times 10^{12}$.
There are two reasons why Sco16 found a significantly higher $\alpha_{850}$
compared to our analysis. The first is simple mathematics, as Sco16 quote a
linear mean for a distribution that has a significant tail to higher values;
in contrast, we quote a log-mean which is less sensitive to tails. This
statistical bias results in a linear mean for $\alpha_{850}$ which is 20 per
cent higher than the log-mean (Behroozi et al., 2013) . Assuming the shape of
our $\alpha_{850}$ distribution is similar to that from Sco16, we adjust their
linear mean down by 20 per cent to approximate our log-mean method. Thus our
log-mean estimate of the Sco16 value is $\mbox{$\alpha_{850}$}^{\rm LM}(\rm
Sco16)=8.4\times 10^{12}$. Secondly, there have been changes in the versions
of the Herschel pipeline data used as the basis for the local portion of the
Sco16 sample (see Appendix A). When we re-fitted the local galaxy SEDs to
estimate $T_{\rm mw}$ and $L_{850}$ using the most recent Herschel flux
densities (Chu et al., 2017; Clark et al., 2018) we found an increase in
$L_{850}$ of $\sim 0.1$ dex compared to that reported in Sco16282828This
difference is not just the photometry change; Sco16 used a different method to
estimate $L_{850}$ from the Herschel 500-$\mu$m fluxes.. Once these factors
are accounted for, the Sco16 result is comparable to ours.

Table 11: Summary of our empirical dust continuum–$M_{\rm mol}$ calibration
factor, $\alpha_{850}$, compared to literature values referenced to
$\mbox{$\alpha_{\rm CO}$}=4.3$ (i.e. Galactic $\alpha_{\rm CO}$, including
He).

.

$\alpha_{850}$($\times 10^{12}$) Sample $N_{\rm gal}$ Notes Reference W Hz-1
M${}_{\odot}^{-1}$ $6.4\pm 0.1$ all 328 log-mean opt this work $7.2\pm 0.2$
SMGs 144 log-mean opt this work $5.9\pm 0.1$ MS 184 log-mean opt this work
$10.1\pm 0.3$ local galaxies and SMGs 72 linear mean Scoville et al. (2016)
8.3 MS 30 linear mean Scoville et al. (2016) 12.7 SMGs 30 linear mean Scoville
et al. (2016) $3.6^{+3.6}_{-1.9}$ local galaxies 136 median H i $+$
4.3$L^{\prime}_{\rm CO}$ Orellana et al. (2017) $6.1\pm 0.14$ $z<0.4$
160-$\mu$m selected 41 log-mean (ex lovales) Hughes et al. (2017) $8.4\pm 1.0$
$z=1.6$–2.9 unlensed SMGs† 9 log-mean Kaasinen et al. (2019) $11.6\pm 1.2$
$z=1.6$–2.9 unlensed SMGs† 9 linear mean Kaasinen et al. (2019)

#### Atomic carbon empirical factor: $\alpha_{\rm C\,I}$

We compare our optimised $\alpha_{\rm C\,I}$ estimates with others from the
literature in Table LABEL:acilitT, and for reference we also compile the
literature values for $X_{\rm C\,I}$ in Table LABEL:XCIlitT, the two are
related simply by the excitation factor $Q_{10}$, as given in Eqn. 11. Our
values are a weighted average of the three samples that contain C i, where we
find $\langle\mbox{$\alpha_{\rm C\,I}$}\rangle=17.3\pm 0.3$ (standard error on
the mean). This compares well with the only truly independent measure,
$\mbox{$\alpha_{\rm C\,I}$}=21.4^{+13}_{-8}$, from absorber systems across a
range of redshift by Heintz & Watson (2020), but is considerably higher than
reported by many literature studies, e.g. $\mbox{$\alpha_{\rm
C\,I}$}=4.9$–10.3 for a large study of (U)LIRGs by Jiao et al., 2017 and
$\mbox{$\alpha_{\rm C\,I}$}=7.3^{+6.9}_{-3.6}$ for local disks by Crocker et
al. (2019). Many literature studies assume a fixed value for either $X_{\rm
C\,I}$ or $\alpha_{\rm CO}$ in order to derive $\alpha_{\rm C\,I}$ (typically
$\mbox{$X_{\rm C\,I}$}=3\times 10^{-5}$ or $\mbox{$\alpha_{\rm CO}$}=1$ for
high-$z$ SMG or local (U)LIRGs). These assumed values are very different from
those we have derived here under our minimal assumption that metal rich
galaxies have similar dust properties. Table LABEL:acilitT describes the
assumptions made for each literature source.

The MS galaxies in this work have a value of $\alpha_{\rm C\,I}$ of $19.1\pm
0.6$, again significantly higher than that found in Crocker et al. (2019).
However, J19, using a largely overlapping sample, derived a value of
$\mbox{$\alpha_{\rm C\,I}$}=19.9\pm 1.9$ using the same Herschel FTS C i
mapping observations. Crocker et al. use $L^{\prime}_{\rm[CO](2-1)}$ images
and spatially resolved $\alpha_{\rm CO}$ estimates from Sandstrom et al.
(2013), which were derived using a robust method which minimises the scatter
in the gas-to-dust ratio (see also Eales et al., 2012). While there are no a
priori assumptions about $\alpha_{\rm CO}$ in Crocker et al. (2019), implicit
assumptions are required for the average CO $r_{21}$ excitation. The limited
sensitivity of the FTS instrument to [C i](1–0) meant that C i was primarily
detected in the brighter nuclear regions, where $\alpha_{\rm CO}$ tends to be
lower than is typical in spiral disks ($\sim 1$ compared to $\sim 3$–4), (e.g.
Sandstrom et al., 2013) – as noted by Crocker et al. Their measured
$L^{\prime}_{\rm C\,I}$/$L^{\prime}_{\rm CO}$ ratios in the resolved regions
are compatible with other MS galaxies in our sample (though still higher than
the ratios derived for the same set of sources by J19, see Fig. 5), meaning
that the $\alpha_{\rm C\,I}$/$\alpha_{\rm CO}$ ratios are also similar. As
$\alpha_{\rm CO}$ in these regions is determined to be low in the Sandstrom et
al. (2013) analysis, the $\alpha_{\rm C\,I}$ inferred by Crocker et al. is
correspondingly lower as well.

Jiao et al. (2021) (henceforth J21) use CO(1–0), H i, [C i](1–0), dust
continuum and metallicity maps to investigate the variation of $\alpha_{\rm
C\,I}$ and $\alpha_{\rm CO}$ across the disks of six well-resolved local
galaxies from their J19 study. They use the FIR/submm dust maps from Spitzer
and Herschel and the method of Draine et al. (2007) to model the dust mass
across the galaxy and relate this to a gas mass via a relationship between
dust-to-gas and metallicity (Muñoz-Mateos et al., 2009; Draine et al., 2007;
Sandstrom et al., 2013). As they explicitly use the dust mass together with a
model of the $\delta_{\rm GDR}$ dependence on metallicity, their results are
normalised to the dust properties of the Draine & Li (2007) model (hereafter
DL07), which assumes $\mbox{$\delta_{\rm GDR}$}=100$ and
$\mbox{$\kappa_{850}$}=0.034$ m2 kg-1 at solar metallicity. Their self-
consistent DGR(ii) method derives weighted mean values of
$\langle\mbox{$\alpha_{\rm C\,I}$}\rangle=19.9\pm 1.9$ (including the
information from lower limits) and $\langle\mbox{$\alpha_{\rm
CO}$}\rangle=2.0\pm 0.3$ ($2.6\pm 0.4$) over the same area as the C i
observations (the entire CO detection region)292929We have multiplied the
values in J21 by 1.36 to include He for consistency with our convention.. In
the central region, the $\alpha_{\rm CO}$ values are significantly lower at
$\mbox{$\alpha_{\rm CO}$}^{\rm C}=1.5\pm 0.3$, while $\alpha_{\rm C\,I}$ is
not found to be significantly different with $\mbox{$\alpha_{\rm C\,I}$}^{\rm
C}=21.8\pm 0.5$. These values are comparable to to our average of
$\mbox{$\alpha_{\rm C\,I}$}=19.1\pm 0.6$ for MS galaxies, and to the average
value for $Z_{\odot}$ derived independently by Heintz & Watson (2020) of
$\mbox{$\alpha_{\rm C\,I}$}(\rm HW20)=21.4^{+13.3}_{-8.2}$303030While there
are differences in the normalisation for the dust mass model chosen by J21 and
ourselves, the introduction of a metallicity dependence for the dust-to-gas
ratio by J21 means that there is no simple way to scale their results to our
method. However, we can calculate their average ‘effective $\kappa_{\rm H}$’,
$\kappa_{\rm eff}=1300$, which indicates that J21 derive a lower H2 for a
given $L_{850}$ compared to our normalisation (and hence a lower value of
$\alpha_{\rm C\,I}$ and $\alpha_{\rm CO}$). However, the six galaxies in J21
are a subset of the C icor objects, at the low luminosity end where there are
potential decreases in $\kappa_{\rm H}$ and increases in $\alpha_{\rm C\,I}$..

Table 12: Summary of our empirical $\alpha_{\rm C\,I}$ calibration compared to
work from the literature; $\alpha_{\rm C\,I}$ is quoted including He.

.

$\alpha_{\rm C\,I}$ Sample $N_{\rm gal}$ Notes Reference M⊙ (k km s-1 pc2)-1
$17.0\pm 0.3$ log $L_{\rm IR}>11$ weighted average this work $19.1\pm 0.6$ MS
weighted average this work $16.2\pm 0.4$ SMGs weighted average this work
$10.3\pm 0.3$ (U)LIRGs 71 assuming $\mbox{$X_{\rm C\,I}$}=3\times 10^{-5}$
Jiao et al. (2017) $4.9\pm 0.3$ (U)LIRGs 71 CO(1–0) with $\mbox{$\alpha_{\rm
CO}$}=1.1$ Jiao et al. (2017) $7.3^{+6.9}_{-3.6}$ resolved local disks 18
CO(2–1) and resolved $\alpha_{\rm CO}$ from S13 Crocker et al. (2019) $19.9\pm
1.9$ resolved local disks 6 H i, CO(1–0), C i, dust, $Z$ ($\kappa_{\rm
eff}\sim 1300$) Jiao et al. (2021) $16.2\pm 7.9$ $z=3$ lensed SMGs 16
multi-$J$ CO, C i dust modelling Harrington et al. (2021)
$21.4^{+13.3}_{-8.2}$ GRB/QSO absorbers 19 $H_{2}$ and C i absorption lines at
$\rm{Z_{\odot}}$ Heintz & Watson (2020) 17.6 theory for $\mbox{$\zeta_{\rm
CR}$}=5\times 10^{-17}\rm{s^{-1}}$ Offner et al. (2014)

Table 13: Summary of our $X_{\rm C\,I}$ calibrations compared to other work in
the literature.

$X_{\rm C\,I}$($\times 10^{-5}$) Sample $N_{\rm gal}$ Notes Reference
$1.6^{+0.5}_{-0.4}$ $z=0$–5 log $L_{\rm IR}>11$ 90 $L_{850}$, CO, C i with
$\mbox{$\kappa_{\rm H}$}^{\rm N}=1884$ kg m-2 this work $2.5\pm 1.0$ local SF
11 CO(1–0) and $\mbox{$\alpha_{\rm CO}$}=1$ Jiao et al. (2019) $1.3\pm$ local
SF 9 CO(1–0) and $\alpha_{\rm CO}$ from S13 Jiao et al. (2019) $1.6\pm 0.7$
$z\sim 1.2$ MS 11 CO(2–1) and $\alpha_{\rm CO}$(Z) ($\langle\mbox{$\alpha_{\rm
CO}$}\rangle=3$) Valentino et al. (2018)† $2.0\pm 0.5$ $z\sim 1.2$ MS 11 dust
and $\delta_{\rm GDR}$(Z) ($\langle\mbox{$\delta_{\rm GDR}$}\rangle=134$)
Valentino et al. (2018)† $3.9\pm 0.4$ $z=2-3$ SMGs 14 CO(4–3), CO(1–0) and
$\mbox{$\alpha_{\rm CO}$}=1$ Alaghband-Zadeh et al. (2013)† $8.4\pm 3.5$
SMGs/QSOs 10 CO(3–2) and $\mbox{$\alpha_{\rm CO}$}=0.8$ Walter et al. (2011)
$8.3\pm 3.0$ local (U)LIRGS 23 CO(1–0) and $\mbox{$\alpha_{\rm CO}$}=0.8$ Jiao
et al. (2017); Jiao et al. (2019)† $0.9\pm 0.3$ $z=1$ ISM selected 2 CO(2–1),
C i and $\mbox{$\alpha_{\rm CO}$}=2.6$ Boogaard et al. (2020) $2.0\pm 0.4$
$z=1$ ISM selected 3 1.2 mm, C i and $\mbox{$\alpha_{850}$}=6.7\times 10^{12}$
from Sco16 Boogaard et al. (2020) ${}^{\ast}1.6^{+1.3}_{-0.7}$ $z=2$–4 GRB/QSO
absorbers 19 H2 and $\rm{C^{0}}$ absorption lines for $\rm{Z_{\odot}}$ Heintz
& Watson (2020) ${}^{\ast}7^{+7}_{-3.5}$ NGC 7469 (CND) 1 AGN, dynamical mass,
$L^{\prime}_{\rm C\,I}$ and $L^{\prime}_{\rm CO}$ Izumi et al. (2020)
${}^{\ast}1.4-5$ NGC 6240 1 $\alpha_{\rm CO}$ from CO–SLED, high-density
tracers Cicone et al. (2018) and two-phase LVG modelling Papadopoulos et al.
(2014)

#### CO empirical factor: $\alpha_{\rm CO}$

We compare our optimised $\alpha_{\rm CO}$ estimates with others from the
literature in Table LABEL:acolitT. Sophisticated LVG modelling with very large
datasets which include high-density gas tracers, optically thin CO
isotopologues, full CO SLEDs, and sometimes the C i lines and dust emission
(e.g. Weiß et al., 2007; Papadopoulos et al., 2012b; Papadopoulos et al.,
2014; Israel, 2020; Harrington et al., 2021) can break some of the model
degeneracies of the optically thick CO lines, though the method is still
reliant on assumptions for [CO/H2], isotopologue ratios, the number of
components allowed (single components give very different results to multiple
components) and the allowed range of velocity gradients in the models.

The best examples are NGC 6240 (Papadopoulos et al., 2014) and the Planck
lensed galaxies (Harrington et al., 2021) where detailed LVG modelling and
comprehensive datasets have sufficient constraints to break the degeneracies
which usually be-devil this method. The two-component LVG result for NGC 6240
is $\mbox{$\alpha_{\rm CO}$}=2$–4 (Papadopoulos et al., 2014) (cf.
$\alpha_{\rm CO}$=0.6 when using a single component LVG model Papadopoulos et
al., 2012b) and we can further use the ratio of $\mbox{$L^{\prime}_{\rm
C\,I}$}/\mbox{$L^{\prime}_{\rm CO}$}$ measured by Cicone et al. (2018) and our
relationship, $\mbox{$L^{\prime}_{\rm C\,I}$}/\mbox{$L^{\prime}_{\rm
CO}$}=\mbox{$\alpha_{\rm CO}$}/\mbox{$\alpha_{\rm
C\,I}$}=3324\,\mbox{$\alpha_{\rm CO}$}\mbox{$X_{\rm C\,I}$}$, to infer that
$\mbox{$X_{\rm C\,I}$}=1.4$–$2.9\times 10^{-5}$ in the starburst region. In
fact, our optimised values for this galaxy using global fluxes, are
$\mbox{$\alpha_{\rm CO}$}({\rm daX})=2.9\pm 0.6$, $\mbox{$X_{\rm C\,I}$}({\rm
daX})=(2.4\pm 0.5)\times 10^{-5}$, $\mbox{$\kappa_{\rm H}$}({\rm daX})=2800\pm
700$ ($\mbox{$\delta_{\rm GDR}$}=200$), in excellent agreement. The Planck
lensed galaxies analysed by Harrington et al. (2021) do not have the same
degeneracy-breaking lines used by (Papadopoulos et al., 2014) in their
analysis, but they do have multi-$J$ CO coverage and incorporate the C i lines
and the dust continuum emission in their model fitting, based on Weiß et al.
(2007). They assume similar dust parameters as we do for their normalisations
($\mbox{$\delta_{\rm GDR}$}=120$–150 with $\mbox{$\kappa_{850}$}=0.08$ m2
kg-1). With this, they infer an average $\mbox{$\alpha_{\rm CO}$}=3$–4 and an
average $\mbox{$\alpha_{\rm C\,I}$}=16.2\pm 7.9$ (incl. He), remarkably
consistent with our results, given our very simple approach.

Table 14: Summary of our empirical $\alpha_{\rm CO}$ calibrations compared to
work in the literature, $\alpha_{\rm CO}$ is quoted including a factor of 1.36
for He.

$\alpha_{\rm CO}$ Sample $N_{\rm gal}$ Notes Reference M⊙ (k km s-1 pc2)-1
$3.6^{+1.3}_{-1.0}$ $z=0$–5, log $L_{\rm IR}>11$ 90 $L_{850}$, CO, C i with
$\mbox{$\kappa_{\rm H}$}^{\rm N}=1884$ kg m-2 this work $4.2^{+1.8}_{-1.1}$
$z=0$–5, log $L_{\rm IR}>11$ 240 $L_{850}$, CO with $\mbox{$\kappa_{\rm
H}$}^{\rm N}=1884\,\mbox{kg\,m${}^{-2}$}$ this work $4.8^{+1.3}_{-1.1}$ local
MS, log $L_{\rm IR}<11$ 88 $L_{850}$, CO with $\mbox{$\kappa_{\rm H}$}^{\rm
N}=1884$ kg m-2 this work ${}^{\ast}3.1^{+3.1}_{-1.5}$ local disks 26 CO(2–1),
$r_{21}=0.7$, H i, dust Sandstrom et al. (2013) ${}^{\ast}4.2$ (3.5–5.4) MW
large scale $\gamma$-ray various Remy et al. (2017) ${}^{\ast}3.4\pm 2.1$
Planck lensed SMGs 24 LVG: multi-$J$ CO, C i and dust Harrington et al. (2021)
${}^{\ast}4.1^{+4}_{-2}$ NGC 7469 (CND) 1 AGN, dynamical mass,
$L^{\prime}_{\rm CO}$ Izumi et al. (2020) ${}^{\ast}2-4$ NGC 6240 1 LVG:
multi-$J$ CO, dense gas tracers Papadopoulos et al. (2014)
${}^{\ast}3.8^{+1.0}_{-0.7}$ $z=0$–5 22 CO, C i, Z and absorber based
$\alpha_{\rm C\,I}$ Heintz & Watson (2020) ${}^{\ast}4.4^{+2.0}_{-1.4}$ local
galaxies 24 C ii, CO(1–0) and modelling at $\rm{Z_{\odot}}$ Accurso et al.
(2017) ${}^{\ast}0.6\pm 0.2$ (U)LIRGs 28 LVG: Single component, multi-$J$ CO
Papadopoulos et al. (2012b) ${}^{\ast}2-6$ (U)LIRGs 28 LVG: Two-comp, free
d$V$/d$R$, dense-gas tracers Papadopoulos et al. (2012b) ${}^{\ast}3.9\pm 1.1$
local disks 9 CO(1–0), H i, dust Eales et al. (2012) $1.8\pm 0.5$ MW local
clouds 6 H i, CO(1–0), $\gamma$-ray Remy et al. (2017) $2.9\pm 0.5$ MW local
clouds 6 H i, CO(1–0), 850-$\mu$m dust Remy et al. (2017) 2.9 Taurus 1 H i,
CO(1–0) and extinction/reddening Chen et al. (2015) $2.4\pm 0.4$ local disks 7
CO(1–0), H i, dust Cormier et al. (2018) $1.9\pm 0.3$ resolved local disks 6 H
i, CO(1–0), [C i](1–0), dust, $Z$, $\kappa_{\rm eff}\sim 1300$ Jiao et al.
(2021) $3.2\pm 1.0$ $z=4$ lensed SMGs 9 CO(2–1), [C i](1–0) $\mbox{$X_{\rm
C\,I}$}=3\times 10^{-5}$ Bothwell et al. (2017)

###  5.3 Lack of bi-modality in the conversion factors

Our sample contains normal star forming galaxies – those obeying the
SFR–$M_{\star}$ correlation that forms as a result of the more intimate
relationship between SFR and H2 – as well as many extreme star-forming
systems, which belong to the (U)LIRG and high-$z$ submillimeter selected
samples. Here we remind the reader that we refer to the extreme SF group –
those that supposedly require a lower $\alpha_{\rm CO}$ – as ‘SMGs’, and the
normal star forming sources as ‘MS galaxies’, or sometimes just ‘MS’. As
mentioned in Section 2, the assignment of the galaxies to either category is
by nature of the data rather ‘fuzzy’ as we do not have a measure of SFR or
stellar mass for all sources, nor any homogeneous way to estimate them. We
thus rely on the categories used by previous authors where possible,
especially for high-$z$ sources. The $z=1$ galaxies from the samples of Bourne
et al. (2019); Valentino et al. (2018, 2020) are deemed to be ‘MS’, as are the
sources from ASPECs (Boogaard et al., 2020). Most low redshift sources with
log $L_{\rm IR}$<12 are classed as ‘MS’ though there are some exceptional LIRG
class sources in the local Universe which have extreme properties as evidenced
by their FIR, MIR lines and vibrational HCN (Díaz-Santos et al., 2017; Falstad
et al., 2021). We note that using a more conservative separation when
assigning galaxies into MS and extreme starburst categories does not change
any of the results. We therefore conclude that while our assignment of sources
into the two SF categories is not perfect, this categorisation is not capable
of masking any strong bi-modality in the observable ratios.

Fig. 11 and Table 8 detail the distributions of conversion factors for each
sample, split into MS galaxies and SMGs. While formally there are significant
differences in the parameters for some samples, these are very small – around
10–20 per cent in the mean, rather than the factor $\sim 3$–$4\times$ often
assumed for $\alpha_{\rm CO}$ (e.g. Downes & Solomon, 1998, $\alpha_{\rm
CO}$=0.8, derived for four ULIRGs). In fact, only the ad sample shows any
difference in $\alpha_{\rm CO}$ between the MS galaxies and SMGs, while the
estimates based on C i and CO, or on all three tracers, show no significant
difference. This is partially explained by the larger luminosity range in the
ad sample, combined with the previously noted negative correlation between
$\alpha_{\rm CO}$ and luminosity (Figs 9 and 10), with a factor $\sim 2\times$
reduction in $\alpha_{\rm CO}$ for a factor $\sim 100\times$ increase in
$L^{\prime}_{\rm CO}$. We cannot rule out that the correlation of $\alpha_{\rm
CO}$ with luminosity is the true reason that the ad sample shows a significant
difference between MS galaxies and SMGs313131$\mbox{$L_{\rm
IR}$}{\rm(ad)}{\rm(MS)}=10.95$ ($\mbox{$L^{\prime}_{\rm
CO}$}{\rm(ad)}{\rm(MS)}=9.31$) while $\mbox{$L_{\rm
IR}$}{\rm(Xa)}{\rm(MS)}=11.22$ ($\mbox{$L^{\prime}_{\rm
CO}$}{\rm(Xa)}{\rm(MS)}=9.52$). Using the relation in Table LABEL:calfitsT,
the expected $\Delta\mbox{$\alpha_{\rm CO}$}=\mbox{$\alpha_{\rm
CO}$}{\rm(MS)}-\mbox{$\alpha_{\rm CO}$}{\rm(SMG)}=0.36$ for the Xa sample and
– due to the lower numbers in the Xa and daX samples – such a difference would
not be detected at a significant level, if it existed..

This is not the first time323232However our current dataset is more
homogeneous, using only CO(1–0) or CO(2–1) and a consistent approach to
modelling the dust with our empirical relations for $T_{\rm mw}$. that lack of
bi-modality in $\alpha_{\rm CO}$ has been reported when compared to dust-based
determinations (e.g. Magdis et al., 2012; Rowlands et al., 2014; Genzel et
al., 2015). The range of $\alpha_{\rm CO}$ we find for SMGs (see Fig. 9) is
well within the framework set out by Papadopoulos et al. (2012a), who noted
that galaxies with a highly turbulent ISM (e.g. ULIRGs and SMGs) can have
$\alpha_{\rm CO}$ similar to galaxies with a much more quiescent ISM, the only
difference being that in a turbulent ISM, the distribution of gas mass as a
function of density is weighted to higher densities than in a less-turbulent
ISM.

Recent joint SLED/SED modelling of an exquisite dataset that includes CO, C i
and dust continuum for lensed SMGs (Harrington et al., 2021) finds a mean
$\mbox{$\alpha_{\rm CO}$}=3.4$–4.2 for these highly turbulent galaxies (albeit
with a large dispersion). The Harrington et al. radiative transfer models
employ a continuum distribution of molecular gas mass as a function of average
Mach number (and average density of the molecular cloud ensemble), making them
better equipped to ‘capture’ any re-distribution of the underlying molecular
gas mass towards higher densities. While important for other issues (e.g. the
initial conditions of star formation in SMGs/ULIRGs), such a re-distribution
in a highly turbulent ISM may actually leave $\alpha_{\rm CO}$ statistically
unaffected. The initial reports of a bimodal $\alpha_{\rm CO}$ factor in the
local Universe, with $\sim$4-5$\times$ lower values for ULIRGs than LIRGs and
ordinary spirals, can possibly be explained by a CO-luminous, strongly
unbound, low-density molecular gas component found preferentially in ULIRGs.
Such a component can dominate the global CO(1–0) line luminosities of
ULIRGs/SMGs (even if containing only small fractions of their total molecular
gas), while its large $\rm K_{vir}$ values will yield systematically low
$\alpha_{\rm CO}$ factors, under one-component LVG modelling (Equation
9).333333Also we must consider the size of the original sample – four ULIRGs
in the first study by Downes & Solomon (1998).

For individual galaxies, only multi-component models of SLED/SED (that also
include molecules/transitions tracing the dense gas) can properly account for
this effect (e.g. Papadopoulos et al., 2014; Harrington et al., 2021), while
for large galaxy samples, our cross-calibration of $\alpha_{\rm CO}$ against
the other two molecular gas mass tracers, is the most economical method. In
that regard it is worth noting that dust continuum is immune to the gas-
dynamics effects described above, i.e. a diffuse low-density, unbound, $\rm
H_{2}$ gas component will contribute very little to the total dust continuum
if its gas/dust mass is indeed low. The optically thin C i line emission will
also be much less sensitive than CO(1–0) to such gas-dynamics effects exactly
because of its low optical depths. These are perhaps the reasons why our
cross-calibration of $\alpha_{\rm CO}$ against dust and C i emission has not
uncovered any obvious bimodality of its values in MS galaxies compared to
SMGs.

The range of values we find for $\alpha_{\rm CO}$ is consistent with expected
values for $Z>0.5\,Z_{\odot}$ galaxies (Accurso et al. 2017, based on
calibrating $\alpha_{\rm CO}$ using Cii). Using their predictions, we would
expect $2.7<\mbox{$\alpha_{\rm CO}$}<15.2$ for the likely range of metallicity
and offset from the MS in our sample.

Figure 12: Results assuming a bi-modal behaviour for $\alpha_{\rm CO}$. We run
the optimisation process with the normalisation for SMGs set to
$\mbox{$\alpha_{\rm CO}$}^{\rm N}(\rm SMG)=0.8$, keeping MS galaxies at
$\mbox{$\kappa_{\rm H}$}^{\rm N}(\rm MS)=1884$ kg m-2(equivalent to
$\alpha_{\rm CO}$=2.8). This creates the strong bi-modality in $\alpha_{\rm
CO}$ but also induces a similarly strong or even stronger effect in
$\kappa_{\rm H}$ and $X_{\rm C\,I}$, which would be mirrored in $\alpha_{850}$
and $\alpha_{\rm C\,I}$.

Our underlying assumption: that the dust–gas properties of MS galaxies and
SMGs can be described as a uni-modal distribution with well defined mean and
scatter, is based on our finding that the luminosity ratios (Fig. 2) – the
most basic observables used in deriving the empirical conversion factors –
have such a distribution. They show no evidence for the strong bi-modality
advocated for $\alpha_{\rm CO}$ in some of the literature. That the
distribution of the observed luminosity ratio is, to first order, similar to
the distribution of the conversion parameters, is the simplest ‘Occam’s Razor’
assumption we can make.