# COLDSIM predictions of [C II] emission in primordial galaxies

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

► To cite this version:

Benedetta Casavecchia, Umberto Maio, Céline Péroux, Benedetta Ciardi. COLDSIM predictions of [C II] emission in primordial galaxies. Astronomy & Astrophysics - A&A, 2024, 689, 10.1051/0004-6361/202450332. insu-04726414

HAL Id: insu-04726414

https://insu.hal.science/insu-04726414v1

Submitted on 9 Oct 2024

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 Attribution 4.0 International License
# COLD SIM predictions of [C II] emission in primordial galaxies

Benedetta Casavecchia¹,*, D., Umberto Maio²,³, Céline Péroux⁴,⁵ D., and Benedetta Ciardi¹

¹ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, 85748 Garching bei München, Germany

² INAF-Italian National Institute of Astrophysics, Observatory of Trieste, Via G. Tiepolo 11, 34143 Trieste, Italy

³ IFPU-Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

⁴ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany

⁵ Aix Marseille Université, CNRS, Laboratoire d'Astrophysique de Marseille (LAM) UMR 7326, 13388 Marseille, France

Received 11 April 2024 / Accepted 1 June 2024

## ABSTRACT

**Context.** A powerful tool with which to probe the gas content at high redshift is the [C II] 158 µm submillimetre emission line, which, due to its low excitation potential and luminous emission, is considered a possible direct tracer of star forming gas.

**Aims.** In this work, we investigate the origin, evolution, and environmental dependencies of the [C II] 158 µm emission line, as well as its expected correlation with the stellar mass and star formation activity of the high-redshift galaxies observed by JWST.

**Methods.** We use a set of state-of-the-art cold-gas hydrodynamic simulations (COLD SIM) with fully coupled time-dependent atomic and molecular non-equilibrium chemistry and self-consistent [C II] emission from metal-enriched gas. We accurately track the evolution of H I, H II, and H$_2$ in a cosmological context and predict both global and galaxy-based [C II] properties.

**Results.** For the first time, we predict the cosmic mass density evolution of [C II] and find that it is in good agreement with new measurements at redshift $z = 6$ from high-resolution optical quasar spectroscopy. We find a correlation between [C II] luminosity, $L_{[CII]}$, and stellar mass, which is consistent with results from ALMA high-redshift large programs. We predict a redshift evolution in the relation between $L_{[CII]}$ and the star formation rate (SFR), and provide a fit to relate $L_{[CII]}$ to SFR, which can be adopted as a more accurate alternative to the currently used linear relation.

**Conclusions.** Our findings provide physical grounds on which to interpret high-redshift detections in contemporary and future observations, such as the ones performed by ALMA and JWST, and to advance our knowledge of structure formation at early times.

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

## 1. Introduction

The last few years have seen the beginning of a golden era in the search for galaxies in the epoch of reionisation. This has been made possible by a combination of space facilities, such as the Hubble Space Telescope (HST) and 8–10 metre diameter ground-based telescopes (e.g. the ESO Very Large Telescope, VLT, and Keck). Dust and metals have been identified in some high-redshift ($z$) objects with the Atacama Large Millimeter/submillimeter Array (ALMA). These datasets are now being supplemented by cutting-edge observations with the James Webb Space Telescope (JWST). This observational progress has naturally given rise to a plethora of theoretical models, ranging from analytic and semi-analytic approaches to numerical simulations.

Therefore, it is now possible to detect and interpret the properties of early galaxies during the epoch of reionisation. As star-forming galaxies are typically bright in the UV due to the emission from young and massive stars, UV luminosity can be used as a primary tracer to quantify their star formation rate (SFR). However, part of this radiation is absorbed by interstellar dust and re-emitted at far-infrared (FIR) wavelengths (Fixsen et al. 1998; Dole et al. 2006; Burgarella et al. 2013; Whitaker et al. 2017; Salim & Narayanan 2020). For this reason, accurate estimates of the SFR depend on a multi-wavelength analysis ranging from the rest-UV, to the far-infrared (FIR) and redshifted to millimetre (mm) wavelengths, which thus includes both com-

ponents of star formation; that is, the unobscured emission and the thermal emission from dust (Béthermin et al. 2020; Khusanova et al. 2021; Ferrara et al. 2022; Bowler et al. 2024).

As stars are also responsible for gas metal enrichment (Sargent et al. 1988; Songaila 2001; Becker et al. 2015; Codoreanu et al. 2017; Tumlinson et al. 2017; Péroux & Howk 2020), an additional source of information on star formation comes from observations of the metallicity ($Z$) of the gas. The global evolution of the metal mass in the Universe can be captured by the mass density parameter $\Omega$ of single species. Among the possible tracers of $Z$, carbon is one of the most observed at all redshifts, given its high luminosity (Simcoe et al. 2006; D’Odorico et al. 2013, 2022; Davies et al. 2023). However, it remains challenging to detect the total amount of carbon mass, because this would require that each ionisation state (C I, C II, C III, C IV, etc.) be observed individually given the different ionisation potentials. Moreover, some of these ions are known to be weak due to their low oscillation strength. Historically, observational campaigns have therefore focused on estimates of C IV, which is the most common ion in hot gas and is expected to be the dominant ionisation state at redshifts of $z < 5$. Bright quasars have been used as background sources to probe the incidence of metal-rich absorbers along the line of sight, and to quantify the cosmological C IV mass density parameter, $\Omega_{CIV}$ (Songaila 2001; Pettini et al. 2003; Boksgaard et al. 2003; Simcoe et al. 2006, 2020; D’Odorico et al. 2013, 2022). Various authors have used large samples of C IV absorbers and have reported a significant increase in $\Omega_{CIV}$ towards lower redshifts (see also Becker et al.

<footnote>* Corresponding author; benecasa@mpa-garching.mpg.de</footnote>
2009; Ryan-Weber et al. 2009; Simcoe et al. 2011). Recently, Davies et al. (2023) employed a large number of VLT/X-shooter high-resolution quasar spectra to measure $\Omega_{CIV}$, finding that it decreases by a factor of $4.8 \pm 2.0$ over the ~300 Myr interval between $z \sim 4.7$ and ~5.8. These authors interpret this feature as an indication of rapid changes in the gas ionisation state driven by the evolution of the UV background at the end of hydrogen reionisation. They also show that $C^+$ is the dominant state at $z > 5$ and for the first time report observations of $\Omega_{C^+}$. The strength of probing one ionisation state in particular as opposed to global carbon is that, in this way, the incident ionising background can be directly constrained.

Beyond providing information on the global metal enrichment of the Universe, the [C II] emission line in particular has been very valuable for studying high-redshift galaxies. Ultra-deep observations with HST and JWST are now detecting the UV stellar light in galaxies up to $z \gtrsim 10$ (Bouwens et al. 2019; Bowler et al. 2020; Leethochawalit et al. 2023a,b; Donnan et al. 2023, 2024; Harikane et al. 2023). In this context, redshifted mm observations provide an important means to both confirm the high-redshift candidates and determine their physical properties (Fujimoto et al. 2019; Schaerer et al. 2020; Ferrara et al. 2022; Sommovigo et al. 2021, 2022; Pozzi et al. 2024; Romano et al. 2024). To this end, the [C II] 158 µm (1900.5 GHz) fine-structure transition ($^2P_{3/2} \rightarrow ^2P_{1/2}$) has proven to be key, as it is the brightest ionic emission line in star-forming galaxies. Due to its low excitation potential, [C II] emission can be produced in photodissociation regions, neutral atomic gas, and molecular clouds (Sargsyan et al. 2012; Pineda et al. 2013; Rigopoulou et al. 2014; Velusamy & Langer 2014; Cormier et al. 2015; Glover & Clark 2016; Nordon & Sternberg 2016; Díaz-Santos et al. 2017; Fahrion et al. 2017). Moreover, this line is not strongly affected by dust absorption and its relation with the SFR has already been reported and calibrated in the nearby Universe (Stacey et al. 1991; Leech et al. 1999; Boselli et al. 2002; Sargsyan et al. 2012; De Looze et al. 2011, 2014; Herrera-Camus et al. 2015; Zanella et al. 2018; Madden et al. 2020). For these reasons, [C II] emission is an excellent tool for investigating the properties of high-z structures, including their SFR and kinematics, as well as their neutral atomic and molecular gas content (Cicone et al. 2015; Hernandez-Monteagudo et al. 2017; Padmanabhan 2019; Dessauges-Zavadsky et al. 2020; Schaerer et al. 2020; Ginolfi et al. 2020; Aravena et al. 2024; Jones et al. 2024; Kaasinen et al. 2024; Posses et al. 2024; Rowland et al. 2024). Recently, this bright line has gained in popularity as ALMA is detecting an increasing number of objects with [C II] emission at $z \gtrsim 4$. While ALMA is expected to observe the redshifted [C II] emission line from $z > 1$, in practice its higher bands (i.e. 8, 9 and 10) at frequencies of >385 GHz are more challenging to schedule, because they require lower precipitable water vapour conditions, which rarely occur. In addition, observations with these frequencies offer a smaller primary beam or field of view, which limits the number of objects that can be detected in a given exposure time and therefore decreases the survey efficiency at the telescope. For this reason, most of the available surveys are done at $z > 4$, making ALMA an ideal instrument to observe high-z cold gas. Recently, two major ALMA campaigns observed selected UV-bright galaxies at $z \sim 4-8$. The ALPINE survey (Le Fèvre et al. 2020) reported over 100 galaxies in the range of $z = 4.5-6$, while the REBELS survey (Bouwens et al. 2022) identified dozens of galaxies at $z = 6.5-7.7$. Some works found a clear correlation between $L_{[C II]}$ and SFR (Dessauges-Zavadsky et al. 2020; Schaerer et al.

2020; Fujimoto et al. 2021; Ferrara et al. 2022; Schouws et al. 2022), while others found a [C II] luminosity below the $L_{[C II]}$-SFR linear relation (Maiolino et al. 2015; Pentericci et al. 2016; Laporte et al. 2019; Carniani et al. 2020; Fujimoto et al. 2024). However, the uncertainties are large at such high redshifts and the $L_{[C II]}$ versus SFR relationship (Harikane et al. 2020; Romano et al. 2022) is calibrated only for the most massive galaxies, making it inaccurate in the lower SFR regime. The CONCERTO instrument (CONCERTO Collaboration 2020; Gkogkou et al. 2023) on the Atacama Pathfinder Experiment telescope (APEX) and purpose-built facilities such as the Atacama Large Aperture Submillimeter Telescope (AtLAST; Klaassen et al. 2020) will provide further samples of [C II] detections in the future.

Theoretical models play a central role in elucidating which mechanisms influence the emission of $C^+$. Three main approaches have been proposed to derive the [C II] luminosity. Semi-analytical models (see e.g. Lagache et al. 2018; Popping et al. 2019) use empirical relations to quantify $L_{[C II]}$ in a cosmological context. The strength of these methods is the fast runtime of the calculations, which leads to the possibility of exploring a large parameter space. These models are also implemented in order to explain single-halo properties, such as the extended (~10 kpc) [C II] emission around galaxies at $z = 6$ generated by supernova-driven cooling outflows (Pizzati et al. 2020, 2023). Zoom-in simulations of a limited number of halos are able to capture the small-scale physics at play within a cosmological framework (see e.g. Katz et al. 2019; Bisbas et al. 2022; Lahén et al. 2024; Schimek et al. 2024). Finally, simulations run in cosmological boxes offer a complementary avenue, which provides much larger statistics while remaining somewhat limited in terms of mass resolution, in particular for the cold gas (see e.g. Leung et al. 2020; Kannan et al. 2022; Garcia et al. 2023; Katz et al. 2023; Liang et al. 2024). Importantly, we stress that modelling below the hydrodynamic simulation resolution ('subgrid' modelling) is required to address this phase of the gas in both zoom-in and cosmological simulations. In this respect, simple analytic techniques, for example those based on either gas metallicity or pressure (Blitz & Rosolowsky 2006; Krumholz 2013), are often used to split neutral hydrogen into its atomic and molecular components and infer possible [C II] properties (Lagos et al. 2015; Popping et al. 2019; Popping & Péroux 2022; Szakacs et al. 2022; Vizgan et al. 2022a,b). Alternatively, different authors have suggested either idealised models of [C II] line emission (Olsen et al. 2017), or post-processing of zoom-in simulations; these latter, although lacking the [C II] treatment at runtime, have provided some interesting information on the resolved properties of specific individual haloes (Pallottini et al. 2019). Simulating the atomic and molecular phases of the cold gas in a full cosmological context, coupled to a consistent implementation for [C II] emission, will therefore be one of the most crucial and challenging objectives of studies of galaxy formation and evolution in the decades to come.

The goal of this work is to make a step forward in capturing the complex physical processes at play by running and analysing a set of the state-of-the-art hydrodynamic cosmological simulations, COLDSIM (Maio et al. 2022), which include fully coupled time-dependent atomic and molecular non-equilibrium chemistry, and lead to self-consistent calculations of [C II] emission. By extending our previous work on primordial molecular gas, we focus here on the [C II] line emission in the early Universe and its correlations with the physical properties of the host halo. Specifically, we study the evolution of the $C^+$ mass density parameter,
$\Omega_{C^+}$, and the dependence of [C II] luminosity, $L_{[CII]}$, on stellar mass, $M_\star$, and SFR in the redshift range $z \sim 6-12$.

The paper is organised as follows: Section 2 presents the hydrodynamical cosmological simulations used in this study, as well as a detailed description of the methodology to compute the [C II] luminosity. Section 3 details our analysis and associated results. We discuss the impact of our findings in Section 4, and present our conclusions in Section 5. We adopt a standard cosmology with cold dark matter and cosmological constant $\Lambda$. The expansion parameter is $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$ while baryon, matter, and $\Lambda$ density parameters are $\Omega_{b,0} = 0.045$, $\Omega_{m,0} = 0.274$, and $\Omega_{\Lambda,0} = 0.726$, respectively. We adopt the notation cMpc to indicate lengths in comoving megaparsecs and, unless otherwise stated, express masses and metallicities in solar units ($M_\odot$ and $Z_\odot$, respectively).

## 2. Methods

In this section, we introduce the simulations used for our analysis, as well as the methodology adopted to calculate the amount of [C II] produced and the associated luminosity.

### 2.1. The COLDSIM numerical simulations

The set of cosmological hydrodynamical simulations used throughout this work is based on the larger COLDSIM project (Maio et al. 2022). The latter explores the formation of galaxies at high redshift by employing a physics- and chemical-rich implementation of the evolution of atomic and molecular gas in primordial environments, and represents a remarkable extension of our previous works (Maio et al. 2007, 2010, 2011, 2016; Biffi & Maio 2013; Maio & Tescari 2015; Ma et al. 2017a,b). In addition, the ColdSIM simulations have been successfully validated against available observational constraints and recent determinations of the molecular content of $z > 7$ galaxies (see Feruglio et al. 2023). Their implementation relies on an extended version of the parallel code P-Gadget3 (Springel 2005), which, in addition to gravity and smoothed particle hydrodynamics (SPH), includes an ad hoc time-dependent non-equilibrium network of first-order differential equations to solve the processes of ionisation, dissociation, and recombination (Abel et al. 1997; Yoshida et al. 2003; Maio et al. 2007) for the following species: $e^-$, $H$, $H^+$, $H^-$, He, $He^+$, $He^{++}$, $H_2$, $H_2^+$, D, $D^+$, HD, and He$H^+$. In addition to pristine-gas chemistry, our implementation also includes metal radiative losses from resonant and fine-structure transitions in metal-enriched gas, such as [C II] 158 µm emission. Here, we use the same chemical and physical implementation as in Maio et al. (2022) with slight differences related to assumptions about some model parameters, box size, and the resolution of the runs (see details below). In the following, we briefly summarise the characteristics of the COLD SIM runs analysed in this work and refer the reader to the original paper for more details.

The primordial chemistry network includes $H_2$ formation through $H^-$ catalysis ($H + e^- \rightarrow H^- + \gamma$ and $H^- + H \rightarrow H_2 + e^-$), $H_2^+$ catalysis ($H + H^+ \rightarrow H_2^+ + \gamma$ and $H_2^+ + H \rightarrow H_2 + H^+$), and three-body interaction ($3H \rightarrow H_2 + H$). In gas enriched by stellar feedback, dust-grain catalysis is also accounted for. In detail, $H_2$ dust-grain catalysis is coupled to the non-equilibrium network and followed for different gas temperatures and metallicities by assuming a Z-dependent dust-to-metal ratio and a grey-body dust emission with index $\beta = 2$. This allows us to estimate the redshift evolution expected for dust-grain temperature in agreement with recent ALMA constraints (da Cunha et al. 2021).

<table><caption>**Table 1.** Simulations characteristics.</caption><thead><tr><th>Name</th><th><i>L</i><br>[cMpc/h]</th><th><i>N</i><sub>part</sub></th><th><i>m</i><sub>DM</sub><br>[M<sub>☉</sub>/h]</th><th><i>m</i><sub>gas</sub><br>[M<sub>☉</sub>/h]</th></tr></thead><tbody><tr><td>CDM Ref</td><td>10</td><td>2 × 512<sup>3</sup></td><td>9.5 × 10<sup>5</sup></td><td>4.7 × 10<sup>4</sup></td></tr><tr><td>CDM HR</td><td>10</td><td>2 × 1000<sup>3</sup></td><td>6.3 × 10<sup>4</sup></td><td>1.3 × 10<sup>4</sup></td></tr><tr><td>CDM LB</td><td>50</td><td>2 × 1000<sup>3</sup></td><td>7.9 × 10<sup>6</sup></td><td>1.6 × 10<sup>6</sup></td></tr></tbody></table>

**Notes.** From left to right, the columns refer to the name of the simulation, the box side $L$, the number of particles $N_{part}$, the mass of DM particles $m_{DM}$, and the initial mass of gas particles $m_{gas}$.

The chemical network is coupled with the evolution of the heavy elements He, C, N, O, Ne, Mg, Si, S, Ca, and Fe, which are traced individually and spread by stellar feedback through asymptotic giant branch winds and Type Ia and Type II supernova explosions (Maio et al. 2010, 2016). Collisionless star particles are formed stochastically at each integration time according to the star formation model by Springel & Hernquist (2003), which has been slightly modified to follow molecular runaway cooling. Stellar particles are considered as a single stellar population with a Salpeter initial mass function (IMF). Chemical enrichment from stellar particles is implemented through SPH kernel smoothing. A UV background is turned on at $z \approx 6$ using the standard Haardt & Madau (1996) prescription. Additional physical processes, such as HI and $H_2$ self-shielding, the photoelectric effect, and cosmic-ray heating are also implemented, as they are likely to play an important role in the cooling of molecular gas (Habing 1968; Draine & Sutin 1987; Bakes & Tielens 1994).

As mentioned, we considered three simulations; these explore different box sizes and mass resolutions to investigate how the physical and chemical processes described above affect the properties of galaxies, particularly their atomic and molecular budgets and the resulting [C II] emission at various scales. The reference run, dubbed CDM Ref, has a box size of $10 \text{ cMpc/h}$ with $2 \times 512^3$ particles evenly divided between gas and dark matter. The simulation with higher resolution (CDM HR) has the same box size, but the cosmic field is sampled with $2 \times 1000^3$ particles. Finally, we employed a larger box (CDM LB) of $50 \text{ cMpc/h}$ with $2 \times 1000^3$ particles to account for rarer, larger structures that would be missing in the other runs. The characteristics of these simulations are listed in Table 1. In each cosmological volume, cosmic structures have been identified via a friends-of-friends algorithm and a substructure finder. Only objects with at least 100 total particles were considered. The detailed chemical properties of each simulated halo were retrieved in post-processing by matching the halo-constituting particles and the snapshot particles.

As an example, the top row of Fig. 1 displays gas column density maps at $z = 12$ and 6 obtained by projecting a slice of the high-resolution box along its orthogonal dimension. Cosmological evolution shapes the cosmic medium in dense and void regions (visible in the maps) through the balance between heating and cooling processes. This latter leads to gas collapse and galaxy formation in the densest filamentary environments, as well as metal pollution on stellar timescales. The injection of carbon atoms – synthesized by stellar interiors and tracked by our modelling – into the surrounding medium is then responsible for the [C II] emission signal displayed in the bottom row (see further discussion below).

In Figure 2 we show the SFR density, $\Psi$, obtained from the simulations, together with observations based on near-
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** CDM HR column density maps and [C II] surface-brightness. Top row: Column density maps at z = 12 and z = 6 obtained by projecting a slice passing through the centre of the CDM HR simulation box along a width equal to one-tenth of the box side (10 cMpc/h) and discretising it on a grid of 512 × 512 pixels. This results in a pixel linear size of about 19.5 ckpc/h. Bottom row: [C II] surface-brightness maps at z = 12 and z = 6 obtained by projecting the [C II] emission from the same slice passing through the centre of the CDM HR simulation box.</figcaption></figure>

infrared (NIR; Merlin et al. 2019; Bhatawdekar et al. 2019), dust-corrected UV (Harikane et al. 2022; Donnan et al. 2023) data, as well as recent JWST spectroscopic data (Harikane et al. 2024). We note that observational samples vary in the way the SFR is calculated, adding important uncertainties to the comparison with simulated data. Additionally, NIR estimates are more accurate than UV measurements because they are less prone to dust effects. At the highest redshifts, the results from CDM HR (purple line) are above those from the other simulations. The higher resolution of CDM HR leads to an earlier collapse of gas and consequent star formation. Similarly, stellar feedback becomes effective in quenching star formation at earlier times, meaning that the trend is reversed at lower redshifts. We find that our predictions reproduce the observed trend at z ≥ 6, although they lie above the UV estimates. Nevertheless, we note that the

more reliable estimates from NIR data (e.g. Merlin et al. 2019; Bhatawdekar et al. 2019) give higher values of Ψ, albeit in some cases with large error bars.

Figure 3 shows the galaxy SFR as a function of stellar mass $M_*$ at different redshifts. We find a strong correlation between the stellar mass of galaxies and their SFR for sources on the galaxy main sequence, and so studying its evolution at different redshifts informs us about the efficiency with which galaxies transform gas into stars. We compare our results to a compilation of data from the ALMA large programs ALPINE (Béthermin et al. 2020; Le Fèvre et al. 2020) and REBELS (Bouwens et al. 2022; Ferrara et al. 2022), as well as to JWST results from Nakajima et al. (2023) regarding the ERO, GLASS, and CEERS surveys, from Curti et al. (2023) for the JADES survey, and to the few detections at $z > 10$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Redshift evolution of the SFR density, Ψ, for CDM HR (purple line), CDM Ref (blue), and CDM LB (green). Data points are obtained from NIR (Merlin et al. 2019, triangles, and Bhatawdekar et al. 2019, pentagons), dust-corrected UV (Donnan et al. 2023; circles), and JWST spectroscopic data (Harikane et al. 2024; diamonds). The latter are lower limits.</figcaption></figure>

reported by Harikane et al. (2023). Globally, our SFR–M* relation agrees well with the trends of most observations. We note that the predicted SFR of individual haloes is sometimes lower than observed, while in Figure 2 the SFR density indicates that our predictions are slightly higher than observations. This could indicate that the number of haloes per unit volume in the simulations is higher than the observed number, which is likely because we simulate haloes with masses of lower than the present detection limit. The JWST observations reported by Curti et al. (2023) and Nakajima et al. (2023) reveal a systematically higher SFR for a given stellar mass. Also, the ALMA-based REBELS observations report SFR values that are higher than predicted, albeit with large error bars. We note that the data from Nakajima et al. (2023) are averaged over galaxies spanning the range z = 3–10, and hence might be biased towards high SFRs, and it is possible that other observational studies are similarly biased. Interestingly, the ALMA-based ALPINE results at z ≈ 6 are in good agreement with our predictions from CDM LB.

Finally, in Figure 4 we show the gas metallicity, $Z_{\text{gas}}$, as a function of stellar mass at various redshifts. Our results indicate a strong relation among these two quantities. ColdSIM implementation features a non-equilibrium chemical network that includes the individual evolution of heavy elements. Therefore, $Z_{\text{gas}}$ is computed for each galaxy as the sum of the masses of the elements heavier than helium within the gas particles divided by the total gas mass. Metal yields are calculated at each time step self-consistently with stellar IMF, lifetimes, and feedback mechanisms. The metal production from SNe and AGB stars comes naturally from the simulations and is consistent with the stellar evolution model. Non-equilibrium chemistry is tracked while accounting for creation and destruction processes, as well as heating and cooling in the primordial Universe, as described in for example Maio et al. (2007, 2016, 2022). In particular, Maio et al. (2010) have already shown that carbon and oxygen roughly trace each other at early times. Therefore, the comparisons made in Figure 4 with oxygen-based JWST metallicities should be valid (Curti et al. 2023; Nakajima et al. 2023; Tacchella et al. 2022). The mass–metallicity relation seems to be recovered quite well, with theoretical predictions matching both

the trend and the scatter of the observed determinations. We note that observational estimates of gas metallicities are extremely challenging at these early times and the error bars are very large. Beyond the relations of the physical galaxy properties reported here, we remind the reader that Maio & Viel (2023) present additional diagnostics, including UV luminosity functions, halo mass functions, and stellar mass fractions.

## 2.2. Galaxy [CII] luminosity

Here we describe how we estimate the galaxy [C II] emission at 158 µm, which is associated with [C II] ions excited to the $^2$P$_{3/2}$ level. The [C II] luminosity of a galaxy, $L_{[C II]}$, is evaluated as the sum of the luminosities of each gas particle, $L_{[C II],p}$, contained within it. This latter is estimated as $L_{[C II],p} = \Lambda_{[C II]}V_p$, where $V_p$ is the volume of a sphere with a radius equal to the gas particle smoothing length, and $\Lambda_{[C II]}$ is the power emitted per unit volume as a consequence of the atomic fine-structure transition. In the specific case of the [C II] 158 µm emission modelled as a two-level system (Hollenbach & McKee 1989; Goldsmith et al. 2012; Maio et al. 2007), we can write

$$ \Lambda_{[C II]} \equiv n_{C^+,2} A_{21} \Delta E_{21}, \quad (1) $$

where $n_{C^+,2}$ is the number density of the C$^+$ ions excited to the upper state, $A_{21}$ is the Einstein coefficient for spontaneous emission, and $\Delta E_{21}$ is the energy level separation. The equation above can be rewritten in terms of the total C$^+$ density $n_{C^+}$ as (see e.g. Maio et al. 2007)

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

where $\gamma_{12}^{H_2}$, $\gamma_{12}^{H}$, and $\gamma_{12}^{e}$ represent the H$_2$-impact, H-impact, and e$^-$-impact collisional excitation rates, $\gamma_{21}^{H_2}$, $\gamma_{21}^{H}$, and $\gamma_{21}^{e}$ denote the rates of de-excitation processes (see the appendix for specific values), while $n_{C^+}$, $n_{H_2}$, $n_H$, and $n_e$ are the number densities of C$^+$, H$_2$, atomic hydrogen, and electrons, respectively.

We note that the resulting [C II] emission is a function of temperature, density, and metallicity, and that the emission flux will be especially influenced by the thermal state of the gas below ~10⁴ K, where the coefficient rates of two-body reactions are valid Maio et al. (2007), McElroy et al. (2013). Above this temperature, carbon is excited at higher excitation levels and does not contribute to [C II] emission.

The bottom row of Fig. 1 shows a pictorial view of the expected [C II] surface brightness stemming from a slice of the CDM HR box. For each gas particle, [C II] emission is computed according to the actual content, density, and temperature of carbon according to Eq. (2). The derived [C II] signal in each pixel is given by the integrated emitted power along the line of sight of the pixel divided by the pixel area. By comparing with the upper panels of Figure 1, we observe that [C II] emission arises from regions with column densities larger than roughly 10¹⁹ cm⁻² and traces metal-enriched parts of the filamentary structures constituting the cosmic medium. This is in broad agreement with C II absorption from damped Lyman-α or subdamped Lyman-α systems Sebastian et al. (2024) as well as with preliminary estimates of bright extended [C II] haloes at redshift $z \approx 6-7$ Pizzati et al. (2023), Bischetti et al. (2024). Tenuous, diffuse gas that lies far away from galaxy build-up sites
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Star formation rate as a function of stellar mass at z = 12, 10, 8, and 6. The coloured dots refer to predictions from CDM HR (purple), CDM Ref (blue), and CDM LB (green). Grey symbols are data from ALPINE (squares; Béthermin et al. 2020; Le Fèvre et al. 2020), REBELS (crosses; Bouwens et al. 2022; Ferrara et al. 2022), and JWST results described in Nakajima et al. (2023, stars), Curti et al. (2023, pentagons), and Harikane et al. (2023, diamonds). Globally, our results agree well with observations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Gas metallicity $Z_{\text{gas}}$ as a function of stellar mass $M_{\star}$ at z = 12, 10, 8, and 6. The solid lines represent the mean values for CDM Ref (blue), CDM HR (purple), and CDM LB (green), while the darker (lighter) shaded regions represent the $1\sigma$ (2$\sigma$) standard deviation. Grey symbols refer to JWST results described in Nakajima et al. (2023, stars), Curti et al. (2023, pentagons), and Tacchella et al. (2022, triangles). Globally we find that our simulations are in good agreement with JWST observations.</figcaption></figure>

is either little enriched with metals (including carbon) or too rarefied to produce significant amounts of [C II] surface brightness. More quantitative findings about features and dependencies of the [C II] 158 µm signal are presented in the following section.

## 3. Results

Here we discuss the results of our simulations in terms of evolution of the cosmological carbon mass density and [C II] luminosity functions.

### 3.1. Cosmological carbon mass density

In this section, we analyse the total carbon content and C+ content in the full cosmological boxes using the corresponding cosmic mass density parameter $\Omega_X \equiv \rho_X/\rho_{\text{crit},0}$, where $\rho_X$ is the comoving mass density of element X and $\rho_{\text{crit},0}$ is the present-day critical density.

While past simulations made predictions as to the cosmological evolution of metal mass densities (e.g. Maio & Tescari 2015; Bird et al. 2016) or C IV at low redshift (Tescari et al. 2011), here

we investigate the evolution of $\Omega_C$ and $\Omega_{C^+}$ during the epoch of reionisation, as shown in Figure 5. The upper panel shows that, in all simulations, $\Omega_C$ increases with cosmic time as a natural consequence of the build up of metals. At the higher redshifts, the predictions from CDM HR (purple line) lay above those from the other models because its higher resolution leads to an earlier star formation, and hence metal enrichment. For the same reason, in CDM HR, stellar feedback becomes effective in quenching star formation and metal pollution at earlier times, and so the trend is reversed at z ~ 9.5 (6.5) with respect to CDM Ref (CDM LB).

As mentioned in Sect. 1, observational estimates of the total C mass density parameter are not available, and thus our results can only be compared to estimates from C II and C IV detections (D'Odorico et al. 2022; Davies et al. 2023), which should be considered as lower limits for $\Omega_C$. As our simulated values lay above the data, we conclude that our estimates are consistent with observations.

The bottom panel of Figure 5 shows the evolution of the C+ mass density parameter, which presents features similar to those of $\Omega_C$. Indeed, the flattening of $\Omega_{C^+}$ here is also related to the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Redshift evolution of the total C mass density, $\Omega_C$ (top panel), and the $C^+$ mass density, $\Omega_{C^+}$ (bottom panel). The lines refer to results from CDM Ref (light blue), CDM HR (purple), and CDM LB (green), while symbols refer to observations of C IV (diamonds; D'Odorico et al. 2022), C II + C IV (stars; Davies et al. 2023), and C II (triangles; Davies et al. 2023).</figcaption></figure>

earlier effect of stellar feedback in the CDM HR simulation. In this case, a direct comparison to observational data (Davies et al. 2023) is possible and we find excellent agreement with our predictions from CDM LB and CDM Ref, while results from CDM HR are smaller by about a factor of 2. The fact that $\Omega_C$ is consistent with observations and the differences in $\Omega_{C^+}$ among simulations are larger than those in $\Omega_C$ at $z < 7$ suggests that such differences are likely not due to chemical evolution but rather to discrepancies in the ionisation state possibly associated to stellar feedback.

## 3.2. [CII] luminosity function

In this section, we discuss the [C II] luminosity functions obtained from our simulations and compare them to available high-$z$ data.

### 3.2.1. Dependence on stellar mass

In Figure 6, we show $L_{[CII]}$ as a function of stellar mass, $M_\star$, at different redshifts. We observe that, typically, more massive galaxies are more luminous in [C II], and that for a given stellar mass, $L_{[CII]}$ increases with decreasing redshift. This is expected from the build up of stars in galaxies and the ongoing metal enrichment. We also note that in CDM LB we are able to sample bigger objects because of the larger box size. The scatter in these relations is of the order of 1 dex, and increases in the low-SFR regime, especially at stellar masses below $10^{6.5} M_\odot$ for CDM HR, which is likely related to the fact that these low-mass haloes have fewer particles. We note a slightly enlarged scatter of 2 dex in the CDM LB, especially for galaxies with stellar masses of $10^{8-9} M_\odot$. We observe that the curves present an irregular trend, with a turnover that becomes visible in the high-mass end, mostly at high redshift. This is discussed further in the following section.

We compare our results with the ALPINE (Béthermin et al. 2020; Le Fèvre et al. 2020) and REBELS (Bouwens et al. 2022; Ferrara et al. 2022) ALMA large programs, which target [C II] emission at $z = 4.5-6$ and $z = 6.5-7.7$, respectively. These observations probe rare bright haloes with stellar masses of ~$10^9-10^{10} M_\odot$ that might even resemble early galaxy groups rather than normal galaxies. While at $z = 6$ ALPINE detections align well with the trends inferred from the simulated samples, REBELS galaxies have typical luminosities for a given stellar mass that are higher than both the simulated haloes and the ALPINE data. At such redshifts, there are only a handful of measurements and these are probably biased towards the brightest $L_{[CII]}$ emitters. More precise comparisons will only be possible with future, deeper ALMA observations targeting fainter galaxies.

### 3.2.2. Dependence on star formation rate

The top row of Figure 7 shows the estimated $L_{[CII]}$ as a function of SFR at $z=12, 10, 8$, and 6. Consistently with the previous section, objects with higher SFR are typically more luminous, with a global scatter around the mean value of 1 dex. We note that CDM LB has a slightly larger scatter than the other runs, likely because it contains a greater number of haloes. We stress that there is a redshift evolution in the amplitude of the relation, which increases by about 1 dex from $z=12$ to $z=6$. An important consequence of this result is that, using the observed $L_{[CII]}$ to estimate the SFR of a galaxy by means of the $L_{[CII]}$–SFR relation calibrated in the local Universe (e.g. De Looze et al. 2014) can lead to inaccurate estimates of high-$z$ SFRs.

At each redshift, we see a turnover in the [C II] luminosity for galaxies with higher SFR, as also found in the $L_{[CII]}-M_\star$ behaviour of Section 3.2.1. We believe that this feature is mainly determined by feedback effects, in particular starbursts, which strongly affect the physical properties of the [C II]-emitting gas. A similar feature has indeed been found by Katz et al. (2023) in simulated ‘bursty leakers’, with SFR averaged over 10 Myr, and also in the radiation hydrodynamic simulations SPICE by Bhagwat et al. (in prep.) in their ‘bursty-sn’ model. A less prominent role is additionally played by numerical effects, as suggested by the fact that the feature is more pronounced in CDM LB, because a higher resolution is required to properly model the gas phase that gives rise to [C II] emission.

Although data available at these epochs are sparse and probably biased towards large, bright structures, we highlight that the $L_{[CII]}$–SFR trend of the simulated galaxies aligns very
<figure><img src="image_7.png" /><figcaption>Fig. 6. [C II] luminosity as a function of halo stellar mass. The panels from left to right display predictions at z = 12, 10, 8, and 6. The solid lines represent the mean values and the darker (lighter) shaded regions represent the 1σ (2σ) standard deviation for CDM Ref (blue), CDM HR (purple), and CDM LB (green). The symbols refer to the ALMA large observational programs ALPINE (squares; Béthermin et al. 2020; Le Fèvre et al. 2020) and REBELS (crosses; Bouwens et al. 2022; Ferrara et al. 2022), which target [C II] emission at z = 4.5–6 and z = 6.5–7.7, respectively.</figcaption></figure>

well with observations. For a more quantitative comparison, we performed a linear regression fit on the CDM HR sample across the four redshifts. From the fit, we exclude galaxies with SFR > 10⁻¹.⁵ M☉ yr⁻¹, that is, those affected by the turnover in [C II] luminosity. Our final relation is

$$ \log_{10}(L_{[C II]}/[L_{\odot}]) = 1.19 \log_{10}(SFR/[M_{\odot} \text{yr}^{-1}]) - 0.12 z + 7.58. \quad (3) $$

This is displayed as a purple dashed line in the top panels of Figure 7, and confirms that the extrapolation of our results to higher SFR is in good agreement with available high-z observations. We emphasise that, using a single linear relation to infer SFR from L_{[C II]}, as is often done in the literature, is a relatively crude approximation given the redshift dependence of the relation. As a reference, in the top row of Figure 7 we show the fit to observations of galaxies in the local Universe with SFR = (10⁻³–10²) M☉ yr⁻¹ (De Looze et al. 2014). We also show more recent fits reported in Harikane et al. (2020) and Romano et al. (2022), derived from galaxies with SFR > 10⁰.⁵ M☉. While at z = 6 the fit by De Looze et al. (2014) is very similar to ours, as we move to higher z the behaviour of COLDSIM galaxies departs from the L_{[C II]}–SFR relation found locally, as a consequence of the mentioned redshift evolution. The fit by Harikane et al. (2020) at z ≈ 7–8 is steeper than both the fit by De Looze et al. (2014) and the fit to our simulations, while the fit by Romano et al. (2022) obtained from the ALPINE sample – which also includes [C II] non-detections – is in good agreement with our results.

## 4. Discussion

In this paper, we present results from the COLDSIM suite of cosmological simulations, focusing in particular on modelling the [C II] fine-structure line emission at 158 µm. Previous works have evaluated [C II] luminosities in the context of semi-analytical models of structure formation, zoom-in simulations of galaxies, and post-processing of cosmological hydrodynamical boxes. Among the latter, COLDSIM includes a complex time-dependent non-equilibrium chemical network particularly suitable for calculating [C II] luminosity, which is expected to arise from cold (T < 10⁴ K) and dense multi-phase gas in its

atomic and molecular components. One of the strengths of the COLDSIM simulations is that the evolution of these components, as well as that of C, are computed separately for each gas particle. Thus, the [C II] emissivity can be evaluated for each particle by taking into account collisions between carbon and molecular, atomic, and ionised gas. Furthermore, the [C II] radiative power enters the balance between the thermal cooling and heating of each gas parcel at runtime, providing consistent results for gas temperature and chemical abundances.

Here we have focused on the analysis of three COLDSIM simulations, which include the same modelling and parameters regulating the various physical processes, and differ only in their box size and resolution. Indeed, the larger box size employed in CDM LB allows us to model galaxies with SFRs, stellar masses, and luminosities that are closer to the observed ones, while CDM HR is more appropriate for capturing the small-scale physics regulating the properties of the dense cold gas. We note that despite some differences observed in the evolution of the mass density estimates, the three simulations give consistent results in terms of L_{[C II]}. This is due to the fact that while the global amounts of metals in boxes with different sizes and resolutions may vary from one run to another, the [C II] luminosities converge with physical halo properties, such as stellar mass and SFR.

Our findings highlight that, while in the literature L_{[C II]} is often assumed to correlate linearly with SFR, there is an evident redshift evolution of its amplitude that is independent of resolution and box size. Here we provide a more physically motivated L_{[C II]}–SFR relation that can be employed in place of the typically used linear relation.

As the predictions discussed in this paper are sensitive to the details of the adopted feedback models (Casey et al. 2014; Narayanan & Krumholz 2017), in the bottom row of Figure 7 we compare our results for the L_{[C II]}–SFR relation with those of a number of previous theoretical investigations. We stress that the various models differ, among others, in terms of resolution, box size, feedback, subgrid physics implementation, and L_{[C II]} post-processing, meaning that a direct, quantitative comparison is not feasible. Lagache et al. (2018), for example, use a semi-analytical approach that includes a chemodynamical model from metal-free primordial accretion to account for the unresolved physics. Hydrodynamic zoom-in simulations are instead run by
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** [C II] luminosity as a function of star formation rate. The top and bottom panels show comparisons with observations and other numerical works, respectively. The panels from left to right display predictions at $z=12, 10, 8$, and $6$. The solid lines represent the mean values and the darker (lighter) shaded regions represent the $1\sigma$ ($2\sigma$) standard deviation for CDM Ref (blue), CDM HR (purple), and CDM LB (green). The top row also displays the linear fit to the CDM HR results (purple dashed line), the De Looze et al. (2014) fit made with star forming galaxies in the Local Universe (solid grey), the fit of the ALPINE data – which accounts for [C II] non-detections – from Romano et al. (2022, solid black), and the fit by Harikane et al. (2020) to $z=6–9$ ALMA observed galaxies (dotted black). The symbols refer to the ALMA large observational programs ALPINE (squares; Béthermin et al. 2020; Le Fèvre et al. 2020) and REBELS (crosses; Bouwens et al. 2022; Ferrara et al. 2022), which target [C II] emission at $z=4.5–6$ and $z=6.5–7.7$, respectively. The other data are from a compilation from Liang et al. (2024), where downward-pointing triangles refer to upper limits in $L_{[C II]}$ and left-pointing triangles to upper limits in SFR. In the bottom row, the other grey lines refer to fits of results of previous investigations, while the grey symbols to simulations of individual galaxies. Specifically, we show semi-analytical calculations by Lagache et al. (2018, dashed-dotted line), post-processed zoom-in haloes by Pallottini et al. (2019, stars), Katz et al. (2019, diamonds) and Bisbas et al. (2022, dashed line), and post-processed cosmological boxes by Leung et al. (2020, dotted line), Kannan et al. (2022, solid line), Liang et al. (2024, crosses), and Bhagwat et al. (in prep.; black dotted and grey dotted lines for the bursty and smooth stellar feedback, respectively). The comparison indicates that the slope and the amplitude of the relation vary from one model to the other, although we note that our results are in good agreement with most observational and theoretical studies reported in this plot.</figcaption></figure>

Katz et al. (2019), Bisbas et al. (2022), Pallottini et al. (2019, SERRA suite), and Liang et al. (2024, MASSIVEFIRE simulations). These simulations model the evolution of galaxies with different properties, such as stellar mass, metallicity, and resolution, spanning a range of galactic viral masses that goes from $10^{10} M_{\odot}$ in Bisbas et al. (2022) to $10^{13} M_{\odot}$ in Liang et al. (2024). The former simulation is designed to reproduce a single dwarf-galaxy merger with a gas mass particle resolution of $4 M_{\odot}$, while the latter simulates the formation and evolution of giant quiescent galaxies. Great effort was also made to extract [C II] luminosities from cosmological volumes, which cover a range in box size from ~100 cMpc/h in SIMBA (Leung et al. 2020) down to 10 cMpc/h in SPICE (Bhagwat et al. 2024), in which a resolution of ~$10^3 M_{\odot}$ for stellar particles is reached. We stress that all these efforts only follow the evolution of the global metallicity, except for Bisbas et al. (2022) who model individual elements as

in the COLDSIM approach. Despite the very different implementations employed in the various investigations, we note a general overlap and agreement in the trend between other theoretical works and COLDSIM results, except for Kannan et al. (2022), who predict systematically lower values of $L_{[C II]}$.

The simulations presented here, as all others in the literature, rely on a number of assumptions and parameters that have been extensively described and discussed in previous works (see e.g. Maio et al. 2022; Maio & Viel 2023). Specifically relevant to the calculation of $\Omega_{C^+}$ and $L_{[C II]}$ are the uncertainties related to carbon stellar yields, as the two quantities correlate with the mass of carbon injected by stellar processes. In COLDSIM, individual elements, including carbon, are tracked separately and not inferred from a global metallicity. This means that the amount of carbon of gas and star particles is estimated at each simulation time step according to the underlying
ing stellar evolution model (i.e. consistently with the adopted IMF, stellar lifetimes, mass-dependent metal yields, and feedback mechanisms, which may increase or decrease the actual chemical abundances). These features make our method robust and precise. Assumptions as to the initial stellar metallicity and the modelling of the explosive nucleosynthesis may result in differences in the carbon yields by up to a factor of two (Wheeler et al. 1995; Chieffi & Limongi 2004; Kobayashi et al. 2006). The major uncertainty is related to the fact that yield calculations are usually based on one-dimensional stellar models, which is because realistic three-dimensional treatments are much more costly in terms of computing time. The modelling of angular-momentum loss due to stellar winds and of stellar rotation provide additional sources of uncertainty in the computation of C yields, which can amount up to 1 dex for stars of masses of greater than 30 M⊙ (Limongi & Chieffi 2018).

As in our implementation we explicitly include the [C\II] 158 µm line transition within gas-cooling calculations, we are able to quantify the [C\II] luminosity of each gas particle as a function of its hydrodynamic temperature and density in a fully consistent manner. The regime in which [C\II] emission is the brightest spans the temperature range between about $10^2$ and $10^4$ K. However, uncertainties might come from the adopted collisional rates for the interactions between C$^+$ and other species (electrons, H atoms, H$_2$ molecules) and could affect the estimates of [C\II]-emitted power. Collisional rates are computed according to standard atomic physics and are often available in the literature (e.g. Hollenbach & McKee 1989; Santoro & Shull 2006; Maio et al. 2007). Rates of H-impacts (usually the most abundant at temperatures below $10^4$ K) are generally well determined. On the contrary, e$^-$-impact rates, albeit a few thousand times higher, are less robust. As the residual electron abundance in cold neutral gas is negligible, this lack of knowledge should not represent a serious problem for our calculations of [C\II] emission. H$_2$-impact collisional rates are lower than the H-impact ones by more than a factor of two at temperatures of ~$10^2$ K, while the typical reservoir of cold-gas H$_2$ is rarely larger than the neutral H abundance (Saintonge et al. 2017; Calette et al. 2018; Hunt et al. 2020; Zanella et al. 2023; Hagedorn et al. 2024). This means that our results should be unaffected by atomic-physics uncertainties (see e.g. Blum & Pradhan 1992; Wilson & Bell 2002; Goldsmith et al. 2012, and references therein).

## 5. Conclusions

[C\II] 158 µm line emission from high-z galaxies was recently measured for the first time by facilities such as ALMA (Béthermin et al. 2020; Le Fèvre et al. 2020). It has thus become important to have accurate theoretical models of cold low-temperature gas ($T < 10^4$ K) in order to understand the physical conditions that give rise to the [C\II] signal.

In this work, we studied the cold-gas properties of high-z galaxies by analysing COLDSIM, a set of state-of-the-art numerical simulations that include gravity and hydrodynamic calculations, time-dependent atomic and molecular non-equilibrium chemistry, gas cooling and heating, star formation, stellar evolution, and feedback effects in a cosmological context (see Maio et al. 2022, for details). We have thus been able to track at each time step the evolution of H\II, H\I, H$_2$, and several metal species in the regime where [C\II] fine-structure line emission takes place. Our main results can be summarised as follows:

- COLDSIM results in general reproduce the same trends observed by ALMA and JWST with respect to SFR density and SFR-$M_\star$ and $Z_{gas}-M_\star$ relations. COLDSIM SFRs of

individual haloes reach regimes that are lower than the ones observed by Bouwens et al. (2022), Curti et al. (2023), and Nakajima et al. (2023).

- For the first time, we compute the redshift evolution of the global mass densities of total carbon and C$^+$ in cosmological simulations, finding good agreement between our predictions and recent data from the XQR-30 survey (Davies et al. 2023).

- Our estimates of the dependence of [C\II] luminosities on galactic properties, such as SFR and stellar mass, are in good agreement with ALPINE and REBELS observations. It should nevertheless be noted that the range of masses and SFRs analysed and/or observed is not always the same and thus the comparison has to rely on extrapolation of the results.

- The amplitude and the scatter of the relations differ from one redshift to another, indicating a possible evolution in time of the mentioned correlations.

- We provide a fitting formula that relates $L_{[C\II]}$ to SFR and redshift, noting that using a constant relation to link $L_{[C\II]}$ to SFR at $z > 6$, as often done in observational studies, is not a good approximation, because not only does the relation evolve with redshift, but so do its amplitude and scatter.

Further studies are still required to fully capture the cold-gas physics and its dependence on modelling assumptions. Nonetheless, this work provides pivotal physical grounds for the interpretation of high-z [C\II] detection in contemporary and future observations and demonstrates a new way to theoretically investigate the cold gas in primordial times.

**Acknowledgements.** We thank the anonymous referee for the constructive report. We are thankful to Aniket Bhagwat for providing a preview of the results on [C\II] luminosity in the SPICE galaxies and for fruitful discussions on emission processes in high-z galaxies. 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 for Astrophysics. 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 tools.

## References

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

- 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

- Becker, G. D., Rauch, M., & Sargent, W. L. W. 2009, ApJ, 698, 1010

- Becker, G. D., Bolton, J. S., & Lidz, A. 2015, PASA, 32, e045

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

- Bhagwat, A., Costa, T., Ciardi, B., Pakmor, R., & Garaldi, E. 2024, MNRAS, 531, 3406

- Bhatawdekar, R., Conselice, C. J., Margalef-Bentabol, B., & Duncan, K. 2019, MNRAS, 486, 3805

- Biffi, V., & Maio, U. 2013, MNRAS, 436, 1621

- Bird, S., Rubin, K. H. R., Suresh, J., & Hernquist, L. 2016, MNRAS, 462, 307

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

- Bischetti, M., Choi, H., Fiore, F., et al. 2024, ApJ, 970, 9

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

- Blum, R. D., & Pradhan, A. K. 1992, ApJS, 80, 425

- Boksenberg, A., Sargent, W. L. W., & Rauch, M. 2003, ArXiv e-prints [arXiv:astro-ph/0307557]

- Boselli, A., Gavazzi, G., Lequeux, J., & Pierini, D. 2002, A&A, 385, 454

- Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25

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

- Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059

- Bowler, R. A. A., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 5808
Burgarella, D., Buat, V., Gruppioni, C., et al. 2013, A&A, 554, A70

Calette, A. R., Avila-Reese, V., Rodríguez-Puebla, A., Hernández-Toledo, H., & Papastergis, E. 2018, *Rev. Mex. Astron. Astrofis.*, **54**, 443

Carniani, S., Ferrara, A., Maiolino, R., et al. 2020, MNRAS, 499, 5136

Casey, C. M., Narayanan, D., & Cooray, A. 2014, *Phys. Rep.*, **541**, 45

Chieffi, A., & Limongi, M. 2004, ApJ, 608, 405

Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14

Codoreanu, A., Ryan-Weber, E. V., Crighton, N. H. M., et al. 2017, MNRAS, 472, 1023

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

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

Curti, M., D'Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425

da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, ApJ, 919, 30

Davies, R. L., Ryan-Weber, E., D'Odorico, V., et al. 2023, MNRAS, 521, 314

De Looze, I., Baes, M., Bendo, G. J., Cortese, L., & Fritz, J. 2011, MNRAS, 416, 2712

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

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

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

D'Odorico, V., Cupani, G., Cristiani, S., et al. 2013, MNRAS, 435, 1198

D'Odorico, V., Finlator, K., Cristiani, S., et al. 2022, MNRAS, 512, 2389

Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, ArXiv e-prints [arXiv:2403.03171]

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

Fahrion, K., Cormier, D., Bigiel, F., et al. 2017, A&A, 599, A9

Ferrara, A., Sommovigo, L., Dayal, P., et al. 2022, MNRAS, 512, 58

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

Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123

Fujimoto, S., Ouchi, M., Ferrara, A., et al. 2019, ApJ, 887, 107

Fujimoto, S., Oguri, M., Brammer, G., et al. 2021, ApJ, 911, 99

Fujimoto, S., Ouchi, M., Nakajima, K., et al. 2024, ApJ, 964, 146

Garcia, K., Narayanan, D., Popping, G., et al. 2023, ArXiv e-prints [arXiv:2311.01508]

Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 643, A7

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

Glover, S. C. O., & Clark, P. C. 2016, MNRAS, 456, 3596

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

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

Habing, H. J. 1968, Bull. Astron. Inst. Neth., **19**, 421

Hagedorn, B., Cicone, C., Sarzi, M., et al. 2024, A&A, 687, A244

Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93

Harikane, Y., Ono, Y., Ouchi, M., et al. 2022, ApJS, 259, 20

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56

Hernandez-Monteagudo, C., Maio, U., Ciardi, B., & Sunyaev, R. A. 2017, ArXiv e-prints [arXiv:1707.01910]

Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1

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

Hunt, L. K., Tortora, C., Ginolfi, M., & Schneider, R. 2020, A&A, 643, A180

Jones, G. C., Witstok, J., Concas, A., & Laporte, N. 2024, MNRAS, 529, L1

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

Kannan, R., Smith, A., Garaldi, E., et al. 2022, MNRAS, 514, 3857

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

Katz, H., Saxena, A., Rosdahl, J., et al. 2023, MNRAS, 518, 270

Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2021, A&A, 649, A152

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

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145

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

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

Lahén, N., Naab, T., & Szécsi, D. 2024, MNRAS, 530, 645

Laporte, N., Katz, H., Ellis, R. S., et al. 2019, MNRAS, 487, L81

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

Leech, K. J., Völk, H. J., Heinrichsen, I., et al. 1999, MNRAS, 310, 317

Leethochawalit, N., Roberts-Borsani, G., Morishita, T., Trenti, M., & Treu, T. 2023a, MNRAS, 524, 5454

Leethochawalit, N., Trenti, M., Santini, P., et al. 2023b, ApJ, 942, L26

Leung, T. K. D., Olsen, K. P., Somerville, R. S., et al. 2020, ApJ, 905, 102

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

Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13

Ma, Q., Maio, U., Ciardi, B., & Salvaterra, R. 2017a, MNRAS, 466, 1140

Ma, Q., Maio, U., Ciardi, B., & Salvaterra, R. 2017b, MNRAS, 472, 3532

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

Maio, U., & Tescari, E. 2015, MNRAS, 453, 3798

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

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

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

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

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

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

Maiolino, R., Carniani, S., Fontana, A., et al. 2015, MNRAS, 452, 54

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309

Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33

Narayanan, D., & Krumholz, M. R. 2017, MNRAS, 467, 50

Nordon, R., & Sternberg, A. 2016, MNRAS, 462, 2804

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

Padmanabhan, H. 2019, MNRAS, 488, 3014

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

Pentericci, L., Carniani, S., Castellano, M., et al. 2016, ApJ, 829, L11

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

Pettini, M., Madau, P., Bolte, M., et al. 2003, ApJ, 594, 695

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

Pizzati, E., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 495, 160

Pizzati, E., Ferrara, A., Pallottini, A., et al. 2023, MNRAS, 519, 4608

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

Posses, A., Aravena, M., González-López, J., et al. 2024, A&A, submitted [arXiv:2403.03379]

Pozzi, F., Calura, F., D'Amato, Q., et al. 2024, A&A, 686, A187

Rigopoulou, D., Hopwood, R., Magdis, G. E., et al. 2014, ApJ, 781, L15

Romano, M., Morselli, L., Cassata, P., et al. 2022, A&A, 660, A14

Romano, M., Donevski, D., Junais, A., et al. 2023, A&A, 683, L9

Rowland, L. E., Hodge, J., Bouwens, R., et al. 2024, ArXiv e-prints [arXiv:2405.06025]

Ryan-Weber, E. V., Pettini, M., Madau, P., & Zych, B. J. 2009, MNRAS, 395, 1476

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

Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529

Santoro, F., & Shull, J. M. 2006, ApJ, 643, 26

Sargent, W. L. W., Boksenberg, A., & Steidel, C. C. 1988, ApJS, 68, 539

Sargsyan, L., Lebouteiller, V., Weedman, D., et al. 2012, ApJ, 755, 171

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

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

Schouws, S., Stefanon, M., Bouwens, R., et al. 2022, ApJ, 928, 31

Sebastian, A. M., Ryan-Weber, E., Davies, R. L., et al. 2024, MNRAS, 530, 1829

Simcoe, R. A., Sargent, W. L. W., Rauch, M., & Becker, G. 2006, ApJ, 637, 648

Simcoe, R. A., Cooksey, K. L., Matejek, M., et al. 2011, ApJ, 743, 21

Simcoe, R. A., Onoue, M., Eilers, A. C., et al. 2020, ArXiv e-prints [arXiv:2011.10582]

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

Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122

Songaila, A. 2001, ApJ, 561, L153

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. A., et al. 2022, MNRAS, 512, 4736

Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170

Tescari, E., Viel, M., D'Odorico, V., et al. 2011, MNRAS, 411, 826

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

Velusamy, T., & Langer, W. D. 2014, A&A, 572, A45

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

Wheeler, J. C., Harkness, R. P., Khokhlov, A. M., & Hoeflich, P. 1995, Phys. Rep., 256, 211

Whitaker, K. E., Pope, A., Cybulski, R., et al. 2017, ApJ, 850, 208

Wilson, N. J., & Bell, K. L. 2002, MNRAS, 337, 1027

Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976

Zanella, A., Valentino, F., Gallazzi, A., et al. 2023, MNRAS, 524, 923
## Appendix A: [C II] Cooling

In this Appendix, we provide the atomic data adopted in Section 2.2 to compute $\Lambda_{[C II]}$. We model [C II] as a two-level system considering the fine structure transition $(2p)[^2P_{3/2} -^2P_{1/2}]$ between the quantum number $J = 3/2$ and $J = 1/2$ states. The following data are taken from Hollenbach & McKee (1989) and Goldsmith et al. (2012) and refer to e⁻-impact, H-impact and H₂-impact collisional rates, spontaneous transition rate, and energy level separation:

$$ \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.1 \cdot 10^{-4} \text{ s}^{-1}, $$

$$ \Delta E_{21} = 5.71 \cdot 10^{-14} \text{ erg.} $$

where $T_{100} = T/(100 \text{ K})$ and the validity regime is for gas with temperature $T < 10^4 \text{ K}$.