# Atomic and molecular gas as traced by [C II] emission

Benedetta Casavecchia, Umberto Maio, Céline Péroux, Benedetta Ciardi

► To cite this version:

Benedetta Casavecchia, Umberto Maio, Céline Péroux, Benedetta Ciardi. Atomic and molecular gas as traced by [C II] emission. Astronomy & Astrophysics - A&A, 2025, 693, pp.A119. 10.1051/0004-6361/202452282. hal-05469669

HAL Id: hal-05469669

https://hal.science/hal-05469669v1

Submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Atomic and molecular gas as traced by [C II] emission

Benedetta Casavecchia<fnref n="1" />,*, D, Umberto Maio<fnref n="2" />,<fnref n="3" />, Céline Péroux<fnref n="4" />,<fnref n="5" /> D, and Benedetta Ciardi<fnref n="1" />

Received 17 September 2024 / Accepted 2 December 2024

## ABSTRACT

The latest ALMA and JWST observations provide new information on the birth and evolution of galaxies in the early Universe at the epoch of reionization. Measurements at redshift $z > 5$ of their cold-gas budget are particularly important because this budget is known to be the main fuel for star formation. A powerful tool for probing the physics characterising galaxies at high redshift is the [C II] 158 µm emission line. Due to its low excitation potential, [C II] emission can be produced in photodissociation regions, neutral atomic gas, and molecular clouds. To properly capture the cold-gas processes taking place in these environments (molecule formation, self-shielding, dust grain catalysis, and photoelectric and cosmic-ray heating), we made use of a new set of dedicated hydrodynamic simulations (COLDSIM) including time-dependent non-equilibrium chemistry, star formation, stellar evolution, metal spreading, and feedback mechanisms. We were able to accurately track the evolution of H I, H II and H$_2$ in a cosmological context and predict the contribution of each gas phase to [C II] luminosity. We provide formulas that can be used to estimate the mass of molecular and atomic gas from [C II] detections. Furthermore, we analysed the evolution of conversion factors with galactic properties, such as stellar metallicity, star formation rate, and stellar mass. We demonstrate that [C II] emission is dominated by H I gas and that most of the [C II] luminosity is generated in warm, dense, star-forming regions. We conclude that although [C II] predominantly traces atomic rather than molecular gas, the [C II] luminosity remains a robust indicator of the H$_2$ mass.

**Key words.** galaxies: abundances – galaxies: evolution – galaxies: formation – galaxies: high-redshift – galaxies: star formation

## 1. Introduction

The molecular phase is mainly composed of H$_2$ and is most important for the physical processes that convert gas into stars. This is traced by the history of the star formation rate (SFR). Neutral hydrogen provides the essential fuel, but this fuel has to cool and transform into the molecular phase for star formation. These processes are often referred to as the baryon cycle (Péroux & Howk 2020).

H$_2$ detections are extremely challenging because H$_2$ lacks a permanent dipole moment and because its lowest rotational transition requires temperatures of 510 K at least for excitation (Saslaw & Zipoy 1967; Hollenbach & McKee 1979; Lepp & Shull 1983; Galli & Palla 1998; Roussel et al. 2007; Fukui & Kawamura 2010; Dobbs & Pringle 2013; Krumholz 2014; Togi & Smith 2016). Hence, direct detections of molecular hydrogen (see e.g. recent JWST determinations by Hernandez et al. 2023) trace a small portion of the warmer molecular gas, but not the majority of the cold gas (at temperatures $T < 10^2$ K). In this context, most estimates of molecular gas masses relied on CO rotational transitions, which enable strong and easily observable emission lines.

Through the detection of CO emission, it is possible to constrain the molecular content of large samples of galaxies (Saintonge et al. 2011; Bolatto et al. 2013; Andreani et al. 2020; Tacconi et al. 2020; Boogaard et al. 2023; Salvestrini et al. 2024). Unbiased surveys for molecular gas have been con-

ducted with the Atacama Large Millimeter/submillimeter Array (ALMA) (Decarli et al. 2019, 2020). These observations provided an assessment of the cosmic evolution of molecular gas up to redshift $z \sim 4$ (Riechers et al. 2019; Decarli et al. 2019; Walter et al. 2022; Aravena et al. 2024). The Very Large Array (VLA) based CO Luminosity Density at High-$z$ (COLDz) survey made complementary measurements at higher $z$ (Riechers et al. 2019), while the IRAM 30 m-based xCOLD GASS provided robust $z = 0$ measurements (Saintonge et al. 2017) and ALMACAL low-redshift measurements over a larger volume (Hamanowicz et al. 2023).

However, CO is often an unreliable tracer of H$_2$ in diffuse and low-metallicity environments (see e.g. Ebagezio et al. 2023, for recent studies). In these cases, the reduced dust content allows far-ultraviolet (UV) radiation from star-forming regions to penetrate deeply into molecular clouds and to photodissociate CO molecules. This creates a layer of ionised carbon (C$^+$) around a central CO core. In contrast, H$_2$ molecules become self-shielded from photodissociation and can survive within C$^+$ layers and in the CO core (Krumholz & Gnedin 2011; Glover & Clark 2012; Gnedin & Draine 2014). Moreover, at higher redshift, CO transitions become progressively more challenging to observe at higher frequencies (see e.g. Feruglio et al. 2023 for a broader discussion). Moreover, converting the observed CO transitions with higher rotational quantum numbers into the reference CO(1-0) introduces an uncertainty that arises from variations in the poorly understood spectral line energy distribution across different galaxy types (Carilli & Walter 2013; Klitsch et al. 2019). Furthermore, the

<footnote n="1">¹ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching b. München, Germany</footnote>

<footnote n="2">² INAF-Italian National Institute of Astrophysics, Observatory of Trieste, via G. Tiepolo 11, 34143 Trieste, Italy</footnote>

<footnote n="3">³ IFPU-Institute for Fundamental Physics of the Universe, via Beirut 2, 34014 Trieste, Italy</footnote>

<footnote n="4">⁴ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching b. München, Germany</footnote>

<footnote n="5">⁵ Aix Marseille Université, CNRS, Laboratoire d'Astrophysique de Marseille (LAM) UMR 7326, 13388 Marseille, France</footnote>

<footnote>* Corresponding author; benecasa@mpa-garching.mpg.de</footnote>
$L_{\text{CO}}$-to-$M_{\text{H}_2}$ conversion factor introduces a significant systematic uncertainty, and CO may no longer serve as a reliable tracer of $\text{H}_2$ mass in extreme environments. Specifically, this conversion factor increases as the metallicity decreases because CO photodissociates to a larger depth in each cloud. Consequently, it varies not only between objects and with redshift (Bolatto et al. 2013), but also in different regions. Alternative indicators of molecular mass are often used. Dunne et al. (2022) exhaustively compared the molecular gas-mass estimates from various tracers, namely the dust continuum detection and observations of CO and CI emission lines. The authors concluded that the various indicators agree well within the uncertainties. The question remains whether some or all of these indicators trace the whole molecular gas, including the more diffuse component of the gas. While $\text{H}_2$ is typically traced by CO observations, Madden et al. (2020) demonstrated that in low-metallicity galaxies, gas is not self-shielded from UV photons, and this causes the formation of a layer of [C II]-emitting gas that may not be detectable through CO observations, but might still host $\text{H}_2$. The authors used a photoionisation CLOUDY modelling to emphasise that these CO-dark $\text{H}_2$ regions can be traced by [C II] or C I in principle. Their findings imply that the [C II] 158 µm line emission is an excellent tracer of the total molecular gas mass in galaxies, and this suggests that a significant reservoir of CO-dark molecular mass might have been overlooked by previous empirical determinations. This highlights the importance of the tracer used to accurately measure molecular gas in galaxies.

The fine-structure emission line of $\text{C}^+$ at a wavelength of ~158 µm ([C II] 158 µm) has proven to be a robust tracer of gas in the early Universe, in part because of its coverage by ALMA at $z > 4$. Because its excitation potential is low, [C II] emission can be produced in gas with temperatures below ~$10^4$ K, such as neutral atomic gas and molecular gas clouds. Three ALMA large programs in particular have observed selected UV-bright galaxies at $z \sim 4-8$. The ALPINE survey (Le Fèvre et al. 2020) has reported 118 normal star-forming galaxies at $z = 4.5-6$, and 20 of these galaxies and 6 others were re-observed within the CRISTAL program (Mitsuhashi et al. 2024; Ikeda et al. 2024; Li et al. 2024; Telikova et al. 2024, Herrera-Camus et al. in prep.) in order to spatially resolve the [C II] line and the far-infrared dust continuum emissions. The REBELS survey (Bouwens et al. 2022) characterised dozens of the most luminous galaxies at $z = 6.5-7.7$. At these redshifts, direct detections of molecular gas are challenging, and $\text{H}_2$ mass ($M_{\text{H}_2}$) estimates typically rely on empirical relations that link molecular gas to the SFR (Kennicutt 1998; Tacconi et al. 2013), infrared luminosity (Leroy et al. 2008; Aravena et al. 2016), dust content (Rémy-Ruyer et al. 2014; Genzel et al. 2015; Tacconi et al. 2018), and dynamical mass (Daddi et al. 2010).

All these relations are empirically well calibrated in the local Universe, but are not reliable for high-redshift systems. For this reason, a physically motivated conversion factor, $\alpha_{[C\text{II}]}$, that connects the [C II] luminosity ($L_{[C\text{II}]}$) to the $\text{H}_2$ mass is extremely important for an independent estimate of the molecular gas mass at the epoch of reionisation. Recent determinations of $\alpha_{[C\text{II}]}$ were provided for dwarf galaxies in the local Universe (Madden et al. 2020) and for star-forming galaxies up to $z = 2$ (Zanella et al. 2018). The $\alpha_{[C\text{II}]}$ from Zanella et al. (2018) was frequently applied to $z = 4-8$ galaxies (Béthermin et al. 2020; Schaerer et al. 2020; Tacconi et al. 2020; Dessauges-Zavadsky et al. 2020; Aravena et al. 2024; Kaasinen et al. 2024), but it is unclear whether [C II] detections can be solid tracers of molecular gas at these high redshifts. A well-calibrated method for inferring $M_{\text{H}_2}$ from $L_{[C\text{II}]}$ is crucial

for current and future observations. This bright emission line is the target of new facilities, including the CONCERTO instrument (CONCERTO Collaboration 2020; Gkogkou et al. 2023) on the Atacama Pathfinder Experiment telescope (APEX), and future projects, such as the proposed Atacama Large Aperture Submillimeter Telescope (AtLAST) (Klaassen et al. 2020).

Additionally, the physical conditions of the gas giving rise to $L_{[C\text{II}]}$ is still poorly understood. It is instrumental for a full interpretation of this emission line to establish the dependences on the star-forming processes and the density of the gas traced by [C II]. Most importantly, $L_{[C\text{II}]}$ can be produced by ionised, atomic, and molecular gas, and some observational efforts attempted to disentangle the [C II] emission arising from different phases in the Local Universe (e.g. Croxall et al. 2017; Cormier et al. 2019; Tarantino et al. 2021). Nevertheless, it is challenging to interpret the wealth of new observational data without a deeper understanding of the dominant phase of the gas it traces.

It is rather demanding to compute the cold phases of the gas in a cosmological context in a simulation. On one hand, the physics involved is complex and has many parameters; and on the other hand, the resulting computing times are extremely long. This challenge is often referred to as “molecular cosmology”. Numerical works have investigated the physical mechanisms that emit [C II] in galaxies in the local Universe and at the epoch of reionisation. Three main methods were commonly adopted to model the [C II] luminosity. Lagache et al. (2018) and Popping et al. (2019) implemented a semi-analytical model (SAM) based on empirical relations to determine the [C II] luminosity in galaxies from the local Universe up to $z \sim 7$ by linking their $L_{[C\text{II}]}$ to the SFR. The advantage of SAMs such as those lies in the fast calculation runtime, which enables the exploration of a large parameter space. However, the simplifications inherent to SAMs can limit the accuracy of the results. An additional method relies on zoom-in simulations of single galaxies. These were typically used to resolve small-scale processes. For instance, Vallini et al. (2015), Pallottini et al. (2017), Katz et al. (2019), Bisbas et al. (2022), and Schimek et al. (2024) investigated a few individual galaxies at $z \sim 6$ with different properties, such as virial mass, metallicity, presence or absence of mergers, and resolution. While zoom-in simulations provide high-resolution insights into specific galaxies, they are computationally more expensive and generally limited to a small number of cases. This prevents broad statistical conclusions. Finally, simulations that are run in cosmological boxes, although with lower resolution, can provide statistics on much larger samples of haloes. Ramos Padilla et al. (2021, 2023), Vizgan et al. (2022a,b) and Khatri et al. (2024b) focused on modelling the [C II] emission in post-processing and the contribution from the ionised atomic and molecular gas phases in the EAGLE, SIMBA and MARIGOLD simulations, respectively.

Sub-grid prescriptions are crucial for modelling the cold-gas phase in zoom-in and cosmological simulations. Simple empirical approaches that are based on gas metallicity or pressure (Blitz & Rosolowsky 2006; Krumholz 2013) are often used to separate neutral hydrogen into its atomic and molecular components (Lagos et al. 2015; Popping et al. 2019; Popping & Péroux 2022; Szakacs et al. 2022). Some other works employ idealised models for line emission (Olsen et al. 2018) or zoom-in simulations. Although these did not directly include [C II], they yielded valuable insights into the properties of individual haloes (Pallottini et al. 2019). Therefore, it remains one of the most critical and challenging objectives for advancing our understanding of galaxy formation and evolution at the epoch of reionisation to
accurately simulate the neutral and molecular phases of cold gas within a cosmological context.

In this context, the COLDSIM simulation suite (see Maio et al. 2022; Maio & Viel 2023; Casavecchia et al. 2024, from now on referred to as C24) proposed a novel numerical approach to investigate the physical processes underlying the $L_{[CII]}$ emission that is detected in galaxies at $z > 6$ within a cosmological framework. The suite consists of a set of cosmological hydrodynamic simulations incorporating a time-dependent atomic and molecular non-equilibrium chemical network coupled with star formation, stellar evolution, feedback effects, and cooling processes. All these features enable capturing the mechanisms that cause [C II] emission in galaxies. In C24 we investigated the amount of $C^+$ mass density that is produced in cosmic cold gas at different epochs, as well as the relations between $L_{[CII]}$, SFR, and stellar mass.

The aim of this paper is to show that there are physical relations between $L_{[CII]}$, the atomic hydrogen mass ($M_{HI}$), and the molecular hydrogen mass ($M_{H_2}$). Moreover, we investigate the evolution of the conversion factors $L_{[CII]}$ to $M_{H_2}$ ($\alpha_{[CII]}$) and $L_{[CII]}$ to $M_{HI}$ ($\beta_{[CII]}$) with galactic properties, such as stellar metallicity ($Z_*$), SFR, and stellar mass ($M_*$).

The paper is organised as follows: In Sect. 2 we outline our method, and Sect. 3 reports our findings and relative comparisons with observational works. We discuss the impacts of our results in Sect. 4 and present our main conclusions in Sect. 5.

## 2. Method

In the following, we present the numerical simulations we adopted throughout this work. We describe how we extracted the expected [C II] signal from simulated galaxies and how it can be related to H I and $H_2$ gas through the $\alpha_{[CII]}$ and $\beta_{[CII]}$ conversion factors.

### 2.1. The ColdSIM numerical simulations

For our analysis, we employed the COLDSIM High Resolution (CDM HR) and COLDSIM Large Box (CDM LB) simulations described in Maio & Viel (2023) and C24, which are part of the COLDSIM suite that was reported in detail by Maio et al. (2022). The COLDSIM simulations have an ad hoc implementation to reproduce the formation of the first structures in the Universe. Different from other approaches that were calibrated to reproduce the properties of the local Universe, COLDSIM simulations were performed with chemical and physical processes that are typical of the early Universe. This ensured a more appropriate modelling of the pristine gas in its atomic and molecular form, and of its collapse into stars and galaxies. Gravity and smoothed-particle-hydrodynamics (SPH) computations were performed with a modified version of P-GADGET3 (Springel 2005), which includes a unique time-dependent non-equilibrium chemical network that resolves ionisation, dissociation, and recombination processes of the primordial gas via first-order differential equations (Abel et al. 1997; Yoshida et al. 2003; Maio et al. 2007). In this way, we tracked the evolution of the following atomic and molecular species: $e^-$, H, $H^+$, $H^-$, He, $He^+$, $He^{++}$, $H_2$, $H_2^+$, D, $D^+$, HD, and $HeH^+$. The chemical network was evolved self-consistently with stellar feedback processes, cooling, and heating (Maio et al. 2022). At each time step, collisionless stellar particles were formed in a stochastic way following a modified version of the Springel & Hernquist (2003) model, in which we included molecular runaway cooling (Maio et al. 2009). Each of these particles represents a simple stellar population with a

Salpeter initial mass function (IMF). Population III (Pop III) stars were distributed with a Salpeter IMF, similarly to population II (Pop II) and population I (Pop I). We followed their stellar evolution according to mass and metallicity-dependent stellar yields and lifetimes. In particular, stellar particles produce heavy elements through asymptotic giant branch (AGB) winds and Type II and Type Ia supernova (SN) explosions (Tornatore et al. 2007; Maio et al. 2010, 2016), including He, C, N, O, Ne, Mg, Si, S, Ca, and Fe. Each of these elements was tracked individually and mixed through the SPH kernel smoothing. As detailed in Maio et al. (2022), radiative losses of metals from resonant and fine-structure transitions, such as the [C II] emission at 158 µm, were included in the chemical network. Primordial gas cools through $H_2$ formation via $H^-$, $H_2^+$ catalysis, and three-body interactions in pristine gas, and additionally, via dust catalysis in metal-enriched gas particles. The $H_2$ dust-grain catalysis was performed for different gas temperatures and metallicities, assuming a Z-dependent dust-to-metal ratio and an energy balance between the cosmic microwave background radiation and dust-grain emission with a power law of slope $\gamma = 2$ (in agreement with recent ALMA determinations). The standard Haardt & Madau (1996) UV background generated by the formation of the first structures was turned on at $z = 6$. Photoelectric and cosmic-ray heating were implemented in the chemical network and complemented with H I and $H_2$ shielding (Habing 1968; Draine & Sutin 1987; Bakes & Tielens 1994; Draine & Bertoldi 1996).

The simulations were performed within the standard cosmological framework based on cold dark matter and the cosmological constant $\Lambda$. The present-day expansion parameter was assumed to be $H_0 = 100 \text{ h km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.7$, while the baryon, matter, and $\Lambda$ density parameters were $\Omega_{b,0} = 0.045$, $\Omega_{m,0} = 0.274$, and $\Omega_{\Lambda,0} = 0.726$, respectively. The adopted power spectrum had a spectral index of 0.96 and a normalisation over an 8 Mpc $h^{-1}$ sphere of $\sigma_8 = 0.8$. We adopted the notation cMpc to indicate lengths in comoving megaparsecs, and unless otherwise stated, we expressed masses, luminosities, and metallicities in solar units ($M_\odot$, $L_\odot$, and $Z_\odot$, respectively).

The two simulated volumes for which we report results in this paper are CDM LB and CDM HR. The former has the largest volume and is able to capture the more massive galaxies. It has a size of $50 \text{ cMpc} \ h^{-1}$ and a number of particles of $2 \times 1000^3$ that are evenly distributed between dark matter and gas at the initial redshift $z = 100$. In CDM LB, dark-matter and gas particles have a mass of $m_{\text{DM}} = 7.9 \times 10^6 \text{ M}_\odot$ and $m_{\text{gas}} = 1.6 \times 10^6 \text{ M}_\odot$, respectively. CDM HR instead has the highest resolution because it has a size of $10 \text{ cMpc} \ h^{-1}$ and the same number of particles, resulting in $m_{\text{DM}} = 6.3 \times 10^4 \text{ M}_\odot$ and $m_{\text{gas}} = 1.3 \times 10^4 \text{ M}_\odot$. We show the results at $z = 6, 7, 8$, and 10 in order to make comparisons with ALMA [C II] observations at the epoch of reionisation and provide predictions for higher redshifts.

### 2.2. Computing [C II] luminosities and conversion factors

In this section, we describe how we estimated the [C II] emission at 158 µm for haloes in COLDSIM. The calculations for [C II] luminosities were directly included in the simulations at runtime and contributed to cosmic gas cooling. This guaranteed that thermodynamical variables were evaluated correctly and the [C II] emission was derived consistently. More detailed explanations were provided in C24, and we limit our presentation here to the most relevant methodological steps.

We considered haloes composed of at least 100 particles, and for each halo, $L_{[CII]}$ was calculated as the sum of the [C II]
<figure><img src="image_2.png" /><figcaption>Fig. 1. [C II] luminosity functions at z = 6. The coloured histograms represent predictions from CDM HR (purple) and CDM LB (green). The solid grey line shows the results from the G.A.S. SAM (Lagache et al. 2018). The grey symbols show ALMA data for galaxies at z = 5.5 (dots, Capak et al. 2015), the results from the ALPINE survey (stars, Yan et al. 2020), and lower limits from lensed [C II] line emitters at z = 6 (crosses, Yamaguchi et al. 2017).</figcaption></figure>

luminosities of each individual particle. Specifically, we modelled [C II] as a two-level system considering the quantum number J = 3/2 and J = 1/2 states for the fine-structure transition (2p)[2P$_{3/2}$ - 2P$_{1/2}$]. Following for example Maio et al. (2007) and C24, we can write for each gas particle

$$L_{[CII]} = \Lambda_{[CII]} V = n_{C^+,2} A_{21} \Delta E_{21} V, \quad (1)$$

where V is the volume of a sphere with a radius equal to the smoothing length of the gas particle, and $\Lambda_{[CII]}$ is the emission rate of the [C II] fine-structure transition. The quantity $\Lambda_{[CII]}$ is expressed as a function of the number density of C$^+$ ions excited to the upper state $n_{C^+,2}$, the Einstein coefficient for spontaneous emission $A_{21}$, and the energy level separation $\Delta E_{21}$ (see Appendix A for further details).

In Figure 1 we show the [C II] luminosity function $\Phi$ for galaxies from CDM HR and CDM LB. CDM LB has a lower mass resolution than CDM HR, but because its cosmic volume is larger, we were also able to model more massive haloes. This is clear from the higher [C II] luminosities that were reached in this latter case (~$10^{7.5}$ L$_{\odot}$). Because its resolution is, CDM LB was not able to fully capture the trend of the luminosity function $\Phi$ at the low-luminosity end, where the statistics of faint structures is underestimated by almost one dex in comparison to CDM HR. The faint end of [C II] luminosity functions seems to be sensitive to the model parameters (e.g. metallicity, stellar yields, numerical resolution, and star formation stocasticity, and even to the way [C II] calculations are performed) and to other popular diagnostics, such as UV luminosity functions Maio & Viel (2023), Williams et al. (2024), Kravtsov & Belokurov (2024), Gelli et al. (2024). The two simulations converge well around $10^7$ L$_{\odot}$. The G.A.S. semi-analytical model by Lagache et al. (2018) and CDM LB agree well because they rely on a similar stellar-mass resolution. On the other hand, the CDM HR results are more robust in characterising galaxies with $L_{[CII]} < 10^7$ L$_{\odot}$

than G.A.S. because the resolution of our simulation is high. Thus, CDM HR captures more galaxies with fainter [C II] emission.

A few observational determinations at z = 6 can help us to constrain the [C II] luminosity function at these redshifts (Capak et al. 2015; Yamaguchi et al. 2017; Yan et al. 2020). These data were obtained with ALMA and cover luminosities that are slightly higher than those of COLD SIM galaxies. They have $L_{[CII]} \ge 10^{7.5}$ L$_{\odot}$, and the resulting $\Phi$ features a large scatter of a few dex. In this respect, our results provide predictions for a luminosity range that is still inaccessible with current instruments. We note that Hemmati et al. (2017) reported the first local [C II] emission line luminosity function from 500 galaxies from the Revised Bright Galaxy Sample, where the [C II] luminosities were measured from the Herschel PACS observations of luminous infrared galaxies. The authors reported an evolution in the [C II] luminosity function with cosmic time similar to the evolution trend of the cosmic SFR density.

As discussed in the Introduction, $L_{[CII]}$ is often used as tracer of atomic and molecular gas in high-z galaxies. To investigate how the [C II] emission is related to different gas phases, we introduced the $L_{[CII]}$-to-$M_{H_2}$ conversion factor $\alpha_{[CII]}$ and the $L_{[CII]}$-to-$M_{HI}$ conversion factor $\beta_{[CII]}$. They are defined as $\alpha_{[CII]} = M_{H_2}/L_{[CII]}$ and $\beta_{[CII]} = M_{HI}/L_{[CII]}$, where $M_{HI}$ and $M_{H_2}$ are the atomic and molecular hydrogen masses of each galaxy. In practice, these conversion factors are mass-to-light ratios relating [C II] to neutral gas. They are commonly reported in $M_{\odot}/L_{\odot}$ units and allow us to link the measured $L_{[CII]}$ values to the underlying $H_2$ or H I gas mass.

The results from different numerical simulations in the literature do not converge yet when the molecular gas content in galaxies is estimated. In the case of COLD SIM, different resolutions can lead to variations in $H_2$ masses by a factor of a few in the most extreme scenarios (see Maio & Viel 2023). We aim to link [C II] emission to precise atomic and molecular mass determinations in a cosmological context. We therefore use in Sect. 3 the most accurate high-resolution simulation box, CDM HR, if not otherwise stated.

## 3. Results

In this section, we present our results about the relations between $L_{[CII]}$ and the different gas phases, and we explore their dependence on various galactic properties, such as stellar mass, metallicity, and star formation rate.

### 3.1. Relations of $L_{[CII]}$ with atomic and molecular masses

The [C II] luminosity is often used as tracer of cold gas (with temperatures below $10^4$ K), which is typically composed of atomic (H I) and molecular ($H_2$) hydrogen. In this regime, H I gas is expected to be dominant in mass with respect to $H_2$. However, because of the density and temperature dependences of the [C II] emission signal (see previous section), this does not necessarily mean that most of the $L_{[CII]}$ is produced by neutral atomic gas.

In Figure 2 we plot $L_{[CII]}$ as a function of $M_{H_2}$, $M_{HI}$, and $M_{HII}$ at z = 10 and 6. As expected, the H II phase contributes least to the total mass of hydrogen in these galaxies. Concerning the neutral gas, we observe that galaxies with higher molecular and atomic masses are more luminous in [C II]. This is expected from the cooling of atomic and molecular gas and the consequent collapse into stars, which is consistent with the linear trend found in C24 between the logarithms of $L_{[CII]}$ and SFR. A similar corre-
lation is shown in Figure 2 between $L_{[CII]}$ versus $M_{H_2}$ and $M_{HI}$ at both redshifts. The $L_{[CII]}$ versus $M_{H_2}$ relation extends to higher molecular masses and $L_{[CII]}$ at lower redshifts as a consequence of the galaxy build-up through gas streams and mergers. It shows a $1\sigma$ scatter of up to ~1.5 dex, suggesting that molecular hydrogen is more sensitive to feedback processes, and it is thus more challenging to constrain than atomic gas. An analogous trend is also observed for the $L_{[CII]}$ versus $M_{HI}$ relation.

The COLDSIM galaxies with higher $L_{[CII]}$ can retain more cold gas, atomic and molecular, and fuel the global star formation activity. At all redshifts and for each $L_{[CII]}$ value, $M_{HI}$ is always ~2 dex higher than $M_{H_2}$ and roughly ~3 dex higher than $M_{HII}$. Only a small fraction (locally up to $\approx 10^{-2} - 10^{-1}$) of cold hydrogen gas becomes molecular, while the majority is atomic. Moreover, the scatter exhibited by the $L_{[CII]}$ vs. $M_{HI}$ relation is slightly smaller than that from the $L_{[CII]}$ versus $M_{H_2}$ relation. Below, we provide linear fits, displayed in Figure 2, for the logarithms of $L_{[CII]}$ versus $M_{H_2}$, $M_{HI}$, and $M_{HII}$ gas masses that also include redshift evolution,

$$ \log_{10} (L_{[CII]}/L_{\odot}) = 1.36(\pm 0.01) \log_{10} (M_{H_2}/M_{\odot}) - 0.14(\pm 0.01) z - 4.13(\pm 0.09), \quad (2) $$

$$ \log_{10} (L_{[CII]}/L_{\odot}) = 1.94(\pm 0.01) \log_{10} (M_{HI}/M_{\odot}) - 0.12(\pm 0.01) z - 10.7(\pm 0.1), \quad (3) $$

$$ \log_{10} (L_{[CII]}/L_{\odot}) = 1.45(\pm 0.06) \log_{10} (M_{HII}/M_{\odot}) + 0.06(\pm 0.04) z - 5.58(\pm 0.43). \quad (4) $$

These fits were estimated by including all the snapshots presented in this paper ($z = 10, 8, 7$, and $6$), and the uncertainty of the various parameters corresponds to one standard deviation of the residuals. Equations (2), (3), and (4) can be used to estimate the mass of molecular, neutral, and ionised hydrogen from the detected [CII] emission in high-redshift galaxies. These expressions clearly point out that $L_{[CII]}$ increases with $H_2$, $H$ I and $H$ II gas masses over-linearly. Furthermore, the redshift correction shows that for any fixed mass, earlier objects have suppressed $L_{[CII]}$ by a factor of $\approx 10^{-0.1} z$ for $H_2$, $H$ I and $H$ II. This is a consequence of the typically lower metallicities that are encountered in primordial epochs and of their evolution towards lower $z$.

Our formulae for molecular and atomic hydrogen can be directly compared to literature studies, such as the relations derived from the SIMBA simulations at $z = 6$ (Vizgan et al. 2022a,b). The authors combined the output of simulations of different box sizes to cover a wider mass range, and they estimated $L_{[CII]}$, $M_{H_2}$, and $M_{HI}$ in post-processing. Their fits were computed for galaxies with a molecular hydrogen mass $M_{H_2} > 10^6 M_{\odot}$ and an atomic hydrogen mass $M_{HI} > 10^7 M_{\odot}$. In the figure, we show their results and extrapolate them to lower $H_2$ and $H$ I masses. COLDSIM galaxies are dominated in mass by atomic gas, and the [CII] emission is mostly due to collisions with atomic hydrogen and not molecular hydrogen. In Vizgan et al. (2022a,b), however, molecular gas is the dominant phase, and the majority of [CII] emission is expected to originate by interactions with $H_2$. The authors recognised that in their model, atomic gas does not contribute by construction to $L_{[CII]}$, however. Some of the $H$ I gas might be hidden within the ionised gas, so that the true neutral atomic mass is likely higher than indicated by their results.

In Figure 3 we show the relation between $L_{[CII]}$ and $M_{H_2}$ of $z = 10, 8, 7$, and $6$ combined, and we compare our results with observational determinations. We report a

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** [CII] luminosity $L_{[CII]}$ as a function of $H_2$ mass $M_{H_2}$ (dark pink), $H$ I mass $M_{HI}$ (orange), and $H$ II mass $M_{HII}$ (green) at $z = 10$ (top) and $6$ (bottom). The dashed lines show the median value for CDM HR, and the darker (lighter) shaded regions represent the $1\sigma$ ($2\sigma$) standard deviation. The solid dark pink line represents the fit to the CDM HR results for $M_{H_2}$, the solid orange line shows the fit for $M_{HI}$, and the green line shows $M_{HII}$ (see the text). The solid red and black lines show the linear fit to post-processed SIMBA simulations for $M_{H_2}$ and $M_{HI}$ in their validity range (see Vizgan et al. 2022a,b), and the extrapolations of the fit to lower masses is highlighted with dotted lines. We highlight the redshift evolution of the linear relation between $\log_{10}(L_{[CII]})$ and $\log_{10}(M_{H_2})$, $\log_{10}(M_{HI})$, and $\log_{10}(M_{HII})$. Moreover, we find that $H$ I is always the dominant phase for COLD SIM galaxies.</figcaption></figure>

correlation between $L_{[CII]}$ and $M_{H_2}$ that agrees well with a compilation of observations from the literature (Stacey et al. 1991; Díaz-Santos et al. 2013, 2017; Ferkinhoff et al. 2014; Huynh et al. 2014; Magdis et al. 2014; Capak et al. 2015; Cormier et al. 2015; Gullberg et al. 2015; Schaerer et al. 2015; Accurso et al. 2017; Contursi et al. 2017; Hughes et al. 2017; Zanella et al. 2018; Madden et al. 2020; Kaasinen et al. 2024). These observations included different types of galaxies: local
<figure><img src="image_4.png" /></figure>

**Fig. 3.** [C II] luminosity $L_{[CII]}$ as a function of H$_2$ mass $M_{H_2}$ at z = 10, 8, 7, and 6 combined. The dashed line refers to the mean value for CDM HR (dark pink), and the darker (lighter) shaded region represents the 1σ (2σ) standard deviation. The solid dark pink line represents the linear fit to the CDM HR results at z = 6. The solid black and grey lines refer to the fits to observations by Zanella et al. (2018) and Madden et al. (2020), and the extrapolations of the fits to lower luminosities are highlighted with dotted lines. The grey dots represent observations of local dwarf galaxies (Cormier et al. 2015; Madden et al. 2020), z = 0.5 main-sequence and starburst galaxies (Stacey et al. 1991; Díaz-Santos et al. 2013, 2017; Magdis et al. 2014; Accurso et al. 2017; Contursi et al. 2017; Hughes et al. 2017), and star-forming galaxies in the range z = 2–6 (Ferkinhoff et al. 2014; Huynh et al. 2014; Capak et al. 2015; Gullberg et al. 2015; Schaerer et al. 2015; Zanella et al. 2018; Kaasinen et al. 2024). Globally, we find that our simulations agree well with observations.

dwarves, and main-sequence and star-forming galaxies that covered a wide range of redshift from z = 0 to 7. Although these observations covered molecular gas masses and $L_{[CII]}$ slightly higher than those of the COLDSIM simulations, our results are consistent with the data of local dwarf galaxies from Cormier et al. (2015) and Madden et al. (2020) and with the data of main-sequence galaxies from Accurso et al. (2017). The extrapolation to higher masses of our fit in equation (2) agrees with most observations. In addition, we compared it to the empirical fits by Zanella et al. (2018) and the fit to a model by Madden et al. (2020), which were derived from galaxies with $L_{[CII]} > 10^7 L_\odot$ and $L_{[CII]} > 10^5 L_\odot$, respectively. The fit by Zanella et al. (2018) applies to lensed galaxies at z = 2, although it is routinely used to infer $M_{H_2}$ from $L_{[CII]}$ for higher-redshift galaxies as well. This fit systematically lies above COLDSIM haloes. The fit by Madden et al. (2020) was instead derived from local dwarf galaxies, which have lower metallicities and masses than more typical local galaxies and are usually considered a good analogue for higher-redshift objects (Izotov et al. 2021; Schaerer et al. 2022; Telford et al. 2023). This fit agrees better with our predictions and with the conclusion that the total H$_2$ mass is well traced by $L_{[CII]}$.

Finally, we highlight that our physically motivated fits extend the empirical relations between $L_{[CII]}$ and the mass of cold molecular gas to galaxies with a [C II] luminosity $L_{[CII]} < 10^7 L_\odot$,

<figure><img src="image_5.png" /></figure>

**Fig. 4.** $L_{[CII]}$-to-$M_{H_2}$ conversion factor $\alpha_{[CII]}$ (dark pink) and $L_{[CII]}$-to-$M_{HI}$ conversion factor $\beta_{[CII]}$ (orange) as a function of stellar metallicity $Z_\star$, at z = 10, 8, 7, and 6 combined. The solid dark pink and orange lines represent the mean values obtained from CDM HR for $\alpha_{[CII]}$ and $\beta_{[CII]}$, respectively. The darker (lighter) shaded regions represent the 1σ (2σ) standard deviation. The grey points are $\alpha_{[CII]}$ estimates from observations of local galaxies (Cormier et al. 2015; Accurso et al. 2017; Contursi et al. 2017; Hughes et al. 2017) and z = 2–4 galaxies (Huynh et al. 2014; Schaerer et al. 2015). The black squares are $\beta_{[CII]}$ derived by observations of γ-ray burst afterglows in star-forming galaxies at z = 2–5 (Heintz et al. 2021). The COLD SIM results agree well with lower-redshift observations and provide predictions for future observational works. Both $\alpha_{[CII]}$ and $\beta_{[CII]}$ increase by ~3 dex at low stellar metallicities.

for which observations are still missing. Our fits also provide a reference for high-redshift observations.

## 3.2. Conversion factors and galaxy properties

Based on the linear trend of the relations presented in Sect. 3.1, we constrained the conversion factors $L_{[CII]}$-to-$M_{H_2}$, $\alpha_{[CII]}$, and $L_{[CII]}$-to-$M_{HI}$, $\beta_{[CII]}$. In this section, we explore their dependencies on galactic properties such as $Z_\star$, $M_\star$, and SFR. In Figure 4 we show the relation between $\alpha_{[CII]}$ and $\beta_{[CII]}$ and the galactic stellar metallicity $Z_\star$ at z = 10, 8, 7, and 6 combined. $\alpha_{[CII]}$ and $\beta_{[CII]}$ are roughly constant for $Z_\star < 10^{-2} Z_\odot$. In this regime, the scatter around the two relations is large because these low-mass haloes contain fewer gas particles.

In the range $Z_\star = (10^{-2}-10^{-1})Z_\odot$, both factors experience a drop of about 2 dex on average. The reason for this is that these galaxies have experienced star formation activity (typically, SFR > $10^{-2} M_\odot \text{yr}^{-1}$), and stellar feedback has enriched the gas with metals. The haloes are therefore brighter in $L_{[CII]}$. When an increase in [C II] luminosity is not followed by an equally efficient increase in neutral mass, $\alpha_{[CII]}$ and $\beta_{[CII]}$ may drop to even lower values. At $Z_\star > 10^{-1} Z_\odot$, $\alpha_{[CII]}$ and $\beta_{[CII]}$ are again nearly constant, with a 2σ scatter of ~2 dex for $\alpha_{[CII]}$ and slightly larger than 1 dex for $\beta_{[CII]}$. This regime represents more massive galaxies that have experienced star formation episodes and significant metal enrichment. Since they are in a more evolved stage, they have had time to restore their neutral
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Conversion factor $\alpha_{[CII]}$ as a function of stellar metallicity $Z_*$ colour-coded according to SFR (left panel) and stellar mass (right panel) for COLDSIM galaxies. We highlight that the scatter around the relations presented in Figure 4 is mainly affected by the SFR at $Z_* > 10^{-2.2} Z_\odot$.</figcaption></figure>

atomic and molecular gas content. This is reflected in the shallower $\alpha_{[CII]}$ and $\beta_{[CII]}$ trends.

For all the metallicities, a large $2\sigma$ scatter of ~5–3 dex for $\alpha_{[CII]}$ and $\beta_{[CII]}$ indicates that it is difficult to constrain the mass of H$_2$ and H I from $L_{[CII]}$ precisely, as is also displayed in Figure 2. Our estimates of $\alpha_{[CII]}$ agree well with observations of galaxies with $Z_* > 10^{-1} Z_\odot$ in the range $z = 0–4$ (Huynh et al. 2014; Cormier et al. 2015; Schaerer et al. 2015; Accurso et al. 2017; Contursi et al. 2017; Hughes et al. 2017), especially with the local dwarf galaxies described in Cormier et al. (2015). Moreover, we compared our results of $\beta_{[CII]}$, with observations of star-forming galaxies at $z = 2–5$ by Heintz et al. (2021). The metallicities $Z_* > 10^{-1.5} Z_\odot$ agree well! At lower metallicities, the COLD SIM simulations predict higher values of $\beta_{[CII]}$. The reason might be that Heintz et al. (2021) estimated $M_{HI}$ from $\gamma$-ray bursts, which mainly probe the very dense gas ($N_{HI} > 10^{21} \text{ cm}^{-2}$), while our values also take the diffuse H I retained in the galaxies into account. Early work from Vizgan et al. (2022b) reported a trend with metallicity that is consistent with the findings presented here. Our theoretical results can be useful for observational work in order to constrain the role that the metallicity plays in determining the trend of conversion factors. Especially for low-metallicity galaxies, our findings explore a regime in which observations struggle to characterise molecular and atomic gas masses.

To better understand the trend and scatter of $\alpha_{[CII]}$ as a function of $Z_*$, we show in Figure 5 the same relation, colour-coded according to SFR (left panel) and stellar mass (right panel). Galaxies with $Z_* < 10^{-2} Z_\odot$ typically have an SFR $< 10^{-2.2} M_\odot \text{ yr}^{-1}$ and $M_* < 10^6 M_\odot$. These haloes are cold and poor in metals, and their $L_{[CII]}$ is faint. For higher metallicities, however, the drop in $\alpha_{[CII]}$ coincides with an increase in SFR and $M_*$. As also demonstrated in C24, more massive and star-forming galaxies are also richer in metals and thus in carbon. This result agrees with the evolution of the mass-metallicity relation (C24), but while masses vary by one or two orders of magnitude, $L_{[CII]}$ can vary by four or five dex (see e.g. previous Figure 2).

The wide scatter of ~3 dex at $Z_*$ above $10^{-2} Z_\odot$ for a fixed metallicity causes a larger scatter in SFR. For example, galaxies

with a higher SFR exhibit higher values of $\alpha_{[CII]}$ than quiescent galaxies in the same $Z_*$ bin. This is a result of the stellar feedback, so that galaxies that form more stars are also hotter, with a mass-weighted temperature above ~$10^5$ K. Therefore, carbon in these haloes is present in higher ionisation degrees, and $L_{[CII]}$ decreases. This effect implies that $\alpha_{[CII]}$ values for galaxies with SFR ~1 $M_\odot \text{ yr}^{-1}$ are higher by up to 2 dex than those of galaxies with SFR ~$10^{-1.5} M_\odot \text{ yr}^{-1}$ at a given $Z_*$. No such strong dependence of the scatter on the stellar masses is observed at fixed metallicity. This agrees with the tighter mass-metallicity relation discussed by C24 for example. We additionally confirmed that the behaviour of $\alpha_{[CII]}$ with the specific SFR shows no clear trend. This is consistent with the SFR and $M_*$ results shown in Figure 5.

### 3.3. Gas phases traced by [CII] emission

As mentioned in the previous sections, $L_{[CII]}$ is a powerful tracer of cold gas in general. In this regime, different gas phases may coexist, however. Cold gas typically hosts atomic (H I) and molecular (H$_2$) hydrogen and residual ionised H II gas. As mentioned above, although at temperatures below few $10^4$ K the H I mass is expected to be dominant with respect to the H$_2$ and H II masses, it is unclear which gas phase is associated with most of the [CII] emission because $L_{[CII]}$ depends on the density and temperature.

The contribution to the [CII] luminosity from each gas phase (ionised, molecular, and atomic) is shown in Figure 6 as a function of $M_*$, SFR, $Z_*$, and $M_{HI}$. In order to compute the contribution of each gas phase to the total $L_{[CII]}$ of a galaxy, we decomposed the total [CII] emission (Eq. (1)) into its components (see Appendix A) for each gas particle, and we then summed the $L_{[CII]}$ particle contributions in each galaxy. We specify that the contribution in luminosity arising from H II gas is mainly due to interactions with electrons.

We note that at each redshift, most of the [CII] emission is generated by collisions with H I. However, the contribution of each phase varies with $M_*$ and SFR. At all redshifts, for
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Fractional phase contribution to [C II] luminosity as a function of galactic properties. From top to bottom: Contribution to $L_{[CII]}$ from H I (orange), H$_2$ (dark pink), and H II (green) as a function of stellar mass, star formation rate, stellar metallicity, and H I mass at $z = 10$ (left panels) and 6 (right). The solid line displays the median values, and the shaded areas show the $1\sigma$ scatter. The stars represent the contribution from the three phases in bins with fewer than three galaxies. The [C II] emission mostly traces neutral atomic gas.</figcaption></figure>

stellar masses between $10^5$ and $10^{7.5}$ M$_\odot$ and SFR of $10^{-2}$ to $10^{-0.5}$ M$_\odot$ yr$^{-1}$, the contribution to $L_{[CII]}$ from collisions with H I is approximately 90%, and that of H II can reach up to 10%, while the contribution from H$_2$ is often negligible. In practice, stellar feedback makes these galaxies brighter in $L_{[CII]}$ until the gas particles reach temperatures higher than a few $10^4$ K. Higher temperatures are reached in galaxies with $M_\star > 10^{7.5}$ M$_\odot$ and

SFR >$10^{-0.5}$ M$_\odot$ yr$^{-1}$ (right panels in the figure). As a result, $L_{[CII]}$ in these galaxies does not increase with stellar mass or SFR, and it therefore contributes less H I than less massive and star-forming galaxies. This decrease in the H I contribution is compensated for by an increase in the fraction of $L_{[CII]}$ produced by H$_2$ (up to almost 30%) instead of H II (a few percent). These haloes have high values of $M_{H_2}$, as shown in Figure 2. By contrast, in gas particles with higher temperatures that are dominated by ionised gas, carbon is excited at higher ionisation levels. Haloes with $M_\star < 10^6$ M$_\odot$ and SFR < $10^{-2}$ M$_\odot$ yr$^{-1}$ show a remarkable scatter. As already mentioned in previous sections, these small and quiescent galaxies are more susceptible to feedback effects, and this increases the scatter around the median values. However, it is challenging to estimate $M_\star$ and SFR in these galaxies observationally. Therefore, we show in the two lower panels of Figure 6 the trends for the contributions of the three phases to $L_{[CII]}$ versus $Z_\star$ and $M_{HII}$. The trends are the same as described above: H I is the phase that dominates $L_{[CII]}$, with a fractional contribution >80% for $Z_\star > 10^{-2}$ Z$_\odot$ and $M_{HII} > 10^8$ M$_\odot$. The trend with $M_{HII}$ denotes large spreads at low H I masses, which corresponds to haloes with a lower gas particle number.

Our COLDSIM simulations suggest that $L_{[CII]}$ largely traces atomic H I gas. This result is different from some interpretations of observations in the literature, which concluded that ~60–80% of the [C II] emission arises from H$_2$ gas (Pineda et al. 2013; De Breuck et al. 2019). Our findings are a non-trivial consequence of the [C II] emission properties and the fact that H I is the dominant species in gas below ~$10^4$ K, while H$_2$ is predominant in rarer cold and dense regimes. This result is valid at any redshift and stellar metallicity for galaxies with $M_\star > 10^5$ M$_\odot$ and SFR > $10^{-2.5}$ M$_\odot$ yr$^{-1}$, and it agrees with observational estimates (e.g. Croxall et al. 2017; Cormier et al. 2019). In particular, the recent results of Ikeda et al. (2024) from the CRISTAL survey on 34 spatially resolved galaxies in the range $z = 4–6$ agree well with our results. The authors argued that the extended [C II] emission in half of the sample originated in the diffuse medium, either H I or H II. We also emphasise that although $L_{[CII]}$ largely traces the atomic H I gas phase, this does not contradict our conclusions about the $L_{[CII]}$ versus $M_{H_2}$ relation in Sect. 3.1 and previous observational claims by Madden et al. (2020) or Zanella et al. (2018), for instance. $L_{[CII]}$ can be a robust indicator of the mass of molecular gas, notwithstanding the typically smaller amounts of H$_2$ in this temperature regime. [C II] correlates with cooling and clumping of gas, which in turn enhances H$_2$-driven star formation and increases collisional excitation, and consequently, emissivity.

### 3.4. [C II] relation to dense and star-forming gas

In this section, we discuss the contribution of star-forming gas to the [C II] emission and characterise how different regimes in terms of density and temperature contribute to the total $L_{[CII]}$ in the simulated box. This was done for the sake of simplicity by various works that did not explicitly decompose the gas into its molecular or atomic form, but into star-forming or non-star-forming and dense or diffuse gas.

In the first row of table 1 we display the fraction of that $L_{[CII]}$ originated in star-forming parcels with an SFR > 0 M$_\odot$ yr$^{-1}$ ($f_{[CII],SF}$) at $z = 10, 8, 7$, and 6. Moreover, we divide gas particles into three different phases: diffuse (diff, with $n_H < 10$ cm$^{-3}$), cold dense (CD, with $n_H > 10$ cm$^{-3}$ and $T < 10^4$ K), and warm dense (WD, with $n_H > 10$ cm$^{-3}$ and $10^4$ K < $T < 4 \times 10^4$ K). We chose $4 \times 10^4$ K as the upper limit for the WD phase since
<table><caption>**Table 1.** Fraction of [C II] luminosity.</caption><thead><tr><th></th><th>z = 10</th><th>z = 8</th><th>z = 7</th><th>z = 6</th></tr></thead><tbody><tr><td>f<sub>[CII],SF</sub></td><td>0.72</td><td>0.92</td><td>0.96</td><td>0.97</td></tr><tr><td>f<sub>[CII],diff</sub></td><td>0.36</td><td>0.04</td><td>0.04</td><td>0.04</td></tr><tr><td>f<sub>[CII],CD</sub></td><td>0.12</td><td>0.04</td><td>0.04</td><td>0.02</td></tr><tr><td>f<sub>[CII],WD</sub></td><td>0.52</td><td>0.92</td><td>0.92</td><td>0.94</td></tr></tbody></table>

**Notes.** Produced in gas particles with SFR > 0 M$_{\odot}$ yr$^{-1}$ ($f_{[CII],SF}$), in diffuse gas with $n_H < 10 \text{ cm}^{-3}$ ($f_{[CII],diff}$), in cold dense gas with $n_H > 10 \text{ cm}^{-3}$ and $T < 10^4 \text{ K}$ ($f_{[CII],CD}$), and in warm dense gas with $n_H > 10 \text{ cm}^{-3}$ and $10^4 \text{ K} < T < 4 \times 10^4 \text{ K}$ ($f_{[CII],WD}$) at $z = 10, 8, 7$, and 6.

particles hotter than this do not emit in [C II] (see C24 for more details). Based on these definitions, we then computed the fraction of $L_{[C II]}$ from the different regimes (i.e. $f_{[CII],SF}$, $f_{[CII],diff}$, $f_{[CII],CD}$, and $f_{[CII],WD}$). At each redshift, most of the [C II] emission is generated by star-forming and WD gas, as displayed in Table 1. This results is expected because star-forming particles typically have number densities above $10 \text{ cm}^{-3}$ and are also close to WD particles because stellar feedback quickly increases the temperatures to above $10^4 \text{ K}$. This result indicates that stellar feedback is essential for producing [C II] emission in our simulations. On the one hand, gas is enriched with metals, including carbon, by SN explosions and AGB winds. On the other hand, this gas has the physical conditions in terms of density and temperature to excite [C II] for the emission at 158 µm. The increase in the contribution from WD and star-forming gas from $z = 10$ to $z = 6$ is in line with the increasing SFR density reported in C24. This shows that increasingly more gas in the simulation turns into stars. Moreover, in C24, we found that [C II] can be used to trace star formation in the same redshift range.

This result might be sensitive to star formation and feedback prescriptions, metal enrichment details, and particle resolution. Different-resolution runs resolved gas regions that were dense enough to emit radiation via [C II] 158 µm fine-structure line differently. Schimek et al. (2024), for example, studied the fraction of $L_{[C II]}$ that is emitted by dense particles with $n_H > 10 \text{ cm}^{-3}$ and $T < 10^4 \text{ K}$ and by diffuse gas particles at lower densities. They examined a zoom-in simulation of a single galaxy at redshift $z = 6.5$ and computed the [C II] luminosity in post-processing, finding that most of the $L_{[C II]}$ originated from dense particles, although a high percentage of 44% is emitted by more diffuse gas. Recent results by Gurman et al. (2024) showed that most of the [C II] emission is produced by warm medium, which they defined as gas with $n_H > 10 \text{ cm}^{-3}$ and a temperature of $3 \times 10^3 \text{ K} < T < 3 \times 10^4 \text{ K}$. They performed high-resolution hydrodynamical simulations of a patch of the interstellar medium (with a pixel resolution of about 2 pc) including H$_2$ chemistry and Type II SN feedback, and post-processed their simulated gas cells to infer the carbon abundances and $L_{[C II]}$ signal. They reported that most of the [C II] luminosity originated in warm neutral medium and H II regions. This result is not affected by metallicity, but it might vary with the choice of the temperature upper limit and with a more detailed treatment of low-temperature metal cooling and feedback.

## 4. Discussion

The conversion of molecular gas into stars is traced by the depletion time, which is based on the ratio of the H$_2$ gas mass and SFR. From a theoretical perspective, it is difficult to properly model the physical processes involving molecular gas. Various

scales need to be considered, ranging from large-scale cosmological environments to interstellar physics on a parsec scale. In this context, our work is intended to understand the origin of the measured $L_{[C II]}$ signals and to support observers, especially in interpreting the properties of galaxies detected with facilities such as ALMA and JWST.

In the following, we first discuss the assumptions and caveats of our COLD SIM numerical simulations. We then compare our results with other numerical works in the literature, and we finally show their usefulness for observational studies in accurately estimating the molecular gas mass in high-$z$ galaxies.

### 4.1. COLD SIM assumptions and caveats

Our simulations, like all others in the literature, relied on a number of assumptions that affect the quantitative results shown in this paper. We stress that the adoption of a detailed non-equilibrium model for chemical abundance calculations, similar to the one used here, is crucial for capturing the processes that shape atomic and molecular gas, as highlighted by Davé et al. (2020) or Hu et al. (2021), for instance.

Possible caveats caused by the choices in the set-up were extensively described by Maio et al. (2022), Maio & Viel (2023), and C24, for example. There, we already discussed the main assumptions about the initial stellar metallicity, different models of explosive nucleosynthesis, stellar wind variations in the values of the carbon yields by up to 1 dex, and the atomic physics uncertainties in the collisional rates needed to calculate the [C II] luminosity.

Sources of uncertainty may arise from the way in which stellar feedback is implemented. In COLD SIM, each stellar particle is considered as a stellar population with a Salpeter (Salpeter 1955) initial mass function. Other IMF models, such as Chabrier (Chabrier 2003), which are more bottom heavy, would inject more $\alpha$ elements from supernova Type II explosions into the surrounding medium. This can be seen, for example, by considering the supernova Type II fraction of the two IMFs, which is ~0.2% for a Salpeter and ~3% for a Chabrier IMF. In our study, the choice of the IMF might have an impact on $L_{[C II]}$ because the mass of ejected carbon depends on the slope of the IMF. In contrast, we do not expect the choice of a Chabrier over a Salpeter IMF to have strong consequences on the values of $M_{HII}$ and $M_{H\alpha}$. Maio et al. (2022) analysed both IMFs and did not note any significant differences in the primordial H$_2$ and H I cosmological mass density.

Our stellar feedback model assumes that Pop III stars are distributed with a Salpeter IMF, similarly to Pop II and Pop I stars. However, different choices for the early Pop III IMF are possible. For example, for a top-heavy Pop III IMF, primordial stars would have much shorter lifetimes than Pop II and Pop I stars. Consequently, at high redshifts, they would more quickly enrich the interstellar medium with metals (see results in e.g. Maio et al. 2010), with possible implications for early [C II] luminosities.

In addition to Type II SNe, we also considered AGB stars and Type Ia SNe, which are driven by lower-mass stars with longer lifetimes. Their typical timescales are about $10^8-10^9$ yr, and thus, they can enrich the cosmic medium already by $z \approx 6$. This has direct implications for feedback effects and chemical abundances during the epoch of reionisation (as demonstrated by the first spectroscopic evidence of a mature stellar population in a $z \sim 8$ galaxy, A2744-YD4; Witten et al. 2024 and massive quenched galaxy at $z \sim 7$, RUBIES-UDS-QG-z7; Weibel et al. 2024). As an example, AGB stars produce relevant amounts of carbon in only few hundred million years, and therefore, the
chemical gas patterns and the expected [CII] emission may be severely affected by them. The metal yields of Type Ia SNe do not strongly depend on the initial stellar metallicity, but their occurrence requires timescales shorter than $10^8$ to more than $10^9$ yr and assumptions about the realisation probability of the Type Ia SN scenario. This latter is basically unknown for early epochs. We therefore relied on standard prescriptions suggested by stellar evolution models (for usual stellar parameters the delay time distribution peaks well below one Gyr; Tinsley 1980; Greggio 2005; Maoz & Graur 2017). In this way, we were able to take their possible impact on our [CII] luminosity estimates into account as well. We stress that these concerns about AGB enrichment and Type Ia SNe must be kept in mind when focusing on primordial metal signatures because most high-redshift numerical simulations tend to neglect them (although roughly half Type Ia SNe in a single stellar population explode within less than a billion years; Cappellaro et al. 2015).

In our work, galactic winds take place at a fixed velocity of 350 km/s and are temporarily decoupled from cooling and star formation (according to Maio et al. 2011a). This allowed us to obtain a more efficient mass ejection and to mitigate potential overcooling problems, and at the same time, to prevent unphysical behaviour related to the formation of stars in the wind phase. The uncertainties about wind velocities (or mass-loading factors) are large, but we verified that at these early epochs, the exact adopted value has no drastic effects for the H I and H$_2$ phases because the mass of the hosting haloes is relatively low (Maio et al. 2022). At intermediate and low redshifts, the impacts might be more relevant, however Nagamine et al. (2004), while a decrease in outflow velocity might lead to more prominent galactic blue peaks at $z > 6$ (Hassan & Gronke 2021).

Throughout this paper, we did not consider feedback from primordial black holes (BHs). Their presence might cause a displacement in the gas that would create turbulence and heat the surrounding medium. In addition to the huge uncertainties in early-BH models, it is possible that this effect reduces the amount of cold gas, particularly in its molecular form, which would cause the baryons to be in a hotter phase. However, at these redshifts, the masses of BHs are expected to be lower than those in the local Universe, and if the BH mass-bulge mass relation still holds, their masses should be ~$10^{-3}$ times the hosting stellar masses at most (see detailed figures in C24). Severe implications from early-BH seeds would only be found in very extreme cases, such as those discussed in Woods et al. (2024), Patrick et al. (2023), and Maio et al. (2019). Thus, the inclusion of primordial BHs is not expected to have drastic consequences for our findings in general.

In COLDSIM heavy elements produced by stars are spread over the neighbouring particles, and metal mixing is mimicked through SPH kernel smoothing. This means that in a low-density region at larger distances from the point of origin of the metals, the metallicity might be underestimated. This is typical of many SPH simulations, however, and winds are usually helpful to alleviate this problem. Since we focused on cold-gas phases and particles close to star-forming regions, this effect is not expected to create significant consequences in our results.

Variations in the adopted cosmological parameters, higher-order corrections to linear theory, or in the assumed dark-matter scenario are likely to play a minor role because [CII] emission essentially stems from cold gas that is bound within the structure potential wells and is not affected long-distance cosmological phenomena. They might alter the high-redshift statistics, however (for a deeper discussion on the topic see e.g. Yoshida et al.

2003; Maio et al. 2006, 2011b; Maio & Viel 2015, 2023, and references therein).

## 4.2. Comparison with other numerical studies

Numerically, it is one of the greatest challenges to accurately model physical processes that cause the [CII] fine-structure line emission at 158 µm in a cosmological context. COLDSIM has a non-equilibrium chemical network implemented on the fly that is designed to describe the primordial gas chemistry. Specifically, at each time step, the abundances of each atomic and molecular gas species are computed in a consistent way with cooling and heating processes. Moreover, individual elements, including carbon, are tracked separately and are not inferred from a global metallicity. This means that the amount of carbon in each particle is estimated at each simulation time step according to the underlying stellar evolution model. This dedicated set-up enables the calculation of the contribution of each gas phase, ionised, atomic, and molecular, to the [CII] emission.

As the results discussed in this paper are sensitive to the implementation adopted to track the different gas phases, we compared our findings with those from other numerical work in the literature. We emphasise that the various models differ in terms of initial set-up, box size, resolution, feedback implementation, sub-grid physics, and post-processing of $L_{[CII]}$ and gas-phase contributions.

Previous works investigated the contribution of different gas phases to the [CII] emission in galaxies at high $z$ via semi-analytical models (Vallini et al. 2015; Lagache et al. 2018; Popping et al. 2019; Yang et al. 2022), zoom-in simulations of individual galaxies (Olsen et al. 2017, 2018; Pallottini et al. 2019; Lupi & Bovino 2020; Bisbas et al. 2022; Schimek et al. 2024; Khatri et al. 2024a), or a section of the interstellar medium (Gurman et al. 2024), and post-processing of cosmological simulations (Ramos Padilla et al. 2021, 2023; Vizgan et al. 2022a,b; Garcia et al. 2024). All these investigations were based on a number of different assumptions and sub-grid models. This can lead to large discrepancies in the results. For example, with the exception of Olsen et al. (2017, 2018), Lupi & Bovino (2020), and Bisbas et al. (2022), all other works traced a global metallicity and inferred the carbon mass in post-processing. COLDSIM instead self-consistently follows individual elements according to stellar evolution and the associated gas cooling. Additionally, as it is computationally expensive to include H$_2$ processes in the chemical network as here, only a few works included them, and for single haloes, they either implemented [CII] emission as well (Bisbas et al. 2022) or neglected it (Pallottini et al. 2017, 2019).

Only a handful of studies investigated the contribution of different gas phases to $L_{[CII]}$ at the epoch of reionisation through cosmological simulations. All these works estimated the H I and H$_2$ masses in post-processing, however. Here, we mention Vizgan et al. (2022a,b), Garcia et al. (2024) with SIMBA, and Ramos Padilla et al. (2021, 2023) with EAGLE. They covered box sizes ranging from 25 cMpc/h to 100 cMpc/h.

As described in Sect. 3.1, our results robustly indicate that the logarithms of $L_{[CII]}$ and those of $M_{H_2}$ and $M_{HI}$ are linearly correlated. Previously, this relation was only studied via [CII] post-processing modelling by Vizgan et al. (2022a,b). The slope and amplitude of their relations are different from ours, but the range of masses and luminosities from which our fits and the SIMBA fits were derived are different.

Another important result is the evolution with redshift of the different galaxy properties, including $L_{[CII]}$, molecular and atomic gas mass, SFR, metallicity, and stellar mass. In this
respect, Liang et al. (2024) studied the empirical $L_{[CII]}$-SFR relation with the FIRE cosmological zoom-in simulations. They focused on the so-called [CII]-deficit galaxies (Knudsen et al. 2016). The authors reported two distinct physical regimes: H$_2$-rich galaxies, in which the depletion time is the main driver of the [CII] deficit; and H$_2$-poor galaxies, in which gas metallicity dominates. Consistently with C24, their findings suggested that the [CII] deficit is a widespread characteristic of galaxies, highlighting the need for caution for using a constant $L_{[CII]}$-to-SFR conversion factor (derived from local star-forming galaxies) to estimate SFRs at high redshift (see, e.g., Righi et al. 2008; Sommovigo et al. 2021; Roy & Lapi 2024; Khatri et al. 2024b).

In our study, we investigated the dependence of the conversion factors $\alpha_{[CII]}$ and $\beta_{[CII]}$ on galactic properties. We found that the two conversion factors are dependent on the stellar metallicity, $Z_\star$. Additionally, for $Z_\star > 10^{-2} Z_\odot$, galaxies with larger $\alpha_{[CII]}$ are those with a higher SFR, which shows that the conversion factor depends strongly on the star formation activity in the galaxy.

In Sect. 3.3 we showed that HI gas phase dominates [CII] emission, while the H$_2$ gas phase can contribute no more than ~30% to a galaxy $L_{[CII]}$. Thus, some of the previous findings (such as the ones by Vallini et al. 2015; Lupi & Bovino 2020; Ramos Padilla et al. 2021, 2023; Bisbas et al. 2022; Gurman et al. 2024) agree with ours.

For the sake of simplicity, a number of works did not explicitly decompose the gas into molecular or atomic gas, but into star-forming or non-star-forming and dense or diffuse gas (Schimek et al. 2024; Gurman et al. 2024). Di Cesare et al. (2024) recently speculated that in galaxy mergers at $z > 4$, about 50% of the [CII] emission might come from the diffuse medium that originated in the middle of the merger. The authors proposed that even at lower redshifts, [CII] emission from diffuse gas is not negligible. Muñoz-Elgueta et al. (2024) studied the [CII] emission in a sample of 357 galaxies at redshift $z = 5$. Despite the limited simulation resolution and possible shortcomings in the adopted modelling, they suggested that the majority of $L_{[CII]}$ originated from dense star-forming sites (assuming a star formation density of $0.1 \text{ cm}^{-3}$ and a fixed cold-phase gas temperature of 1000 K), while the [CII] extended emission mainly stems from satellite galaxies within the virial radius. We note again that different simulations and post-processing methods provide non-converging results.

In addition to these numerical efforts that attempted to characterise $M_{H_2}$ with $L_{[CII]}$, others used the CO emission. Among them, Hu et al. (2023) presented high-resolution (~0.2 pc) hydrodynamical simulations of an isolated low-metallicity ($0.1 Z_\odot$) dwarf galaxy coupled with a time-dependent chemistry network. This metallicity regime is of particular relevance to the high-redshift objects that are now unrevealed by JWST Curtis-Lake et al. (2023). On a cosmological scale, Keating et al. (2020) used the FIRE zoom-in simulations to reproduce the $L_{CO}$-to-$M_{H_2}$ conversion factor. They relied on an equilibrium solver for estimating H$_2$ and CO abundances in the post-processing and reported little evolution of the $L_{CO}$-to-$M_{H_2}$ conversion factor with SFR. This contradicts our non-equilibrium findings about the $L_{[CII]}$-to-$M_{H_2}$ conversion factor.

### 4.3. Predictions for high-$z$ observations

Observationally, many works reported a linear relation between $L_{[CII]}$ and the mass of H$_2$ in the local Universe and then estimated a conversion factor. Nevertheless, it is complicated to estimate it empirically at high redshift or at low metallicities. In these

regimes, it becomes challenging to observe metal emission lines or to infer the mass of molecular gas without relying on $L_{[CII]}$. The maximum redshift reached to estimate $\alpha_{[CII]}$ from galactic gas is $z \approx 4$ (Huynh et al. 2014), but the hosting objects are usually very luminous ($L_{[CII]} > 10^9 L_\odot$), massive in H$_2$, and are not statistically representative of the majority of high-redshift galaxies, which are less luminous and have lower masses. For galaxies at the epoch of reionisation, when CO cannot be observed, the conversion factor of $\alpha_{[CII]} = 31 M_\odot/L_\odot$ empirically derived by Zanella et al. (2018) from star-forming lensed galaxies at $z = 2$ is often applied. Dessauges-Zavadsky et al. (2020), Aravena et al. (2024), and Kaasinen et al. (2024) used this value to estimate the H$_2$ mass of galaxies in the redshift range $z = 4-8$. Kaasinen et al. (2024) reported that estimates of $M_{H_2}$ can vary by more than one dex depending on whether [CII] or CO emissions are used as tracers. Recently, Salvestrini et al. (2024), demonstrated that conversion factor calibrated in the local Universe overestimate significantly the $M_{H_2}$ measurement at the epoch of reionisation. Moreover, Glazer et al. (2024) reported that low-mass galaxies at $z \approx 7$ present a deficit in [CII] luminosity when their $L_{[CII]}$ is compared to the expectations by the $L_{[CII]}$-SFR relation calibrated on local dwarf galaxies (e.g. De Looze et al. 2014). The authors justified this conclusion by invoking the low metallicities of their sample compared to those of the local Universe. This results therefore indicates that galaxies at high $z$ might have different properties from those at lower redshift and would require a physically motivated $\alpha_{[CII]}$ calibration.

Our findings highlight that $L_{[CII]}$ is dominated by the HI atomic phase of the gas in high-redshift galaxies. Nonetheless, it is important to note that this result does not imply that all available HI interacts with carbon and produces $L_{[CII]}$. The HI gas is clearly expected to be more extended than the [CII] region, which requires ionising photons with energy greater than 11.3 eV, as also observed in nearby galaxies (Péroux et al. 2019; Szakacs et al. 2021). Moreover, pristine gas, such as the gas that is far from star-forming regions, could be rich in HI, but contain no metals (hence would not emit any $L_{[CII]}$ signal). This is particularly true in primordial galaxies, when the amount of pristine material is expected to be greater than in the local Universe. We therefore need a $L_{[CII]}$-to-$M_{HI}$ conversion factor $\beta_{[CII]}$ computed at $z > 6$ that evolves with metallicity.

Our work is useful for estimating the mass of cold gas, atomic and molecular, in galaxies at the epoch of reionisation, such as those of the ALPINE and REBELS surveys. Furthermore, it would be interesting to understand how [CII] emission is spatially distributed with respect to HI and H$_2$. This type of analysis will be very useful in view of the new results from the CRISTAL survey, for example. The CRISTAL project targets 34 galaxies in the range $z = 4-6$ and spatially resolves their [CII] emission. Our simulations so far agree with recent results by Ikeda et al. (2024) and Li et al. (2024), who reported extended [CII] emission that was primarily linked to star formation activity, even at kiloparsec scales, with contributions from atomic gas and mergers that might further expand the [CII] line distribution.

Finally, in order to determine the observability of the simulated objects, we computed the exposure time needed for detecting COLDSIM galaxies with the ALMA 12 m array. For this purpose, we used the CASA sensitivity calculator (McMullin et al. 2007) and determined the time required to detect the most luminous galaxies in our sample at $5\sigma$. We selected [CII] emission at $z = 6$, corresponding to an observed frequency of 271.5 GHz, which is observable in ALMA Band 7. To detect a galaxy with $L_{[CII]} \sim 10^7 L_\odot$, the exposure time needed is seven hours on-source, or shorter for lower redshifts. Our predictions indicate
that although it is still challenging to probe these galaxies, with the physical properties predicted by COLDSIM, more massive haloes are within reach for ALMA.

## 5. Conclusions

It is not yet possible to observationally resolve the spatial distribution of cold gas and to estimate the masses of H I and H$_2$ components in galaxies at the epoch of reionisation. An important tool for tracing atomic and molecular gas at these redshifts is the fine-structure 158 µm [C II] emission line observed by dedicated facilities such as ALMA, and it is within the reach of the upcoming AtLAST project. However, it is still unclear which phase of the gas at these epochs mostly causes the [C II] luminosity, whether ionised, neutral atomic, or molecular, and what the high-z trends of the corresponding mass-to-light ratios are.

In this work, we have extended our previous numerical analysis Casavecchia et al. (2024) and explored the processes that cause [C II] emission by using selected runs from the COLDSIM numerical simulations (Maio et al. 2022; Maio & Viel 2023; Casavecchia et al. 2024). These cosmological simulations include in addition to the usual implementation of gravity and hydrodynamic calculations a physics-rich and accurate non-equilibrium chemical network that follows the evolution of atomic and molecular species, in particular, those we focused on here: H I, H II and H$_2$. Primordial chemistry was computed consistently with star formation and feedback effects and with cooling and heating processes. With COLDSIM, we tracked the time evolution of individual heavy elements, including carbon, and consequently derived the [C II] emission from collision processes with H I, H II, and H$_2$. In this way, we were able to provide physically motivated conversion factors that link the [C II] luminosity, $L_{[CII]}$, to the H$_2$ mass ($\alpha_{[CII]}$) and H I mass ($\beta_{[CII]}$), and their scaling with stellar mass and SFR.

Our main findings are summarised below.

1. The COLDSIM results reproduce the [C II] luminosity function trends observed by ALMA. Our galaxies extend the [C II] luminosity function to $L_{[CII]}$ lower than in accessible to observations at $z = 6$ and provide first predictions in a range in which data are still limited.

2. We found a linear relation between $\log_{10}(L_{[CII]})$ and $\log_{10}(M_{H_2})$ at all redshifts. Our results at $z = 6$ agree with observational estimates from galaxies in the range $z = 0-7$. In particular, the model by Madden et al. (2020), derived from local dwarf galaxies, broadly agrees with our estimates. The range of $M_{H_2}$ and $L_{[CII]}$ that is analysed and/or observed is not always the same, and thus, the comparison must refer to an extrapolation of the results.

3. We also studied the relation between $\log_{10}(L_{[CII]})$ and $\log_{10}(M_{HI})$. We found that for a fixed $L_{[CII]}$ the mass of H I is typically ~2 dex higher than that of H$_2$.

4. The amplitude and scatter of the mentioned relations evolve with redshift. We provided redshift-dependent fitting formulae that link $L_{[CII]}$ to $M_{H_2}$ and $M_{HI}$, and that can be used to estimate the mass of molecular and atomic gas of a galaxy from its (detected) [C II] luminosity.

5. We investigated the dependence of the $\alpha_{[CII]}$ and $\beta_{[CII]}$ conversion factors on the stellar metallicity, $Z_\star$, finding a good agreement with observational estimates at lower redshift. For $Z_\star < 10^{-1} M_\odot$, the lack of carbon causes an increase of about 2 dex in the values of $\alpha_{[CII]}$ and $\beta_{[CII]}$, and this might help explain the common [C II] deficit in galaxies around $z \sim 7$ (consistently with observational findings by Glazer et al. 2024).

6. We explored the scatter in the relation between $\alpha_{[CII]}$ and $Z_\star$, finding that for metallicities $Z_\star > 10^{-2} M_\odot$, the large scatter depends on the galaxy SFR. Highly star-forming galaxies ionise carbon at higher ionisation states, making them less luminous in $L_{[CII]}$ than galaxies with similar metallicity, but lower SFR.

7. We examined the contribution of each gas phase, ionised, neutral atomic, and molecular, to the total emission in $L_{[CII]}$, and we studied the correlation with SFR and $M_\star$. We found that for an SFR > $10^{-2.5} M_\odot \text{yr}^{-1}$ and $M_\star > 10^5 M_\odot$, the phase that dominates $L_{[CII]}$ is the atomic phase, with fractional contributions ranging from 70 to 90%.

8. The contribution from molecular gas varies, with values up to ~30% for galaxies with SFR $\ge 10^{-1} M_\odot \text{yr}^{-1}$ and $M_\star \ge 10^8 M_\star$. Even though $L_{[CII]}$ mainly traces H I gas, it is still extremely useful to estimate the mass of the sub-dominant molecular phase because the fluctuations of H$_2$ the contributions across the considered mass range are relatively minor.

9. We investigated the relation between the three gas-phase contributions to the total $L_{[CII]}$ with $Z_\star$ and $M_{HI}$. In this case, H I is also the dominating phase, with a fractional contribution >80% for $Z_\star > 10^{-2} M_\odot$ and $M_{HI} > 10^8 M_\odot$.

10. In the entire simulated volume, the vast majority of $L_{[CII]}$ is predicted to originate in warm dense and star-forming regions, regardless of the redshift. This agrees with the latest observational results by Ikeda et al. (2024).

The work presented here offers new theoretical insights to study the processes at the origin of [C II] 158 µm fine-structure line emission during the epoch of reionisation and supports recent and upcoming observations of star-forming gas in primordial times. However, to better link theory and observations future investigations are still required. They will be crucial to reach a full understanding of carbon chemical signatures in the first Gyr and interpret correctly the origin of the first galaxies at such primeval epochs.

**Acknowledgements.** We acknowledge the anonymous referee for the constructive report, that allowed us to greatly improve the original manuscript. We are thankful to Aniket Bhagwat for fruitful discussions on emission processes in high-z galaxies. We are grateful to Rodrigo Herrera-Camus for providing a preview of the results on [C II] luminosity in the CRISTAL galaxies. BC thanks Victoria Bollo for the help in the CASA calculations and the interesting discussion about molecular gas emission and absorption features through cosmic time. UM acknowledges financial support from the theory grant no. 1.05.23.06.13 "FIRST - First Galaxies in the Cosmic Dawn and the Epoch of Reionization with High Resolution Numerical Simulations" and the travel grant no. 1.05.23.04.01 awarded by the Italian National Institute of Astrophysics. UM was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP), Santa Barbara (USA). This research was supported by the International Space Science Institute (ISSI) in Bern (Switzerland), through ISSI International Team project no. 564 (The Cosmic Baryon Cycle from Space). The numerical calculations done throughout this work have been performed on the machines of the Max Planck Computing and Data Facility of the Max Planck Society, Germany. We acknowledge the NASA Astrophysics Data System for making available their bibliographic research tool

## References

Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New A, 2, 181

Accurso, G., Saintonge, A., Catinella, B., et al. 2017, MNRAS, 470, 4750

Andreani, P., Miyamoto, Y., Kaneko, H., et al. 2020, A&A, 643, L11

Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 71

Aravena, M., Heintz, K., Dessauges-Zavadsky, M., et al. 2024, A&A, 682, A24

Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

Bisbas, T. G., Walch, S., Naab, T., et al. 2022, ApJ, 934, 115

Blitz, L., & Rosolowsky, E. 2006, ApJ, 650, 933

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
Boogaard, L. A., Decarli, R., Walter, F., et al. 2023, ApJ, 945, 111

Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

Cappellaro, E., Botticella, M. T., Pignata, G., et al. 2015, A&A, 584, A62

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Casavecchia, B., Maio, U., Péroux, C., & Ciardi, B. 2024, A&A, 689, A106

Chabrier, G. 2003, PASP, 115, 763

CONCERTO Collaboration (Ade, P., et al.) 2020, A&A, 642, A60

Contursi, A., Baker, A. J., Berta, S., et al. 2017, A&A, 606, A86

Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53

Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23

Croxall, K. V., Smith, J. D., Pellegrini, E., et al. 2017, ApJ, 845, 96

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nat. Astron., 7, 622

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146

De Breuck, C., Weiß, A., Béthermin, M., et al. 2019, A&A, 631, A167

De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62

Decarli, R., Walter, F., González-López, J., et al. 2019, ApJ, 882, 138

Decarli, R., Aravena, M., Boogaard, L., et al. 2020, ApJ, 902, 110

Dessauges-Zavadsky, M., Ginolfi, M., Pozzi, F., et al. 2020, A&A, 643, A5

Di Cesare, C., Ginolfi, M., Graziani, L., et al. 2024, A&A, 690, A255

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32

Dobbs, C. L., & Pringle, J. E. 2013, MNRAS, 432, 653

Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Dunne, L., Maddox, S. J., Papadopoulos, P. P., Ivison, R. J., & Gomez, H. L. 2022, MNRAS, 517, 962

Ebzagezio, S., Seifried, D., Walch, S., et al. 2023, MNRAS, 525, 5631

Ferkinhoff, C., Brisbin, D., Parshley, S., et al. 2014, ApJ, 780, 142

Feruglio, C., Maio, U., Tripodi, R., et al. 2023, ApJ, 954, L10

Fukui, Y., & Kawamura, A. 2010, ARA&A, 48, 547

Galli, D., & Palla, F. 1998, A&A, 335, 403

Garcia, K., Narayanan, D., Popping, G., et al. 2024, ApJ, 974, 197

Gelli, V., Mason, C., & Hayward, C. C. 2024, ApJ, 975, 192

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Gkogkou, A., Béthermin, M., Lagache, G., et al. 2023, A&A, 670, A16

Glazer, K., Bradăc, M., Sanders, R. L., et al. 2024, MNRAS, 531, 945

Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 421, 9

Gnedin, N. Y., & Draine, B. T. 2014, ApJ, 795, 37

Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13

Greggio, L. 2005, A&A, 441, 1055

Gullberg, B., De Breuck, C., Vieira, J. D., et al. 2015, MNRAS, 449, 2883

Gurman, A., Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2024, ApJ, 965, 179

Haardt, F., & Madau, P. 1996, ApJ, 461, 20

Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421

Hamanowicz, A., Zwaan, M. A., Péroux, C., et al. 2023, MNRAS, 519, 34

Hassan, S., & Gronke, M. 2021, ApJ, 908, 219

Heintz, K. E., Watson, D., Oesch, P. A., Narayanan, D., & Madden, S. C. 2021, ApJ, 922, 147

Hemmati, S., Yan, L., Diaz-Santos, T., et al. 2017, ApJ, 834, 36

Hernandez, S., Jones, L., Smith, L. J., et al. 2023, ApJ, 948, 124

Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555

Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306

Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2021, ApJ, 920, 44

Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2023, ApJ, 952, 140

Hughes, T. M., Ibar, E., Villanueva, V., et al. 2017, A&A, 602, A49

Huynh, M. T., Kimball, A. E., Norris, R. P., et al. 2014, MNRAS, 443, L54

Ikeda, R., Tadaki, K. I., Mitsuhashi, I., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202451811

Izotov, Y. I., Guseva, N. G., Fricke, K. J., et al. 2021, A&A, 646, A138

Kaasinen, M., Venemans, B., Harrington, K. C., et al. 2024, A&A, 684, A33

Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902

Keating, L. C., Richings, A. J., Murray, N., et al. 2020, MNRAS, 499, 837

Kennicutt, R. C., Jr 1998, ApJ, 498, 541

Khatri, P., Porciani, C., Romano-Díaz, E., Seifried, D., & Schäbe, A. 2024a, A&A, 688, A194

Khatri, P., Romano-Díaz, E., & Porciani, C. 2024b, A&A, submitted [arXiv:2411.09755]

Klaassen, P. D., Mroczkowski, T. K., Cicone, C., et al. 2020, SPIE Conf. Ser., 11445, 114452F

Klitsch, A., Péroux, C., Zwaan, M. A., et al. 2019, MNRAS, 490, 1220

Knudsen, K. K., Richard, J., Kneib, J.-P., et al. 2016, MNRAS, 462, L6

Kravtsov, A., & Belokurov, V. 2024, arXiv e-prints [arXiv:2405.04578]

Krumholz, M. R. 2013, MNRAS, 436, 2747

Krumholz, M. R. 2014, Phys. Rep., 539, 49

Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36

Lagache, G., Cousin, M., & Chatzikos, M. 2018, A&A, 609, A130

Lagos, C. d. P., Crain, R. A., Schaye, J., et al. 2015, MNRAS, 452, 3815

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

Lepp, S., & Shull, J. M. 1983, ApJ, 270, 578

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782

Li, J., Da Cunha, E., González-López, J., et al. 2024, ApJ, 976, 70

Liang, L., Feldmann, R., Murray, N., et al. 2024, MNRAS, 528, 499

Lupi, A., & Bovino, S. 2020, MNRAS, 492, 2818

Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141

Magdis, G. E., Rigopoulou, D., Hopwood, R., et al. 2014, ApJ, 796, 63

Maio, U., & Viel, M. 2015, MNRAS, 446, 2760

Maio, U., & Viel, M. 2023, A&A, 672, A71

Maio, U., Dolag, K., Meneghetti, M., et al. 2006, MNRAS, 373, 869

Maio, U., Dolag, K., Ciardi, B., & Tornatore, L. 2007, MNRAS, 379, 963

Maio, U., Ciardi, B., Yoshida, N., Dolag, K., & Tornatore, L. 2009, A&A, 503, 25

Maio, U., Ciardi, B., Dolag, K., Tornatore, L., & Khochfar, S. 2010, MNRAS, 407, 1003

Maio, U., Khochfar, S., Johnson, J. L., & Ciardi, B. 2011a, MNRAS, 414, 1145

Maio, U., Koopmans, L. V. E., & Ciardi, B. 2011b, MNRAS, 412, L40

Maio, U., Petkova, M., De Lucia, G., & Borgani, S. 2016, MNRAS, 460, 3733

Maio, U., Borgani, S., Ciardi, B., & Petkova, M. 2019, PASA, 36, e020

Maio, U., Péroux, C., & Ciardi, B. 2022, A&A, 657, A47

Maoz, D., & Graur, O. 2017, ApJ, 848, 25

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Mitsuhashi, I., Tadaki, K.-I., Ikeda, R., et al. 2024, A&A, 690, A197

Muñoz-Elgueta, N., Arrigoni Battaia, F., Kauffmann, G., et al. 2024, A&A, 690, A392

Nagamine, K., Springel, V., & Hernquist, L. 2004, MNRAS, 348, 421

Olsen, K., Greve, T. R., Narayanan, D., et al. 2017, ApJ, 846, 105

Olsen, K., Greve, T. R., Narayanan, D., et al. 2018, ApJ, 857, 148

Pallottini, A., Ferrara, A., Bovino, S., et al. 2017, MNRAS, 471, 4128

Pallottini, A., Ferrara, A., Decataldo, D., et al. 2019, MNRAS, 487, 1689

Patrick, S. J., Whalen, D. J., Latif, M. A., & Elford, J. S. 2023, MNRAS, 522, 3795

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Péroux, C., Zwaan, M. A., Klitsch, A., et al. 2019, MNRAS, 485, 1595

Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103

Popping, G., & Péroux, C. 2022, MNRAS, 513, 1531

Popping, G., Narayanan, D., Somerville, R. S., Faisst, A. L., & Krumholz, M. R. 2019, MNRAS, 482, 4906

Ramos Padilla, A. F., Wang, L., Ploeckinger, S., van der Tak, F. F. S., & Trager, S. C. 2021, A&A, 645, A133

Ramos Padilla, A. F., Wang, L., van der Tak, F. F. S., & Trager, S. C. 2023, A&A, 679, A131

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7

Righi, M., Hernández-Monteagudo, C., & Sunyaev, R. A. 2008, A&A, 489, 489

Roussel, H., Helou, G., Hollenbach, D. J., et al. 2007, ApJ, 669, 959

Roy, A., & Lapi, A. 2024, arXiv e-prints [arXiv:2407.19007]

Saintonge, A., Kauffmann, G., Kramer, C., et al. 2011, MNRAS, 415, 32

Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22

Salpeter, E. E. 1955, ApJ, 121, 161

Salvestrini, F., Feruglio, C., Tripodi, R., et al. 2024, A&A, submitted [arXiv:2412.02688]

Saslaw, W. C., & Zipoy, D. 1967, Nature, 216, 976

Schaerer, D., Boone, F., Jones, T., et al. 2015, A&A, 576, L2

Schaerer, D., Ginolfi, M., Béthermin, M., et al. 2020, A&A, 643, A3

Schaerer, D., Marques-Chaves, R., Barrufet, L., et al. 2022, A&A, 665, L4

Schimek, A., Decataldo, D., Shen, S., et al. 2024, A&A, 682, A98

Sommovigo, L., Ferrara, A., Carniani, S., et al. 2021, MNRAS, 510, 4878

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423

Szakacs, R., Péroux, C., Zwaan, M., et al. 2021, MNRAS, 505, 4746

Szakacs, R., Péroux, C., Zwaan, M. A., et al. 2022, MNRAS, 512, 4736

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Tarantino, E., Bolatto, A. D., Herrera-Camus, R., et al. 2021, ApJ, 915, 92

Telford, O. G., McQuinn, K. B. W., Chisholm, J., & Berg, D. A. 2023, ApJ, 943, 65

Telikova, K., González-López, J., Aravena, M., et al. 2024, arXiv e-prints [arXiv:2411.09033]

Tinsley, B. M. 1980, Fund. Cosmic Phys., 5, 287
Togi, A., & Smith, J. D. T. 2016, ApJ, 830, 18

Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050

Vallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813, 36

Vizgan, D., Greve, T. R., Olsen, K. P., et al. 2022a, ApJ, 929, 92

Vizgan, D., Heintz, K. E., Greve, T. R., et al. 2022b, ApJ, 939, L1

Walter, F., Neeleman, M., Decarli, R., et al. 2022, ApJ, 927, 21

Weibel, A., de Graaff, A., Setton, D. J., et al. 2024, arXiv e-prints
[arXiv:2409.03829]

Williams, C. E., Lake, W., Naoz, S., et al. 2024, ApJ, 960, L16

Witten, C., McClymont, W., Laporte, N., et al. 2024, A&A, submitted
[arXiv:2407.07937]

Woods, T. E., Patrick, S., Whalen, D. J., & Heger, A. 2024, ApJ, 960, 59

Yamaguchi, Y., Kohno, K., Tamura, Y., et al. 2017, ApJ, 845, 108

Yan, L., Sajina, A., Loiacono, F., et al. 2020, ApJ, 905, 147

Yang, S., Popping, G., Somerville, R. S., et al. 2022, ApJ, 929, 140

Yoshida, N., Sokasian, A., Hernquist, L., & Springel, V. 2003, ApJ, 591, L1

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976
## Appendix A: [C II] 158 µm line emission data

Here we provide the atomic data adopted in this work and in the companion paper C24 to compute the [C II] 158 µm line emission. We model the fine-structure transition (2p)[$^2$P$_{3/2}$ – $^2$P$_{1/2}$] between the quantum number J = 3/2 and J = 1/2 of C$^+$ as a two-level system considering e$^-$ impacts, H impacts and H$_2$ impacts. Collisional rates, spontaneous-transition rate, and energy level separation are taken from Hollenbach & McKee (1989), Maio et al. (2007) and Goldsmith et al. (2012) and listed below:

$$ \gamma_{21}^{e} = 2.8 \cdot 10^{-7} T_{100}^{-0.5} \text{ cm}^{3} \text{ s}^{-1}, $$

$$ \gamma_{21}^{H} = 8 \cdot 10^{-10} T_{100}^{0.07} \text{ cm}^{3} \text{ s}^{-1}, $$

$$ \gamma_{21}^{H_2} = 3.8 \cdot 10^{-10} T_{100}^{0.14} \text{ cm}^{3} \text{ s}^{-1}, $$

$$ A_{21} = 2.4 \cdot 10^{-6} \text{ s}^{-1}, $$

$$ \Delta E_{21} = 1.259 \cdot 10^{-14} \text{ erg}, $$

where the indeces 2 and 1 refer to the transition upper and lower levels, T$_{100}$ = T/(100 K) with T gas temperature, and ΔE$_{21}$ corresponds to the transition wavelength (frequency) of 157.74 µm (1900.5 GHz), while the excitation temperature is about 91.2 K. With these inputs, the [C II] radiative losses can be computed as:

$$ \Lambda_{[C II]} = \frac{(n_{H_2}\gamma_{12}^{H_2} + n_H\gamma_{12}^{H} + n_e\gamma_{12}^{e}) n_{C^+} A_{21} \Delta E_{21}}{n_{H_2}(\gamma_{12}^{H_2} + \gamma_{21}^{H}) + n_H(\gamma_{12}^{H} + \gamma_{21}^{H}) + n_e(\gamma_{12}^{e} + \gamma_{21}^{e}) + A_{21}}, \quad (A.1) $$

which highlights the contributions of H$_2$, H I and e$^-$ (ionised/H II) gas phases for given H$_2$, H I, e$^-$ and C$^+$ number densities (n$_{H_2}$, n$_H$, n$_e$, n$_{C^+}$, respectively). The relation between excitation (1 → 2 transition) and de-excitation (2 → 1 transition) rates satisfies the Boltzmann limit (see e.g. Maio et al. 2007). The [C II] luminosity emitted by a gas volume V can be estimated as L$_{[C II]}$ = Λ$_{[C II]}$V.