# Tracing Molecular Gas Mass in ${z}\simeq 6$ Galaxies with [CII]

D. Vizgan11affiliation: Astronomy Department and Van Vleck Observatory,
Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA
22affiliation: Department of Astronomy, University of Illinois at Urbana-
Champaign, 1002 West Green St., Urbana, IL 61801, USA 33affiliation: The
Cosmic Dawn Center 44affiliation: National Space Institute, DTU Space,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark T. R.
Greve33affiliation: The Cosmic Dawn Center 44affiliation: National Space
Institute, DTU Space, Technical University of Denmark, 2800 Kgs. Lyngby,
Denmark 55affiliation: Dept. of Physics and Astronomy, University College
London, Gower Street, WC1E6BT London K. P. Olsen66affiliation: Department of
Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721,
USA A. Zanella77affiliation: Istituto Nazionale di Astrofisica (INAF), Vicolo
dell’Osservatorio 5, I-35122 Padova, Italy D. Narayanan33affiliation: The
Cosmic Dawn Center 88affiliation: Department of Astronomy, University of
Florida, Gainesville, Florida 32611 R. Davè99affiliation: Institute for
Astronomy, University of Edinburgh, Edinburgh EH9 3HJ, U.K. 1010affiliation:
University of the Western Cape, Bellville, Cape Town 7535, South Africa G. E.
Magdis33affiliation: The Cosmic Dawn Center 44affiliation: National Space
Institute, DTU Space, Technical University of Denmark, 2800 Kgs. Lyngby,
Denmark 1111affiliation: Niels Bohr Institute, University of Copenhagen,
Lyngbyvej 2, DK-2100 Copenhagen, Denmark G. Popping1212affiliation: European
Southern Observatory, 85478 Garching bei München, Germany F.
Valentino33affiliation: The Cosmic Dawn Center 1111affiliation: Niels Bohr
Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen, Denmark
K. E. Heintz33affiliation: The Cosmic Dawn Center 1111affiliation: Niels Bohr
Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen, Denmark
1313affiliation: Centre for Astrophysics and Cosmology, Science Institute,
University of Iceland, Dunhagi 5, 107 Reykjavík, Iceland

###### Abstract

We investigate the fine-structure [CII] line at $158\,\mu$m as a molecular gas
tracer by analyzing the relationship between molecular gas mass ($M_{\rm
mol}$) and [CII] line luminosity ($L_{\rm[CII]}$) in 11,125 $z\simeq 6$ star-
forming, main sequence galaxies from the simba simulations, with line emission
modeled by Sígame. Though most ($\sim 50-100\,\%$) of the gas mass in our
simulations is ionized, the bulk ($>50\,\%$) of the [CII] emission comes from
the molecular phase. We find a sub-linear (slope $0.78\pm 0.01$) $\log
L_{\rm[CII]}-\log M_{\rm mol}$ relation, in contrast with the linear relation
derived from observational samples of more massive, metal-rich galaxies at
$z\mathrel{\raise 1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}6$. We derive a median [CII]-to-$M_{\rm
mol}$ conversion factor of $\alpha_{\rm[CII]}\simeq 18\,{\rm
M_{\rm\odot}/L_{\rm\odot}}$. This is lower than the average value of $\simeq
30\,{\rm M_{\rm\odot}/L_{\rm\odot}}$ derived from observations, which we
attribute to lower gas-phase metallicities in our simulations. Thus, a lower,
luminosity-dependent, conversion factor must be applied when inferring
molecular gas masses from [CII] observations of low-mass galaxies. For our
simulations, [CII] is a better tracer of the molecular gas than CO $J=1-0$,
especially at the lowest metallicities, where much of the gas is ’CO-dark’. We
find that $L_{\rm[CII]}$ is more tightly correlated with $M_{\rm mol}$ than
with star-formation rate (${\rm SFR}$), and both the $\log L_{\rm[CII]}-\log
M_{\rm mol}$ and $\log L_{\rm[CII]}-\log{\rm SFR}$ relations arise from the
Kennicutt-Schmidt relation. Our findings suggest that $L_{\rm[CII]}$ is a
promising tracer of the molecular gas at the earliest cosmic epochs.

###### Subject headings:

galaxies: formation — galaxies: evolution — galaxies: ISM — galaxies: high-
redshift

##  1\. Introduction

Molecular gas and star formation are observed to be tightly correlated (e.g.
Kennicutt, 1998; Kennicutt & Evans, 2012; Bolatto et al., 2013) in our Galaxy
as well as in local and high-redshift galaxies, with the possible exception of
star-formation in extremely metal-poor environments where stars may form
directly from atomic gas inflows (Krumholz, 2012; Michałowski et al., 2015).
Unfortunately, the dominant component of molecular gas, H2, is invisible under
the typical conditions prevailing in the molecular interstellar medium (ISM).
Due to its lack of a permanent dipole moment, H2 only emits from warm
($>100\,{\rm K}$) gas (Draine, 2011), and its two lowest rotational
transitions have upper energy levels of $510\,{\rm K}$ and $1015\,{\rm K}$.
The low ($<10-50\,{\rm K}$) temperatures typical of the giant molecular clouds
(Cox, 2005), which make up the bulk of the star formation fuel, are unable to
excite H2 and are therefore ’$\rm{H_{2}}$-dark’ to observers. Our inability to
directly trace the distribution and amount of molecular gas in galaxies is a
long-standing challenge for extragalactic studies and has led to the
development of a number of indirect H2-tracing methods (Carilli & Walter,
2013; Li et al., 2018), which we briefly summarize below.

The $J=1-0$ rotational transition line of carbon monoxide (CO) is the most
rigorously calibrated and most regularly used tracer of molecular gas mass
(e.g., Dickman et al., 1986; Solomon & Barrett, 1991; Downes & Solomon, 1998;
Narayanan et al., 2012; Bolatto et al., 2013). The utility of the CO line
comes from the fact that it is relatively bright and easily excited in
molecular ISM conditions.

To convert the observed CO(1-0) line luminosity into a molecular gas mass
estimate, a CO-to-$\rm{H_{2}}$ conversion factor, $\alpha_{\rm CO}$, is
needed. The conversion factor is calibrated for giant molecular clouds (GMCs)
in our Galaxy and nearby galaxies (e.g., Bolatto et al., 2013). A galaxy-
averaged value for $\alpha_{\rm CO}$ has been determined for local IR-luminous
galaxies and found to be $\sim 4\times$ lower than the Galactic value (Downes
& Solomon, 1998), albeit with significant scatter (Papadopoulos et al., 2012).
For high-$z$ galaxies, we do not have strong constraints on $\alpha_{\rm CO}$
(Tacconi et al., 2013, 2018; Valentino et al., 2020a); both observations and
simulations suggest that $\alpha_{\rm CO}$ depends on the gas temperature,
velocity dispersion, and metallicity (Tacconi et al., 2008; Ivison et al.,
2011). In regions of high gas temperature and velocity dispersion, the
optically thick CO(1-0) emission line emits more luminosity per unit H2 column
density, which leads to lower conversion factors as seen in the local IR-
luminous galaxies. In low-metallicity galaxies, $\alpha_{\rm CO}$ scales
inversely with metallicity since in such environments, CO is photo-dissociated
in a larger fraction of the molecular gas, which is thus referred to as ’CO-
dark’ gas. (e.g., Wilson, 1995; Arimoto et al., 1996; Wolfire et al., 2010;
Narayanan et al., 2012; Genzel et al., 2012). Intense cosmic-ray ionization
rates that come along with extreme star-formation activity can destroy CO deep
within molecular clouds, thereby further degrading its H2-tracing capabilities
(Bisbas et al., 2015, 2017). Furthermore, at high redshifts, the cosmic
microwave background (CMB) becomes brighter than the low-$J$ CO emission from
cold molecular gas; for example, da Cunha et al. (2013) showed that at
$z\mathrel{\raise 1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}3$, the CMB could significantly affect
both CO(1-0) and CO(2-1) line fluxes, and neglecting these effects can lead to
non-negligible underestimations of the derived luminosities and, therefore,
the gas masses (see also Tunnard & Greve, 2016; Zhang et al., 2016).

The fine-structure transitions of neutral carbon –
[Ci](${}^{3}P_{1}$$-$${}^{3}P_{0}$) and [Ci](${}^{3}P_{2}$$-$${}^{3}P_{1}$)
have also been proposed as tracers of the molecular gas in galaxies (Gerin &
Phillips, 2000; Papadopoulos & Greve, 2004; Narayanan & Krumholz, 2017).
Carbon is abundant ($\sim 10^{-5}$ relative to H2), yet the [Ci] lines have
low to moderate optical depths ($\tau_{\rm\textbf{[C{\sc i}]}}\sim 0.1-1$).
Both lines are easily excited in typical molecular gas conditions, which,
together with its simple 3-level partition function, makes it a promising
tracer of H2. [Ci] surveys of molecular clouds in our Galaxy (Ikeda et al.,
2002) showed a remarkable spatial and kinematical similarity in the [Ci] and
CO emission, suggesting that [Ci] does indeed trace the bulk H2 reservoir.
Simulations of molecular clouds have shown that cloud inhomogeneity (resulting
in deeper UV penetration and thus an increase in the [Ci] layer; Offner et al.
2014) and turbulence (smoothing any initial [Ci]/CO abundance gradient; Glover
et al. 2015) are responsible for the observed [Ci]-CO concomitant. Offner et
al. (2014) further find that over a wide range of column densities [Ci] is a
comparable or superior tracer of the molecular gas, especially when an intense
UV field is impinging on it. Another work on absorption lines in high-redshift
galaxies suggests a universal metallicity-dependent [Ci]-to-H2 conversion
factor (Heintz & Watson, 2020). Recent work (Valentino et al., 2018, 2020b)
has also investigated the utility of [Ci] as a gas mass tracer in the high-$z$
universe.

The far-IR/sub-millimeter dust continuum emission has been used as an indirect
tracer of the total gas content (molecular and atomic) of the Galactic ISM
(Hildebrand, 1983), of nearby galaxies (Guelin et al., 1993, 1995; Eales et
al., 2012), and of high-$z$ galaxies (Magdis et al., 2012; Scoville et al.,
2014, 2016; Liang et al., 2018; Privon et al., 2018; Kaasinen et al., 2019).
One approach has been to derive the dust mass from careful modelling of the
far-IR/sub-millimeter spectral energy distribution (SED), and apply a dust-to-
gas mass ratio. The problem with this method is not only the uncertainties
associated with the SED fitting but a number of other factors, such as the
dust temperature, the grain emissivity, and the dust-to-gas mass conversion
factor. The latter depends on metallicity (Rémy-Ruyer et al., 2014; Li et al.,
2019), which implies a redshift dependence via the mass-metallicity relation.
Furthermore, this method requires flux measurements over a range of far-
IR/(sub-)mm bands. To circumvent these problems, a direct calibration between
the gas mass and the dust continuum emission at a single, broadband (sub-)mm
wavelength has been sought (e.g., Scoville et al., 2014). Such calibrations
have been established for optically thin continuum emission along the
Rayleigh-Jeans tail, e.g., at $850\,{\rm\mu m}$, and successfully applied to
both local and high-$z$ galaxies (Scoville et al., 2014, 2016). This method
has the advantage over CO and [Ci] in that it is easier to detect the
continuum emission in a single broadband than it is to detect a single line.
At redshifts $\mathrel{\raise
1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}4-5$, however, this method too runs into
the limitations set by the increasing CMB temperature, which will not only
couple thermally to the cold dust but also ’drown out’ a significant fraction
of the dust continuum emission (e.g. da Cunha et al., 2013; Zhang et al.,
2016).

In recent years, there has been a growing effort, from both the observation
and simulation sides, to explore the $158\,{\rm\mu m}$ ($1900.5\,{\rm GHz}$)
fine-structure transition of singly-ionized carbon ([Cii]) as a molecular gas
mass tracer (Accurso et al., 2017; Hughes et al., 2017; Zanella et al., 2018;
Dessauges-Zavadsky et al., 2020; Madden et al., 2020) The [Cii]158${\rm\mu m}$
line is one of the strongest cooling lines of the ISM, and can carry up to a
few per-cent of the total far-IR energy emitted from galaxies (e.g., Malhotra
et al., 1997). Carbon is the fourth most abundant element and has an
ionization potential of $11.3\,{\rm eV}$, lower than that of neutral hydrogen.
[CII], therefore, permeates much of the ISM. It is found in photodissociation
regions (PDRs), diffuse ionised and atomic regions, and even molecular gas
(e.g., Kaufman et al., 1999). Because the [Cii] line arises in most phases of
the ISM, it is important to disentangle the different contributions to
interpret galaxy-wide, integrated [Cii] observations correctly. How much of
the [Cii] emission comes from the molecular phase, and to what extent the line
can be used as a molecular gas tracer is of particular interest. Observations
of [CII] in the Milky Way, for instance, show that $\sim$75% of the [CII]
emission comes from dense PDRs and CO-dark H2 gas (Pineda et al., 2014).

Because of its association with PDRs, [CII] has been viewed as a tracer of the
star-formation rate (SFR) (e.g. De Looze et al., 2014; Magdis et al., 2014;
Herrera-Camus et al., 2015; Schaerer et al., 2020). At low redshift, metal-
poor dwarf galaxies and star-forming galaxies display slightly offset, but
tight log-linear relations between the [CII] luminosity and the sum of total
obscured and unobscured SFR (De Looze et al., 2014). However, compact
starbursts with high SFR surface densities and intense UV radiation fields
tend to exhibit a ’[CII]/FIR deficit’ by as much as an order of magnitude
compared to the local [CII]-SFR relations (e.g. Magdis et al., 2014; Díaz-
Santos et al., 2013; Narayanan & Krumholz, 2017). The [CII]/FIR deficit is
typically observed for local ultra-luminous IR galaxies (ULIRGS; $L_{\rm
IR}\geq 10^{12}\,{\rm\mbox{$\rm L_{\odot}\,$}}$) and high-$z$
starbursts/mergers, while normal star-forming galaxies at high-redshifts are
in line with the local [CII]-SFR relationship (e.g., Carniani et al., 2018).
Schaerer et al. (2020) also found that there is little evolution in the
[CII]-SFR relationship for main-sequence galaxies both locally and in the
high-redshift (i.e., $4<z<8$) universe, albeit the scatter is larger at higher
redshift.

By extension, this explains why [CII] traces molecular gas since the
Kennicutt-Schmidt law shows a strong linear relationship between SFR and
molecular gas mass. However, this relationship is not unique; from
observations, normal galaxies seem to differ from starbursts/mergers in that
the latter seem to deplete their gas reservoirs around 10 times faster. As
pointed out by Zanella et al. (2018), if one assumes [CII] traces molecular
gas mass, the [CII]/IR relationship becomes $M_{\rm mol}$/SFR, which is the
gas depletion timescale of a galaxy. In this case, the observed ’FIR deficit’
of starbursts reflects their $10\times$ shorter gas depletion time-scales
compared to normal star-forming galaxies. Theoretically, a strong connection
between [CII] luminosity and molecular gas mass is also expected from the fact
that the low critical density of carbon ($11.3\,$eV) allows for the excitation
of [CII] in dense,molecular regions by collisions with H2 molecules. Indeed,
simulations of high-redshift normal galaxies (e.g. Vallini et al., 2015;
Popping et al., 2019; Leung et al., 2020) found that most [CII] emission comes
from dense PDRs associated with regions of molecular gas; Though Olsen et al.
(2017) found that most ($\sim 64\%$) of [CII] emission arose from ionized gas
clouds in their simulations, they amended this result in an erratum, after
accounting for the CMB in Cloudy, and found that about 50% of [CII] emission
arose from molecular gas clouds (see Olsen et al., 2018).

The use of [CII] as a molecular gas mass tracer has been implemented on a
sample of 10 main-sequence galaxies at $z\sim 2$ by Zanella et al. (2018)
(henceforth; Z18). The authors combine Band 9 ALMA (Atacama Large
Millimeter/submillimeter Array) observations of [CII] emission with existing
multiwavelength observations to derive an average [CII] mass-to-luminosity
conversion factor $\alpha=31M_{\sun}/L_{\sun}$ with an uncertainty of 0.3 dex.
The molecular gas mass was derived from the integrated Schmidt-Kennicutt
relation using the SFR measured from pixel-by-pixel SED fitting, while an
independent and consistent gas mass estimate was derived from the dust
continuum. The resulting conversion factor was found to be mostly invariant
with galaxies’ redshift, depletion time, and gas-phase metallicity. At lower
redshifts [CII] has also been demonstrated to trace molecular gas mass (Madden
et al., 2020). Notably, [CII] has also been shown to trace atomic hydrogen in
galaxies until the epoch of reionization (Heintz et al., 2021).

The aim of this paper is to investigate whether [CII] is still a reliable
tracer of the molecular gas mass for lower-mass and lower-metallicity galaxies
at earlier epochs (i.e. z $\approx$ 6). Repeating the experiment on a much
larger sample of simulated high-redshift galaxies, we aim to derive a
conversion factor for the molecular gas mass of high-redshift galaxies using
[CII] line luminosities, as has been done in the nearby universe ($z\sim
0.03-0.2$) for CO (Hughes et al., 2017), and as done for galaxies in the local
Herschel Dwarf Galaxy Survey for [CII] (see Madden et al., 2020).

In this paper we present Sígame (Simulator of Galaxy Millimeter/Submillieter
Emission) simulations of [CII] line emission from 11,125 galaxies at $z\simeq
6$ in order to examine whether the line can be used to trace the gas content
in normal star-forming galaxies at cosmic dawn. Section 2 of the paper
describes the observation samples that we compare our simulations with in this
work; Section 3 describes the simulation sample and the post-processing done
by Sígame. Section 4 compares the integrated physical properties of the
simulation sample in the context of observations. Section 5 examines [CII] as
a tracer of molecular gas mass and as a tracer of star formation rate, derives
a [CII]-to-H2 conversion factor ($\alpha_{\rm[C{\scriptsize II}]}$), and
investigates the origin of the $M_{\rm mol}-L_{\rm[C{\scriptsize II}]}$
relation. Section 6 summarizes the main conclusions from our work.

##  2\. Observed comparison samples

We employ two observed galaxy samples for comparison with our simulations. The
first is from Z18, who presented six [CII]-detections out of a sample of ten
$z\sim 2$ main sequence galaxies observed with ALMA (Atacama Large
Millimetre/submillimetre Array), and combined these with other [CII]
detections from the literature. The latter spanned the redshift range $0\leq
z\leq 6$ and included diverse galaxy populations, from local dwarfs and normal
star-forming galaxies to luminous starbursts. Independent estimates of their
molecular gas masses for all of these galaxies were given in Z18, mostly from
existing observations of CO but also via more indirect methods such as from
parametrizations of the gas depletion time ($\tau_{\rm depl}=M_{\rm mol}/{\rm
SFR}$) as a function of the specific star-formation rate (${\rm sSFR}={\rm
SFR}/M_{\rm\star}$). In this paper, we split the Z18 sample into $z\geq 4$ and
$z<4$ galaxies.

We additionally compare our simulations to results obtained from the ALPINE-
ALMA [CII] Survey (henceforth referred to as ALPINE; see Le Fèvre et al.,
2020; Dessauges-Zavadsky et al., 2020). ALPINE has observed 118 galaxies in
the redshift range $4<z<6$ for [CII] and dust FIR emission. In 64% of galaxies
observed, [CII] was detected. These 75 galaxies, with which we compare our
simulations, consist mainly of main sequence galaxies, with only a handful of
sources lying either above or below the main sequence at these redshifts. The
molecular gas masses of the ALPINE sample have been estimated from their
restframe $850\,{\rm\mu m}$ continuum luminosity extrapolated from restframe
158 $\mu$m continuum; see Dessauges-Zavadsky et al. (2020) for details. The
SFRs of the ALPINE sample are derived, via Eq. 4 in Bell (2003), from their IR
luminosities, which were inferred from the IRX-$\beta$ relation (see Fudamoto
et al., 2020). These $L_{\rm IR}$-derived SFRs were added to the UV-derived
SFRs (Faisst et al., 2020) in order to get the total SFRs.

We use the Z18 and ALPINE samples in combination with our simulations (see
Section 3) to examine [CII] as a tracer of molecular gas.

##  3\. Simulations

###  3.1. Hydrodynamic simulations

This work builds on the analysis of snapshots taken from the simba suite of
cosmological galaxy formation simulations, which themselves were evolved using
the meshless finite mass hydrodynamics technique of gizmo (Hopkins, 2015,
2017; Davé et al., 2019). The simba simulation set consists of three cubical
volumes of 25, 50 and $100\,{\rm cMpc}\,h^{-1}$ on a side, all of which are
used in this work to search for galaxies at $z\sim 6$. For each volume, a
total of 10243 gas elements and 10243 dark matter particles are evolved from
$z=249$. Compared to its predecessor, mufasa, new features in simba include
the growth and feedback of supermassive black holes as well as a subgrid model
to form and destroy dust during the simulation run; for further details we
refer to Davé et al. (2019) and Li et al. (2018) for these two processes,
respectively. The galaxy properties in simba have been compared to various
observations across cosmic time (Thomas et al., 2019; Appleby et al., 2020),
including the epoch of reionization (Wu et al., 2020; Leung et al., 2020), and
are in reasonable agreement.

Since version 2 of Sígame used here relies heavily on the molecular gas mass
fraction in the simulated galaxies, we briefly describe how this fraction is
derived in simba, referring to Davé et al. (2020) for a more detailed
description. The molecular gas mass content of each fluid element in simba is
calculated at each time step following the subgrid prescription of Krumholz &
Gnedin (2011). This prescription relies on local metallicity and gas column
density, modified to account for variations in resolution (Davé et al., 2016).
The H2 mass function (H2MF) was found be in good agreement with observations
at redshifts $z\sim 0-2$, albeit overpredicting the H2MF at the high mass end
(Davé et al., 2020).

Importantly, when comparing the [CII] luminosity derived with Sígame to the
molecular gas content of the simulated galaxies, we will not be using the
molecular gas mass fractions from simba, but rather the re-calculated
molecular gas masses from Cloudy, as described in the following section. This
is because the molecular gas masses as calculated in simba are rough estimates
based on local gas densities ($\sim$ 1 $\times$ 102 cm3) and metallicities
within the simulations, whereas the post-processed molecular gas masses from
Cloudy account for the much higher densities ($\sim$ 105 – 106 cm3) of the
molecular clouds in the simulations. As a consistency check, we compared the
Sígame vs simba molecular gas masses and found, as expected, a strong
correlation between the two quantities.

From the simba simulations, 11,137 galaxies were selected via the yt (Turk et
al., 2011)-based package caesar; yt is a Python package used to analyze and
visualize volumetric data, and caesar is a six-dimensional friends-of-friends
algorithm which is applied to the simulated gas and stellar particles to
identify and select simba galaxies. Our final sample consists of 11,125
galaxies between a [CII] luminosity range of 103.82 to 108.91 L⊙ and a
molecular gas mass range of 106.25 to 1010.33 M⊙ The galaxy properties derived
with caesar include star formation rate, computed by dividing the stellar mass
formed over a 100 Myr timescale (Leung et al., 2020), star formation rate
surface density, stellar mass (M⋆) and total gas mass. In addition, we derive
a SFR-weighted gas-phase metallicity, $Z$. All of these properties will be
used for the analysis in Section 4.

###  3.2. Sígame post processing

The sample of simba galaxies are post-processed with version 2 of the Sígame
module (Olsen et al., 2017).111https://kpolsen.github.io/SIGAME/index.html At
its core, this version of Sígame uses the spectral synthesis code Cloudy
(v17.01; Ferland et al., 2013, 2017) to model the line emission from the
multi-phased ISM within each simulated galaxy. Using physically motivated
prescriptions, Sígame calculates, throughout the simulated galaxy, the local
interstellar radiation field (ISRF) spectrum, the cosmic ray (CR) ionization
rate, and the gas density distribution of the ionized, atomic and molecular
ISM phases. The local equilibrium gas temperatures are calculated for each ISM
phase by balancing the gas heating and cooling.

Sígame calculates the local ISRF that impinges on a gas particle in the
simulations as the sum of the radiation field from stars within one smoothing
length of that gas particle. The radiation from a stellar particle is the
integrated spectrum from starburst99 modeling of the stellar population
contained in the stellar particle, based on its mass, age and metallicity. The
local CR ionization rate is assumed to scale linearly with the local ISRF far-
UV strength. As for the metal abundances, gas-phase metallicities are carried
over from the simba simulations. Finally, the gas density requires a sub-
gridding of the fluid elements which are themselves too large ($\geq 10^{5}$
M⊙) to follow the intricate processes taking place inside giant molecular
clouds (GMCs) of mass $10^{4}$–$10^{6}$ M⊙. This sub-gridding is carried out
via subdivision of the gas into three ISM phases. First, the gas mass of each
fluid element is divided into a dense and a diffuse part based on the H2 gas
mass fraction of the simulation itself. The diffuse gas is further subdivided
into neutral and ionized diffuse gas, leading to three ISM phases; dense gas
which is later on organized in giant molecular clouds (GMCs), diffuse neutral
gas (DNG) and diffuse ionized gas (DIG). Finally, the level populations and
line emission are derived with Cloudy, which also provides mass fractions of
atomic, ionized and molecular hydrogen. It is from these mass fractions that
we derive the final molecular gas mass of each galaxy.

In general the molecular gas masses re-calculated in this manner do not exceed
those of the original simba simulation, as that sets the maximum mass of all
GMCs, of which Cloudy will deem a certain fraction (close to $100\%$) to be
molecular. For a full description of the main methods and assumptions of
Sígame, we refer to Olsen et al. (2017) and references therein. For this
paper, we have used a slightly updated version of Sígame and we refer to Leung
et al. (2020) for a description of these updates and the creation of the data
set analyzed here.

##  4\. Sample properties

###  4.1. Integrated properties

Figure 1.— Top: The ${\rm SFR}$ vs $M_{\rm\star}$ plane with the distribution
of our simulated galaxies (grey hexbins) shown, along with the ALPINE sample
(blue filled circles) and the $z>4$ galaxies from Z18 (magenta filled
circles). The lines show best-fit power-laws to observations of the galaxy
main sequence at $z\sim 5-6$ from Speagle et al. (2014) (solid black lines,
with 1-$\sigma$ scatter indicated by dotted lines), Schreiber et al. (2015)
(dashed purple), and Salmon et al. (2015) (dot-dashed green). The highlighted
segments of the lines indicate the stellar mass ranges of the galaxy samples
used to infer the main sequence. Middle: The integrated Kennicutt-Schmidt
relation (${\rm SFR}$ vs $M_{\rm mol}$) for our simulated and the ALPINE and
$z>4$ Z18 samples. Where they overlap, the simulations and observations follow
the same relation although the latter exhibit a larger scatter. Bottom: The
molecular fundamental mass-metallicity relation, i.e., metallicity
($12+\log({\rm O/H})$) vs the optimum linear combination of $\log
M_{\rm\star}$ and $\log M_{\rm H_{2}}$ as determined from a PCA analysis of
the local galaxy sample ALLSMOG (see Bothwell et al. (2016) for details). The
solid curve represents the quadratic fit to the ALLSMOG data, here shown using
the revised data values from Cicone et al. (2017). Our simulations as well as
the full Z18 sample are in good agreement with the ALLSMOG results.

In this section we perform a series of sanity checks of our simulated galaxies
by examining their star formation rates, masses (stellar and molecular), and
metallicities in the context of known empirical correlations between these
quantities. Before doing so, we stress that observed integrated properties
such as stellar and molecular gas masses are likely to be uncertain by at
least a factor of two (Lower et al., 2020; Gilda et al., 2021).

Fig. 1a shows the location of our simulated galaxies in the ${\rm
SFR}-M_{\rm\star}$ plane, where they define a sequence that extends over more
than three orders of magnitude in stellar mass. The three “clouds” of
simulations correspond with the three volumes of simba galaxies ($25$, $50$,
$100\,{\rm cMpc}\,h^{-1}$; see §3). The cloud with the highest mass galaxies
corresponds with the largest volume, and the cloud with the lowest mass
galaxies corresponds with the lowest volume. As explained in Leung et al.
(2020), the galaxies in the larger volumes have a larger smoothing length and
lower initial gas mass resolution than in the lowest volume.

Also shown in Fig. 1a, are the ALPINE and $z\geq 4$ Z18 samples. The galaxy
main sequence is usually modelled by a log-linear relation of the form:
$\log{\rm SFR}=\alpha\log M_{\rm\rm\star}+\beta$, where the parameters
$\alpha$ and $\beta$ are functions of cosmic time, or redshift. The
normalisation, $\beta$, in particular, is found to increase at earlier cosmic
epochs. Observational constraints on the galaxy main sequence at $z\simeq 5-6$
come from Speagle et al. (2014), Schreiber et al. (2015), Salmon et al.
(2015), Pearson et al. (2018), and Khusanova et al. (2020). The Speagle et al.
(2014) main sequence is based on observed galaxy samples at $z\simeq 4-6$ with
$M_{\rm\star}=10^{9.5}-10^{11.5}\,{\rm\mbox{$\rm M_{\odot}\,$}}$, while
Schreiber et al. (2015) derived their main sequence from a small sample of
galaxies at $z\simeq 3.5-5$ with $M_{\rm\star}\geq 10^{11}\,{\rm\mbox{$\rm
M_{\odot}\,$}}$, and an extrapolation to $z\simeq 6$. Finally, the Salmon et
al. (2015) main sequence is based on a sample of $z\simeq 6$ galaxies with
stellar masses in the range $M_{\rm\star}=10^{9}-10^{10.5}\,{\rm\mbox{$\rm
M_{\odot}\,$}}$. While the main sequence parametrizations by Speagle et al.
(2014), Salmon et al. (2015), and Khusanova et al. (2020) are in agreement
with the ALPINE and Z18 data, the parametrizations by Schreiber et al. (2015)
and Pearson et al. (2018) appear to have normalizations that fall below that
of the data (albeit the slopes appears to match). This may be because the
Schreiber et al. (2015) fit was made using $z<4$ data, and is thus an
extrapolation towards higher redshift. Nevertheless, the main sequence
parametrizations by Schreiber et al. (2015) and Pearson et al. (2018) are seen
to match our simulations across the full stellar mass range spanned by the
simulations, although these are based on extrapolations to lower masses than
those probed by the observations.

For stellar masses $\geq 2\times 10^{9}\,{\rm\mbox{$\rm M_{\odot}\,$}}$ there
is fair agreement between our simulations, and the ALPINE and $z\geq 4$ Z18
data. At lower stellar masses (down to $\sim 2\times 10^{8}\,{\rm\mbox{$\rm
M_{\odot}\,$}}$) there is an apparently increasing offset between the
simulations and the observational data, although, there are few observational
data available at these low masses to warrant a fair comparison. We note that
both the ALPINE and Z18 data may be somewhat biased toward high SFR-values
simply by virtue of having SFR-estimates at these low masses. Furthermore,
there may be a selection bias in the ALPINE sources due to the selected parent
sample of ALPINE; this is discussed further in Khusanova et al. (2021).

Fig. 1b shows ${\rm SFR}$ vs $M_{\rm mol}$ for our simulated galaxies, along
with $z\geq 4$ galaxies from our comparison samples. Our simulated galaxies
trace out a correlation which is consistent with the galaxy-integrated version
of the Kennicutt-Schmidt (KS) law (Kennicutt, 1998). We show parametrizations
of the redshift-independent galaxy-integrated KS-law for starbursts (yellow
line) and normal disk galaxies (green line) from Sargent et al. (2014). The
galaxies simulated within the $25\,{\rm cMpc}\,h^{-1}$ box are seen to agree
with the KS-law for starbursts, while the simulated galaxies extracted from
the $50\,{\rm cMpc}\,h^{-1}$ box fall slightly below the starburst KS-law.
Galaxies from the $100\,{\rm cMpc}\,h^{-1}$ box– as well as the observed
samples – fall in between the KS-law for starbursts and disks. They agree well
with the Z18 and ALPINE samples, which are seen to also straddle the KS-laws
for starbursts and normal disk galaxies.

Finally, in Fig. 1c we examine the relationship between gas-phase metallicity,
$M_{\rm\star}$, and $M_{\rm mol}$. Bothwell et al. (2016) showed that the SFR-
metallicity relation at a given stellar mass (the so-called fundamental
metallicity relation) is in fact a result of an underlying relation between
metallicity and molecular gas mass. Using the ALLSMOG sample of nearby
galaxies (Bothwell et al., 2016; Cicone et al., 2017), they derived a relation
between metallicity ($12+\log({\rm O/H})$), $M_{\rm\star}$, and $M_{\rm mol}$.
This relation is shown as the solid line in Fig. 1c. It represents the linear
combination in $\log M_{\rm\star}$ and $\log M_{\rm H_{2}}$, which yields the
least scatter in metallicity (see Bothwell et al. (2016) for details).

More recently, however, Cicone et al. (2017) reanalysed the ALLSMOG sample and
derived improved metallicities, as well as stellar and molecular gas masses.
Using these revised and final values, we plot the ALLSMOG sample in Fig. 1c
(orange symbols). We find that the galaxies from Z18 are consistent with the
ALLSMOG sample, independent of their redshifts, and even extends the relation
to the point where it turns over. To overlay the simba simulations, we use the
calibration $12+\log({\rm O/H})=\log(\langle Z\rangle_{\rm mass})+9$ derived
by Ma et al. (2016) using high-resolution cosmological simulations over the
redshift range $z=0-6$. Ma et al. (2016) uses the mass-weighted average gas-
phase metallicity, $\langle Z\rangle_{\rm mass}$, and for that reason we have
used this quantity from our simulations, and not the SFR-weighted average gas-
phase metallicities (see §3.1), to infer $12+\log({\rm O/H})$. Our simulations
are seen to be in good agreement with the observations in the region where
there is overlap, and they show a declining trend towards lower metallicities.

The analysis above demonstrates that within the range of masses where a
comparison with observations is possible, our simulated simba galaxies are in
overall agreement with key scaling relations observed for $z\geq 4$ main
sequence galaxies. This suggest that the simulations are representative, even
at the low stellar mass-end that extends well below the observed range, of the
bulk of the normal star forming galaxies at these epochs.

###  4.2. ISM properties

Before we explore the H2-tracing capabilities of [Cii], it is instructive to
first examine the gas mass fractions of the three ISM phases in our
simulations and their [Cii] emission properties. We note here that the
diffuse, molecular, and atomic gas masses are calculated to be the total mass
of all gas particles associated with each respective phase according to
caesar. However, for observations these phases are typically calculated within
a certain radius. Therefore the gas mass fractions are calculated to be
independent of galactic region or radius for our simulations.

Fig. 2 shows fractions of total [CII] luminosity ($f_{\rm[C{\scriptsize
II}]}$) and ISM gas mass ($f_{\rm mass}$) as a function of stellar mass. It is
seen that the molecular phase dominates the [CII] emission: it contributes
$\sim 50\%$ to the total [CII] emission for galaxies with $M_{\rm\star}\sim
10^{7}\,{\rm\mbox{$\rm M_{\odot}\,$}}$ and $\sim 100\%$ for
$M_{\rm\star}\mathrel{\raise
1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{9}\,{\rm\mbox{$\rm M_{\odot}\,$}}$
galaxies. Over the same stellar mass range, the ionized gas phase contribution
to the total [CII] emission spans $\sim 50\%$ at the low-mass-end to only a
few per-cent at the high-mass-end. The atomic gas hardly contributes anything
to the [CII] emission. Looking at the gas mass fraction, the ionized and
molecular gas phases dominate, with the former accounting for $\sim 50-100\%$
and the latter with $\sim 10-75\%$. Thus, while the molecular gas does not
account for most of the ISM mass in our simulated galaxies, except for in the
most massive galaxies ($M_{\rm star}\mathrel{\raise
1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{9.5}\,{\rm\mbox{$\rm
M_{\odot}\,$}}$), it is responsible for the bulk of the total [CII] emission.

In contrast, the atomic phase makes up $\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10\%$ of the total ISM mass, and even
less of the total [CII] emission. The low mass fraction for the atomic gas is
somewhat surprising, given that the cosmic atomic gas density exceeds that of
molecular gas at redshifts $z\mathrel{\raise
1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}3$ (Walter et al., 2020). We believe that
this is due to the way that gas is classified as ’ionized’ in Sígame, which is
simply when the electron-to-hydrogen fraction is $x_{e^{-}}>0.5$. As a result,
some of the atomic gas is lumped in with the ionized gas, and both the true
atomic gas mass and [CII] emission fraction are thus likely to be higher than
what is indicated in Fig. 2.

The fraction of the [CII] emission coming from the molecular phase in our
simulations is consistent with what is found in local star-forming galaxies
($\sim 60-80\%$; Accurso et al., 2017), as well as in high-resolution
simulations of $z\simeq 6$ Lyman-break galaxies ($\sim 95\%$; Pallottini et
al., 2017). The contribution to the [CII] emission from the ionized phase is
seen to increase as the stellar mass, and thus the metallicity, decreases.
Although, one would expect the [CII] cooling rate to decrease at lower
metallicity, this effect is negligible compared to the increase in CO photo-
dissociation rate, and thus available C+ ions, that comes with lower
metallicities; a similar effect was seen by Accurso et al. (2017). The fact
that the bulk of the [CII] emission is coming from the molecular gas suggests
that we can indeed use the line as a molecular gas tracer.

Figure 2.— Top: The fraction of the total [CII] luminosity emitted by the
molecular (blue), ionised (red), and atomic (green) ISM phases for each galaxy
as a function of its stellar mass. Bottom: the same as in the top panel but
now showing the gas mass fraction of the different ISM phases.

##  5\. Results and Discussion

###  5.1. How well does [CII] trace molecular gas mass?

In this section we examine whether there is a correlation between the [Cii]
emission and the molecular gas masses of our simulated $z\simeq 6$ galaxies,
and to what extent this correlation agrees with the empirical relation found
by observations.

Figure 3.—  $L_{\rm[CII]}$ vs $M_{\rm mol}$ (top) and $M_{\rm
mol}/L_{\rm[CII]}$ vs $M_{\rm mol}$ (bottom) for our simulated $z=6$ galaxies
(grey hexbin) and observed comparison samples from Z18 and Dessauges-Zavadsky
et al. (2020). Log-linear fits to our simulations (solid line), the observed
samples (dotted lined), and to simulations by Madden et al. (2020) of
nearby20warf galaxies (dot-dashed line) are indicated.

Fig. 3 shows the [Cii] luminosity vs the total molecular gas mass and $M_{\rm
mol}/L_{\rm[C{\scriptsize II}]}$ vs $M_{\rm mol}$ for our simulated galaxies,
along with observed $z<4$ and $z\geq 4$ galaxies from Z18, and $z=4-6$
galaxies from the ALPINE survey (Dessauges-Zavadsky et al., 2020). We fit a
log-linear model of the form $\log L_{\rm[C{\scriptsize II}]}=a\log M_{\rm
mol}+b$ to the simulations only, as well as the ALPINE data only. This yields
the following relations:

Simulations:

$\displaystyle\log L_{\rm[C{\scriptsize II}]}$ $\displaystyle=$
$\displaystyle(0.78\pm 0.01)\log M_{\rm mol}+(0.44\pm 0.04),$ (1)

shown as the thick solid line in Fig. 3. The r.m.s. scatter of the simulations
around the fitted line is $0.45\,{\rm dex}$.

ALPINE:

$\displaystyle\log L_{\rm[C{\scriptsize II}]}$ $\displaystyle=$
$\displaystyle(1.00\pm 0.01)\log M_{\rm mol}-(1.35\pm 0.1),$ (2)

shown as the thick dashed line in Fig. 3, with an r.m.s. scatter of the data
around the fitted line of $0.12\,{\rm dex}$.

For comparison, Z18 derived a relation given by $\log L_{\rm[C{\scriptsize
II}]}=(0.98\pm 0.02)\log M_{\rm mol}-(1.28\pm 0.01)$, which is shown as the
thin dotted line in Fig. 3. They reported an r.m.s. scatter of $0.2\,{\rm
dex}$. There is good agreement between the observed data from ALPINE and Z18
(see also Dessauges-Zavadsky et al., 2020), with the fit to the ALPINE data
being very similar to the Z18-relation. The Z18 data exhibit somewhat larger
scatter compared to the ALPINE data. This is likely due to the ALPINE galaxies
being uniformly selected as main-sequence galaxies at $z=4-6$, while the Z18
sample is heterogeneous, consisting of local (U)LIRGs, as well as main-
sequence galaxies and starburst at low and high redshifts.

Our simulations appear to be in broad agreement with the observations, despite
lying in a very different part of the ’galaxy parameter space’ (see §4). Upon
closer inspection, it is seen from Fig. 3 that our simulations contain a
number of galaxies which exhibit higher [CII] luminosities given their gas
masses than would be expected from the Z18 relation. On average, the
simulations have lower $M_{\rm mol}/L_{\rm[CII]}$ values than bulk of the
observed data. This is most pronounced for the lowest gas mass simulations
($M_{\rm mol}\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{7}\,{\rm M_{\odot}}$), while the
simulations with the highest gas masses ($M_{\rm mol}\sim 10^{9}\,{\rm
M_{\odot}}$) have $M_{\rm mol}/L_{\rm[CII]}$ ratios in line with the observed
data. The overall effect of this is a shallower slope in the simulations only
relation, and a higher r.m.s. scatter ($0.45\,{\rm dex}$). The shallower slope
($\sim 0.78$) suggested by our simulations implies a $L_{\rm[C{\scriptsize
II}]}$-dependent, i.e., non-universal $M_{\rm mol}/L_{\rm[C{\scriptsize II}]}$
ratio for $z\sim 6$ galaxies (cf. Z18). The small number of [CII] observations
of galaxies with gas masses $\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{8}\,{\rm M_{\odot}}$ prevents us from
assessing whether a similar trend is seen in observations.

From the $z\sim 0$ simulations presented by Madden et al. (2020), they derive:
$\log L_{\rm[C{\scriptsize II}]}=1.03\log M_{\rm mol}-2.19$, with an r.m.s.
scatter of $0.14\,{\rm dex}$. This relation, shown as the thin dash-dotted
line in Fig. 3, has a slope consistent with unity, i.e., similar to the Z18
and ALPINE relations, but the normalisation is significantly lower. Using this
relation would result in $\sim 3\times$ larger molecular gas masses, which
they ascribe to Z18 not accounting for CO-dark gas and thus using a CO-to-H2
conversion factor that is too low (Madden et al., 2020). This explanation,
however, is not borne out by our simulations, which examine the relationship
between $L_{\rm[C{\scriptsize II}]}$ and the full molecular gas reservoir (CO-
dark or not).

###  5.2. How well does [CII] trace star formation rate?

The traditional interpretation of the [CII] line has been as a diagnostic of
the physical conditions in photon-dominated regions, and as a tracer of the
star formation activity in galaxies (Stacey et al., 1991; De Looze et al.,
2014). In Fig. 4 we show $L_{\rm[C{\scriptsize II}]}$ vs ${\rm SFR}$ for our
simulations, the ALPINE, the Z18 sample, and the Carniani et al. (2018)
(henceforth C18) sample. Also shown is the log-linear fit to the combined Z18,
C18, and ALPINE samples, as well as our simulations.

For comparison the global $L_{\rm[C{\scriptsize II}]}-{\rm SFR}$ relation
inferred for high-$z$ starbursts by De Looze et al. (2014) is also shown,
along with relations derived from observations (Harikane et al., 2020) and
simulations (Lagache et al., 2018) of $z\simeq 6$ star-forming galaxies.

Figure 4.— $L_{\rm[CII]}$ vs. SFR (top) and SFR/$L_{\rm[CII]}$ vs. SFR
(bottom) for our simulated $z=6$ galaxies (grey hexbin), and the observed
comparison samples from Z18, Dessauges-Zavadsky et al. (2020), and Carniani et
al. (2018) (circles and crosses).

The r.m.s. scatter of the $\log L_{\rm[CII]}-\log{\rm SFR}$ relation fitted to
our simulations alone is $0.50\,{\rm dex}$. A fit to the Z18, ALPINE, and C18
data exhibits an even larger scatter, i.e., $0.55\,{\rm dex}$, which is higher
than the scatter reported by De Looze et al. (2014) for high-$z$ starbursts
($0.40\,{\rm dex}$). This is also higher than Schaerer et al. (2020), who
found a scatter of 0.28 dex when plotting $L_{\rm[CII]}$ vs. SFR of the ALPINE
sample after using the IRX-$\beta$ relation (Fudamoto et al., 2020) to add an
$L_{\rm IR}$-derived SFR to the UV-derived SFRs (Faisst et al., 2020, also see
§2).

The $\log L_{\rm[C{\scriptsize II}]}-\log{\rm SFR}$ relation exhibits more
scatter than the $\log M_{\rm mol}-\log L_{\rm[C{\scriptsize II}]}$ relation,
both for our simulations and for observed data, and we thus conclude that the
[CII] luminosity is more tightly correlated with the molecular gas mass than
the star-formation rate. Our simulations are seen to agree well with the
high-$z$ relation derived by De Looze et al. (2014), which also predict a
slightly decreasing ${\rm SFR}/L_{\rm[C{\scriptsize II}]}$ with increasing
${\rm SFR}$. This trend can account for some of the scatter in the $\log
L_{\rm[C{\scriptsize II}]}-\log{\rm SFR}$. The scatter can also arise from
other factors; we note, for instance, that SFR is prone to large observational
systematics – particularly at high-redshift – and this may be the reason why
the scatter in the $\log L_{\rm[C{\scriptsize II}]}-\log{\rm SFR}$ relation is
larger than the $\log M_{\rm mol}-\log L_{\rm[C{\scriptsize II}]}$ relation.
Furthermore, the methods of measuring SFR, as well as uncertainties on SFR
measurements, can affect the total scatter on the $\log L_{\rm[CII]}-\log{\rm
SFR}$ relation for our observations.

Figure 5.— Top panel: The distributions of $\alpha_{\rm[C{\scriptsize II}]}$
values for our simulations (grey horizontal band) and comparison samples.
Bottom panel: $\alpha_{\rm[C{\scriptsize II}]}$ vs redshift for our
simulations and comparison samples. The vertical dashed lines indicate the
median $\alpha_{\rm[CII]}$ values of the different samples. For clarity, our
$z=6$ simulations have been spread out in the $z$-direction slightly. Figure
6.— $\alpha_{\rm[CII]}$ vs offset from the main sequence for our simulations
(grey hexbins) and observed galaxies from Z18 (green and magenta symbols), and
the ALPINE survey (blue filled circles, Dessauges-Zavadsky et al., 2020).

###  5.3. The [CII] conversion factor

The fact that Z18 found the $L_{\rm[C{\scriptsize II}]}-M_{\rm mol}$ relation
to be log-linear with a slope close to unity allowed them to define a [Cii]
luminosity to molecular gas mass conversion factor, $\alpha_{\rm[C{\scriptsize
II}]}=M_{\rm mol}/L_{\rm[C{\scriptsize II}]}$. They inferred a value of
$\alpha_{\rm[C{\scriptsize II}]}\sim 30\,{\rm\mbox{$\rm
M_{\odot}\,$}/\mbox{$\rm L_{\odot}\,$}}$ with an uncertainty of $0.3\,{\rm
dex}$ averaged across their sample, and argued that $\alpha_{\rm[C{\scriptsize
II}]}$ was effectively invariant with redshift, metallicity, and galaxy type
(e.g., starburst vs main sequence galaxy), see also Fig. 3. For the ALPINE
sample, we find a median conversion factor of $\alpha_{\rm[C{\scriptsize
II}]}\sim 27\,{\rm\mbox{$\rm M_{\odot}\,$}/\mbox{$\rm L_{\odot}\,$}}$ and a
median absolute deviation (MAD) of $4\,\rm M_{\sun}/L_{\sun}$. From our
simulated galaxies, we find a median value of $\alpha_{\rm[C{\scriptsize
II}]}\sim 18\,{\rm\mbox{$\rm M_{\odot}\,$}/\mbox{$\rm L_{\odot}\,$}}$, with a
MAD of $10\,\rm M_{\sun}/L_{\sun}$.

In Fig. 5 we plot $\alpha_{\rm[C{\scriptsize II}]}$ vs $z$ for our
simulations, the Z18 and ALPINE samples along with histograms of their overall
$\alpha_{\rm[C{\scriptsize II}]}$ distributions (top panel). The
$\alpha_{\rm[C{\scriptsize II}]}$-distribution for our simulated galaxies
peaks at $15\,\rm M_{\sun}/L_{\sun}$ with an extended tail towards higher
values. The $\alpha_{\rm[C{\scriptsize II}]}$-distribution for the Z18 sub-
sample at $z<4$ also exhibit a tail towards high $\alpha_{\rm[C{\scriptsize
II}]}$-values. In addition to high $\alpha_{\rm[C{\scriptsize II}]}$-values,
the simulated sample also contain many galaxies with
$\alpha_{\rm[C{\scriptsize II}]}$-values $<10\,\rm M_{\sun}/L_{\sun}$. These
low values of $\alpha_{\rm[C{\scriptsize II}]}$ derived from our simulations
are in agreement with high-$z$ [CII] observations by Rizzo et al. (2021), who
inferred a median $\alpha_{\rm[C{\scriptsize II}]}=7^{+4}_{-1}$ M⊙/L⊙ for five
strongly-lensed starburst galaxies between $4\leq z\leq 5$ rizzo2020. Finally,
we note that the Z18 sub-sample at $z\geq 4$ overlaps with the
$\alpha_{\rm[C{\scriptsize II}]}$-values of the ALPINE sample, and both are
tightly distributed around their peak. For the ALPINE sample, in particular,
this might be due to it being a relatively homogeneous sample of main sequence
galaxies. In contrast, the full Z18 sample, as well as our simulations span a
wider range in star formation rate and masses.

In Fig. 6 we show $\alpha_{\rm[C{\scriptsize II}]}$ vs offset from the main
sequence for our simulated galaxies, which for a given galaxy is defined as
the ratio of its specific star formation rate (${\rm sSFR}={\rm
SFR}/M_{\rm\star}$) and the specific star formation of a galaxy with the same
stellar mass on the main sequence relation (at the given redshift). We adopt
the prescription for the main sequence at $z=6$ as given by Schreiber et al.
(2015), since it agrees well with our simulations (§4.1). The Z18 and ALPINE
comparison samples have also been normalized with the Schreiber et al. (2015)
main-sequence parametrisations at the relevant redshifts of each individual
galaxy. Our simulations show a weak trend of increasing $\alpha_{\rm[CII]}$
with main-sequence offset. In contrast, both the Z18 and ALPINE samples show
constant $\alpha_{\rm[CII]}$. This is true even for the sources with large
offsets from the main-sequence.

###  5.4. Is [Cii] a metallicity-invariant tracer of molecular gas?

In Fig. 7 (top) we show $\alpha_{\rm[C{\scriptsize II}]}$ vs gas-phase
metallicity (i.e., $12+\log({\rm O/H})$) for our simulations along with the
Z18 sample (note, we do not have metallicity estimates for the ALPINE sample).
Across the metallicity range probed by our simulations ($12+\log{\rm O/H}\sim
7.0$ to $\sim 9.0$), the $\alpha_{\rm[C{\scriptsize II}]}$ values show no
strong dependence on metallicity. As already noted, the simulated
$\alpha_{\rm[C{\scriptsize II}]}$ values are, on average, lower than observed
values, although with significant overlap between the two distributions. Only
few observations, however, are of galaxies with as low metallicities as the
simulations, and where there is overlap in terms of metallicity, the agreement
is good between the observed and simulated $\alpha_{\rm[C{\scriptsize II}]}$
values (although, the data points are sufficiently sparse that a quantitative
comparison is difficult).

Figure 7.— Top: $\alpha_{\rm[CII]}$ vs gas-phase metallicity for our simulated
galaxies (gray hexbins), along with observed samples from Z18 (green and
magenta symbols). The black solid line indicates the median
$\alpha_{\rm[CII]}$ value ($31\,{\rm M_{\odot}/L_{\odot}}$) found by Z18 and
its $0.2\,{\rm dex}$ scatter (dotted lines), while the red line indicates the
median $\alpha_{\rm[CII]}$ value ($18\,{\rm M_{\odot}/L_{\odot}}$) of our
$z=6$ simulations. Bottom: The CO(1-0)-to-[CII] luminosity ratio (both in
units of solar luminosities) vs the gas-phase metallicity. The luminosity
ratio, which scales with the ratio of the [CII]-to-CO(1-0) conversion factors
($\alpha_{\rm[C{\scriptsize II}]}/\alpha_{\rm CO}$), is seen to strongly
decrease with decreasing metallicity. This is a result of the increasing
photo-dissociation of CO in unshielded environments in metal-poor galaxies.
The blue and red lines are fits from (Zanella et al., 2018) and (Accurso et
al., 2017), respectively.

Fig. 7 suggests that the on average lower $\alpha_{\rm[CII]}$-value for the
simulations is, at least in part, due to their much lower metallicities. Z18
found $\alpha_{\rm[CII]}$ to be invariant with the metallicity. It should be
noted, however, that bulk of their sample fall within a relatively small range
of metallicities, making it hard to tease out any trend. Our simulations
suggest that there is in fact a weak dependency on metallicity, that becomes
apparent when one probes extremely low metallcities ($12+\log({\rm
O/H})\mathrel{\raise 1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}8.0$), i.e., much lower metallicities
than probed by Z18, and as low as those of local metal-poor dwarfs De Looze et
al. (2014).

Naively, one would expect $L_{\rm[C{\scriptsize II}]}$ to increase with
metallicity, due to a higher abundance of carbon. This in turn would imply a
decreasing trend in $\alpha_{\rm[C{\scriptsize II}]}$ with metallicity, which
is not what is seen in Fig. 7(top). However, an increase in metallicity also
implies a higher dust content; thus, fewer UV-photons are available to ionize
neutral carbon (photoionization potential of $11.3\,{\rm eV}$). It also
implies fewer UV-photons capable of photodissociating CO (requires photon
energies $>11.1\,{\rm eV}$) in neutral PDRs, which ultimately means less
carbon available for ionization. A higher dust content also implies more UV-
shielded regions, which would give rise to an increase in $\alpha_{\rm[CII]}$
with gas-phase metallicity.

In Sígame, the interstellar radiation field that impinges on the diffuse and
molecular clouds is attenuated by the dust inside the clouds. The fraction of
ionized carbon thus depends on the radiation field hitting the clouds, as well
as the dust fraction (which scales with the metallicity) of the clouds. CO as
a tracer of molecular gas is known to deteriorate in low-metallicity
environments, where the photo-dissociation of CO leads to a significant
fraction of the molecular ISM being CO-dark, (e.g. Narayanan & Krumholz, 2017;
Li et al., 2018; Madden et al., 2020).

This effect is illustrated in Fig. 7 (bottom), where the CO-to-[CII]
luminosity ratio ($L_{\rm CO}/L_{\rm[CII]}\propto\alpha_{\rm[CII]}/\alpha_{\rm
CO}$) is seen to sharply decrease with decreasing metallicity. This is a
direct result of $\alpha_{\rm CO}$ increasing with decreasing metallicity
(e.g., CO-dark gas), while $\alpha_{\rm[CII]}$ remains largely invariant as is
seen in Fig. 7 (top). Across the metallicity range $7.0\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}12+\log({\rm O/H})\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}9.0$, $\alpha_{\rm[CII]}/\alpha_{\rm CO}$
increases by more than two orders of magnitude, which is in line with a
similar trend found by Z18 and Accurso et al. (2017) (red and blue lines in
Fig. 7 (bottom)).

###  5.5. What is the origin of the $\log M_{\rm mol}-\log L_{\rm[CII]}$
relation?

In this section we examine whether combining the log-linear relations $\log
L_{\rm[CII]}-\log M_{\rm mol}$ and $\log L_{\rm[CII]}-\log{\rm SFR}$, derived
for our simulated galaxies in the previous section, leads to a better
determination of $M_{\rm mol}$. Both $L_{\rm[CII]}$ and ${\rm SFR}$ are
observables that can be measured with some confidence towards $z\sim 6$
galaxies.

Figure 8.— Equation 7 (top) and 8 (bottom) for our simulated galaxies, ALPINE
galaxies, and Z18 galaxies.

The unobscured star formation rate can be measured in a number of ways, from
the luminosity of the rest-frame continuum at $1500\,{\rm\AA}$, for example,
or from the Ly$\alpha$ emission. Obscured star formation rates are typically
inferred from (sub-)millimeter or far-IR continuum measurements. We consider
here the total star formation rates of our simulated galaxies and do not
distinguish between obscured or unobscured star formation rates. We will also
examine whether including the stellar masses might lead to more accurate gas
mass estimates. To this end, we perform a principal component analysis (PCA)
on the parameter space spanned by $\log L_{\rm[CII]}$, $\log M_{\rm mol}$,
$\log{\rm SFR}$, and $\log M_{\rm\star}$. This gives us the following
principal components:

$\displaystyle{\rm PC1}$ $\displaystyle=$ $\displaystyle 0.47\log
L_{\rm[CII]}+0.51\log{\rm SFR}+0.51\log M_{\rm mol}+0.51\log M_{\rm\star}$ (3)
$\displaystyle{\rm PC2}$ $\displaystyle=$ $\displaystyle-0.87\log
L_{\rm[CII]}+0.40\log{\rm SFR}+0.26\log M_{\rm mol}+0.14\log M_{\rm\star}$ (4)
$\displaystyle{\rm PC3}$ $\displaystyle=$ $\displaystyle-0.14\log
L_{\rm[CII]}-0.41\log{\rm SFR}-0.31\log M_{\rm mol}+0.85\log M_{\rm\star}$ (5)
$\displaystyle{\rm PC4}$ $\displaystyle=$ $\displaystyle 0.08\log
L_{\rm[CII]}+0.64\log{\rm SFR}-0.76\log M_{\rm mol}+0.05\log M_{\rm\star}.$
(6)

Of the four principal components, PC1 is responsible for 91% of the variance,
and PC2 and PC3 are responsible for 7 and 2%, respectively. Given its
negligible contribution to the variance, we can set PC4 to zero, which in turn
allows us to derive the following expression for the molecular gas:

$\log M_{\rm mol}=0.10\log L_{\rm[CII]}+0.85\log SFR+0.07\log M_{\rm\star}.$
(7)

The dominant principal component in this prescription is the star formation
rate, as expected given that star formation is fueled by molecular gas. The
[Cii] luminosity and stellar mass also account for some of the correlation
with molecular gas mass, each with about 10% of the contribution. We thus
interpret the $\log L_{\rm[C{\scriptsize II}]}-\log M_{\rm mol}$ relation at
high-$z$ as a second-order result of the Kennicutt-Schmidt Law, where [CII]
luminosity is a coincidental tracer of the molecular gas mass as a result of
the [CII]-SFR relation.

In Fig. 8 we have plotted eq. 7 for our simulations and comparison data, with
(top) and without the $M_{\rm\star}$-term (bottom). The residual scatter in
the two cases is virtually identical, i.e., $0.18$ and $0.19\,{\rm dex}$,
respectively, and we can therefore drop $M_{\rm\star}$, which is the most
difficult observable from eq. 7. Thus, a useful relation to obtain the
molecular gas mass to within a fractional scatter of $0.19\ln(10)\simeq 44\%$
is:

$\log M_{\rm mol}=0.10\log L_{\rm[CII]}+0.85\log\rm{SFR}.$ (8)

The conclusion from our PCA analysis is that the $\log L_{\rm[CII]}-\log
M_{\rm mol}$ relation comes about from a combination of the $\log
L_{\rm[CII]}-\log{\rm SFR}$ relationship and the Kennicutt-Schmidt relation,
which relates ${\rm SFR}$ to $M_{\rm mol}$.

Sommovigo et al. (2021) came to a similar conclusion in that they provided a
plausible physical explanation for the invariant $\alpha_{\rm[CII]}$ by
combing the $\Sigma_{\rm[CII]}-\Sigma_{\rm SFR}$ relation (De Looze et al.,
2014), and the $\Sigma_{\rm SFR}-\Sigma_{\rm mol}$ (Kennicutt, 1998; Kennicutt
& Evans, 2012). They found that $\alpha_{\rm[CII]}=30.3\tau_{\rm
depl}\Sigma_{\rm SFR}^{-0.075}$, where $\tau_{\rm depl}=M_{\rm mol}/{\rm SFR}$
is the molecular gas depletion time-scale in Gyr, and $\Sigma_{\rm SFR}$ is
the star-formation rate surface density in ${\rm\mbox{$\rm
M_{\odot}\,$}\,yr^{-1}\,kpc^{-2}}$. Given the weak dependence on $\Sigma_{\rm
SFR}$ and observational indications that for main sequence galaxies,
$\tau_{\rm depl}$ is approximately constant ($\sim 0.4-0.7\,{\rm Gyr}$) across
the redshift range $z\sim 0-2$ they find a nearly constant value of
$\alpha_{\rm[CII]}=(12-21)\Sigma_{\rm SFR}^{-0.075}$, which is in good
agreement with the Z18 values ($\simeq 30\,{\rm M_{\rm\odot}/L_{\rm\odot}}$)
as well as the $\alpha_{\rm[CII]}$ values for our $z=6$ simulations ($\simeq
18\,{\rm M_{\rm\odot}/L_{\rm\odot}}$).

Following Sommovigo et al. (2021), we derive an integrated version of the
above relation. To this end, we adopt ${\rm
SFR}=10^{-8.52}L_{\rm[C{\scriptsize II}]}^{1.18}$, which De Looze et al.
(2014) found applied to $z>0.5$ starburst galaxies. For the integrated
Kennicutt-Schmidt relation we adopt $M_{\rm mol}=10^{8.00}{\rm SFR^{0.83}}$,
which Sargent et al. (2014) found applied to starbursts. Combining the two, we
derive:

$\alpha_{\rm[C{\scriptsize II}]}=60.3\ \tau_{\rm depl}\ {\rm SFR}^{0.16},$ (9)

where $\tau_{\rm depl}$ is in Gyr. This relation is similar to that of
Sommovigo et al. (2021), both in terms of its normalization and the weak ${\rm
SFR}$-dependence, albeit with a positive exponent rather than a negative one.
In Fig. 9 we plot $\alpha_{\rm[C{\scriptsize II}]}$ for our simulations
derived from the Sommovigo et al. (2021) relation (top) and eq. 9 (bottom) vs
their $\alpha_{\rm[C{\scriptsize II}]}=M_{\rm mol}/L_{\rm[C{\scriptsize II}]}$
values.

Figure 9.— $\alpha_{\rm[CII]}$-values calculated for our simulations using the
relation from Sommovigo et al. (2021) (panel a)) and eq. 9 (panel b)) vs the
true $\alpha_{\rm[CII]}=M_{\rm mol}/L_{\rm[CII]}$ values for our simulations.
The red dotted lines denotes the 1:1 relationship.

It is seen that neither the $\alpha_{\rm[C{\scriptsize II}]}$ estimator
provided by Sommovigo et al. (2021) nor its volumetric version derived here
correlate well with the $\alpha_{\rm[C{\scriptsize II}]}=M_{\rm
mol}/L_{\rm[C{\scriptsize II}]}$ values of the simulations, and significant
scatter around the 1:1 relation, is seen along the horizontal direction.
However, both estimators agree in the mean with the simulation values. The
apparent bimodality seen in Fig. 9b is due to the fact that our simulations,
which come from three simulation boxes (see §3) span a wide range in both SFR
and $\tau_{\rm depl}$ (Fig. 1 top and bottom panels). The bottom cloud of
points in Fig. 9b is predominantly from the low-SFR and low-$\tau_{\rm depl}$
galaxies extracted from the $25\,{\rm cMpc}\,h^{-1}$ box, while the top cloud
consists of galaxies from the $50$ and $100\,{\rm cMpc}\,h^{-1}$ box. This
bimodality is not seen in Fig. 9a since i) the dependency on $\Sigma_{\rm
SFR}$ in the relation by Sommovigo et al. (2021)is significantly weaker and
ii) the simulations span a relatively narrow range in $\Sigma_{\rm SFR}$.

Extreme starburst galaxies have molecular gas depletion time-scales of $\sim
100\,{\rm Myr}$ (e.g., Aravena et al., 2016), i.e., of the order $4-7\times$
shorter than main-sequence galaxies, which would imply correspondingly lower
$\alpha_{\rm[C{\scriptsize II}]}$-values for starbursts (see eq. 9). The $\sim
5-10\times$ higher star-formation rates of starburst galaxies is not
sufficient to fully offset the effect of the shorter gas-depletion time-scales
on $\alpha_{\rm[C{\scriptsize II}]}$. Assuming that the gas depletion time-
scales and star-formation rates for the main-sequence galaxies (MS) and
starburst galaxies (SB) roughly scale as $\tau_{\rm
SB}\sim(0.4-0.7)\times\tau_{\rm MS}$ and ${\rm SFR_{\rm SB}\sim(5-10)\times
SFR_{\rm MS}}$, we find from eq. 9 that $\alpha_{\rm[C{\scriptsize
II}],SB}\sim(0.2-0.4)\times\alpha_{\rm[C{\scriptsize II}],MS}$. Inserting our
median value of $\alpha_{\rm[C{\scriptsize II}]}=18\,{\rm
M_{\rm\odot}/L_{\rm\odot}}$ from our simulations, then yields an expected
conversion factor for starbursts of $\alpha_{\rm[C{\scriptsize
II}],SB}\sim(4-7)\,{\rm M_{\rm\odot}/L_{\rm\odot}}$. This is consistent with
the low values found in high-$z$ dusty starburst galaxies by Rizzo et al.
(2021) (see §5.3).

##  6\. Conclusion

Using SÍGAME postprocessing of galaxies from simba cosmological simulations,
we have simulated and analysed the [CII] emission from 11,129 $z=6$ galaxies
with the goal of examining the viability of [CII] as a tracer of the ISM gas
mass, and in particular the molecular gas mass, in normal star-forming
galaxies at this epoch. The simba simulated galaxies span the stellar mass
range $1.5\times 10^{7}-5\times 10^{10}\,{\rm\mbox{$\rm M_{\odot}\,$}}$ and
define a star-forming main sequence at $z=6$ that at the high-mass-end agrees
with observations, and at the low-mass-end agrees with extrapolations of the
observed main-sequence Schreiber et al. (2015). The simulated galaxies also
follow the empirical Kennicutt-Schmidt law and the mass-metallicity relation.
Our main findings are presented below:

  1. 1.

For the most massive ($M_{\rm\star}\mathrel{\raise
1.50696pt\hbox{$\scriptstyle>$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{9}\,{\rm\mbox{$\rm M_{\odot}\,$}}$)
galaxies in our simulation, the [CII] emission is almost entirely coming from
the molecular gas phase. For galaxies with lower stellar masses ($M_{\rm
star}=10^{7}-10^{9}\,{\rm\mbox{$\rm M_{\odot}\,$}}$) there is a larger spread
in the fraction ($\sim 50-100\%$) of the total [CII] emission coming from the
molecular gas, but the general trend is that the average fraction decreases
towards lower masses. The atomic gas phase contributes a negligible amount
($\mathrel{\raise 1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}5\%$) to the total [CII] emission. As a
result, the contribution from the ionized ISM to the [CII] emission is the
mirror opposite to that of the molecular gas.

  2. 2.

We find and parametrize a log-linear correlation between the [CII] emission
and molecular gas masses of our galaxies: ${\log L_{\rm[C{\scriptsize
II}]}=(0.78\pm 0.01)\log M_{\rm mol}+(0.44\pm 0.04)}$. The slope is shallower
than the near-unity slope derived by Z18, and the scatter in the correlation
is $0.45\,{\rm dex}$, which is larger than the $0.2\,{\rm dex}$ scatter in the
Z18 relation. The shallower slope is due the simulated galaxies at the low-
mass end ($\mathrel{\raise
1.50696pt\hbox{$\scriptstyle<$}\kern-6.00006pt\lower
1.72218pt\hbox{{$\scriptstyle\sim$}}}10^{8}\,{\rm\mbox{$\rm M_{\odot}\,$}}$)
generally having higher [CII] luminosities compared to the extrapolation of
the linear Z18 relation to these low masses. However, the observations that
are available at the low-mass end (with the exception of Madden et al., 2020)
are consistent with the simulation predictions, albeit with a larger scatter.
Thus, the difference in slope might be attributable to a large amount of dwarf
galaxies in our sample which lack comparable observational data at high
redshifts. At the high-mass end, our simulations are in good agreement with
the Z18-relation and observational data from Z18 and ALPINE.

  3. 3.

We derive a [CII]-to-$M_{\rm mol}$ conversion factor,
$\alpha_{\rm[C{\scriptsize II}]}$, based on our simulations, finding a median
value of $\alpha_{\rm[C{\scriptsize II}]}=18$ M⊙/L⊙ with a median absolute
deviation of 10 M⊙/L⊙. This is lower than the average conversion factor
derived by Z18, although there is a significant overlap in the distributions
between our simulated and the observed $\alpha_{\rm[C{\scriptsize
II}]}$-values. We attribute the lower average $\alpha_{\rm[C{\scriptsize
II}]}$-values in our simulations to their much lower gas-phase metallicities,
which allow for a higher CO photo-dissociation rate and for neutral carbon
atoms in the ISM to be more readily ionized.

  4. 4.

Using principal component analysis, we determine a relation between molecular
gas mass, star formation rate, and [CII] emission, ${\log M_{\rm mol}=0.10\log
L_{\rm[CII]}+0.85\log\rm{SFR}}$, which has a scatter of only $0.18\,{\rm
dex}$. Our analysis suggests that the tight relationship between the star-
formation rate and the molecular gas mass is ultimately responsible for [CII]
as a tracer of molecular gas mass.

###  6.1. Acknowledgements

We thank the referee for valuable insight and constructive feedback which
helped to significantly improve the results of this paper. We thank Robert
Thompson for developing caesar, and the yt team for development and support of
yt. This research was made possible by the National Science Foundation
(NSF)-funded DAWN-IRES program and would not be possible without support from
the Cosmic Dawn Center (DAWN) in Copenhagen, Denmark. The Cosmic Dawn Center
(DAWN) is funded by the Danish National Research Foundation under grant No.
140. DV is funded by an Open Study/Research Award from the Fulbright U.S.
Student Program in Denmark. KPO is funded by NASA under award No
80NSSC19K1651. GEM acknowledges the Villum Fonden research grant 37440, ”The
Hidden Cosmos” and the Cosmic Dawn Center of Excellence funded by the Danish
National Research Foundation under then grant No. 140. DN acknowledges support
from the NSF via grant AST 1909153. RD acknowledges support from the Wolfson
Research Merit Award program of the U.K. Royal Society. KEH acknowledges
support by a Postdoctoral Fellowship Grant (217690–051) from The Icelandic
Research Fund and from the Carlsberg Foundation Reintegration Fellowship Grant
CF21-0103. simba was run on the DiRAC@Durham facility managed by the Institute
for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. The
equipment was funded by BEIS capital funding via STFC capital grants
ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC
operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

## References

  * Accurso et al. (2017) Accurso, G., Saintonge, A., Bisbas, T. G., & Viti, S. 2017, MNRAS, 464, 3315
  * Appleby et al. (2020) Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053
  * Aravena et al. (2016) Aravena, M., Spilker, J. S., Bethermin, M., et al. 2016, MNRAS, 457, 4406
  * Arimoto et al. (1996) Arimoto, N., Sofue, Y., & Tsujimoto, T. 1996, PASJ, 48, 275
  * Bell (2003) Bell, E. F. 2003, ApJ, 586, 794
  * Bisbas et al. (2015) Bisbas, T. G., Papadopoulos, P. P., & Viti, S. 2015, ApJ, 803, 37
  * Bisbas et al. (2017) Bisbas, T. G., van Dishoeck, E. F., Papadopoulos, P. P., et al. 2017, ApJ, 839, 90
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Bothwell et al. (2016) Bothwell, M. S., Maiolino, R., Cicone, C., Peng, Y., & Wagg, J. 2016, A&A, 595, A48
  * Carilli & Walter (2013) Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105
  * Carniani et al. (2018) Carniani, S., Maiolino, R., Amorin, R., et al. 2018, MNRAS, 478, 1170
  * Cicone et al. (2017) Cicone, C., Bothwell, M., Wagg, J., et al. 2017, A&A, 604, A53
  * Cox (2005) Cox, D. P. 2005, ARA&A, 43, 337
  * da Cunha et al. (2013) da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827
  * Davé et al. (2020) Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146
  * Davé et al. (2016) Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265
  * De Looze et al. (2014) De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62
  * Dessauges-Zavadsky et al. (2020) Dessauges-Zavadsky, M., Ginolfi, M., Pozzi, F., et al. 2020, A&A, 643, A5
  * Díaz-Santos et al. (2013) Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68
  * Dickman et al. (1986) Dickman, R. L., Snell, R. L., & Schloerb, F. P. 1986, ApJ, 309, 326
  * Downes & Solomon (1998) Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Eales et al. (2012) Eales, S., Smith, M. W. L., Auld, R., et al. 2012, ApJ, 761, 168
  * Faisst et al. (2020) Faisst, A. L., Schaerer, D., Lemaux, B. C., et al. 2020, ApJS, 247, 61
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, RMxAA, 49, 137
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMxAA, 53, 385
  * Fudamoto et al. (2020) Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4
  * Genzel et al. (2012) Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69
  * Gerin & Phillips (2000) Gerin, M., & Phillips, T. G. 2000, ApJ, 537, 644
  * Gilda et al. (2021) Gilda, S., Lower, S., & Narayanan, D. 2021, ApJ, 916, 43
  * Glover et al. (2015) Glover, S. C. O., Clark, P. C., Micic, M., & Molina, F. 2015, MNRAS, 448, 1607
  * Guelin et al. (1995) Guelin, M., Zylka, R., Mezger, P. G., Haslam, C. G. T., & Kreysa, E. 1995, A&A, 298, L29
  * Guelin et al. (1993) Guelin, M., Zylka, R., Mezger, P. G., et al. 1993, A&A, 279, L37
  * Harikane et al. (2020) Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93
  * Heintz & Watson (2020) Heintz, K. E., & Watson, D. 2020, ApJ, 889, L7
  * Heintz et al. (2021) Heintz, K. E., Watson, D., Oesch, P., Narayanan, D., & Madden, S. C. 2021, arXiv e-prints, arXiv:2108.13442
  * Herrera-Camus et al. (2015) Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53
  * Hopkins (2017) —. 2017, arXiv e-prints, arXiv:1712.01294
  * Hughes et al. (2017) Hughes, T. M., Ibar, E., Villanueva, V., et al. 2017, A&A, 602, A49
  * Ikeda et al. (2002) Ikeda, M., Oka, T., Tatematsu, K., Sekimoto, Y., & Yamamoto, S. 2002, ApJS, 139, 467
  * Ivison et al. (2011) Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913
  * Kaasinen et al. (2019) Kaasinen, M., Scoville, N., Walter, F., et al. 2019, ApJ, 880, 15
  * Kaufman et al. (1999) Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795
  * Kennicutt (1998) Kennicutt, R. C. 1998, ARA&A, 36, 189
  * Kennicutt & Evans (2012) Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531
  * Khusanova et al. (2020) Khusanova, Y., Le Fèvre, O., Cassata, P., et al. 2020, A&A, 634, A97
  * Khusanova et al. (2021) Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2021, A&A, 649, A152
  * Krumholz (2012) Krumholz, M. R. 2012, ApJ, 759, 9
  * Krumholz & Gnedin (2011) Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36
  * Lagache et al. (2018) Lagache, G., Cousin, M., & Chatzikos, M. 2018, A&A, 609, A130
  * Le Fèvre et al. (2020) Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1
  * Leung et al. (2020) Leung, T. K. D., Olsen, K. P., Somerville, R. S., et al. 2020, ApJ, 905, 102
  * Li et al. (2019) Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425
  * Li et al. (2018) Li, Q., Narayanan, D., Davè, R., & Krumholz, M. R. 2018, ApJ, 869, 73
  * Liang et al. (2018) Liang, L., Feldmann, R., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 478, L83
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140
  * Madden et al. (2020) Madden, S. C., Cormier, D., Hony, S., et al. 2020, arXiv e-prints, arXiv:2009.00649
  * Magdis et al. (2012) Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6
  * Magdis et al. (2014) Magdis, G. E., Rigopoulou, D., Hopwood, R., et al. 2014, ApJ, 796, 63
  * Malhotra et al. (1997) Malhotra, S., Helou, G., Stacey, G., et al. 1997, ApJ, 491, L27
  * Michałowski et al. (2015) Michałowski, M. J., Gentile, G., Hjorth, J., et al. 2015, A&A, 582, A78
  * Narayanan & Krumholz (2017) Narayanan, D., & Krumholz, M. R. 2017, MNRAS, 467, 50
  * Narayanan et al. (2012) Narayanan, D., Krumholz, M. R., Ostriker, E. C., & Hernquist, L. 2012, MNRAS, 421, 3127
  * Offner et al. (2014) Offner, S. S. R., Bisbas, T. G., Bell, T. A., & Viti, S. 2014, MNRAS, 440, L81
  * Olsen et al. (2017) Olsen, K., Greve, T. R., Narayanan, D., et al. 2017, ApJ, 846, 105
  * Olsen et al. (2018) —. 2018, ApJ, 857, 148
  * Pallottini et al. (2017) Pallottini, A., Ferrara, A., Gallerani, S., et al. 2017, MNRAS, 465, 2540
  * Papadopoulos & Greve (2004) Papadopoulos, P. P., & Greve, T. R. 2004, ApJ, 615, L29
  * Papadopoulos et al. (2012) Papadopoulos, P. P., van der Werf, P., Xilouris, E., Isaak, K. G., & Gao, Y. 2012, ApJ, 751, 10
  * Pearson et al. (2018) Pearson, W. J., Wang, L., Hurley, P. D., et al. 2018, A&A, 615, A146
  * Pineda et al. (2014) Pineda, J. L., Langer, W. D., & Goldsmith, P. F. 2014, A&A, 570, A121
  * Popping et al. (2019) Popping, G., Narayanan, D., Somerville, R. S., Faisst, A. L., & Krumholz, M. R. 2019, MNRAS, 482, 4906
  * Privon et al. (2018) Privon, G. C., Narayanan, D., & Davé, R. 2018, ApJ, 867, 102
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31
  * Rizzo et al. (2021) Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952
  * Salmon et al. (2015) Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183
  * Sargent et al. (2014) Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19
  * Schaerer et al. (2020) Schaerer, D., Ginolfi, M., Béthermin, M., et al. 2020, A&A, 643, A3
  * Schreiber et al. (2015) Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74
  * Scoville et al. (2014) Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84
  * Scoville et al. (2016) Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83
  * Solomon & Barrett (1991) Solomon, P. M., & Barrett, J. W. 1991, in Dynamics of Galaxies and Their Molecular Cloud Distributions, ed. F. Combes & F. Casoli, Vol. 146, 235
  * Sommovigo et al. (2021) Sommovigo, L., Ferrara, A., Carniani, S., et al. 2021, MNRAS, 503, 4878
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15
  * Stacey et al. (1991) Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423
  * Tacconi et al. (2008) Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246
  * Tacconi et al. (2013) Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179
  * Thomas et al. (2019) Thomas, N., Davé, R., Anglés-Alcázar, D., & Jarvis, M. 2019, MNRAS, 487, 5764
  * Tunnard & Greve (2016) Tunnard, R., & Greve, T. R. 2016, ApJ, 819, 161
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
  * Valentino et al. (2018) Valentino, F., Magdis, G. E., Daddi, E., et al. 2018, ApJ, 869, 27
  * Valentino et al. (2020a) Valentino, F., Daddi, E., Puglisi, A., et al. 2020a, A&A, 641, A155
  * Valentino et al. (2020b) Valentino, F., Magdis, G. E., Daddi, E., et al. 2020b, ApJ, 890, 24
  * Vallini et al. (2015) Vallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813, 36
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111
  * Wilson (1995) Wilson, C. D. 1995, ApJ, 448, L97
  * Wolfire et al. (2010) Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  * Wu et al. (2020) Wu, X., Davé, R., Tacchella, S., & Lotz, J. 2020, MNRAS, 494, 5636
  * Zanella et al. (2018) Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976
  * Zhang et al. (2016) Zhang, Z.-Y., Papadopoulos, P. P., Ivison, R. J., et al. 2016, Royal Society Open Science, 3, 160025