#  $\rm[C_{II}]$ 158 $\rm\mu m$ emission as an indicator of galaxy star
formation rate

Lichen Liang1,2, Robert Feldmann2, Norman Murray1,3, Desika Narayanan4,5,6,
Christopher C. Hayward7, Daniel Anglés-Alcázar8,7, Luigi Bassini2, Alexander
J. Richings9,10 Claude-André Faucher-Giguère11, Dongwoo T. Chung1,12, Jennifer
Y. H. Chan1,12, Onur Çatmabacak2, Dušan Kereš13, Philip F. Hopkins14
1Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto, ON M5S 3H8, Canada
2Institute for Computational Science, University of Zurich,
Winterthurerstrasse 190, Zurich CH-8057, Switzerland
3Canada Research Chair in Theoretical Astrophysics
4Department of Astronomy, University of Florida, 211 Bryant Space Sciences
Center, Gainesville, FL 32611, USA
5University of Florida Informatics Institute, 432 Newell Drive, CISE Bldg
E251, Gainesville, FL 32611, USA
6Cosmic Dawn Center at the Niels Bohr Institute, University of Copenhagen and
DTU-Space, Technical University of Denmark, Denmark
7Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA
8Department of Physics, University of Connecticut, 196 Auditorium Road,
U-3046, Storrs, CT 06269-3046, USA
9E. A. Milne Centre for Astrophysics, University of Hull, Cottingham Road,
Hull, HU6 7RX, UK
10DAIM, University of Hull, Cottingham Road, Hull, HU6 7RX, UK
11Department of Physics and Astronomy and CIERA, Northwestern University,
Evanston, IL 60208, USA
12Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto, ON M5S 3H4, Canada
13Department of Physics, Center for Astrophysics and Space Sciences,
University of California at San Diego, La Jolla, CA 92093, USA
14TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA Email<EMAIL_ADDRESS>0000-0001-9422-0095
0000-0002-1109-1919 0000-0002-8659-3729 0000-0002-7064-4309
0000-0003-4073-3236 0000-0001-5769-4945 0000-0002-6864-7762
0000-0003-0502-9235 0000-0002-4900-6628 0000-0003-2618-6504
0000-0003-0314-7027 0000-0003-4067-1434 0000-0002-1666-7067
0000-0003-3729-1684

(Accepted 2022. Received 2022; in original form 2022)

###### Abstract

Observations of local star-forming galaxies (SFGs) show a tight correlation
between their singly ionized carbon line luminosity ($L_{\rm[C_{II}]}$) and
star formation rate (SFR), suggesting that $L_{\rm[C_{II}]}$ may be a useful
SFR tracer for galaxies. Some other galaxy populations, however, are found to
have lower $L_{\rm[C_{II}]}{}/{}\rm SFR$ than the local SFGs, including the
infrared-luminous, starburst galaxies at low and high redshifts, as well as
some moderately star-forming galaxies at the epoch of re-ionization (EoR). The
origin of this ‘$\rm[C_{II}]$ deficit’ is unclear. In this work, we study the
$L_{\rm[C_{II}]}$-SFR relation of galaxies using a sample of $z=0-8$ galaxies
with $M_{*}\approx 10^{7}-5\times 10^{11}\,M_{\odot}$ extracted from
cosmological volume and zoom-in simulations from the Feedback in Realistic
Environments (FIRE) project. We find a simple analytic expression for
$L_{\rm[C_{II}]}$/SFR of galaxies in terms of the following parameters: mass
fraction of $\rm[C_{II}]$-emitting gas ($f_{\rm[C_{II}]}$), gas metallicity
($Z_{\rm gas}$), gas density ($n_{\rm gas}$) and gas depletion time ($t_{\rm
dep}{}={}M_{\rm gas}{}/{}\rm SFR$). We find two distinct physical regimes,
where $t_{\rm dep}$ ($Z_{\rm gas}$) is the main driver of the $\rm[C_{II}]$
deficit in $\rm H_{2}$-rich ($\rm H_{2}$-poor) galaxies. The observed
$\rm[C_{II}]$ deficit of IR-luminous galaxies and early EoR galaxies,
corresponding to the two different regimes, is due to short gas depletion time
and low gas metallicity, respectively. Our result indicates that $\rm[C_{II}]$
deficit is a common phenomenon of galaxies, and caution needs to be taken when
applying a constant $L_{\rm[C_{II}]}$-to-SFR conversion factor derived from
local SFGs to estimate cosmic SFR density at high redshifts and interpret data
from upcoming $\rm[C_{II}]$ line intensity mapping experiments.

###### keywords:

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

††pubyear: 2022††pagerange: $\rm[C_{II}]$ 158 $\rm\mu m$ emission as an
indicator of galaxy star formation rate–25

##  1 Introduction

The census of cosmic star formation from the present day to the highest
redshifts imposes a key constraint on galaxy evolution theory and physical
cosmology (see e.g., Madau & Dickinson, 2014; Dayal & Ferrara, 2018, and
references therein). The rest-frame ultra-violet (UV) luminosity ($L_{\rm
UV}$) of galaxies, tracing the young, massive stars, is a common star
formation rate (SFR) indicator of galaxies. However, a large fraction of the
UV light from galaxies in the Universe is absorbed by interstellar dust and
gets re-emitted as thermal radiation at far-infrared (far-IR) wavelength (e.g.
Fixsen et al., 1998; Takeuchi et al., 2005; Dole et al., 2006; Magnelli et
al., 2009; Gruppioni et al., 2013; Burgarella et al., 2013; Whitaker et al.,
2017; Salim & Narayanan, 2020). Therefore, an accurate estimate of the cosmic
SF history depends on a multi-wavelength, UV-to-millimetre (mm) analysis that
accounts for both the direct, unobscured stellar light and the dust thermal
emission of galaxies over cosmic time.

In practice, however, our capability of constraining the two components of
stellar radiation is largely imbalanced (e.g. Casey et al., 2018a). While the
rest-frame UV-based, unobscured component has been constrained to as early as
redshift $z{}\sim{}15$ (e.g. Bouwens et al., 2007, 2011, 2015, 2019; Oesch et
al., 2012, 2016; Oesch et al., 2018; Ellis et al., 2013; McLure et al., 2013;
Finkelstein et al., 2015; McLeod et al., 2015; Bowler et al., 2020;
Leethochawalit et al., 2022a, b; Naidu et al., 2022; Donnan et al., 2022;
Harikane et al., 2022) through deep imaging with the Hubble and the James Webb
Space Telescopes, the obscured component is still not well constrained beyond
$z{}\sim{}3$ due to the lack of statistically-representative, unbiased galaxy
samples in that regime (Casey et al., 2014, 2018a; Dayal & Ferrara, 2018;
Zavala et al., 2021). It is therefore important to have other SFR diagnostics
in addition to UV+IR for early galaxies (see e.g. Khusanova et al., 2021, and
references therein).

The $\rm 158\,\mu m$ (1900.5 GHz) fine structure transition
(${}^{2}P_{3/2}\rightarrow{}^{2}P_{1/2}$) of singly ionized carbon
($\rm[C_{II}]$) has been thought as a promising alternative SFR indicator,
particularly for high-$z$ galaxies (Hodge & da Cunha, 2020). It is a major
coolant of the neutral atomic gas of the interstellar medium (ISM) and the
strongest emission line of star-forming galaxies at rest-frame far-IR
wavelength (Carilli & Walter, 2013). The $\rm[C_{II}]$ line of galaxies is
usually not much affected by dust extinction.

To first order, a correlation between $L_{\rm[C_{II}]}$ and global SFR of
galaxies is expected. Much of the $\rm[C_{II}]$ emission of galaxy originates
from the neutral atomic gas regions (Hollenbach & Tielens, 1999; Wolfire et
al., 2003; Ferrara et al., 2019), where the far-UV (FUV) photons produced by
the young O and B-type stars heat the gas via the photoelectric (PE) effect on
small dust grains and polycyclic aromatic hydrocarbon (PAH) molecules (Tielens
& Hollenbach, 1985; Hollenbach et al., 1991; Weingartner & Draine, 2001a;
Helou et al., 2001). The photo-electrons ejected from the dust grains/PAH
molecules collisionally couple to and heat the gas. Since the PE heating rate
($\dot{E}_{\rm PE}$) traces galaxy SFR, and $L_{\rm[C_{II}]}$ balances
$\dot{E}_{\rm PE}$ given that $\rm[C_{II}]$ line is the dominant coolant in
those regions (assuming a thermal equilibrium), $L_{\rm[C_{II}]}$ should
therefore be correlated to SFR. Observations of local star-forming galaxies
(SFGs) have indeed found a linear correlation between $L_{\rm[C_{II}]}$ and
SFR over the broad SFR range of $\approx{}10^{-4}-10\,M_{\odot}\,\rm yr^{-1}$
(e.g. Stacey et al., 1991; Leech et al., 1999; Boselli et al., 2002; De Looze
et al., 2011; De Looze et al., 2014; Herrera-Camus et al., 2015). These
observations suggest that the $\rm[C_{II}]$ line can be a useful SFR indicator
for galaxies.

There is evidence, however, showing that this scaling relationship does not
hold on all occasions. For instance, observations find that local ultra-
luminous infrared galaxies (ULIRGs, galaxies having $L_{\rm
IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10^{12}\,L_{\odot}$) show a
significant lower $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ($\sim
L_{\rm[C_{II}]}{}/{}\rm SFR$) ratio than normal star-forming galaxies by up to
an order of magnitude (Malhotra et al., 1997, 2001; Luhman et al., 1998;
Luhman et al., 2003; Brauher et al., 2008; Farrah et al., 2013; Magdis et al.,
2014), the so-called ‘$\rm[C_{II}]$ deficit’ problem. This result was at first
revealed with the Infrared Space Observatory (ISO; Kessler et al., 1996) and
later confirmed by observations with the Herschel Space Observatory (hereafter
Herschel; Pilbratt et al., 2010) that has improved far-IR observing
capabilities. Subsequent observations with Herschel also show that the
$\rm[C_{II}]$ deficit extends to lower $L_{\rm IR}$ and that the
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio of galaxies exhibits a continuous
decrease with increasing $L_{\rm IR}$ at $L_{\rm IR}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}10^{11}\,L_{\odot}$ (e.g. Graciá-Carpio et al., 2011; Sargsyan
et al., 2012; Díaz-Santos et al., 2013; Cormier et al., 2015; Herrera-Camus et
al., 2015, 2018; Díaz-Santos et al., 2017; Hughes et al., 2017; Contursi et
al., 2017).

Studies have investigated the $L_{\rm[C_{II}]}$-SFR relation of galaxies at
higher redshifts (e.g. Stacey et al., 2010; Gullberg et al., 2015, 2018;
Brisbin et al., 2015; Spilker et al., 2016; Zanella et al., 2018; Cooke et
al., 2018; Rybak et al., 2019; McKinney et al., 2020). At $z{}\approx{}1-5$,
the selected galaxies are mostly uncovered by sub-mm surveys, which are
traditionally classified as ‘sub-millimetre-bright galaxies (SMGs111In the
literature, ‘SMGs’ typically refer to the galaxies detectable by single-dish
sub-mm telescopes, of which the observed sub-mm flux density is above
$\sim{}1$ mJy (Casey et al., 2014; Hodge & da Cunha, 2020).)’. These are
heavily dust-obscured systems having $L_{\rm IR}{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{12}\,L_{\odot}$ (corresponding to SFR${}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}100\,M_{\odot}\,\rm yr^{-1}$; Kennicutt 1998). In general, it
is found that $\rm[C_{II}]$ deficit persists at high $L_{\rm IR}$ at high
redshifts, although the high-$z$ populations appear to show larger scatter of
$L_{\rm[C_{II}]}{}/{}\rm SFR$ at given $L_{\rm IR}$ than the local ones.

The advent of the Atacama Large Millimetre/submillimetre Array (ALMA)
Telescope (e.g. Wootten & Thompson, 2009) has triggered particular interest in
searching for $\rm[C_{II}]$ emitters at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$,
and accumulating efforts have been made to constrain the $L_{\rm[C_{II}]}$-SFR
relation of galaxies at this epoch (e.g. Ouchi et al., 2013; Ota et al., 2014;
Maiolino et al., 2015; Capak et al., 2015; Willott et al., 2015b; Pentericci
et al., 2016; Matthee et al., 2017, 2019; Carniani et al., 2018a; Smit et al.,
2018; Schaerer et al., 2020; Fujimoto et al., 2021; Ferrara et al., 2022;
Schouws et al., 2022a). The ALMA observational programs are often designed to
target the Lyman-$\alpha$ emitters (LAEs), Lyman-break galaxies (LBGs) and the
quasar host galaxies (hereafter quasar hosts for simplicity) having pre-
determined redshift (Hodge & da Cunha, 2020). Though the earliest attempts
targeting the bright LAEs were mostly unsuccessful (e.g. Maiolino et al.,
2005; Ouchi et al., 2013; Ota et al., 2014; Inoue et al., 2016), follow-up
programs targeting the LBGs and quasar hosts generally have had much higher
success rate of $\rm[C_{II}]$ line detection. Overall, there have been
$>{}200$ galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ that have
confirmed detection of $\rm[C_{II}]$ line to date. While the quasar hosts are
typically very luminous and have substantial SFR (e.g. Bañados et al., 2016;
Decarli et al., 2018; Venemans et al., 2020), many of the selected LBGs/LAEs
at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}5$ are normal star-forming galaxies
having moderate SFR ($\approx 10\,M_{\odot}\,\rm yr^{-1}$). In particular, the
ALMA Large Program to INvestigate $\rm[C_{II}]$ at Early times (ALPINE) survey
(Le Fèvre et al., 2020; Béthermin et al., 2020; Faisst et al., 2020a) in
Cycle-5, targeting a sample of 118 star-forming galaxies at $z{}\approx{}5-6$,
has contributed more than a third ($\sim 75/200$) of the total number of
successful detections at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ (Schaerer et
al., 2020). More recently, the ALMA Reionization Era Bright Emission Line
Survey (REBELS; Bouwens et al., 2022) in Cycle-7 has targeted a sample of 40
UV-bright, star-forming galaxies at $z{}\approx{}7$ and confirmed
$\rm[C_{II}]$ line detection for 18 galaxies in their sample (Ferrara et al.,
2022).

Observations have drawn divergent conclusions on the $L_{\rm[C_{II}]}$-SFR
relation at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$. While some have
argued a clear $\rm[C_{II}]$ deficit of galaxies at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ with respect to the local normal SFGs (e.g. Ouchi et al.,
2013; Ota et al., 2014; Maiolino et al., 2015; Inoue et al., 2016; Knudsen et
al., 2016; Pentericci et al., 2016; Bradač et al., 2017; Ferrara et al., 2019;
Laporte et al., 2019; Carniani et al., 2020; Fujimoto et al., 2022), others
have argued that they follow the same linear scaling relation (e.g. Matthee et
al., 2017; Carniani et al., 2018a; Schaerer et al., 2020; Fujimoto et al.,
2021; Ferrara et al., 2022; Schouws et al., 2022a). It should be noted,
however, that the SFR estimates at such high redshifts can be highly
uncertain. Galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ typically have
very few reliable photometric data points in the dust thermal continuum that
are measured with ALMA (at band 6 or 7). A number of recent studies, both
observational (Capak et al., 2015; Bouwens et al., 2016; Casey et al., 2018a;
Faisst et al., 2020b) and theoretical (Liang et al., 2019, 2021; Ma et al.,
2019; Sommovigo et al., 2020, 2021), have pointed out that based on the ALMA
broad-band flux(es) alone, $L_{\rm IR}$ (and hence the obscured SFR) of
galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ is likely to be
poorly constrained due to the large variation in the shape of the spectral
energy distribution (SED) of their dust emission. The reported
(in)consistencies of the $L_{\rm[C_{II}]}$-SFR relation at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ with the local SFGs by the observations therefore need to
be more carefully assessed.

Much effort has been made to model $\rm[C_{II}]$ emission of galaxies and
explain the origins of the observed $\rm[C_{II}]$ deficit over the last two
decades. A broad variety of different methods are used by different studies,
including pure analytic approaches (e.g. Muñoz & Oh, 2016; Ferrara et al.,
2019), numerical models of idealized gas clouds (e.g. Abel et al., 2009;
Narayanan & Krumholz, 2017), semi-analytic galaxy models (SAMs, e.g. Popping
et al., 2014, 2016, 2019; Lagache et al., 2018; Yang et al., 2021, 2022) and
cosmological hydrodynamic galaxy simulations (e.g. Vallini et al., 2013;
Vallini et al., 2015; Olsen et al., 2015, 2017; Pallottini et al., 2017;
Pallottini et al., 2019; Katz et al., 2019; Leung et al., 2020; Lupi et al.,
2020; Lupi & Bovino, 2020; Kannan et al., 2022; Richings et al., 2022; Bisbas
et al., 2022). A pure analytic approach and/or a simplified cloud model can
capture the key physical mechanisms that determine $L_{\rm[C_{II}]}$ of
galaxies and provide useful insights at low computational cost, but does not
provide the necessary galaxy statistics. SAMs can produce statistically
significant galaxy samples probing a very wide dynamic range (in stellar mass,
SFR, redshift and etc.) and are computationally efficient (Somerville & Davé,
2015), but they do not provide any information of structures on the sub-
galactic scales. Hydrodynamic simulations, in contrast, can calculate the
detailed sub-galactic structures and thus provide more accurate prediction for
the $\rm[C_{II}]$ emission properties of galaxies, at the cost of more
computational expense.

Different explanations for the $\rm[C_{II}]$ deficit in the high $L_{\rm IR}$
regime have been proposed by the theory groups (see also e.g. Casey et al.,
2014; Narayanan & Krumholz, 2017, for a summary). For instance, some studies
argue that the deficit is due to a strong UV radiative intensity ($G$) in the
IR luminous galaxies (e.g. Malhotra et al., 1997; Luhman et al., 1998; Genzel
& Cesarsky, 2000; Helou et al., 2001; Luhman et al., 2003; Abel et al., 2009;
Stacey et al., 2010; Graciá-Carpio et al., 2011; Lagache et al., 2018). This
can have two important effects on the thermal balance of
$\rm[C_{II}]$-emitting gas. First of all, a high $G$ leads to large positive
grain charges, thereby reducing the kinetic energy of the ejected photo-
electrons and hence the rate of PE heating ($\dot{E}_{\rm PE}$) of gas
(Tielens & Hollenbach, 1985; Kaufman et al., 1999). As a result, $\rm[C_{II}]$
cooling rate drops. Besides, $\rm H_{II}$ regions in those galaxies may become
‘dust bounded’ rather than ‘ionization bounded’ (e.g. Bottorff et al. 1998;
Abel et al. 2009; see also Ferrara et al. 2019). In this scenario, most of the
UV radiation from young stars is absorbed by dust in the $\rm H_{II}$ regions,
leading to both an excess of IR emission in the $\rm H_{II}$ regions and a
reduced $\dot{E}_{\rm PE}$ (and hence $L_{\rm[C_{II}]}$) in gas outside the
$\rm H_{II}$ regions due to a starvation of UV photons there.

Alternatively, Narayanan & Krumholz (2017) suggest that a high gas density can
lead to a $\rm[C_{II}]$ deficit of galaxy in addition to having a high $G$.
Using a stratified gas cloud model, the authors demonstrate that with
increasing gas density, a larger fraction of carbon in gas turns into neutral
(i.e. in CO and $\rm C_{I}$) and $L_{\rm[C_{II}]}$ decreases due to a reduced
mass fraction of $\rm[C_{II}]$-emitting gas.

Apart from these studies, Muñoz & Oh (2016) posit an analytic model where
$\rm[C_{II}]$ deficit is due to thermal saturation of the upper fine structure
transition state (${}^{2}P_{3/2}$) of $\rm C_{II}$ ions222Throughout this
paper, we use ‘$\rm[C_{II}]$’ when referring to the observable emission line,
and ‘$\rm C_{II}$’ when discussing ionized carbon under the context of
chemical abundances of gas.. At above 91.8 K (note: $T_{*}{}={}91.8$ K is the
equivalent temperature of the $\rm[C_{II}]$ transition), $L_{\rm[C_{II}]}$
does not increase much with gas kinetic temperature and this has been
suggested to be the reason for $L_{\rm[C_{II}]}$ not increasing much with SFR
at high $L_{\rm IR}$ ($\rm\sim{}SFR$). Note, however, that the Muñoz & Oh
(2016) model assumes that the bulk of the $\rm[C_{II}]$ emission of galaxies
originates from the gas having density in excess of the critical density for
the $\rm[C_{II}]$ transition (Goldsmith et al., 2012).

With the recent success of the ALMA programs in searching for
$\rm[C_{II}]$-emitters, there has been an increasing amount of effort to
predict $\rm[C_{II}]$ emission properties of galaxies at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ by coupling cosmological hydrodynamic simulations or SAMs
with photo-ionization codes (e.g. CLOUDY, Ferland et al. 1998; Ferland et al.
2013; DESPOTIC, Krumholz 2014; RADMC-3D, Dullemond et al. 2012). The predicted
$L_{\rm[C_{II}]}$-SFR relation for galaxies, however, shows non-trivial
discrepancy between different groups in both normalization and slope (see
summary in e.g. Katz et al., 2019; Leung et al., 2020), which can be ascribed
to the differences in the simulation methodology and $\rm[C_{II}]$ modelling
techniques adopted by the different groups. Despite the discrepancy, many have
predicted a $\rm[C_{II}]$ deficit of galaxies at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ with respect to the local normal star-forming galaxies. For
instance, Lagache et al. (2018) couple a sample of $\sim{}20$ K SAM galaxies
at $4{}\leq{}z{}\leq{}8$ with CLOUDY and report a $\rm[C_{II}]$ deficit of
$>{}0.5$ dex and a trend of decreasing normalization of the relation with
redshift. Olsen et al. (2017) post-process 30 star-forming galaxies at
$z{}={}6$ extracted from the MUFASA ‘zoom-in’ simulations (Davé et al., 2016)
using CLOUDY and predict a $\rm[C_{II}]$ deficit of about one decade. A
similarly strong $\rm[C_{II}]$ deficit is reported by Pallottini et al.
(2017); Pallottini et al. (2019) using the SERRA ‘zoom-in’ simulations that
include more sophisticated chemical networks. More recently, Kannan et al.
(2022) predict an even more prominent $\rm[C_{II}]$ deficit at $z{}\geq{}5$
than the above-mentioned earlier studies, especially at low SFR, using a
galaxy sample produced by the THESAN ‘zoom-in’ suite, which includes the
Illustris-TNG galaxy formation model (Pillepich et al., 2018a, b).

It has been generally thought that gas metallicity ($Z_{\rm gas}$) is the key
factor in determining the $\rm[C_{II}]$ luminosity of the early galaxies (e.g.
Vallini et al., 2015; Olsen et al., 2017; Ferrara et al., 2019) since
$\rm[C_{II}]$ emissivity is linearly scaled with $Z_{\rm gas}$. The early work
by Vallini et al. (2015) shows that the $L_{\rm[C_{II}]}$-SFR relation of EoR
galaxies depends sensitively on $Z_{\rm gas}$, and the significant
$\rm[C_{II}]$ deficit found with the LAEs at $z{}\approx{}5-7$, such as Himiko
(Ouchi et al., 2013; Ota et al., 2014) and IOK-1 (Ota et al., 2014), can be
well accounted for by assigning a very low gas metallicity ($Z_{\rm
gas}{}<{}0.05\,Z_{\odot}$) to the simulated galaxy in an ad hoc manner. The
$\rm[C_{II}]$ deficit of galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$
commonly found in the recent simulations, as mentioned above, is likely due to
the much lower $Z_{\rm gas}$ of the early galaxies than the $z{}={}0$ ones
predicted by these simulations. Observationally, however, direct measurement
of $Z_{\rm gas}$ at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ is still very
challenging, though some preliminary attempts have been made recently (e.g.
Rigopoulou et al., 2018; Curti et al., 2022; Heintz et al., 2022a, b; Rhoads
et al., 2022; Schaerer et al., 2022; Trump et al., 2022).

A few recent studies have predicted $\rm[C_{II}]$ emission of galaxies at
lower redshifts using simulations. For instance, Popping et al. (2019) and
Yang et al. (2021) predict the $L_{\rm[C_{II}]}$-SFR relation for the catalog
derived from the ‘Santa Cruz’ semi-analytic models (Somerville & Primack,
1999; Somerville et al., 2015) using DESPOTIC. Their result is in good
agreement with the observational data at $z{}\approx{}2$, except that at high
SFR (i.e. ${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}10\,M_{\odot}\,\rm
yr^{-1}$), they produce a noticeably weaker $\rm[C_{II}]$ deficit than is
observed. More recently, Richings et al. (2022) ran a set of hydrodynamic
simulations of isolated (dwarf and Milky Way-mass) galaxies implemented with
the CHIMES non-equilibrium chemistry module (Richings et al., 2014a, b)
(including a dust-depletion model) and predict the $\rm[C_{II}]$ emission of
their galaxy sample using RADMC-3D. Despite having a small sample size, the
predicted $L_{\rm[C_{II}]}$ of their galaxies appears to be in agreement with
the observational result of local galaxies (e.g. De Looze et al., 2011; De
Looze et al., 2014; Herrera-Camus et al., 2015) at similar SFR (see also
another recent work by Bisbas et al. 2022 using isolated dwarf simulations).

Apart from these studies, there has been limited effort to predict the
$L_{\rm[C_{II}]}$-SFR relation of galaxies at $z{}={}0-5$ using statistically
representative galaxy samples and compare the result to the fruitful
observational data in this regime. In particular, the origin of the
$\rm[C_{II}]$ deficit of the IR-luminous galaxies has not yet been studied in
detail using cosmological hydrodynamic simulations. This is largely because
producing a statistically representative sample in this regime with well-
resolved ISM is computationally demanding, which is possible only for a few
large simulation consortiums. It is, however, of critical importance that a
robust $\rm[C_{II}]$ model should be able to simultaneously reproduce the data
of different galaxy populations over the entire SFR and redshift ranges.

In this study, we conduct a comprehensive analysis of the galaxy
$L_{\rm[C_{II}]}$-SFR relation using a simulated sample spanning an
unprecedentedly broad redshift range of $z{}={}0-8$ extracted from the
MassiveFIRE (Feldmann et al., 2016; Feldmann et al., 2017; Anglés-Alcázar et
al., 2017) and FIREbox (Feldmann et al., 2022) cosmological hydrodynamic
simulations from the Feedback in Realistic Environments (FIRE) project333FIRE
project website: http://fire.northwestern.edu (Hopkins et al., 2014, 2018,
2022). The sample covers a very broad range of galaxy stellar mass and SFR,
allowing us to make direct comparison with the observational data in different
regimes. In particular, the sample includes local normal SFGs (having
SFR${}\approx{}0.1-10\,M_{\odot}\,\rm yr^{-1}$) that can be compared with the
observations where a linear $L_{\rm[C_{II}]}$-SFR correlation has been found
by the observations. It also includes IR-luminous ($L_{\rm
IR}>10^{11}\,L_{\odot}$) galaxies at $z{}={}0-5$ that are candidates for
(U)LIRGs and SMGs, where observations have shown to have $\rm[C_{II}]$
deficit. Moreover, the sample includes early galaxies at above $z{}={}5$
spanning a broad SFR range. Many of these galaxies have similar mass and SFR
to the samples of the ALPINE and REBELS projects and therefore can be used to
provide useful interpretations for a variety of their recent observational
results (e.g. Fujimoto et al., 2020; Ginolfi et al., 2020; Schaerer et al.,
2020; Fudamoto et al., 2021, 2022; Ferrara et al., 2022; Sommovigo et al.,
2022).

The main goal of this work is to predict the $L_{\rm[C_{II}]}$-SFR relation
for the FIRE galaxy sample (spanning $z{}={}0-8$ and ${\rm
SFR}{}\approx{}0.1-10^{3}\,M_{\odot}\,\rm yr^{-1}$) and to understand what
physical parameters of galaxies determine their overall $L_{\rm[C_{II}]}$-to-
SFR ratio. This will then help us find the origin of the observed
$\rm[C_{II}]$ deficit of galaxies at both high $L_{\rm IR}$ and high
redshifts.

Note that the results from this work will be useful for interpreting the data
of several upcoming $\rm[C_{II}]$ line intensity mapping (LIM) experiments
(see e.g. Kovetz et al., 2017, 2019; Bernal & Kovetz, 2022, and references
therein), such as TIME444https://cosmology.caltech.edu/projects/TIME (Crites
et al., 2014; Sun et al., 2021), CCAT-prime555http://www.ccatobservatory.org
(CCAT-Prime collaboration et al., 2021), CONCERTO666https://www.apex-
telescope.org/ns/concerto/ (CONCERTO Collaboration et al., 2020; Gkogkou et
al., 2022) and EXCLAIM (Ade et al., 2020). The LIM experiments have been
designed to measure the emission from spectral lines originating from galaxies
at all luminosities, including the ones that cannot be resolved by the current
surveys (e.g. with ALMA). The experiments that will target $\rm[C_{II}]$
emission, in particular, will be useful for constraining the cosmic star-
formation history (see e.g. Gong et al., 2012; Silva et al., 2015; Serra et
al., 2016; Fonseca et al., 2017; Padmanabhan, 2019; Yue & Ferrara, 2019; Chung
et al., 2020; Padmanabhan et al., 2022; Sun et al., 2022; Karoumpis_2022). It
is, however, not yet certain whether the $\rm[C_{II}]$ line always acts as a
reliable SFR tracer for galaxies of all types and at all redshifts.

This paper is structured as follows. We describe in Section 2 the simulation
methodology and in Section 3, the method used to simulate $\rm[C_{II}]$
emission. In Section 4, we compare the predicted $L_{\rm[C_{II}]}$-SFR
relation of the FIRE galaxy sample with the observational data at different
redshifts. In Section 5, we investigate the origin of the tight
$L_{\rm[C_{II}]}$-SFR linear scaling relation of normal star-forming galaxies
at $z{}={}0$ and the causes of the $\rm[C_{II}]$ deficit of galaxies. We
discuss our results in Section 6 and finally summarize and conclude this study
in Section 7. Throughout this paper, we adopt the cosmological parameters of
the Planck 2015 Cosmology (Planck Collaboration et al., 2016), specifically
$\Omega_{\rm m}=0.309$, $\Omega_{\Lambda}=0.691$, $\Omega_{\rm b}=0.049$,
$\sigma_{8}=0.816$, and $H_{0}=67.74\,\rm km\,s^{-1}\,Mpc^{-1}$.

##  2 Simulation methodology

In this section, we introduce the simulation suites (FIREbox and MassiveFIRE)
from which we extract the galaxy sample used for this study.

###  2.1 Simulation set-up and galaxy catalogue

We adopt a sample that spans the wide redshift range $z{}={}0-8$, stellar mass
($M_{*}$) range $M_{*}{}\approx{}10^{7}-5\times 10^{11}\,M_{\odot}$ and SFR
range ${\rm SFR}{}\approx{}0.1-10^{3}\,M_{\odot}\,\rm yr^{-1}$. The sample
consists primarily of galaxies at $z{}={}0-8$ produced by FIREbox (Feldmann et
al., 2022), the new-generation simulation suite of FIRE run with full
cosmological volume boxes. It is supplemented by a number of high-$z$
($z{}={}1-8$) massive galaxies ($M_{*}{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{10}\,M_{\odot}$) extracted from the ‘zoom-in’ suite,
MassiveFIRE (Feldmann et al., 2016; Feldmann et al., 2017), rerun with FIRE-2
physics (Anglés-Alcázar et al., 2017; Çatmabacak et al., 2022; Bassini et al.,
2022). Many of the MassiveFIRE galaxies have the $L_{\rm IR}$ close to that of
the SMGs (Liang et al., 2018; Cochrane et al., 2019) that are used by the
observational studies on the $L_{\rm[C_{II}]}$-SFR relation at high redshifts.
All simulations used for this study are run with the same FIRE-2 physics and
numerics (Hopkins et al., 2018).

##### FIREbox simulations

Figure 1: Histograms of the stellar mass distribution of the FIREbox sample at
different redshifts. Violet, green, magenta, blue, red, yellow and cyan
histograms correspond to $z{}={}8$, $z{}={}6$, $z{}={}4$, $z{}={}3$,
$z{}={}2$, $z{}={}1$ and $z{}={}0$, respectively. For each redshift, the
unfilled histograms indicate the result of the entire galaxy sample, whereas
the filled histograms indicate specifically the result of the galaxies having
SFR${}\geq{}10\,M_{\odot}\,\rm yr^{-1}$ (corresponding to $L_{\rm
IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10^{11}\,L_{\odot}$ based on the
Kennicutt 1998 relation. Note: $\rm[C_{II}]$ deficit is observed at $L_{\rm
IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10^{11}\,L_{\odot}$.). For clarity of
presentation, we separately show the result of the 7 snapshots in 3 separate
panels (top panel for $z{}={}8$ and $z{}={}6$, middle panel for $z{}={}3$ and
$z{}={}4$ and bottom panel for $z{}={}0$, $z{}={}1$ and $z{}={}2$).

FIREbox (Feldmann et al., 2022) is a new-generation simulation suite using
FIRE physics. Different from all previous simulations of FIRE, FIREbox
simulates full cosmological volumes instead of using ‘zoom-in’ set-up to study
galaxy evolution. FIREbox simulations are run in cubic boxes with periodic
boundary conditions, and with initial conditions at redshift $z{}={}120$
generated using the MUSIC (Multi-Scale Initial Conditions) code (Hahn & Abel,
2011). The simulations use the Planck 2015 Cosmology (Planck Collaboration et
al., 2016).

All FIREbox simulations use the same initial conditions and cosmology but
differ in numerical resolution. For this study, we extract galaxies from the
fiducial FIREbox hydrodynamic simulation, which is run with a box length of
$15\,h^{-1}$ cMpc and with the following number of dark matter (DM) and
baryonic particles: $N_{\rm DM}{}={}1024^{3}$ and $N_{\rm b}{}={}1024^{3}$.
The mass resolution of DM and baryon particles are $m_{\rm DM}{}={}3.3\times
10^{5}$ and $m_{\rm b}{}={}6.3\times 10^{4}\,M_{\odot}$. The gravitational
softening lengths are kept fixed in proper (comoving) coordinates at
$z{}\leq{}9$ ($z{}\geq{}9$) and are set to $h_{\rm DM}{}={}80$ pc for DM
particles and $h_{*}{}={}12$ pc for star particles. The softening length for
gas particles ($h_{\rm gas}$) is fully adaptive and is set equal to their
kernel smoothing length down to a minimum of 1.5 proper pc, which is reached
in the densest parts of the ISM. FIREbox is evolved down to $z{}={}0$.

We identify galaxies in different snapshots of the FIREbox simulation using
the AMIGA Halo Finder777Code available at: popia.ft.uam.es/AHF/Download.html
(AHF; Gill et al., 2004; Knollmann & Knebe, 2009). We use the galaxies
extracted from 7 snapshots corresponding to redshift $z{}={}0$, 1, 2, 3, 4, 6
and 8. For each snapshot, we include the central galaxy of the 30 most massive
halos identified by AHF. To enlarge our sample, we also include the central
galaxy of a number of additional, randomly chosen halos having ${\rm
log}\,(M_{\rm vir}/M_{\odot}){}>{}10$. We show in Fig. 1 the histograms of the
$M_{*}$ distribution of the selected FIREbox galaxies at different redshifts.

The number of galaxies selected at $z{}={}0$, 1, 2, 3, 4, 6, and 8 are 113,
84, 80, 75, 64, 61 and 30, respectively. As is shown in Fig. 1, all but a few
selected galaxies have stellar mass greater than $10^{7}\,M_{\odot}$
(corresponding to $\sim 160$ times of the mass resolution). The most massive
galaxy of the FIREbox sample has $M_{*}{}={}4.8\times 10^{11}\,M_{\odot}$ (at
$z{}={}0$).

In the same figure, we also show the $M_{*}$ distribution of the galaxies
having ${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}10\,M_{\odot}\,\rm
yr^{-1}$ (filled histograms). These galaxies have $L_{\rm
IR}{}\geq{}10^{11}\,L_{\odot}$, the regime where a $\rm[C_{II}]$ deficit is
observed (see Section 3). They apparently are more massive than the galaxies
having SFR${}<{}10\,M_{\odot}\,\rm yr^{-1}$. In our catalogue, we find most
galaxies with ${\rm SFR}{}\geq{}10\,M_{\odot}\,\rm yr^{-1}$ at $z{}={}2$ (red
histogram, $N{}={}29$) and $z{}={}3$ (blue histogram, $N{}={}28$). These
redshifts are at the ‘cosmic noon’, where massive galaxies start to form and
they are more gas-rich and actively star-forming than galaxies at lower
redshifts.

Since the FIREbox simulation is run with a volume of $(15\,h^{-1}\rm
cMpc)^{3}$, it does not produce enough galaxies at high redshifts that are as
massive and luminous as the galaxy samples selected by the observational
studies. We therefore supplement our sample with a handful of more massive
galaxies ($M_{*}{}\approx{}10^{9}-5\times 10^{11}\,M_{\odot}$) extracted from
the MassiveFIRE ‘zoom-in’ simulations (see below).

##### MassiveFIRE simulations

Table 1: List of MassiveFIRE simulations used for this work.

Sim ID $\dagger$ | Box Size | $z_{\rm final}$ |  $M_{\rm vir}$ $\ddagger$ |  $M_{*}$ ($M_{\odot}$)
---|---|---|---|---
| ($h^{-1}$ Mpc) |  | ($10^{12}\,M_{\odot}$) | $z=1$ | $z=2$ | $z=3$ | $z=4$ | $z=6$ | $z=8$
A1 | 100 | 1 | 2.4 |  $5.4\times 10^{11}$ |  $5.1\times 10^{10}$ |  $9.6\times 10^{9}$ |  $1.2\times 10^{9}$ | / | /
A2 | 100 | 1 | 3.0 |  $4.1\times 10^{11}$ |  $2.9\times 10^{11}$ |  $1.3\times 10^{11}$ |  $2.7\times 10^{10}$ | / | /
A4 | 100 | 1 | 2.9 |  $2.3\times 10^{11}$ |  $1.3\times 10^{11}$ |  $2.2\times 10^{10}$ |  $6.5\times 10^{9}$ | / | /
A8 | 100 | 1 | 3.6 |  $2.8\times 10^{11}$ |  $1.8\times 10^{11}$ |  $9.8\times 10^{10}$ |  $5.1\times 10^{10}$ | / | /
D3 | 400 | 6 | 4.5 | / | / | / | / |  $3.9\times 10^{11}$ |  $7.0\times 10^{10}$
D7$\mathsection$ | 400 | 6 | 2.5 | / | / | / | / | / |  $5.8\times 10^{10}$
D9 | 400 | 6 | 1.0 | / | / | / | / |  $3.9\times 10^{10}$ |  $1.3\times 10^{9}$
E1 | 762 | 6 | 6.8 | / | / | / | / |  $1.6\times 10^{10}$ |  $3.2\times 10^{9}$
E2 | 762 | 6 | 6.5 | / | / | / | / |  $7.2\times 10^{9}$ |  $5.3\times 10^{9}$
E3 | 762 | 6 | 6.1 | / | / | / | / |  $8.6\times 10^{9}$ |  $2.7\times 10^{9}$

  * $\dagger$

The A (D and E) Series of MassiveFIRE were published in Anglés-Alcázar et al.
(2017) (Çatmabacak et al. 2022) for the first time.

  * $\ddagger$

Virial mass at $z_{\rm final}$.

  * $\mathsection$

The HR simulation of D7 has been run only down to $z=7.2$.

MassiveFIRE (Feldmann et al., 2016; Feldmann et al., 2017) is a set of
simulations of massive galaxies at high redshifts using the ‘zoom-in’ method.
A number of low-resolution (LR) DM-only simulations were run with the initial
conditions generated using the MUSIC code within periodic boxes. From the
outputs of these LR DM-runs, we then select a number of model haloes to re-
simulate at much higher resolution and with baryons included (HR runs). The
selected haloes have a variety of masses, accretion history, and environmental
over-densities.

For this study, we use the galaxies produced by 10 MassiveFIRE simulations,
which are from the A (Anglés-Alcázar et al., 2017), D and E Series (Çatmabacak
et al., 2022; Bassini et al., 2022). The A, D and E Series were run in the
periodic boxes with size of $(100\,h^{-1}\,\rm Mpc)^{3}$, $(400\,h^{-1}\,\rm
Mpc)^{3}$ and $(762\,h^{-1}\,\rm Mpc)^{3}$, respectively. The model haloes of
the A Series are selected from the snapshot of $z_{\rm final}{}={}1$, those of
the D and E Series are selected from the snapshot of $z_{\rm final}{}={}6$.
All the HR runs were run down to $z_{\rm final}$ except D7, where the HR run
is evolved to only $z{}={}7.2$. This is because part of the ISM in D7 became
too compact so that the gas particles with the highest densities were evolved
at extremely small time-steps and it became infeasible to run the simulation
down to the target redshift.

Initial conditions for the HR runs are set up using a convex hull surrounding
all particles within $3R_{\rm vir}$ at $z_{\rm final}$ of the chosen halo
defining the Lagrangian HR region following the method of Hahn & Abel (2011).
The mass resolutions and force softening lengths of the HR runs are similar to
those of the FIREbox simulation. Specifically, $m_{\rm DM}$ and $m_{\rm b}$
are set to $1.9\times 10^{5}\,M_{\odot}$, $3.6\times 10^{4}\,M_{\odot}$,
respectively. Both $h_{\rm DM}$ and $h_{*}$ are fixed in proper (comoving)
coordinates at $z{}\leq{}9$ ($z{}\geq{}9$) and are set equal to 57 pc and 7
pc, respectively. $h_{\rm gas}$ is set equal to the smoothing length of the
gas particles down to a minimum of 0.7 proper pc.

We include the central galaxy of the chosen haloes at $z_{\rm final}$ except
for that of the D7 run. In addition, we also include the most massive
progenitors (MMPs) of the central galaxies at higher redshifts. Specifically,
for the 4 A Series runs, we include the MMPs at $z{}={}2$, $z{}={}3$ and
$z{}={}4$, while for the D and E Series, we include the MMPs at $z{}={}8$. The
galaxies are identified in the simulation snapshots using AHF (Gill et al.,
2004; Knollmann & Knebe, 2009). In Table 1, we summarize the
information888Physical properties, including e.g. $M_{*}$ , SFR, $L_{\rm IR}$
and $L_{\rm[C_{II}]}$, of the FIRE galaxies reported in this paper are
estimated using a radial kernel of $0.1R_{\rm vir}$ around the DM halo centre,
i.e. the maximum density centre provided by AHF. of the 10 MassiveFIRE
simulations used for this study.

Both the MassiveFIRE and FIREbox simulations used in this work are run using
the N-body+hydrodynamics code GIZMO (FIRE-2 version) in the Meshless-Finite-
Mass (MFM) mode (Hopkins et al., 2018). The simulations incorporate various
gas cooling processes (free-free, photoionization/recombination, Compton,
photoelectric, metal-line, molecular and fine structure processes) and a
uniform UV background following the FG09 prescription (Faucher-Giguère et al.,
2009), Star formation occurs in dense, self-gravitating and self-shielding
molecular gas based on a sink-particle prescription. The simulations
explicitly incorporate several different stellar feedback channels (but not
feedback from supermassive black holes) including 1) local and long-range
momentum flux from radiative pressure, 2) energy, momentum, mass and metal
injection from supernovae (Types Ia and II), 3) stellar mass loss (both OB and
AGB stars) and 4) photo-ionization and photo-electric heating processes. We
refer the reader to Hopkins et al. (2014, 2018) for details of the star
formation and feedback prescriptions of FIRE.

FIRE has demonstrated success at reproducing a variety of key galaxy
properties that are relevant to this work, such as the stellar-to-halo mass
relation (Hopkins et al., 2014; Feldmann et al., 2017), the specific SFR
(sSFR) of galaxies at the cosmic noon ($z{}\sim{}2$) (Hopkins et al., 2014;
Feldmann et al., 2016; Sparre et al., 2017; Feldmann et al., 2022), the galaxy
molecular (atomic) hydrogen gas mass and stellar mass relations at $z{}={}0$
(Feldmann et al., 2022), the gas-phase and stellar mass-metallicity relation
at $z{}={}0-2$ (Ma et al., 2016; Feldmann et al., 2022), the observational
effective dust temperatures at $z{}={}2-4$ (Liang et al., 2019) as well as the
UV luminosity functions and UV-based cosmic star formation rate density
(CSFRD) at $z{}>{}5$ (Ma et al., 2019).

##  3 Simulating observational properties

In this section, we describe the method used to predict the observational
properties for the FIRE galaxy sample, which we compare to the observational
data. In Section 3.1, we describe our $\rm[C_{II}]$ emission model. In Section
3.2, we describe the prescription for the dust RT modelling of the FIRE
galaxies using SKIRT code, based on which we derive the multi-wavelength SED
and the distribution of the interstellar radiation field (ISRF) for the
galaxies. The ISRF distribution is essential for predicting the $\rm[C_{II}]$
emission properties of the galaxies.

###  3.1 Predicting $\rm[C_{II}]$ emission using CLOUDY

We predict the $\rm[C_{\rm II}]$ line luminosity for the FIRE sample using the
spectral synthesis code CLOUDY version 17.01 (Ferland et al., 2017). CLOUDY is
a plasma simulation code designed to simulate the ionization, level
populations, molecular state and thermal state of gas over a wide range of
density and temperature in different astrophysical environments (e.g. black
hole accretion disks, PDRs, molecular clouds, etc). It solves for the
ionization structure for all stages of ionization for the lightest 30 elements
(Abel et al., 2008).

We treat each gas particle of the galaxies as an idealized spherical uniform
‘gas cloud’. The $\rm[C_{II}]$ luminosity of each ‘cloud’ is calculated based
on its physical conditions, including ‘cloud’ (or gas particle) mass ($M_{\rm
cl}$), gas density999‘Gas density’ here refers to H nuclei number density of
gas, $n_{\rm H}\equiv X(\rho_{\rm gas}/m_{\rm H})$, where $X$ and $m_{\rm H}$
represent the mass fraction of hydrogen in gas and the proton mass,
respectively. CLOUDY uses $n_{\rm H}$ as input instead of mass density
$\rho_{\rm gas}$. In this paper, we constantly use ‘gas density’ to refer to
$n_{\rm H}$ for simplicity. ($n_{\rm H}$), gas metallicity ($Z_{\rm gas}$),
gas turbulent velocity dispersion ($\sigma$) and local UV ISRF strength
($G$101010Conventionally, $G$ is used to denote the mean ISRF in the Habing
band ($6.0-13.6$ eV). It is indicated in units of $G_{0}=1.6\times
10^{-3}\,\rm erg\,s^{-1}\,cm^{-2}$, the observed value in the solar
neighbourhood (Habing, 1968).). $M_{\rm cl}$, $n_{\rm H}$, $Z_{\rm gas}$ of
each ‘cloud’ are known directly from the FIRE simulations. $\sigma$ is the
mass-weighted standard deviation of the velocities in gas at the location of
the ‘cloud’, which is calculated in post-processing. Finally, $G$ at the
location of each ‘cloud’ in the galaxy is calculated using the dust RT code
SKIRT (Baes et al., 2011; Baes & Camps, 2015; Camps & Baes, 2015) in post-
processing (see Section 3.2 for the details).

We calculate the $\rm[C_{II}]$ luminosity for each ‘cloud’
($L_{\rm[C_{II}],\,cl}$) by integrating the $\rm[C_{II}]$ line cooling rate,
$\Lambda_{\rm[C_{II}]}$ ($\rm erg\,s^{-1}\,cm^{-3}$; see Appendix A for its
analytic expression), obtained from the output of the CLOUDY simulations, over
the volume of the cloud:

$L_{\rm[C_{II}],\,cl}=4\pi\int^{R_{\rm
cl}}_{0}\Lambda_{\rm[C_{II}]}(x)\;x^{2}{\rm d}x,$ (1)

where

$R_{\rm cl}=\left(\frac{3}{4\pi}\right)^{1/3}\left(\frac{M_{\rm cl}}{\mu_{\rm
H}m_{\rm H}n_{\rm H}}\right)^{1/3}$ (2)

represents the radius of the cloud111111Note that we do not derive
$L_{\rm[C_{II}],\,cl}$ using the ‘emergent intensity’ ($I_{\rm em}$, with
physical unit $\rm erg\,s^{-1}\,cm^{-2}$) output by CLOUDY because $I_{\rm
em}$ is calculated for a plane-parallel geometry instead of a spherical
geometry. The conversion factor between the two geometries is not simply a
constant but depends on the profile of $\rm[C_{II}]$ emissivity (Olsen et al.,
2017, 2018). and $\mu_{\rm H}$ in equation (2) represents the mean molecular
weight of the gas. The $\rm[C_{II}]$ luminosity of the galaxy
($L_{\rm[C_{II}]}$) is then derived by summing over $L_{\rm[C_{II}],\,cl}$ of
all gas clouds calculated using equation (1). We treat the $\rm[C_{II}]$
emission of our galaxy sample as being optically thin.

In practice, to run CLOUDY simulations for every gas particle for the whole
FIRE sample ($>{}400$ galaxies in total) is computationally formidable: a
CLOUDY simulation is typically completed (i.e. when iterative convergence is
reached) in $0.1{}-{}0.5$ CPU hour, depending on the gas column density, and
hence to analyze one single galaxy snapshot that contains $\sim{}1$ million
gas particles would cost $100{}-{}500$ K CPU hours in total. We therefore use
a lookup-table method similar to the previous studies (e.g. Vallini et al.,
2015, 2018, 2021; Katz et al., 2017; Katz et al., 2019; Olsen et al., 2017;
Lagache et al., 2018; Li et al., 2018; Pallottini et al., 2019; Leung et al.,
2020; Lupi et al., 2020; Yang et al., 2021; Lupi & Bovino, 2020).
Specifically, for each of the 7 snapshots, we build a grid of CLOUDY models
that covers a gas density range $-1{}<{}{\rm log}\,(n_{\rm H}/{\rm
cm^{-3}}){}<{}5$, a gas metallicity range $-2{}<{}{\rm log}\,(Z_{\rm
gas}/Z_{\odot}){}<{}0.8$, a turbulent velocity dispersion range $0{}<{}{\rm
log}\,(\sigma/\rm km\,s^{-1}){}<{}2.4$ and a UV ISRF range $-1{}<{}{\rm
log}\,(G/G_{0}){}<{}4$. The grid spacing is set 0.5 dex for $n_{\rm H}$ and
$G$, and 0.4 dex for $Z_{\rm gas}$ and $\sigma$. In total, the default look-up
table that we use for calculating the $\rm[C_{II}]$ luminosity of our galaxy
sample consists of 8,008 ($13{}\times 8{}\times{}7{}\times{}11$) models for
each redshift. We include the CMB background in the CLOUDY simulations for
each redshift and the predicted $\rm[C_{II}]$ luminosity is corrected for the
CMB attenuation effect (da Cunha et al., 2013). Cosmic-ray (CR) hydrogen
ionization rate in these models is fixed to the fiducial value of
$2{}\times{}10^{-16}\,\rm s^{-1}$, the observed value in the Milky Way
(Indriolo et al., 2007; Indriolo & McCall, 2012; Neufeld & Wolfire, 2017). We
assume a constant dust-to-metal mass ratio $\delta_{\rm dzr}{}={}0.4$ (Dwek,
1998; Draine et al., 2007; Watson, 2011; Li et al., 2019) and adopt the
default interstellar medium metal abundances (abundance ISM) stored in CLOUDY.
The simulations are run till sufficiently large distance from the surface of
the slab is reached. Given $n_{\rm H}$, $Z_{\rm gas}$, $G$ and $N_{\rm H}$ of
each gas cloud, we interpolate $\rm[C_{II}]$ luminosity of the cloud from the
values found in the computed grid.

Figure 2: Top and middle panels: Ionization structures of a plane-parallel gas
slab ($n_{\rm H}=50\,\rm cm^{-3}$) irradiated by an external radiation field
($G=200\,G_{0}$) incident from the left in the figure predicted by the CLOUDY
code. Dashed red, solid green and dotted blue lines in the top (middle) panels
represent the abundance profiles for $\rm H_{II}$ ($\rm C_{III}$), $\rm H_{I}$
($\rm C_{II}$) and $\rm H_{2}$ ($\rm C_{I}$), respectively. Solid black line
in the top and middle panels shows the profile of gas kinetic temperature
(normalized by $10^{4}$ K). Solid magenta line in the middle panels indicates
the profile of $\rm[C_{II}]$ cooling rate (normalized by $10^{-22}\,\rm
erg\,s^{-1}\,cm^{-3}$). Bottom panels: Cumulative fraction of $\rm[C_{II}]$
luminosity (thick orange line) and volume (thin blue) as a function of gas
column density (from the surface) of a spherical gas cloud ($M_{\rm
cl}=10^{5}\,M_{\odot}$, $n_{\rm H}=50\,\rm cm^{-3}$) irradiated by an external
radiation field ($G{}={}200\,G_{0}$). Red dotted line marks the surface-to-
centre column density of the cloud ($N_{\rm H}=4\times 10^{21}\,\rm cm^{-2}$).
The left and right columns correspond to the metal-rich and metal-poor models
where gas metallicity of the slab (cloud) is set to $Z_{\odot}$ and
$1/10\,Z_{\odot}$. For the metal-poor model, the dust-to-gas mass ratio
($\delta_{\rm dgr}$) becomes lower and therefore Lyman-Werner photons can
penetrate deeper into the slab (cloud), resulting in larger
$\rm[C_{II}]$-emitting region (Zone I + Zone II).

##### CLOUDY simulation: an example

Figure 3: The UVJ false-colour image (left), ${\rm[C_{II}]}$ surface
brightness (middle) and the distribution of UV ISRF strength ($G$) (right) of
selected FIRE galaxies. The upper panels show the results of a $z{}={}0$ disc
galaxy from FIREbox (c.f. Fig.3 of Feldmann et al. 2022), while the lower
panels correspond to a galaxy undergoing multiple mergers at $z{}={}6$
extracted from the MassiveFIRE ‘zoom-in’ suite.

Here we show the conditions of a plane-parallel gas slab calculated by CLOUDY
(Fig. 2). The slab has a uniform gas density $n_{\rm H}{}={}50\,\rm cm^{-3}$
and is illuminated by an external radiation field having $G{}={}200\,G_{0}$.
We present CLOUDY simulations for two different models, where $Z_{\rm gas}$ is
set to $Z_{\odot}$ and $1/10\,Z_{\odot}$. We include the $z{}={}0$ CMB
background and the CR hydrogen ionization rate is set to the default value. We
show the results of the dust-rich and dust-poor models in the left and right
panels of Fig. 2, respectively.

The slab is characterized by three distinct zones based on the ionization
state of hydrogen gas. In the upper panels, we show the abundance profiles for
ionized hydrogen ($\rm H_{II}$; dashed red line), atomic hydrogen ($\rm
H_{I}$; solid green line) and molecular hydrogen ($\rm H_{2}$; dotted blue
line), as well as the profile for gas temperature (solid black line). We can
see that a $\rm H_{II}$ region (Zone I) is created near the surface of the
slab by the ionizing photons ($E_{\gamma}{}>{}13.6$ eV) of the incident
radiation field. Gas in this region is heated to high temperature
($T{}\approx{}10^{4}$ K). The slab then transits to a $\rm H_{I}$-dominating
region (Zone II) at a distance where ionizing radiation gets fully absorbed.
The photons in the Lyman-Werner (LW) band ($11.2{}<{}E_{\gamma}{}<{}13.6$ eV)
dissociate $\rm H_{2}$ in this region, while maintaining gas temperature at
about $10^{2}$ K. Finally, the slab transits to a $\rm H_{2}$-dominating
region (Zone III) at some larger distance, beyond which the LW radiation
becomes sufficiently absorbed and the majority of hydrogen turns into $\rm
H_{2}$.

Like hydrogen, carbon has a very different ionization state in the three
zones. This can be seen from the middle panels of Fig. 2, where we explicitly
show the abundance profiles for atomic carbon ($\rm C_{I}$; dotted blue line),
singly ionized carbon ($\rm C_{II}$; solid green line) and doubly ionized
carbon ($\rm C_{III}$; dashed red line) for the two models. Carbon is mostly
ionized in Zone I and II. Specifically, in Zone I, it gets excited to $\rm
C_{II}$ level as well as higher ionization levels (e.g. $\rm C_{III}$). In
Zone II, on the contrary, carbon is singly ionized by LW photons but not
excited to higher levels since ionizing photons are shielded from the
region121212Note: the ionization energy of $\rm C_{III}$ is 24.39 eV, which is
above the ionization energy of hydrogen atom (13.6 eV).. Finally, in Zone III,
carbon turns into $\rm C_{I}$ since the region is UV-dark131313Note: the first
ionization energy of carbon is 11.26 eV. This coincides with the lower
frequency limit of the LW band (11.2 eV). Hence, carbon is neutral in the $\rm
H_{2}$ region..

$\rm[C_{II}]$ emission originates mostly from the ionized (Zone I) and atomic
hydrogen (Zone II) phases. We show in the middle panels the profile for
$\rm[C_{II}]$ cooling rate ($\rm erg\,s^{-1}\,cm^{-3}$),
$\Lambda_{\rm[C_{II}]}$, for the two models (solid magenta line). It is clear
that $\Lambda_{\rm[C_{II}]}$ drops sharply in Zone III, which is due to the
very low abundance of $\rm C_{II}$ ions (solid green line) in this region
(note: most carbon is in $\rm C_{I}$ state in Zone III). For the chosen
models, $\Lambda_{\rm[C_{II}]}$ appears to be similar in the ionized and
atomic hydrogen phases, varying by less than a factor of few. Comparing the
metal-rich (left panel) and metal-poor (right panel) models, it can be seen
that $\Lambda_{\rm[C_{II}]}$ of the metal-rich model is about a factor of ten
higher. This is due to the fact that $\Lambda_{\rm[C_{II}]}$ is linearly
scaled to $Z_{\rm gas}$ and $Z_{\rm gas}$ of the metal-rich model is set as
ten times that of the metal-poor model.

Using the $\Lambda_{\rm[C_{II}]}$ profile output by CLOUDY, we subsequently
derive the $\rm[C_{II}]$ luminosity profile (cumulative $\rm[C_{II}]$
luminosity as a function of column depth from the surface) for a uniform
spherical cloud having $n_{\rm H}{}={}50\,\rm cm^{-3}$ (same as the gas slab)
and $M_{\rm cl}{}={}10^{5}\,M_{\odot}$ that is irradiated by an external field
having $G{}={}200\,G_{0}$ (same as the gas slab) following equation (1). We
calculate the result for the metal-rich ($Z_{\rm gas}{}={}Z_{\odot}$) and
metal-poor ($Z_{\rm gas}{}={}0.1Z_{\odot}$) models, which are shown in the
lower left and lower right panels of the figure, respectively. It can be seen
that about $30\%$ ($20\%$) of the total $\rm[C_{II}]$ luminosity of the cloud
is produced by the $\rm H_{II}$ region for the metal-rich (poor) model, while
the remainder originates almost totally from the $\rm H_{I}$ region. The $\rm
H_{2}$ region contributes very limited fraction of the $\rm[C_{II}]$
luminosity. Note that the $\Lambda_{\rm[C_{II}]}$ profile, the size of the
different zones, and their relatively contribution to the total $\rm[C_{II}]$
luminosity of the cloud depends on $G$, $n_{\rm H}$ and $Z_{\rm gas}$ (see
Section 5.1 for a detailed discussion).

One major difference between the two models (metal-rich vs. metal-poor) is
whether or not the gas cloud has an $\rm H_{2}$ region in the core, as can be
seen from the bottom panels. For the metal-poor model (bottom right panel),
because dust column density is small, LW photons are able to penetrate the
entire cloud, making it $\rm H_{2}$-free. The metal-rich model (bottom left
panel), in contrast, has an $\rm H_{2}$ core owing to the high dust column
density, which accounts for nearly half of $M_{\rm cl}$. The two cloud models
correspond to the two distinct regimes where $L_{\rm[C_{II}]\,,cl}$ has
different scaling with $Z_{\rm gas}$. When the cloud has no $\rm H_{2}$ core,
$L_{\rm[C_{II}]\,,cl}$ scales linearly with $Z_{\rm gas}$. As $Z_{\rm gas}$
(and hence the dust-to-gas mass ratio, $\delta_{\rm dgr}$) increases, the
depth of Zone I+Zone II decreases (Ferrara et al., 2019). When $Z_{\rm gas}$
is high enough that $\rm H_{2}$ becomes abundant (i.e. , Zone III forms),
$L_{\rm[C_{II}]\,,cl}$ saturates and no longer depends sensitively on $Z_{\rm
gas}$. In Section 5, we will discuss in detail how the
$L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio of the FIRE galaxies depends on gas
metallicity, and interpret the results using the insights obtained from the
toy models presented here.

###  3.2 Calculating ISRF distribution and multi-wavelength SEDs of galaxies
using SKIRT

To predict the $\rm[C_{II}]$ luminosity of the ISM, it is essential to know
the local UV ISRF strength. We calculate the ISRF distribution for the FIRE
galaxies using the open-source141414Code repository:
https://skirt.ugent.be/version8/ 3D Monte Carlo dust RT code SKIRT (Baes et
al., 2011; Baes & Camps, 2015; Camps & Baes, 2015) (version 8). SKIRT provides
full treatment of absorption and anisotropic scattering by dust, and self-
consistently computes dust thermal re-emission and dust temperature
distribution for various astrophysical systems.

To prepare the galaxy snapshots as RT input models for SKIRT, we follow the
prescription of Camps et al. (2016) (see also Trayford et al. 2017; Camps et
al. 2018). We summarize the key points of the prescription here, and refer
interested readers to the above-mentioned papers for the details.

For the analysis, each star particle of the galaxy is treated as a ‘single
stellar population’ (SSP), and a spectrum of stellar emission is assigned to
each particle using the STARBURST99 (Leitherer et al., 1999; Vazquez &
Leitherer, 2005) SED libraries according to the age, metallicity and initial
mass of the particle. The RT calculations are performed on an equally spaced
logarithmic wavelength grid consisting of 250 wavelength points spanning the
wavelength range $\lambda{}={}0.05-1000\,\rm\mu m$. We launch $10^{6}$ photon
packages for each of the 250 point in the wavelength grid and for each of the
stellar emission and following dust emission stages. The calculation iterates
until convergence. To produce mock images and SEDs for the galaxies, we place
mock detectors at an arbitrary ‘local’ distance of 10 Mpc from galaxy along
multiple viewing angles to accumulate both spatially resolved as well as
integrated fluxes at each wavelength grid point.

We assume that dust mass traces metal mass in galaxies (Hayward et al., 2011;
Narayanan et al., 2015; Camps et al., 2016; Trayford et al., 2017; Liang et
al., 2018, 2019, 2021; Ma et al., 2019; Cochrane et al., 2019; Cochrane et
al., 2022; Vogelsberger et al., 2020; Shen et al., 2022) and adopt a constant
dust-to-metal mass ratio $\delta_{\rm dzr}=0.4$ in gas cooler than $10^{6}$ K.
Hotter gas is assumed to be dust-free due to thermal sputtering (Draine &
Salpeter, 1979; Tielens et al., 1994). We adopt the Weingartner & Draine
(2001b) dust model with Milky-Way size distribution for the case of $R_{\rm
V}=3.1$. We discretize the spatial domain using an octree grid and keep
subdividing grid cells until the cell contains less than $f{}={}3\times
10^{-6}$ of the total dust mass and the $V$-band (0.55 $\rm\mu m$) optical
depth in each cell is less than unity. The highest grid level corresponds to a
cell width of $\sim 20$ pc, i.e. about twice the minimal SPH smoothing length.
We self-consistently calculate the self-absorption of dust emission and
include the transient heating function to calculate non-local thermal
equilibrium dust emission by transiently heated small grains and PAH molecules
(Baes et al., 2011; Camps & Baes, 2015). To account for the heating of dust by
the cosmic microwave background, we adopt a correction to the dust temperature
using equation (12) of da Cunha et al. (2013).

The final output of the SKIRT simulations includes the ISRF, $J_{\lambda}$
($\rm W\,cm^{-3}\,sr^{-1}$), of each adaptive grid cell. We calculate the UV
ISRF strength ($G$) for each cell by integrating $J_{\lambda}$ over the Habing
band ($6-13.6$ eV) and solid angle ($\Omega$). $G$ is assigned to every gas
particle (‘cloud’) inside the cell for predicting its $\rm[C_{II}]$
luminosity.

In Fig. 3, we show the UVJ composite image (left panels), $\rm[C_{II}]$
surface brightness (middle panels), and $G$ distribution (right panels) for
the two selected FIRE galaxies calculated using CLOUDY and SKIRT. The upper
panels show the results of a disc galaxy at $z{}={}0$ extracted from FIREbox,
whilst the lower panels show the results of a galaxy undergoing multiple
mergers at $z{}={}6$ extracted from the MassiveFIRE simulation (Sim ID: D9).
The $z{}={}6$ galaxy system has much stronger strength of ISRF (right panels)
due to its higher SFR ($220\,M_{\odot}\,\rm yr^{-1}$ vs.
$4.5\,M_{\odot}\rm\,yr^{-1}$) and shows higher $\rm[C_{II}]$ surface
brightness. $L_{\rm[C_{II}]}$ of the $z{}={}6$ system and the $z{}={}0$ galaxy
are $4.8\times 10^{8}\,L_{\odot}$ and $2.8\times 10^{8}\,L_{\odot}$,
respectively.

Figure 4: The $L_{\rm[C_{II}]}$ vs. SFR (left panel) and $L_{\rm[C_{II}]}{}/{}\rm SFR$ vs. SFR (right panel) relations of the $z{}={}0$ galaxies. The filled cyan stars in the two panels show the result of the FIRE galaxies. Black triangles and diamonds show the observational data of Herrera-Camus et al. (2015) (H15) and De Looze et al. (2011) (L11), and the orange and green lines indicate the best-fit linear relation of the H15 and L11 samples, respectively. The coloured shaded regions indicate the $1\sigma$ scatter of the data around the best-fit linear relation of the observed samples. Purple dashed line in the two panels represents the best-fit linear relation to the low-metallicity dwarf galaxy sample of De Looze et al. (2014) (L14). The result of the FIRE galaxies at $z{}={}0$ is in good agreement with the observational data.  Table 2: Observed scaling relations between SFR and $L_{\rm[C_{II}]}$ of local galaxies, i.e. $L_{\rm[C_{II}]}{}/{}L_{\odot}{}={}A{}({\rm SFR}{}/{}M_{\odot}\,{\rm yr}^{-1})^{B}$.  Galaxy sample |  SFR range ($M_{\odot}\,\rm yr^{-1}$) |  Median SFR ($M_{\odot}\,\rm yr^{-1}$) | $A$ | $B$ |  $1\sigma$ scatter
---|---|---|---|---|---
De Looze et al. (2011) | $0.02-88$ | 1.75 |  $7.31\pm 0.06$ |  $0.93\pm 0.06$ | 0.26 dex
De Looze et al. (2014) | $6\times 10^{-4}-56$ | 0.12 |  $7.10\pm 0.11$ |  $1.05\pm 0.07$ | 0.43 dex
Herrera-Camus et al. (2015) | $10^{-3}-9.6$ | 0.34 |  $7.63\pm 0.03$ |  $0.97\pm 0.03$ | 0.21 dex
Figure 5: The $L_{\rm IR}$ vs. SFR relation of FIRE galaxies at different
redshifts (cyan stars for $z{}={}0$, yellow hexagons for $z{}={}1$, red
triangles for $z{}={}2$, blue squares for $z{}={}3$, magenta circles for
$z{}={}4$, green diamonds for $z{}={}6$ and purple downward diamonds for
$z{}={}8$). The diagonal solid black line indicates the K98 relation, i.e.
$L_{\rm IR}\,(L_{\odot}){}={}1.36\times 10^{10}\,{\rm SFR}\,(M_{\odot}\,{\rm
yr^{-1}})$. The dashed and dotted lines indicate the modified K98 relations
where the normalization is lower than the solid black line by a factor of 2
and 10, respectively. The K98 relation (solid black line) fits well to the
galaxies at high SFR. Figure 6: The $L_{\rm[C_{II}]}$ vs. $L_{\rm IR}$ (left
panel) and the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ vs. $L_{\rm IR}$ (right panel)
relations of $z=0$ galaxies. In the two panels, cyan stars show the result of
the FIRE galaxies, whereas black symbols indicate the observational data from
different studies, including Malhotra et al. (2001) (diamond), Brauher et al.
(2008) (vertical crosses), Sargsyan et al. (2012) (filled squares), Farrah et
al. (2013) (empty squares), Díaz-Santos et al. (2013) (filled circles), Magdis
et al. (2014) (diagonal crosses), Cormier et al. (2015) (empty stars),
Herrera-Camus et al. (2015) (asterisks), Hughes et al. (2017) (triangles) and
Contursi et al. (2017) (empty circles). Observations show that
$L_{\rm[C_{II}]}/L_{\rm IR}$ ratio of galaxies is nearly a constant at
$10^{9}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
60\relax$}\hss}\mathchar 29208\relax$}}{}L_{\rm IR}{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}10^{11}\,L_{\odot}$, but declines with $L_{\rm IR}$ at $L_{\rm
IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10^{11}\,L_{\odot}$. In the two
panels, black line (solid at $L_{\rm IR}{}<{}10^{11}\,L_{\odot}$ and dotted at
$L_{\rm IR}{}\geq{}10^{11}\,L_{\odot}$) indicates the median
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio ($\approx{}2\times 10^{-3}$) of the
galaxies having $L_{\rm IR}{}<{}10^{11}\,L_{\odot}$ and grey shaded bar (dark
grey at $L_{\rm IR}{}<{}10^{11}\,L_{\odot}$ and light grey at $L_{\rm
IR}{}\geq{}10^{11}\,L_{\odot}$) indicates the $1\sigma$ scatter of the
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio of the same galaxies. FIREbox
successfully reproduces the observed $L_{\rm[C_{II}]}$ vs. $L_{\rm IR}$ (and
the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ vs. $L_{\rm IR}$) relation at $z{}={}0$.

##  4 Comparison with observations

In this section, we compare the $L_{\rm[C_{II}]}$-SFR relation of the FIRE
galaxies predicted by our model with the observational data at various
redshifts. We separately discuss the results for three redshift regimes,
$z{}={}0$ (Section 4.1), $1{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}5$ (Section 4.2) and $z{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$
(Section 4.3). We make this distinction because observations use different
sample selection methods and the SFR of galaxies is estimated by different
means of calibration in the three different regimes.

###  4.1 Local Universe (redshift $z{}={}0$)

Observations of the $L_{\rm[C_{II}]}$-SFR relation at $z{}={}0$ probe a very
wide SFR range across several orders of magnitude. The selected samples
include low-SFR systems such as dwarf galaxies as well as the extreme IR-
luminous starbursts.

There are three main samples of nearby galaxies that have been used for
calibrating the $L_{\rm[C_{II}]}$-SFR relation of normal star-forming galaxies
(${\rm SFR}{}\approx{}10^{-5}-10\,M_{\odot}\,\rm yr^{-1}$): De Looze et al.
(2011, hereafter L11), De Looze et al. (2014, hereafter L14) and Herrera-Camus
et al. (2015, hereafter H15). The L11 sample consists of 24 star-forming
galaxies selected from the early compilation by Brauher et al. (2008) that
have measurements at both the Galaxy Evolution Explorer (GALEX) FUV and the
Multiband Imaging Photometer for Spitzer (MIPS) $\rm 24\,\mu m$ bands. The
sample of L14 includes 48 nearby low-metallicity ($Z_{\rm gas}\approx
0.03-0.55\,Z_{\odot}$) dwarf galaxies extracted from the Dwarf Galaxy Survey
(DGS, Madden et al., 2013) catalogue. Lastly, H15 study a sample consisting of
46 local star-forming galaxies chosen from the KINGFISH catalogue (Kennicutt
et al., 2011), having very diverse integrated galaxy properties and ISM
environments. All these studies have excluded the sources with AGN features.

Both L11 and L14 derive the SFR of their sample using GALEX FUV and MIPS
$24\,\rm\mu m$ fluxes (i.e. ${\rm SFR}=\beta\,(L_{\rm FUV,\,obs}+\alpha\times
L_{\rm 24\,\mu m})$) but with different calibration. Specifically, L11 and L14
use the calibration by Zhu et al. (2008) ($\alpha{}={}6.31$) and Hao et al.
(2011) ($\alpha{}={}3.89$), respectively. H15, on the other hand, derive the
SFR of their sample using a hybrid of different methods: for 27 galaxies in
their sample, SFR is derived using the $\rm H_{\alpha}$+$24\,\rm\mu m$
calibration by Calzetti et al. (2007) (equation 7). For the other 8 galaxies,
they use the FUV+$24\,\rm\mu m$ calibration by Leroy et al. (2008) (equations
D10 and D11). And lastly, for the remaining 11 galaxies having no measurement
of either $\rm H_{\alpha}$ nor FUV flux, SFR is derived based solely on their
$24\,\rm\mu m$ flux using the calibration by Calzetti et al. (2007) (equation
6). In Table 2, we show the SFR range as well as the median SFR of the three
samples (L11, L14 and H15). We also show in the table the best-fit parameter
values for the scaling relation

${\rm log}{}(L_{\rm[C_{II}]}{}/{}L_{\odot}){}={}A+B\;{\rm log}{}({\rm
SFR}{}/{}M_{\odot}\,\rm yr^{-1})$ (3)

for the three samples as well as the $1\sigma$ scatter (in dex) of the data
around the best-fit relation. Note that for the galaxies of the L11 and H15
samples whose SFR is derived using the FUV+$24\,\rm\mu m$ fluxes, we have re-
calibrated their SFR following Hao et al. (2011) as has been done by L14 for a
fair comparison. All the SFR calibrations are based on the Kroupa (2002)
initial mass function (IMF).

From Table 2, we can see that the three samples all exhibit an almost linear
correlation between $L_{\rm[C_{II}]}$ and SFR, though having noticeable
difference in the normalization. The H15 sample has the highest normalization
among the three samples. It is higher than that of the L11 sample by 0.32 dex.
This offset may partly be due to the difference in sample selection. Another
potential cause is that H15 adopt different SFR indicators and calibration
methods compared with L11 for a large fraction of the galaxies in their
sample. The offset between the L11 and L14 samples (0.21 dex), on the other
hand, is mainly due to the difference in sample selection since L11 and L14
adopt the same SFR indicators (FUV+$\rm 24\,\mu m$ fluxes) for their entire
samples and we have re-calibrated their results following the same method of
Hao et al. (2011). The lower normalization of the L14 relation is very likely
due to the relatively lower $Z_{\rm gas}$ of the dwarf galaxies they use for
the study, as has been explicitly stated in L14.

In Fig. 4, we show the $L_{\rm[C_{II}]}$-SFR relation of the three samples
(L11, L14 and H15) in the left panel. To more clearly show the difference in
the normalization of these scaling relations, we present the
$L_{\rm[C_{II}]}{}/{}\rm SFR$ vs. SFR relation of the same samples in the
right panel. In both panels, we also present the results for the FIRE
sample151515We calculate the SFR of the FIRE galaxies by averaging over a
timescale of the last 100 Myrs. $z{}={}0$ (filled cyan stars) for comparison
with the observational data. Note that for the L11 and H15 samples, we show
both the data of the individual sources as well as the best-fit scaling
relation for each sample, whereas for the L14 sample, we only present the
best-fit scaling relation (purple dashed line) for reference. The L14 sample
has systematically lower gas metallicity than the other two observational
samples as well as the FIRE galaxy sample at $z{}={}0$.

The FIRE simulations, combined with our line model, produce the
$L_{\rm[C_{II}]}$-SFR relation at $z{}={}0$ (cyan stars) that is in good
agreement with the local star-forming samples of L11 (black diamonds) and H15
(black triangles). The best-fit parameter values for the FIRE galaxies (over
the SFR range of $0.05-100\,M_{\odot}\,\rm yr^{-1}$) are $A=7.70\pm 0.02$ and
$B=0.80\pm 0.03$, and the $1\sigma$ scatter of the data points around the
best-fit relation is 0.21 dex, similar to the L11 and H15 samples (see Table
2). Note that we have excluded the galaxies having
SFR${}<{}0.05\,M_{\odot}\,\rm yr^{-1}$ for the fitting to avoid the regime
where galaxy statistics can be contaminated by the shot noise due to the
resolution limit of the simulation (Feldmann, 2017).

The FIRE galaxies exhibit a sub-linear relation between $L_{\rm[C_{II}]}$ and
SFR, i.e. $L_{\rm[C_{II}]}\propto{\rm SFR}^{0.80\pm 0.03}$. The sub-linearity
is due to the galaxies having SFR${}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}3\,M_{\odot}\,\rm yr^{-1}$, which have lower
$L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio on average than the galaxies having lower
SFR (see the right panel of Fig. 4). Such a trend of reduced
$L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio at high SFR ($\rm[C_{II}]$ deficit) is not
clearly present in any of the three (L11, L14 and H15) observational samples.
We note, however, that these samples do not contain statistically large number
of galaxies at SFR${}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}3\,M_{\odot}\,\rm
yr^{-1}$. Some other studies probing the local LIRGs and ULIRGs have found
clear evidence of $\rm[C_{II}]$ deficit at high $L_{\rm IR}$ ($\rm\sim SFR$)
(see below).

##### The $L_{\rm[C_{II}]}$ vs. $L_{\rm IR}$ relation of $z{}=0{}$ galaxies

A number of observational studies have probed the relation between
$L_{\rm[C_{II}]}$ and $L_{\rm IR}$ (or $L_{\rm FIR}$161616In the literature,
‘$L_{\rm IR}$’ is used to denote the bolometric IR luminosity of galaxy that
is integrated over the wavelength range $\rm 8-1000\,\mu m$, whereas ‘$L_{\rm
FIR}$’ represents the FIR luminosity of galaxy ($\rm 42.5-122.5\,\mu m$). Both
$L_{\rm IR}$ and $L_{\rm FIR}$ are commonly adopted as SFR indicators for
heavily dust-obscured galaxies.) of local galaxies.

$L_{\rm IR}$ (or $L_{\rm FIR}$) can be a good proxy for galaxy SFR when the
stellar light of a galaxy is heavily absorbed by dust (e.g. Kennicutt, 1998;
Salim & Narayanan, 2020). Galaxies having higher SFR tend to be more gas/dust-
rich and have higher gas density. Therefore, they tend to have higher dust
opacity (e.g. Whitaker et al., 2017). We show in Fig. 5 the $L_{\rm IR}$ vs.
SFR relation of the FIRE galaxies at different redshifts, where $L_{\rm IR}$
is calculated using their SEDs produced by SKIRT. It can be seen that at
$z{}={}0$, the FIRE galaxies (cyan stars) well follow the Kennicutt (1998,
hereafter K98) relation171717We adopt the K98 relation for the Kroupa (2002)
IMF using the stellar population synthesis (SPS) model STARBURST99, assuming a
constant star formation history lasting for 1 Gyr (see Hao et al. 2011 for the
details). The original relation (i.e., $L_{\rm IR}{}/{}L_{\odot}=5.8\times
10^{9}\;{\rm SFR}{}/{}(M_{\odot}\,\rm yr^{-1}$)) was derived for the Salpeter
IMF based on the older SPS model of Leitherer & Heckman (1995), and for a
shorter starburst period ($t_{\rm age}=10-100\,\rm Myrs$). , i.e.

$L_{\rm IR}\;(L_{\odot})=1.36\times 10^{10}\;{\rm SFR}\;(M_{\odot}\,\rm
yr^{-1})$ (4)

at $\rm SFR{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}1\,M_{\odot}\,\rm
yr^{-1}$ (or $L_{\rm IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{10}\,L_{\odot}$). The K98 relation is derived assuming
that all radiative energy of the young stars is absorbed and re-emitted by
dust and AGN radiation does not contribute to dust heating. At $\rm
SFR{}<{}1\,M_{\odot}\,\rm yr^{-1}$, however, the $z{}={}0$ FIRE galaxies show
larger scatter. Some of these galaxies are below the K98 relation by over 0.3
dex (indicating that less than half of the radiative energy of the young stars
gets re-emitted at FIR by dust). These are the galaxies having relatively low
dust opacity181818It can be seen from Fig. 5 that some of the simulated
galaxies (particularly those having low SFR) lie above the K98 relation, which
seem to break the energy conservation law. These are in fact the galaxies that
are recently quenched after a strong starburst whose dust is heated mainly by
the stars older than $100\,\rm Myrs$ (see e.g. Hayward et al., 2014)..
Nonetheless, $L_{\rm IR}$ appears to be a good SFR tracer for the $z{}={}0$
galaxies at $\rm SFR{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}1\,M_{\odot}\,\rm
yr^{-1}$ in the FIRE simulations.

We show in Fig. 6 the observed $L_{\rm[C_{II}]}$ vs. $L_{\rm IR}$ (left panel)
and the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ vs. $L_{\rm IR}$ (right panel)
relations of the local galaxy samples of Malhotra et al. (2001), Brauher et
al. (2008), Sargsyan et al. (2012), Farrah et al. (2013), Díaz-Santos et al.
(2013), Magdis et al. (2014), Cormier et al. (2015), Herrera-Camus et al.
(2015), Hughes et al. (2017) and Contursi et al. (2017). Note that for those
having used $L_{\rm FIR}$ as SFR indicator, we convert the reported $L_{\rm
FIR}$ of the galaxies to $L_{\rm IR}$ by multiplying it with 1.6 (Sanders et
al., 2003). For comparison, we also show in the same figure the data of the
$z{}={}0$ FIRE galaxies. $L_{\rm IR}$ of the FIRE galaxies is computed by
integrating the SKIRT-produced SED over the wavelength range $\rm 8-1000\,\mu
m$.

The observed samples contain a large number of galaxies that are IR-luminous
($L_{\rm IR}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}10^{11}\,L_{\odot}$, corresponding to
${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10\,M_{\odot}\,{\rm yr}^{-1}$
following equation 4). With these statistically large samples, the
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ($\sim L_{\rm[C_{II}]}{}/{}{\rm SFR}$) ratio
of the $z{}={}0$ galaxies appear to show a clear decline with $L_{\rm IR}$ at
$L_{\rm IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}10^{11}\,L_{\odot}$
($\rm[C_{II}]$ deficit), albeit with a large scatter ($1\sigma{}={}0.3$ dex)
at given $L_{\rm IR}$. From $L_{\rm IR}{}={}10^{11}$ to $10^{13}\,L_{\odot}$,
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ decreases from $2{}\times{}10^{-3}$ to
$10^{-4}$, over a factor of ten. At $L_{\rm IR}{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}10^{11}\,L_{\odot}$, on the other hand,
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ of the observed galaxies is a constant.
Overall, the observational and the simulated data agree well with each other
(on both the mean value and level of scatter). In particular, the FIRE sample
exhibits a mild $\rm[C_{II}]$ deficit at $L_{\rm IR}{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{11}\,L_{\odot}$ at $z{}={}0$, which is in agreement with
the observational data. Note, however, that our FIRE sample at $z{}={}0$ does
not include any ULIRGs (i.e. $L_{\rm IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{12}\,L_{\odot}$) at $z{}={}0$.

###  4.2 High redshifts ($1{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}5$)

Figure 7: The $L_{\rm[C_{II}]}$ vs. SFR (left panel) and the
$L_{\rm[C_{II}]}/L_{\rm IR}$ vs. $L_{\rm IR}$ (right panel) relations of
galaxies at $z{}={}0$ and high redshifts. In both panels, filled coloured
symbols represent the FIRE galaxies (cyan stars for $z{}={}0$, yellow hexagons
for $z{}={}1$, red triangles for $z{}={}2$, blue squares for $z{}={}3$ and
magenta circles for $z{}={}4$). Black symbols (filled and empty) show the
observational data of galaxies at $1{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}5$ (see Table 3 for
the details). Specifically, black circles and black triangles correspond to
SMGs and other star-forming galaxies, respectively. For the gravitationally-
lensed galaxies, their $\rm[C_{II}]$ and IR luminosities have been corrected
by the lensing magnification factor $\mu$ reported in the literature. Those
having direct measurement of $\mu$ as well as the un-lensed galaxies are
marked by filled symbols (triangles and circles), whereas the 16 lensed SPT
galaxies whose $\mu$ is extrapolated ($\mu$ has been assumed to be 14.1 by
Gullberg et al. 2015) are shown by empty circles. The two grey empty squares
represent the stacked result of the galaxy samples of Stacey et al. (2010) and
Brisbin et al. (2015). The $\rm[C_{II}]$ line of the two samples is measured
with the redshift and Early Universe Spectrometer (ZEUS) and their data
systematically offsets from that of the other galaxy samples. For reference,
we also show in the left (right) panel the observational results of the local
galaxy samples as shown in Fig. 4 (Fig. 6). Both observations and FIRE
simulations show that high-$z$ ($1{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}5$) galaxies
exhibit a $\rm[C_{II}]$ deficit at high $L_{\rm IR}$ similar to local
galaxies.

Observational studies have investigated the $L_{\rm[C_{II}]}$-SFR relation of
galaxies at $1{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}5$, including e.g. Ivison et al. (2010); Stacey et al. (2010);
Valtchanov et al. (2011); Brisbin et al. (2015); Gullberg et al. (2015, 2018);
Schaerer et al. (2015b); Umehata et al. (2017); Zanella et al. (2018);
Hashimoto et al. (2019b); McKinney et al. (2020). Their samples consist of
roughly 80 galaxies in total (see Table 3 for the details). Most of these
galaxies have substantial SFR (${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}100\,M_{\odot}\,\rm yr^{-1}$) and are IR-luminous ($L_{\rm
IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}10^{12}\,L_{\odot}$). This is in
stark contrast with the local observations (see Section 4.1), which probe the
galaxies having much lower SFR (see Table 2). Note that a large fraction of
the selected galaxies in this redshift regime are uncovered by wide-field sub-
mm galaxy surveys, e.g. the South Pole Telescope (SPT; Vieira et al., 2010;
Carlstrom et al., 2011) survey (Weiß et al., 2013; Gullberg et al., 2015).

We derive the SFR of the selected galaxies from their measured $L_{\rm IR}$
(see Table 3) using the K98 relation (equation 4) assuming that the galaxies
are heavily dust-obscured. Note that at high redshifts, the K98 relation may
only apply to the more massive and starburst galaxies. High-$z$ galaxies are
metal and dust-poorer than the $z{}={}0$ galaxies at given mass (or SFR), and
therefore only the more massive and gas-rich systems have high enough dust
opacity leading to total obscuration of stellar light. We can see from Fig. 5
that the K98 relation (solid black line) fits well the high-$z$ FIRE galaxies
at ${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}100\,M_{\odot}\,\rm
yr^{-1}$ (or $L_{\rm IR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}10^{12}\,L_{\odot}$. Note: For the $z{}={}1$ galaxies, the K98
relation fits well to the data down to $L_{\rm
IR}{}\approx{}10^{11}\,L_{\odot}$). At lower SFR, the high-$z$ galaxies
exhibit larger scatter and they, on the average, have lower $L_{\rm IR}$ at
given SFR than the $z{}={}0$ galaxies due to their reduced dust opacity.

The galaxies selected at $1{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}5$ typically have a good sampling of photometric data points
in the dust continuum, which are obtained by observations with multiple IR and
millimetre instruments (Spitzer, Herschel, ALMA and etc). The shape of the
dust SED of these galaxies is therefore well constrained. This results in
relatively small uncertainty in the estimate of their $L_{\rm IR}$.

The $\rm[C_{II}]$ line of these galaxies is measured with different
instruments (see Table 3). For instance, Stacey et al. (2010) and Brisbin et
al. (2015) measure the $\rm[C_{II}]$ line of the 20 galaxies at
$z{}\approx{}1-2$ of their samples using the redshift ($z$) and Early Universe
Spectrometer (ZEUS; Stacey et al., 2007; Hailey-Dunsheath, 2009) on the 10.4 m
Caltech Submillimeter Observatory (CSO). Gullberg et al. (2015) measure the
$\rm[C_{II}]$ line of the 16 SMGs selected from the SPT catalogue (Weiß et
al., 2013) using the SPIRE Fourier Transform Spectrometer (FTS; Griffin et
al., 2010) onboard Herschel (for the galaxies at $z{}<{}3$) and the First
Light APEX Sub-millimetre Heterodyne receiver (FLASH; Heyminck et al., 2006)
(for the galaxies at $z{}>{}3$). For the remaining galaxies ($\sim 40$), their
$\rm[C_{II}]$ line is measured with ALMA (at Band 7, 8 and 9 for the galaxies
at $z\sim 4$, $z\sim 3$ and $z\sim 2$, respectively). ALMA observations often
marginally resolve a galaxy spatially in $\rm[C_{II}]$, whereas observations
with ZEUS, APEX/FLASH and SPIRE FTS do not.

Table 3: The observed $L_{\rm[C_{II}]}$-SFR relation of galaxies at high
redshifts.

Name$\dagger$ | $z$ |  ${\rm log}\,(L_{\rm IR}{}/{}L_{\odot})$$\mathsection$ |  ${\rm log}\,(L_{\rm[C_{II}]}{}/{}L_{\odot})$$\ddagger$, $\mathsection$, $\parallel$ | Galaxy type$\\#$ | AGN | $\mu$ |  References$\ast$
---|---|---|---|---|---|---|---
ID 7118 |  1.7290 |  $12.06\pm 0.01$ |  $<9.70$ (ALMA 9) | MS | No | $-$ |  [1, 2]
GS IRS61 |  1.759 |  $12.46\pm 0.13$ |  $<8.31$ (ALMA 9) | SB | No | $-$ |  [3, 4]
ID 9834 |  1.7644 |  $11.99\pm 0.02$ |  $9.11\pm 0.07$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 2910 |  1.7686 |  $11.76\pm 0.08$ |  $<9.08$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 2861 |  1.8102 |  $12.00\pm 0.03$ |  $<9.58$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 6515 |  1.8438 |  $11.68\pm 0.04$ |  $9.09\pm 0.12$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 9347 |  1.8505 |  $11.80\pm 0.05$ |  $8.98\pm 0.14$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 9681 |  1.8852 |  $11.84\pm 0.04$ |  $9.26\pm 0.20$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 8490 |  1.9056 |  $11.54\pm 0.06$ |  $8.85\pm 0.20$ (ALMA 9) | MS | No | $-$ |  [1, 2]
ID 10049 |  1.9200 |  $11.60\pm 0.06$ |  $<8.78$ (ALMA 9) | MS | Yes | $-$ |  [1, 2]
GS IRS20 |  1.923 |  $13.06\pm 0.12$ |  $9.17\pm 0.01$ (ALMA 9) | SB | Yes | $-$ |  [3, 4]
ID 10076 |  1.9462 |  $11.91\pm 0.03$ |  $9.38\pm 0.14$ (ALMA 9) | MS | No | $-$ |  [1, 2]
MACS J0451+0006 |  2.013 |  $11.08\pm 0.04$ |  $8.08\pm 0.04$ (ALMA 9) | MS | No | $49\pm 5$ |  [5, 6, 7]
GRB 080207 |  2.0865 |  $12.26\pm 0.05$ |  $8.89\pm 0.12$ (ALMA 9) | MS | No | $-$ |  [8]
SPT 0551-50 |  2.123 |  $11.89\pm 0.05$ |  $<9.33$ (SPIRE FTS) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0512-59 |  2.234 |  $12.29\pm 0.04$ |  $9.45\pm 0.09$ (SPIRE FTS) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SMM J2135 |  2.3259 |  $12.08\pm 0.07$ |  $8.25\pm 0.11$ (SPIRE FTS) | SMG | No |  $32.5\pm 4.5$ |  [12, 13]
SDP.130 |  2.625 |  $12.40\pm 0.02$ |  $<10.14$ (SPIRE FTS) | SMG | No | $6\pm 1$ |  [14, 15]
SPT 0538-50 |  2.782 |  $12.44\pm 0.03$ |  $<9.95$ (SPIRE FTS) | SMG | No |  $20.9\pm 4.2$ |  [9, 10]
ALESS 49.1 |  2.943 |  $12.85\pm 0.06$ |  $9.48\pm 0.12$ (ALMA 8) | SMG | No | $-$ |  [16, 17, 18]
ALESS 57.1 |  2.943 |  $12.87\pm 0.06$ |  $9.04\pm 0.17$ (ALMA 8) | SMG | No | $-$ |  [16, 17, 18]
SDP.81 |  3.042 |  $12.32\pm 0.08$ |  $10.06\pm 0.01$ (SPIRE FTS) | SMG | No | $25\pm 7$ |  [14, 15]
SPT 0103-45 |  3.090 |  $12.38\pm 0.02$ |  $9.41\pm 0.06$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
LAB1-ALMA3 |  3.0993 |  11.76 |  $9.41\pm 0.06$ (ALMA 8) | MS | No | $-$ |  [19, 20]
LAB1-ALMA1 |  3.1 |  11.54 |  $<8.9$ (ALMA 8) | MS | No | $-$ |  [19, 20]
LAB1-ALMA2 |  3.1 |  11.60 |  $<8.9$ (ALMA 8) | MS | No | $-$ |  [19, 20]
SPT 0550-53 |  3.129 |  $12.08\pm 0.09$ |  $9.46\pm 0.09$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0529-54 |  3.369 |  $12.36\pm 0.04$ |  $9.74\pm 0.04$ (APEX/FLASH) | SMG | No | $9.4\pm 1.0$ |  [9, 10]
SPT 0532-50 |  3.399 |  $12.69\pm 0.07$ |  $9.46\pm 0.08$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0300-46 |  3.596 |  $12.40\pm 0.11$ |  $9.05\pm 0.11$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 2147-50 |  3.761 |  $12.39\pm 0.06$ |  $9.38\pm 0.06$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0418-47 |  4.224 |  $12.48\pm 0.03$ |  $9.49\pm 0.03$ (APEX/FLASH) | SMG | No | $21.0\pm 3.5$ |  [9, 10]
SPT 0113-46 |  4.232 |  $12.20\pm 0.09$ |  $9.51\pm 0.10$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SDP.141 |  4.24 |  $12.52\pm 0.12$ |  $9.48\pm 0.07$ (APEX/FLASH) | SMG | No | $10-30$ |  [11]
SPT 2311-54 |  4.281 |  $12.40\pm 0.04$ |  $9.23\pm 0.06$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0345-47 |  4.296 |  $12.84\pm 0.04$ |  $9.37\pm 0.04$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
COSMOS-AzTEC-1 |  4.342 |  $13.21\pm 0.09$ |  $9.80\pm 0.04$ (ALMA 7) | SMG | No | $-$ |  [21, 22]
AS2UDS.0568.0 |  4.404 |  $13.30\pm 0.08$ |  $9.20\pm 0.08$ (ALMA 7) | SMG | No | $-$ |  [23, 24]
ALESS 61.1 |  4.4189 |  $12.49\pm 0.03$ |  $9.18\pm 0.17$ (ALMA 7) | SMG | No | $-$ |  [24, 25, 26]
UDS 47.0 |  4.4201 |  $12.50\pm 0.06$ |  $9.42\pm 0.12$ (ALMA 7) | SMG | No | $-$ |  [24, 26]
AS2UDS.0051.0 |  4.421 |  $12.85\pm 0.20$ |  $9.38\pm 0.05$ (ALMA 7) | SMG | No | $-$ |  [23, 24]
AS2UDS.0104.0 |  4.423 |  $12.85\pm 0.20$ |  $9.46\pm 0.05$ (ALMA 7) | SMG | No | $-$ |  [23, 24]
SGP 38326 (SMG1) |  4.4237 |  $13.20\pm 0.09$ |  $9.92\pm 0.05$ (ALMA 7) | SMG (SB) | No | $-$ |  [27]
SGP 38326 (SMG2) |  4.4289 |  $12.90\pm 0.09$ |  $9.46\pm 0.05$ (ALMA 7) | SMG (SB) | No | $-$ |  [27]
BRI 0952-0115 |  4.4337 |  $12.40\pm 0.25$ |  $9.66\pm 0.25$ (APEX/FLASH) | SMG (SB) | No | $4.5\pm 2.8$ |  [28, 29, 30]
SPT 2103-60 |  4.435 |  $12.41\pm 0.03$ |  $9.70\pm 0.06$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
AS2UDS.0232.0 |  4.443 |  $13.26\pm 0.15$ |  $8.70\pm 0.09$ (ALMA 7) | SMG | No | $-$ |  [23, 24]
ALESS 65.1 |  4.4445 |  $12.49\pm 0.03$ |  $9.51\pm 0.09$ (ALMA 7) | SMG | No | $-$ |  [24, 25, 26]
AS2UDS.0109.0 |  4.450 |  $12.90\pm 0.06$ |  $9.42\pm 0.03$ (ALMA 7) | SMG | No | $-$ |  [23, 24]
AS2UDS.0002.1 |  4.4611 |  $13.38\pm 0.08$ |  $8.90\pm 0.11$ (ALMA 7) | SMG | No | $\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}1.5-2$ |  [23, 24]
AS2UDS.0643.0 |  4.4614 |  $13.11\pm 0.22$ |  $8.95\pm 0.15$ (ALMA 7) | SMG | No | $\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}1.5-2$ |  [23, 24]
AS2UDS.0208.0 |  4.4615 |  $12.89\pm 0.01$ |  $9.42\pm 0.06$ (ALMA 7) | SMG | No | $\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}1.5-2$ |  [23, 24]
SPT 0441-46 |  4.477 |  $12.45\pm 0.02$ |  $9.13\pm 0.11$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 2146-55 |  4.567 |  $12.31\pm 0.05$ |  $9.19\pm 0.10$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
W2246–0526 |  4.601 |  $14.34\pm 0.08$ |  $9.79\pm 0.03$ (ALMA 7) | DOG | Yes | $-$ |  [31]
ALESS 73.1 |  4.7555 |  $12.46\pm 0.03$ |  $9.69\pm 0.14$ (ALMA 7) | SMG (SB) | Yes | $-$ |  [24, 26, 32]
SPT 2132-58 |  4.768 |  $12.37\pm 0.04$ |  $9.17\pm 0.08$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
HDF850.1 |  5.185 |  $12.58\pm 0.07$ |  $9.38\pm 0.05$ (IRAM/PdBI) | SMG | No | $1.5-1.7$ |  [33, 34]
HLSJ091828.6+514223 |  5.24 |  $13.04\pm 0.10$ |  $9.98\pm 0.01$ (SMA) | SMG | No | $8.9\pm 1.9$ |  [35]
SPT 2319-55 |  5.293 |  $12.28\pm 0.03$ |  $9.00\pm 0.06$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]
SPT 0346-52 |  5.656 |  $13.39\pm 0.02$ |  $9.97\pm 0.06$ (APEX/FLASH) | SMG | No | $5.4\pm 0.2$ |  [9, 10]
SPT 0243-49 |  5.699 |  $12.40\pm 0.04$ |  $<9.40$ (APEX/FLASH) | SMG | No |  ($14.1\pm 7.8$) |  [9, 10]

(Continue on next page)

Table 4: continued

The observed $L_{\rm[C_{II}]}$-SFR relation of galaxies at high redshifts.
Name$\dagger$ $z$ ${\rm log}\,(L_{\rm IR}{}/{}L_{\odot})$$\mathsection$ ${\rm
log}\,(L_{\rm[C_{II}]}{}/{}L_{\odot})$$\ddagger$, $\mathsection$, $\parallel$
Galaxy type$\\#$ AGN $\mu$ References$\ast$ HerMESFLS3 6.3369 $13.34\pm 0.05$
$9.83\pm 0.10$ (CARMA) SMG No $2.2\pm 0.3$ [36, 37] SPT 0311-58-E 6.900
$12.66\pm 0.12$ $9.62\pm 0.06$ (ALMA 6) SMG No $1.3$ [38] SPT 0311-58-W 6.900
$13.52\pm 0.09$ $9.66\pm 0.06$ (ALMA 6) SMG No $2.2$ [38]

  * $\dagger$

The table does not include the 20 galaxies ($z\approx 2$) in the samples of
Stacey et al. (2010) and Brisbin et al. (2015), of which the $\rm[C_{II}]$
line is measured by ZEUS. The $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ vs $L_{\rm IR}$
relation of these two samples systematically offsets from the others that use
different instrument to measure $\rm[C_{II}]$ line (see Fig. 7).

  * $\mathsection$

For the gravitationally-lensed galaxies, $L_{\rm[C_{II}]}$ and $L_{\rm IR}$
have been de-magnified by the reported lensing magnification factor $\mu$. For
those SPT galaxies having no direct measurement of $\mu$ (galaxies are not
spatially resolved by any observation), we adopt a constant $\mu{}={}14.1$ as
is done by Gullberg et al. (2015), which is the mean of the four galaxies (SPT
0538-50, SPT 0529-54, SPT 0418-47 and SPT 0346-52) in the same sample that is
observationally determined via lensing modelling.

  * $\ddagger$

For the $\rm[C_{II}]$-undetected galaxies, we show the $3\sigma$ upper
confidence limit.

  * $\parallel$

IRAM/PdBI: the IRAM Plateau de Bure Interferometer (Guilloteau et al., 1992);
SMA: the Submillimeter Array (Ho et al., 2004); CARMA: the Combined Array for
Research in Millimeter-wave Astronomy (Woody et al., 2004). Note that the
three telescopes have produced spatially resolved line emission maps of
$\rm[C_{II}]$ for high-$z$ SMGs (HDF850.1, HLSJ091828.6+514223 and HerMESFLS3)
as ALMA does.

  * $\\#$

SMG: sub-mm galaxies; MS: ‘main-sequence’ galaxies; SB: starburst galaxies;
DOG: hot dust-obscured galaxies (galaxies uncovered by surveys at near-
infrared wavelengths, which have strong IR emission from warm dust, e.g. Dey
et al., 2008; Eisenhardt et al., 2012).

  * $\ast$

References: (1): Zanella et al. (2018), [2]: Elbaz et al. (2011), [3]:
McKinney et al. (2020), [4]: Kirkpatrick et al. (2015), [5]:Schaerer et al.
(2015b), [6]: Sklias et al. (2014), [7]: Jones et al. (2010), [8]: Hashimoto
et al. (2019b), [9]: Gullberg et al. (2015), [10]: Weiß et al. (2013), [11]:
Cox et al. (2011), [12]: Ivison et al. (2010), [13]: Swinbank et al. (2010),
[14]: Valtchanov et al. (2011), [15]: Hopwood et al. (2011), [16]: Rybak et
al. (2019), [17]: Wardlow et al. (2018), [18]: da Cunha et al. (2021), [19]:
Umehata et al. (2017), [20]: Geach et al. (2016), [21]: Tadaki et al. (2018),
[22]: Tadaki et al. (2020), [23]: Cooke et al. (2018), [24]: Swinbank et al.
(2014), [25]: Swinbank et al. (2012), [26]: Gullberg et al. (2018), [27]: Oteo
et al. (2016), [28]: Maiolino et al. (2009), [29]: Priddey & McMahon (2001),
[30]: Lehar et al. (2000), [31]: Díaz-Santos et al. (2015), [32]: Breuck et
al. (2014), [33]: Neri et al. (2014), [34]: Walter et al. (2012), [35]: Rawle
et al. (2014), [36]: Riechers et al. (2013), [37]: Cooray et al. (2014), [38]:
Marrone et al. (2018).

It should be particularly noted that a large number (26) of the selected
galaxies (mostly SMGs) in this regime are gravitationally-lensed systems (see
Table 3). Hence, one important source of uncertainty in the estimates of their
intrinsic $L_{\rm[C_{II}]}$ and $L_{\rm IR}$ ($\sim$SFR) is the lensing
magnification factor $\mu$. To observationally determine $\mu$ of a lensed
source requires spatially resolved imaging. Note that 16 of the selected SPT
galaxies in this regime, however, are not spatially resolved by the
observations and their $\mu$ is unknown. Gullberg et al. (2015) adopt a
constant $\mu{}={}14.1$ to de-magnify the luminosities of all the 16 galaxies.
This is the mean of the $\mu$ of the only 4 galaxies in their selected SPT
sample, which is determined using the spatially resolved ALMA $\rm 860\,\mu m$
broadband imaging of dust continuum by Hezaveh et al. (2013).

In Fig. 7, we show the $L_{\rm[C_{II}]}$-SFR relation (left panel) of the
observed samples at $1{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}5$, where we have converted the SFR of all galaxies from their
$L_{\rm IR}$ using the K98 relation following the observational studies. We
show the stacked result for the samples of Stacey et al. (2010) and Brisbin et
al. (2015) by grey empty squares. Both studies measure $\rm[C_{II}]$ line with
ZEUS, and both obtain systematically higher $L_{\rm[C_{II}]}{}/{}\rm SFR$
ratio of galaxies than the other studies using different instruments (by about
one dex) at similar SFR. For the other studies, we explicitly show the data of
each individual source in their samples. Specifically, we show the result of
the SMGs by black circles (empty and filled), whilst the other star-forming
galaxies are denoted by black triangles. For all the lensed galaxies, both
$L_{\rm[C_{II}]}$ and $L_{\rm IR}$ are de-magnified by the observationally
determined $\mu$ when available. For the 16 SPT galaxies having no determined
$\mu$ (indicated by empty black circles in Fig. 7), we correct their
luminosities by an assumed $\mu{}={}14.1$ following Gullberg et al. (2015).
For reference, we also show the $L_{\rm[C_{II}]}$-SFR relation of local
galaxies by L11, L14 and H15 in the same (left) panel.

The bulk of the selected samples at $1{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}{}z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}5$ have higher SFR
than the local samples of L11, L14 and H15. Only the few galaxies at
$z{}\approx{}1-2$ of the Zanella et al. (2018) sample overlap with the SFR
range of the most actively star-forming galaxies of the L11 sample, and they
appear to follow the same $L_{\rm[C_{II}]}$-SFR relation. At higher SFR (i.e.
${\rm SFR}{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}100\,M_{\odot}\,\rm yr^{-1}$), the
high-$z$ galaxy samples show a larger scatter in the $L_{\rm[C_{II}]}$-SFR
relation compared to the local samples (L11, L14 and H15) . Apart from that,
the high-$z$ samples show a decline of $L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio
with increasing SFR at above $100\,M_{\odot}\,\rm yr^{-1}$ (corresponding to
$L_{\rm IR}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}10^{12}\,L_{\odot}$). This trend can be
more clearly seen in the right panel, where we show the
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ($\approx{}L_{\rm[C_{II}]}{}/{}\rm SFR$ at
high SFR) ratio of the same high-$z$ galaxy samples as a function of their
$L_{\rm IR}$ ($\sim$SFR). From $L_{\rm IR}=10^{12}\,L_{\odot}$ to
$10^{13}\,L_{\odot}$, the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ (or
$L_{\rm[C_{II}]}{}/{}\rm SFR$) ratio of the high-$z$ samples decreases by
roughly a factor of 50 (excluding the Stacey et al. 2010 and Brisbin et al.
2015 samples). This $\rm[C_{II}]$ deficit at high $L_{\rm IR}$ is similar to
what has been found with the local galaxy samples (indicated by the filled
grey symbols in Fig. 7).

In the same figure, we also show the results of the FIRE galaxies at high
redshifts. Specifically, we show the $L_{\rm[C_{II}]}$-SFR (left panel) and
the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$-$L_{\rm IR}$ (right panel) relations of
the FIRE galaxies at $z{}={}1$ (yellow hexagons), $z{}={}2$ (red triangles),
$z{}={}3$ (blue squares) and $z{}={}4$ (magenta circles). For reference, we
also show in the two panels the results of the FIRE sample at $z{}={}0$ (cyan
stars).

The FIRE galaxies follow a roughly linear $L_{\rm[C_{II}]}$-SFR scaling
relation over the SFR range of $\approx{}0.05-100\,M_{\odot}\,\rm yr^{-1}$ at
each redshift (left panel), though having considerable scatter
($1\sigma{}\approx{}0.2-0.35$ dex). The normalization of the relation,
however, shows clear redshift evolution. From $z{}={}0$ to $z{}={}4$, the mean
$L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio of the FIRE sample declines by about one
dex (see the left panel of Fig. 7). This indicates that using the
$L_{\rm[C_{II}]}$-SFR relation derived by L11 or H15 will lead to a systematic
underestimate of SFR of galaxies at high redshifts.

On the other hand, the $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio of the FIRE
galaxies does not evolve as much with redshift between $z{}={}0-4$ (right
panel). From $z{}={}0$ to $z{}={}4$, the mean $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$
ratio of the FIRE galaxies decreases by 0.5 dex, which is less than the
decrease of the $L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio ($\sim 1$ dex). Obviously,
the reason for the discrepancy in the redshift evolution of the two ratios
($L_{\rm[C_{II}]}{}/{}\rm SFR$ and $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$) is the
redshift evolution of the $L_{\rm IR}$-SFR relation of the galaxies (see Fig.
5 for the result of the FIRE galaxies, and also the observational data of e.g.
Whitaker et al. 2017) — at fixed SFR, galaxies at higher redshift have on
average lower dust opacity and thus a smaller fraction of stellar radiation is
absorbed and re-emitted at far-IR. The mean $L_{\rm IR}{}/{}\rm SFR$ ratio of
galaxies therefore decreases with redshift.

Apart from that, it is clear from the right panel that the FIRE galaxies at
$z{}={}1-4$ show a similar decrease of $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio
with $L_{\rm IR}$ like the local $z{}={}0$ FIRE galaxies (cyan stars), and the
decrease appears to be more significant at $L_{\rm IR}{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5\times 10^{11}\,L_{\odot}$. The sharp decrease of
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ at the high $L_{\rm IR}$ end is in line with
the trend in the observational data at similar redshifts. In Section 5, we
will examine in detail the origin of this ‘$[\rm C_{II}]$ deficit’ at high
$L_{\rm IR}$ and we will show that it is mainly driven by the decrease of gas
mass per unit SFR, or depletion timescale ($t_{\rm dep}{}\equiv{}M_{\rm
gas}{}/{}\rm SFR$), of galaxies with SFR.

Note that at $L_{\rm IR}{}\approx{}10^{12}\,L_{\odot}$, the observed
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio of the galaxies at high redshifts
(black symbols) appears to be higher than that of the observed $z{}={}0$
galaxy samples (grey symbols) as well as the FIRE galaxies (coloured symbols).
The mean $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ ratio is roughly in agreement with
the upper bound of the FIRE galaxies at similar $L_{\rm IR}$. This can
possibly be due to selection effect. Those galaxies at $L_{\rm
IR}{}\approx{}10^{12}\,L_{\odot}$ are mostly the ‘main-sequence’ (MS) galaxies
at $z{}\approx{}1.5-2$ selected by Zanella et al. (2018), which are expected
to have longer $t_{\rm dep}$ (i.e. gas mass per unit SFR) than starburst
galaxies at the same redshift (e.g. Genzel et al., 2015; Aravena et al., 2016;
Miettinen et al., 2017; Tacconi et al., 2018; Feldmann, 2020) and hence higher
$L_{\rm[C_{II}]}{}/{}L_{\rm IR}$ (note: $L_{\rm[C_{II}]}{}/{}{\rm
SFR}{}\propto{}t^{0.7}_{\rm dep}$, equation 30). The FIRE sample as well as
the local observed galaxy samples, on the contrary, consist of galaxies across
the star-forming MS as well as starburst galaxies, exhibiting a wide range of
$t_{\rm dep}$.

Finally, we note that the observational data in this redshift regime has large
uncertainties due to the large fraction of gravitationally-lensed galaxies
included in the samples (see Table 3). First of all, as mentioned above, many
of the lensed galaxies do not have determined magnification factor $\mu$
(marked by empty circles in Fig. 7). Even for those whose $\mu$ is derived
from either the rest-UV (with Hubble Space Telescope) or dust continuum
imaging (with ALMA), it is not yet certain whether their $\rm[C_{II}]$
luminosity is magnified by the same level, given that the $\rm[C_{II}]$ and
stellar/dust emission of galaxies may have different spatial configuration
(e.g. Cochrane et al., 2019; Fujimoto et al., 2019; Matthee et al., 2019;
Novak et al., 2020; Fudamoto et al., 2022) and thus the different emission
components may have different $\mu$ due to the effect of differential lensing
(e.g. Blain, 1999; Hezaveh et al., 2012; Serjeant, 2012; Cañameras et al.,
2018; Yang et al., 2019a; Harrington et al., 2021). Hence, it is important to
obtain spatially resolved imaging of both $\rm[C_{II}]$ and dust emission for
lensed galaxies and re-examine the intrinsic $L_{\rm[C_{II}]}{}/{}L_{\rm IR}$
ratio of these galaxies (note: most of the lensed SMGs do not have spatially
resolved $\rm[C_{II}]$ imaging, see Table 3).

Figure 8: Comparison of the $L_{\rm[C_{II}]}$-SFR relation of the FIRE
galaxies with the observational data at high redshifts. In the two panels, we
show the result of the FIRE galaxies at $z{}={}4$, $z{}={}6$ and $z{}={}8$ by
magenta circles, green diamonds and purple downward triangles, respectively.
We also show in the two panels the observational data of the rest-UV-selected
star-forming galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$, including the
ones targeted by the ALPINE (blue symbols) and REBELS (red symbols) ALMA
surveys as well as the others targeted by the other observations (black
symbols) (see Table 5 for the details). The galaxies having both confirmed
$\rm[C_{II}]$ and dust continuum detection are indicated by vertical (REBELS)
and diagonal (red for ALPINE and black for others) crosses. The galaxies
having no $\rm[C_{II}]$ detection are shown by downward arrows in both panels.
The location of the arrows indicate the $3\sigma$ upper limit of their
$L_{\rm[C_{II}]}$. For the ones having $\rm[C_{II}]$ but without dust
detection (meaning that their $\rm SFR_{IR}$ is uncertain), we show the
relation between their $L_{\rm[C_{II}]}$ and the lower (upper) SFR limit in
the left (right) panel by rightward (leftward) triangles. For reference, we
also show the result of local ($z{}={}0$) observations of normal star-forming
galaxies (SFGs) by L11, L14 and H15 in the two panels. The FIRE sample at
$z{}={}4-8$ shows systematically lower $L_{\rm[C_{II}]}{}/{}\rm SFR$ ratio
than the local SFGs ($\rm[C_{II}]$ deficit). The observed galaxy samples at
$z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}5$ show similar $\rm[C_{II}]$ deficit
if $T_{\rm eqv}$ follows equation (7) (assuming $\delta_{\rm dzr}{}={}0.4$).
Table 5: Properties of the star-forming galaxies at $z\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}5$ targeted for search for $\rm[C_{II}]$ emission.

Name$\dagger$ | $z$ |  $\rm SFR_{\rm UV}$$\S$, $\\#$ ($M_{\odot}\,\rm yr^{-1}$) |  $S$ ($\rm\mu Jy$)$\ddagger$, $\P$, $\\#$ |  ${\rm log}\,(L_{\rm IR}/L_{\odot})$$\parallel$ |  $\rm SFR$$\dagger\dagger$ ($M_{\odot}\,\rm yr^{-1}$) |  ${\rm log}\,(L_{\rm[C_{II}]}/L_{\odot})$$\P$, $\\#$ | $\mu$ |  References $\ast$
---|---|---|---|---|---|---|---|---
HZ7 |  5.253 |  31.2 |  $<108$ (ALMA 7) |  $<11.6$ |  $<62.7$ |  8.74 (ALMA 7) | $-$ |  [1, 2, 3]
HZ9 |  5.541 |  22.1 |  516 (ALMA 7) |  11.9 |  174.5 |  9.21 (ALMA 7) | $-$ |  [1, 2, 3]
HZ10 |  5.657 |  58.2 |  1261 (ALMA 7) |  12.7 |  432.8 |  9.13 (ALMA 7) | $-$ |  [1, 2, 3]
NB816-S-61269 |  5.684 |  19.9 |  $<66$ (ALMA 7) |  $<11.4$ |  $<39.6$ |  8.32 (ALMA 7) | $-$ |  [4, 5]
WMH13 |  5.985 |  87.1 |  $<48$ (ALMA 6) |  $<11.7$ |  $<131.3$ |  8.56 (ALMA 6) | $-$ |  [4, 5]
A383-5.1 |  6.029 |  3.5 |  $<2.9$ (ALMA 6) |  $<10.5$ |  $<6.2$ |  6.95 (ALMA 6) |  $11.4\pm 1.9$ |  [6]
J1211-0118 |  6.029 |  55.2 |  220 (ALMA 6) |  12.4 |  257.3 |  9.15 (ALMA 6) | $-$ |  [7]
WMH5 |  6.070 |  63.2 |  218 (ALMA 6) |  12.4 |  263.5 |  8.82 (ALMA 6) | $-$ |  [9, 10]
NTTDF2313 |  6.07 |  18.4 |  $<54$ (ALMA 6) |  $<11.8$ |  $<68.0$ |  $<7.7$ (ALMA 6) | $-$ |  [8]
RXCJ0600-z6 |  6.0719 |  2.8 |  $9.5$ (ALMA 6) |  11.0 |  11.5 |  8.04 (ALMA 6) | $21\pm 10$ |  [11]
J0235-0532 |  6.089 |  58.4 |  $<101$ (ALMA 6) |  $<12.1$ |  $<150.5$ |  8.63 (ALMA 6) | $-$ |  [7]
BDF2203 |  6.12 |  24.2 |  $<69$ (ALMA 6) |  $<11.9$ |  $<87.6$ |  8.1 (ALMA 6) | $-$ |  [8]
CLM1 |  6.166 |  56.0 |  40 (ALMA 6) |  11.7 |  92.9 |  8.33 (ALMA 6) | $1.13$ |  [4, 9]
J0217-0208 |  6.203 |  86.6 |  239 (ALMA 6) |  12.4 |  307.3 |  9.15 (ALMA 6) | $-$ |  [7]
GOODS3203 |  6.27 |  27.2 |  $<123$ (ALMA 6) |  $<12.2$ |  $<140.4$ |  $<8.1$ (ALMA 6) | $-$ |  [8]
COSMOS20521 |  6.36 |  20.2 |  $<60$ (ALMA 6) |  $<11.8$ |  $<75.5$ |  $<7.7$ (ALMA 6) | $-$ |  [8]
VR7 |  6.529 |  58.2 |  $<31.8$ (ALMA 6) |  $<11.6$ |  $<87.5$ |  8.68 (ALMA 6) | $-$ |  [12]
MASOSA |  6.543 |  13.0 |  $<27.6$ (ALMA 6) |  $<11.5$ |  $<35.5$ |  $<7.34$ (ALMA 6) | $-$ |  [12]
HCM6A |  6.56 |  5.9 |  $<680$ (PdBI) |  $<12.9$ |  $<631.1$ |  $<7.81$ (PdBI) | $4.5$ |  [13, 14]
UDS4812 |  6.561 |  19.3 |  $<72$ (ALMA 6) |  $<11.9$ |  $<85.7$ |  $<7.8$ (ALMA 6) | $-$ |  [8]
Himiko |  6.591 |  19.8 |  $<27$ (ALMA 6) |  $<11.5$ |  $<44.8$ |  8.08 (ALMA 6) | $-$ |  [15, 16]
CR7 |  6.600 |  41.7 |  $<21$ (ALMA 6) |  $<11.4$ |  $<61.1$ |  8.34 (ALMA 6) | $-$ |  [17, 18]
COSMOS24108 |  6.629 |  25.6 |  $<54$ (ALMA 6) |  $<11.8$ |  $<68.2$ |  8.04 (ALMA 6) | $-$ |  [19]
UDS16291 |  6.638 |  13.4 |  $<60$ (ALMA 6) |  $<11.8$ |  $<65.4$ |  7.85 (ALMA 6) | $-$ |  [19]
NTTDF6345 |  6.701 |  21.2 |  $<48$ (ALMA 6) |  $<11.7$ |  $<60.2$ |  8.26 (ALMA 6) | $-$ |  [19]
MS0451-H |  6.703 |  0.4 |  $<0.33$ (ALMA 6) |  $<9.6$ |  $<0.7$ |  $<5.48$ (ALMA 6) | $100\pm 20$ |  [6]
UVISTA-Z-007 |  6.7496 |  23.7 |  $<52.2$ (ALMA 6) |  $<11.8$ |  $<72.0$ |  8.75 (ALMA 6) | $-$ |  [20, 21]
UVISTA-Z-019 |  6.7534 |  15.8 |  66 (ALMA 6) |  11.9 |  74.1 |  8.94 (ALMA 6) | $-$ |  [20, 21]
RXJ1347-1216 |  6.766 |  2.4 |  $<45$ (ALMA 6) |  $<11.7$ |  $<44.8$ |  7.18 (ALMA 6) | $5.0\pm 0.3$ |  [22]
COS-2987030247 |  6.808 |  24.6 |  $<75$ (ALMA 6) |  $<11.9$ |  $<94.3$ |  8.56 (ALMA 6) | $-$ |  [23]
A1703-zD1 |  6.827 |  10.1 |  $<24.5$ (NOEMA) |  $<11.5$ |  $<32.8$ |  7.54 (NOEMA) | $9.0\pm 2.7$ |  [24]
SDF-46975 |  6.844 |  15.4 |  $<57.6$ (ALMA 6) |  $<11.8$ |  $<68.7$ |  $<7.75$ (ALMA 6) | $-$ |  [25]
COS-3018555981 |  6.854 |  20.8 |  $<87$ (ALMA 6) |  $<12.0$ |  $<101.3$ |  8.67 (ALMA 6) | $-$ |  [23]
UVISTA-Z-009 |  6.86 |  16.9 |  $<38.0$ (ALMA 6) |  $<11.6$ |  $<52.1$ |  $<8.12$ (ALMA 6) | $\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}{}1.5$ |  [20, 21]
IOK-1 |  6.965 |  20.0 |  $<63$ (ALMA 6) |  $<11.9$ |  $<78.4$ |  $<7.53$ (ALMA 6) | $-$ |  [26]
BDF-512 |  7.008 |  6.0 |  $<55.2$ (ALMA 6) |  $<11.8$ |  $<54.2$ |  $<7.78$ (ALMA 6) | $-$ |  [25]
UVISTA-Z-013 |  7.02 |  22.1 |  $<45.0$ (ALMA 6) |  $<11.7$ |  $<63.8$ |  $<8.30$ (ALMA 6) | $-$ |  [20, 21]
UVISTA-Z-001 |  7.0599 |  45.8 |  104 (ALMA 6) |  12.1 |  137.8 |  8.83 (ALMA 6) | $-$ |  [20, 21]
UVISTA-Z-010 |  7.06 |  17.4 |  $<44.1$ (ALMA 6) |  $<11.7$ |  $<58.3$ |  $<8.30$ (ALMA 6) | $-$ |  [20, 21]
BDF-3299 |  7.109 |  5.7 |  $<23.4$ (ALMA 6) |  $<11.4$ |  $<27.4$ |  7.83 (ALMA 6) | $-$ |  [25, 27, 28]
A1689-zD1 |  7.137 |  4.7 |  60.2 (ALMA 6) |  11.9 |  67.5 |  7.87 (ALMA 6) | $9.3$ |  [29, 30, 31]
COSMOS13679 |  7.145 |  21.1 |  $<42$ (ALMA 6) |  $<11.7$ |  $<60.1$ |  7.85 (ALMA 6) | $-$ |  [19]
B14-65666 |  7.152 |  50.2 |  130 (ALMA 6) |  12.2 |  170.2 |  9.12 (ALMA 6) | $-$ |  [32, 33]
SXDF-NB1006-2 |  7.212 |  21.6 |  $<42$ (ALMA 6) |  $<11.7$ |  $<60.6$ |  $<7.45$ (ALMA 6) | $-$ |  [34]
z8-GND-5296 |  7.508 |  16.6 |  $<480$ (PdBI) |  $<12.7$ |  $<464.1$ |  $<8.55$ (PdBI) | $-$ |  [35, 36]
MACS0416-Y1 |  8.311 |  11.7 |  137 (ALMA 7) |  11.8 |  56.8 |  8.15 (ALMA 5) |  $1.43\pm 0.04$ |  [37, 38, 39]
A2744-YD4 |  8.380 |  11.2 |  99 (ALMA 7) |  11.6 |  43.8 |  7.26 (ALMA 5) | $1.8\pm 0.3$ |  [28, 40, 41]
S04590 |  8.4931 |  0.5 |  $<{}4.81$ (ALMA 7) |  $<10.3$ |  $<2.0$ |  7.22 (ALMA 5) | $8.69\pm 2.5$ |  [42, 43]
MACS1149-JD1 |  9.110 |  4.5 |  $<5.3$ (ALMA 7) |  $<10.4$ |  $<6.5$ |  7.08 (ALMA 5) | $10$ |  [28, 41, 44]
REBELS$\ddagger\ddagger$
REBELS-05 |  6.496 |  15.1 |  67.2 (ALMA 6) |  11.9 |  77.2 |  8.84 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-38 |  6.577 |  19.5 |  163.0 (ALMA 6) |  12.3 |  170.2 |  9.23 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-29 |  6.685 |  27.0 |  56.1 (ALMA 6) |  11.8 |  78.9 |  8.74 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-32 |  6.729 |  15.1 |  60.4 (ALMA 6) |  11.8 |  71.0 |  8.89 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-08 |  6.749 |  17.3 |  101.4 (ALMA 6) |  12.1 |  111.2 |  8.87 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-39 |  6.847 |  40.0 |  79.6 (ALMA 6) |  12.0 |  113.7 |  8.90 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-14 |  7.084 |  37.9 |  60.0 (ALMA 6) |  11.8 |  93.6 |  8.57 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-27 |  7.090 |  21.6 |  50.6 (ALMA 6) |  11.8 |  68.5 |  8.79 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-25 |  7.306 |  16.2 |  56.1 (ALMA 6) |  11.8 |  68.3 |  9.20 (ALMA 6) | $-$ |  [45, 46, 47]
REBELS-12 |  7.349 |  32.5 |  86.8 (ALMA 6) |  12.0 |  113.2 |  9.00 (ALMA 6) | $-$ |  [45, 46, 47]

(Continue on next page)

Table 6: continued

Name$\dagger$ | $z$ |  $\rm SFR_{\rm UV}$$\S$, $\\#$ ($M_{\odot}\,\rm yr^{-1}$) |  $S$ ($\rm\mu Jy$)$\ddagger$, $\P$, $\\#$ |  ${\rm log}\,(L_{\rm IR}/L_{\odot})$$\parallel$ |  $\rm SFR$$\dagger\dagger$ ($M_{\odot}\,\rm yr^{-1}$) |  ${\rm log}\,(L_{\rm[C_{II}]}/L_{\odot})$$\P$, $\\#$ | $\mu$ |  References $\ast$
---|---|---|---|---|---|---|---|---
REBELS-40 |  7.365 |  18.4 |  48.3 (ALMA 6) |  11.8 |  64.5 |  8.69 (ALMA 6) | $-$ |  [43, 44, 45]
REBELS-19 |  7.369 |  15.1 |  71.2 (ALMA 6) |  11.9 |  81.3 |  8.94 (ALMA 6) | $-$ |  [43, 44, 45]
REBELS-18 |  7.675 |  33.5 |  52.9 (ALMA 6) |  11.8 |  82.8 |  9.03 (ALMA 6) | $-$ |  [43, 44, 45]

  * $\dagger$

The table does not include the 118 galaxies ($4.5\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 60\relax$}\hss}\mathchar
29208\relax$}}z\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 60\relax$}\hss}\mathchar 29208\relax$}}6$) selected by the
ALPINE project. The information of the ALPINE galaxies can be downloaded from
the official webpage of the project:
https://cesam.lam.fr/a2c2s/data_release.php. The ALPINE galaxies are unlensed.

  * $\S$

$\rm SFR_{\rm UV}$ is converted from $L_{\rm UV}$ via ${\rm
SFR_{UV}}\,(M_{\odot}\,{\rm yr^{-1})}{}={}1.58\times 10^{-10}\,L_{\rm
UV}\,(L_{\odot})$ following Hao et al. (2011) (see Table 3) for the Kroupa
(2002) IMF.

  * $\ddagger$

The number in the brackets indicates the specific ALMA band at which dust
continuum is measured.

  * $\P$

For the galaxies having no detection of dust thermal continuum ($\rm[C_{II}]$
emission), we show the $3\sigma$ upper confidence limit of $S$
($L_{\rm[C_{II}]}$).

  * $\\#$

For the gravitationally-lensed galaxies, $L_{\rm UV}$ (and hence $\rm
SFR_{UV}$), $S$, $L_{\rm IR}$ and $L_{\rm[C_{II}]}$ are de-magnified by $\mu$.

  * $\parallel$

$L_{\rm IR}$ (or the upper limit of $L_{\rm IR}$ for the dust-undetected
sources) is converted from $S$ (the $3\sigma$ upper limit of $S$) via the
standard MBB function with $T_{\rm eqv}$ calculated by equation (4) (assuming
$\beta_{\rm dust}{}={}2.0$ and $\delta_{\rm dzr}{}={}0.4$).

  * $\dagger\dagger$

SFR is derived using ${\rm SFR}\,(M_{\odot}\,{\rm yr^{-1})}{}={}{\rm
SFR_{UV}}+{\rm SFR_{IR}}{}={}1.58\times 10^{-10}\,(L_{\rm UV}+0.46L_{\rm
IR})\,(L_{\odot})$ following Hao et al. (2011) (see Table 3) for the Kroupa
(2002) IMF.

  * $\ddagger\ddagger$

We only list here the 13 galaxies of the REBELS sample that have confirmed
detection of both $\rm[C_{II}]$ and dust continuum. The information of the
other 5 galaxies having $\rm[C_{II}]$ but no dust detection is not yet
publicly available.

  * $\ast$

References: (1): Capak et al. (2015), [2]: Barisic et al. (2017), [3]: Faisst
et al. (2017), [4]: Fujimoto et al. (2019), [5]: Fujimoto et al. (2016), [6]:
Knudsen et al. (2016), [7]: Harikane et al. (2020), [8]: Carniani et al.
(2018a), [9]: Willott et al. (2015b), [10]: Willott et al. (2013a), [11]:
Fujimoto et al. (2021), [12]: Matthee et al. (2019), [13]: Kanekar et al.
(2013), [14]: Hu et al. (2002), [15]: Ouchi et al. (2013), [16]: Carniani et
al. (2018b), [17]: Sobral et al. (2015), [18]: Matthee et al. (2017), [19]:
Pentericci et al. (2016), [20]: Schouws et al. (2022a), [21]: Schouws et al.
(2022b), [22]: Bradač et al. (2017), [23]: Smit et al. (2018), [24]: Molyneux
et al. (2022), [25]: Maiolino et al. (2015), [26]: Ota et al. (2014), [27]:
Carniani et al. (2017), [28]: Carniani et al. (2020), [29]: Watson et al.
(2015), [30]: Knudsen et al. (2017), [31]: Wong et al. (2022), [32]: Hashimoto
et al. (2019a), [33]: Bowler et al. (2018), [34]: Inoue et al. (2016), [35]:
Schaerer et al. (2015a), [36]: Finkelstein et al. (2013), [37]: Tamura et al.
(2019), [38]: Bakx et al. (2020), [39]: Kawamata et al. (2016), [40]: Laporte
et al. (2017), [41]: Laporte et al. (2019), [42]: Fujimoto et al. (2022),
[43]: Heintz et al. (2022b), [44]: Hashimoto et al. (2018), [45]: Ferrara et
al. (2022), [46]: Sommovigo et al. (2022), [47]: Bouwens et al. (2022).

###  4.3 Early galaxies (redshift $z{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$)

Table 7: Comparison between the mean ‘equivalent dust temperature’ (<$T_{\rm
eqv}$>) assumed by the ALPINE and REBELS projects and by this work.

Project name | Reference | No. of galaxies | <$z$> |  <$T_{\rm eqv}/\rm K$> |  <$T_{\rm eqv}/\rm K$>$\dagger$ |  $\Delta$<${\rm log}(\frac{L_{\rm[C_{II}]}}{\rm SFR}$)>$\ddagger$
---|---|---|---|---|---|---
|  |  |  |  (literature) |  (this work) |  (dex)
ALPINE |  Schaerer et al. (2020) | 118 | 4.58 | 42 | 51.3 | -0.27
REBELS |  Ferrara et al. (2022) | 40 | 7.08 | 55 | 57.4 | -0.12

  * $\dagger$

Calculated using equation (7) with $\delta_{\rm dzr}{}={}0.4$. Note that with
a lower $\delta_{\rm dzr}$, $T_{\rm eqv}$ is higher than the listed value in
this column.

  * $\ddagger$

The resulting difference in the derived mean $L_{\rm[C_{II}]}{}/{}\rm SFR$
ratio (in dex) of the galaxy samples due to the difference in $T_{\rm eqv}$
used by the previous studies (Schaerer et al. 2020 and Ferrara et al. 2019)
and this work.

Observational studies on the $L_{\rm[C_{II}]}$-SFR relation at
$z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char
62\relax$}\hss}\mathchar 29208\relax$}}{}5$ depend mainly on the rest-frame
UV-selected galaxies whose redshift has previously been confirmed either
spectroscopically or via the Lyman break ‘drop-out’ technique (Hodge & da
Cunha, 2020). Their $\rm[C_{II}]$ and dust emission are constrained in follow-
up observational campaigns with ALMA, which has the power to spatially resolve
the distant galaxies down to the scale of $\sim 1$ physical kpc. The majority
of the UV-selected galaxies at this epoch are unlensed.

There have been two major observational campaigns for searching for
$\rm[C_{II}]$ line of galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox
to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$.
The ALPINE ALMA Large Program (Le Fèvre et al., 2020; Béthermin et al., 2020)
in cycle-5 targeted a sample of 118 UV-selected star-forming galaxies at
$4.5{}<{}z{}<{}6$ with $M_{\rm UV,\,AB}{}<{}-20.2$ and identified
$\rm[C_{II}]$ emission (at $>{}3.5\sigma$ level) in 75 galaxies of them
(Schaerer et al., 2020). More recently, the REBELS Large Program (Bouwens et
al., 2022) in Cycle-7 studied a sample of 40 UV-bright ($M_{\rm
UV,\,AB}<-21.4$) galaxies at $6.5<{}z{}<{}7.7$ and confirmed $\rm[C_{II}]$
detection (at $>{}7\sigma$) in 18 galaxies in their sample (Ferrara et al.,
2022). Other observations targeting the LBGs/LAEs at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ have identified $\rm[C_{II}]$ emission in another $>{}35$
sources in total. The most distant galaxy that has a $\rm[C_{II}]$ detection
to date is MACS1149-JD1 (Hashimoto et al., 2018), a gravitationally-lensed
($\mu{}={}10$) galaxy at $z{}={}9.11$ (Carniani et al. 2020; see also Inoue et
al. 2016 and Laporte et al. 2019). We provide a summary of the star-forming
galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$ having confirmed
$\rm[C_{II}]$ detection in Table 5 (excluding quasar host galaxies).

The SFR of these UV-selected galaxies has been derived based on their $L_{\rm
UV}$ and $L_{\rm IR}$. Because the galaxies at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$ typically do not have good photometric sampling of the dust
continuum (e.g. Casey et al., 2018b; Liang et al., 2019; Faisst et al.,
2020b), $L_{\rm IR}$ has frequently been converted from the ALMA broad-band
flux density (measured at band 6 or 7 for galaxies at $z{}\mathbin{\lower
3.0pt\hbox{$\hbox to0.0pt{\raise 4.0pt\hbox{$\char 62\relax$}\hss}\mathchar
29208\relax$}}{}5$) using the standard modified-blackbody (MBB) function of
the form (e.g. Hildebrand, 1983; Hayward et al., 2011)

$S_{\nu_{0}}=\frac{(1+z)}{d^{2}_{\rm L}}\kappa_{\nu}M_{\rm dust}B_{\nu}(T),$
(5)

where $\nu_{0}$ is the observing frequency (note: $\nu_{0}{}={}345$ GHz for
ALMA band 7 and $\nu_{0}{}={}230$ GHz for ALMA band 6), $S_{\nu_{0}}$ is the
broad-band flux density at $\nu_{0}$, $\nu{}={}(1+z)\nu_{0}$ is the rest-frame
frequency, $\kappa_{\nu}$ is the dust opacity (per unit dust mass) at $\nu$,
$M_{\rm dust}$ is the dust mass of galaxy, $T$ is the ‘dust temperature’,
$B_{\nu}(T)$ is the Planck function and $d_{\rm L}$ is the luminosity
distance. $L_{\rm IR}$ is then converted from $S_{\nu_{0}}$ using (see Section
3.1.3 of Liang et al. 2019 for the details)

$L_{\rm IR}=\frac{\mathcal{D}d^{2}_{\rm L}T^{4+\beta_{\rm
dust}}}{(1+z)\kappa_{\nu}B_{\nu}(T)}S_{\nu_{0}},$ (6)

where $\beta_{\rm dust}{}\approx{}2.0$ is the dust emissivity spectral index
(e.g. Dunne et al., 2000; Draine et al., 2007) and $\mathcal{D}$ is a
parameter that depends on the shape of the dust opacity curve. The derived
$L_{\rm IR}$ (and hence the obscured SFR) therefore depends mainly on the
assumed ‘dust temperature’. It should be noted that recent cosmological
simulations show that the true SED of high-$z$ galaxies may significantly
differ from the standard MBB function (e.g. Liang et al. 2019; Ma et al. 2019,
and also Casey 2012, Casey et al. 2018b) and $T$ does not faithfully reflect
the physical temperature of dust in galaxies (e.g. Behrens et al., 2018; Liang
et al., 2019). Liang et al. (2019) defines the ‘dust temperature’ that one
would need to obtain the correct $L_{\rm IR}$ and match the observed
$S_{\nu_{0}}$ under the assumption that the SED has the shape of a standard
MBB function (equation 5) to be the ‘equivalent dust temperature’ ($T_{\rm
eqv}$).

Using a sample of high-$z$ galaxies produced by the MassiveFIRE suite
(Feldmann et al., 2016; Feldmann et al., 2017), Liang et al. (2019) derived
the best-fitting formula for $T_{\rm eqv}$ using redshift and dust-to-gas mass
ratio ($\delta_{\rm dzr}$) as variables, i.e.

$T_{\rm eqv}=T_{0}\,(1+z)^{\alpha}(\delta_{\rm
dzr}/0.4)^{\gamma}.\;\;\rm(L19)$ (7)

For ALMA band 7 (6) fluxes, the best-fitting parameter values are
$T_{0}{}={}26.9$ (24.5) K, $\alpha{}={}0.31$ (0.36) and $\gamma{}={}-0.13$
($-0.15$). The increase of $T_{\rm eqv}$ with redshift is related to the
enhanced level of star formation activity in galaxies (i.e. higher specific
SFR) (Safarzadeh et al., 2016; Ma et al., 2019; Liang et al., 2019; Sommovigo
et al., 2020). The anti-correlation with $\delta_{\rm dzr}$, on the other
hand, is due to the fact that an increase of $\delta_{\rm dzr}$ leads to a
higher dust opacity, which in turn results in a ‘colder’ dust SED shape of
galaxies (Scoville, 2013; Faisst et al., 2017; Liang et al., 2019).
Observationally, $\delta_{\rm dzr}$ of high-$z$ galaxies has not yet been
constrained.

Often, it is easier to detect the $\rm[C_{II}]$ line than the dust continuum
of galaxies at $z{}\mathbin{\lower 3.0pt\hbox{$\hbox to0.0pt{\raise
4.0pt\hbox{$\char 62\relax$}\hss}\mathchar 29208\relax$}}{}5$. For example, 75
out of the 118 ($63.6\%$) galaxies in the ALPINE sample have confirmed
detection of $\rm[C_{II}]$ emission, whilst only 21 ($17.8\%$) of them have
confirmed detection of dust continuum. Almost all dust-detected galaxies have
detection of $\rm[C_{II}]$ line. The detection limit of $\rm[C_{II}]$ of the
current ALMA observations is about $10^{8}\,L_{\odot}$.

We convert the sub-mm broad-band flux density ($S_{\nu_{\rm 0}}$) of the dust-
detected galaxies (or the $3\sigma$ upper limit of $S_{\nu_{0}}$ for the dust-
undetected galaxies) to $L_{\rm IR}$ (the upper limit of $L_{\rm IR}$)
consistently using $T_{\rm eqv}$ that follows equation (7) (assuming
$\delta_{\rm dzr}{}={}0.4$) to make a fair comparison between different
observed samples and our theoretical predictions using FIRE galaxies. We
compute the SFR of the observed galaxies using their measured $L_{\rm UV}$ and
the derived $L_{\rm IR}$ following Hao et al. (2011), i.e. ${\rm
SFR}\,(M_{\odot}\,{\rm yr^{-1})}{}={}1.58\times 10^{-10}\,(L_{\rm
UV}+0.46L_{\rm IR})\,(L_{\odot})$, for the Kroupa (2002) IMF. For the dust-
undetected galaxies, we estimate the lower and upper bounds of their SFR,
where the former is converted from their $L_{\rm UV}$ assuming no dust
emission (i.e. $L_{\rm IR}{}={}0$), whilst the latter accounts for the upper
limit of $L_{\rm IR}$ converted from the $3\sigma$ upper limit of
$S_{\nu_{0}}$.