© The Authors 2024

1. Introduction

The gas bound to a galaxy dark-matter halo, known as the circumgalactic medium (CGM), plays an essential role in understanding galaxy formation and evolution. It is a reservoir of gas that interacts with the galactic environment, regulating key processes such as star formation, accretion, and feedback. Investigating the CGM provides crucial insights into the complex interplay between galaxies and their surrounding gas reservoirs (e.g., Tumlinson et al. 2017). Studying it is challenging as it is predominantly composed of diffuse and tenuous gas, making direct observations difficult. Recent advancements in observational capabilities – for example, the Atacama Large Millimeter/Submillimeter Array (ALMA; Wootten & Thompson 2009)–, combined with state-of-the-art simulations (e.g., IllustrisTNG, Nelson et al. 2019a, EAGLE, Crain et al. 2015; Schaye et al. 2015), have opened new avenues for exploring the CGM and unveiling its properties. However, the limited resolution of the simulations and the requirement of deep observations are obstacles that remain.

The [CII] 158 μm fine-structure line has emerged as a valuable tracer of the CGM. This line is one of the brightest emission lines in the infrared spectra of star-forming galaxies (e.g., Stacey et al. 1991; Brauher et al. 2008) and is an important tracer of cold gas in both the interstellar medium (ISM) and CGM. It is emitted from the upper fine-structure level J = 3/2 to the lower level J = 1/2 of the C⁺ ion. The higher energy level exhibits an equivalent temperature of 91 K, while its critical densities are 16, 2.4 × 10³, and 4.8 × 10³ cm³ for collisions with electrons, hydrogen atoms, and H₂ molecules, respectively (Goldsmith et al. 2012). Since the ionization potential of atomic carbon is only 11.3 eV (versus 13.6 eV of neutral hydrogen), [CII] can originate from atomic ISM, diffuse molecular, and ionized gas (e.g., Abel 2006; Zanella et al. 2018; Heintz et al. 2021). At high redshifts, both simulations (e.g., Olsen et al. 2015; Vallini et al. 2015) and observations (e.g., Stacey et al. 2010; Gullberg et al. 2015) have found that the [CII] line primarily emanates from photodissociation regions (PDRs).

Observations of [CII] emission in extragalactic sources notably increased with the advent of the Herschel Space Observatory (Pilbratt et al. 2010). Studies have focused on local starburst, luminous infrared galaxies, and local dwarf galaxies, investigating the properties of the ISM and its link with the [CII] properties (e.g., Sargsyan et al. 2012; Díaz-Santos et al. 2014, 2017; Herrera-Camus et al. 2015; Samsonyan 2022; Romano et al. 2024). In particular, there has been progress in investigating [CII] as a star formation rate (SFR) tracer. While different works found a tight correlation between the [CII] luminosity and SFR (e.g., De Looze et al. 2014; Herrera-Camus et al. 2015), other Herschel observations showed that the L_[CII]/L_IR (or equivalently, L_[CII]/SFR) ratio decreases continuously with increasing L_IR (at L_IR ≳ 10¹¹ L_⊙, e.g., Sargsyan et al. 2012; Díaz-Santos et al. 2017; Contursi et al. 2017). This could be an indication that the L_[CII]/SFR ratio is influenced by other galaxy properties, such as metallicity and dust temperature (e.g., Bolatto et al. 2013; Díaz-Santos et al. 2014; Bisbas et al. 2023).

While locally the [CII] line can only be observed from space due to atmospheric limitations, at high-z the frequency of the [CII] line is accessible with ground-based observatories such as ALMA. In particular, the Earth’s atmosphere allows for optimal access to the [CII] transition in the frequency range 345.5−211.2 GHz, which corresponds to the redshift range 4.5 < z < 8.5. Therefore, [CII] has become a key target for observational studies of high-redshift galaxies, where it can provide insights into the early stages of galaxy formation and evolution. Extended [CII] emission (known as [CII] halos) has now been regularly observed in several high-redshift objects such as dusty starburst and normal star-forming galaxies (e.g., Fujimoto et al. 2019, 2020; Ginolfi et al. 2020a; Carniani et al. 2020; Rybak et al. 2020; Herrera-Camus et al. 2015, 2021; Meyer et al. 2022; Akins et al. 2022; Solimano et al. 2024). These halos can extend up to tens of kiloparsecs from the galaxies and their physical properties and origins are still poorly understood. As a reference, the typical effective radius of star-forming galaxies at 4 < z < 6 is expected to be ≲1.5 kpc (e.g., Shibuya et al. 2015). In the context of star-forming galaxies, Fujimoto et al. (2019) reported the first identification of [CII] halos (∼10 kpc scales) around z ∼ 5 − 7 star-forming galaxies using deep ALMA data, through uv-visibility plane stacking. Subsequently, with the development of the ALMA Large Program to Investigate C⁺ at Early Times (ALPINE) survey (Le Fèvre et al. 2020), further observations of [CII] halos around high-z star-forming galaxies have been reported.

Ginolfi et al. (2020a) (hereafter G20) studied a sample of 50 [CII] emitting galaxies at z ∼ 4 − 6 from ALPINE. By stacking the [CII] spectra and cubes, they find outflow signatures together with extended [CII] emission (∼15 kpc) in their galaxy sample with the higher SFR (SFR > 25 M_⊙ year⁻¹). In particular, they obtained an average [CII] radial surface-brigthness profile for the high-SFR galaxies, which is characterized by a compact inner component and an extended component (radial scales > 10 kpc). Later, Fujimoto et al. (2020) performed individual measurements for the ALPINE galaxies, finding that ∼30% of them present [CII] halos (∼10 kpc scales). Recently, an additional ALMA effort, the ALMA Reionization Era Bright Emission Line Survey (REBELS) (Bouwens et al. 2022), has been conducted. This study consists of ALMA [CII] observations of a sample of 40 star-forming galaxies at z ∼ 6.5 − 9. Based on this survey, Fudamoto et al. (2022) stacked 28 [CII] emitting galaxies at z ∼ 7. Their results indicate that the average [CII] surface-brightness profile is well fit with a single component (i.e., no extended component), with an effective radius of ∼2 kpc. Possible reasons for this absence of an extended component could be related to the fact that they stacked moment maps rather than the data cubes, and that they included both high- and low-SFR galaxies in the analysis. Despite the aforementioned result, this [CII] emission is spatially more extended compared to the dust continuum and the rest-frame UV emission, in agreement with previous works (e.g., G20; Fujimoto et al. 2020).

To fully understand these [CII] halos, it is essential to perform theoretical investigations. Currently, cosmological simulations are not able to reproduce extended [CII] emission at high-z (e.g., Pallottini et al. 2017, 2019). For instance, Fujimoto et al. (2019) compared their observational results with two sets of numerical simulations of massive galaxies (M_halo ∼ 10¹¹ − 10¹² M_⊙) at z ∼ 6 (Pallottini et al. 2017; Arata et al. 2019, m_baryon = 1.2 × 10⁴ and 1.8 × 10⁵ M_⊙ h⁻¹, respectively). By using mock data from both simulations, they find that the [CII] emission is not as extended as the observational data. Specifically, they obtained surface-brightness profiles that decline rapidly beyond r ∼ 3 kpc, which is significantly smaller than the reported [CII] extension of 10 kpc observed by Fujimoto et al. (2019). These profiles are at least one order of magnitude lower than the observed flux at inner-CGM scales. This disagreement suggests that certain crucial physical processes for the production of [CII], such as feedback or metal enrichment, are not well reproduced by the simulations. It is also possible that there are limitations or uncertainties in the modeling of [CII] emissivities. This discrepancy between simulations and observations highlights the need for additional investigations and theoretical refinement.

The origin of this extended [CII] emission remains unclear. Different scenarios have been proposed, including the presence of satellite galaxies, extended PDRs, extended HII regions, inflows of cold streams, outflows produced by stellar feedback, or active galactic nuclei (AGN) feedback (see Fujimoto et al. 2019 for a detailed description). Currently, the scenario in which star formation driven outflows (from the central galaxy) are responsible for the extended [CII] emission seems to be the most plausible according to observations (e.g., G20; Fujimoto et al. 2020; Romano et al. 2024). The observed [CII] halos are indeed only evident around more star-forming galaxies (SFR > 25 M_⊙ year⁻¹) and the line-emission stacking of these sources seems to highlight the presence of broad wings (G20). This scenario has also been supported by semi-analytical models (e.g., Pizzati et al. 2020, 2023). These models primarily focus on simulating the cooling outflows and their resulting [CII] emission, suggesting that the observed [CII] halos could be attributed to gas that is cooling after being heated and expelled in supernova-driven outflows. The effectiveness of these models is largely attributed to their ability to accurately capture the catastrophic cooling process within central kiloparsec scales. Specifically, they find that supernova-driven outflows expand into the CGM at velocities of 200−500 km s⁻¹. The outflowing gas undergoes rapid cooling (approximately a few hundred Kelvin) within the first kiloparsec, then it is slowly heated (T ∼ 10³ K) by the cosmic UV background. These conditions favor the formation and survival of C⁺ ions.

In this framework, we post-processed a sample of TNG50 halos hosting z ∼ 4 − 6 star-forming galaxies to model their [CII] emission. In particular, we produced mock observations and compared them with the ALPINE observations (G20), ultimately discussing the possible origin for the observed extended [CII] emission in the inner CGM of high-z galaxies. This paper is structured as follows. In Sect. 2 we present a summary of the observations of G20. In Sect. 3 we describe the methods used in this work. In Sect. 4 we present our main results. In Sect. 5 we discuss our main results and explore the different factors involved in the production of [CII] halos and caveats. Finally, Sect. 6 summarizes our findings.

2. Observations

The ALPINE survey is an ALMA large program that targets [CII] and FIR-continuum emission of a sample of 118 normal (i.e., main sequence) star-forming galaxies. These galaxies are at 4.4 < z < 5.8, and have M_* ∼ 10^8.5 − 10¹¹ M_⊙ and SFR ∼ 3−650 M_⊙ year⁻¹. The average resolution of these observations is 0.9 arcsec, that is, ranging from ∼5.2 to 6 kpc for the different redshifts considered. This coarse resolution implies that the galaxies are not resolved in these ALMA observations. The spectral resolution is Δν_channel = 25 − 35 km s⁻¹, enabling the resolution of the typical [CII] emission line width (FWHM ∼ 235 km s⁻¹). All the details regarding the observational setup and data reduction are described in Béthermin et al. (2020). In total, 75 out of 118 galaxies were robustly detected in [CII], with signal-to-noise ratio (S/N) > 3.5. G20 studied these [CII]-emitting galaxies, excluding ∼30% of those objects that present signs of ongoing mergers (major or minor). Thus, their final sample¹ consists of 50 normal star-forming galaxies at 4.4 < z < 5.8, with M_* ∼ 10^8.5 − 10¹¹ M_⊙ and SFR ∼ 5−650 M_⊙ year⁻¹.

In order to investigate the efficiency of galactic feedback in these galaxies and the spatial distribution of the [CII] emission, G20 investigated the [CII] emission data through a stacking analysis. For the purposes of this paper, in the following paragraphs we summarize only their stacking of the cubes. This procedure consists of a combination of the [CII] cubes of each galaxy C_i with i slices (previously aligned), based on a vector variance stacking technique, where the stacked cube C${}_i^{\mathrm{stack}}$ is estimated as:

$$\begin{array}{c}\hfill C_i^{\mathrm{stack}}=\frac{{\sum }_{k=1}^N{C}_{i,k}·w_{i,k}}{{\sum }_{k=1}^N{w}_{i,k}}\end{array}$$

where C_i, k corresponds to the cube of the kth galaxy. The weigthing factor w_i, k is the inverse variance, equal to 1/${\sigma }_{i,k}^2$, where σ_i, k corresponds to the spatial rms associated with each slice of each galaxy. In this way, noisier cubes are down-weighted.

This stacking was performed on two subsamples separated by a threshold in SFR: high-SFR (SFR > 25 M_⊙ yr⁻¹) galaxies and low-SFR (SFR ≤ 25 M_⊙ yr⁻¹) galaxies. The SFRs were computed through spectral energy distribution (SED) fitting, assuming a constant star-formation history over 100 Myr (Faisst et al. 2020). In order to obtain velocity-integrated flux maps of the core of the [CII] emission, the spectral slices of the stacked cubes were collapsed in the velocity range of [−200 : +200] km s⁻¹. The final stacked synthetized beam (or point spread function, PSF) was estimated by stacking the PSF of each input galaxy cube. Then, the channels of this stacked PSF cube are collapsed over the same velocity range mentioned above. This final PSF has a major axis full-width-half-maximum (FWHM) of 0.98 arcsec, a minor axis FWHM of 0.89 arcsec, and a position angle of −30 deg. From the resulting [CII] flux maps of the subsamples, G20 finally computed circularly averaged radial surface-brigthness (SB) profiles, considering radial bins of 0.15 arcsec. For most of our comparative analysis, we present a single observed SB profile, derived by averaging the two SB profiles obtained from G20 for the high and low-SFR galaxy subsamples.

3. Methods

3.1. IllustrisTNG simulations

In this work, we explore different high-z halos from the IllustrisTNG cosmological magnetohydrodynamical simulations (Nelson et al. 2019a), in particular the TNG50 simulation (Nelson et al. 2019b; Pillepich et al. 2019). Using the moving-mesh code AREPO (Springel 2010), these simulations follow the coupled dynamics of the dark-matter and gas. This method is based on a spatial discretization consisting of a Voronoi tessellation of the simulation box. The physical processes of baryonic TNG runs, among which are star-formation, stellar evolution and feedback, supermassive black hole formation and evolution, follow the “TNG galaxy formation model” (methods described in Weinberger et al. 2017; Pillepich et al. 2018). We briefly describe here only some aspects of these simulations relevant for this work.

The simulations include gas radiative processes, including metal-line cooling and heating from the ultraviolet background (UVB). Due to limitations in the resolution, the gas does not cool radiatively below 10⁴ K and the processes of star formation and pressurization in the multiphase ISM are modeled based on the approach of Springel & Hernquist (2003). In this model, gas with densities above a threshold of ≃0.1 cm⁻³ undergoes stochastic star formation according to the Kennicutt–Schmidt relation, assuming a Chabrier initial mass function (Chabrier 2003). To account for the effects of unresolved supernovae, a two-phase effective equation of state model is employed to incorporate the pressurization of star-forming gas. Chemical enrichment is driven from supernovae type II, Ia, and AGB stars. The simulations also include local radiation effects from AGN, and AGN feedback in two modes: a thermal “quasar” state at high accretion rates, and a kinetic “wind” state at low accretion rates (Weinberger et al. 2017).

This TNG model has been calibrated in order to reproduce several of the observed galaxy properties mainly at z = 0. The cosmology assumed has parameters Ω_Λ, 0 = 0.6911, Ω_m, 0 = 0.3089, Ω_b, 0 = 0.0486, σ₈ = 0.8159, n_s = 0.9667 and h = 0.6774 (Planck Collaboration XIII 2016).

Specifically, in this work we used the TNG50-1 box, which presents the smallest volume (51.7³ cMpc³) but the highest resolution (m_baryon = 5.6 × 10⁴ M_⊙ h⁻¹, galaxies with M_* ≳ 10⁷ M_⊙) of the TNG runs. The latter is a key factor to perform our [CII] emission predictions, especially considering that high densities are needed for [CII] production (see critical densities in Sect. 1). Compared to previous works, the TNG50 simulations offer key advantages. They strike a balance between computational efficiency and resolution, allowing for a larger sample size and a more comprehensive representation of the high-z halo population. Additionally, these simulations incorporate a self-consistent treatment of various physical processes, which are uniformly implemented across the simulation volume. This ensures a consistent modeling approach for the entire sample of halos.

3.2. Selection of simulated galaxies

We selected a sample of TNG50-1 star-forming galaxies² with similar physical properties to the galaxies of the full ALPINE survey (Le Fèvre et al. 2020). Specifically, they were chosen to match the ranges in stellar mass (log(M_*) ∼ 8.5–11.0 M_⊙) and star-formation rate (SFR ∼ 3–623 M_⊙ year⁻¹)³. In total, we found 357 TNG50 galaxies matching these properties. We then optimized our selection by constraining our distribution of SFR and M_* to conform to those exhibited by the ALPINE [CII] emitting galaxies (G20). Throughout this work, we adopt for the simulated galaxies: (i) stellar masses within twice the stellar half mass radius, and (ii) SFRs averaged across the last 100 Myr within twice the stellar half mass radius (see Donnari et al. 2019; Pillepich et al. 2019), comparable to the observationally derived SFRs (see Sect. 2). These galaxies were also selected to have similar redshifts to those of the observations (z ∼ 4.4 − 5.7). However, it was not possible to precisely match them because snapshots with intermediate redshifts (the so-called mini snapshots) do not contain information on the metal abundances which are needed to estimate the [CII] emission. Thus, the final selected halos were extracted from the three full snapshots (13, 17 and 21) corresponding to the redshifts 6, 5 and 4, respectively.

The individual central galaxies were cut from the full snapshot (at ∼3 R₂₀₀) in order to facilitate their analysis. In total, we selected 75 galaxies that match the distribution of SFR and M_* of the sample in G20. Analogous to the observations, we excluded from our analysis those systems that would have exhibited observable merging signatures. Specifically, following the observational limitations, we did not consider systems with companion galaxies with SFR > 1.5 M_⊙ year⁻¹, within a 3D radial distance of 1.1 R₂₀₀. We also note that our selected galaxies contain AGNs, while the observed sample excludes type I AGNs. However, these AGNs of our simulated sample have low bolometric luminosities⁴ (L_BOL ≲ 10^45.5 erg s⁻¹), which would not be detected in currently available observations. Figure 1 shows the distribution of the selected simulated galaxies (red) in the SFR and M_* plane, in comparison to the sample of G20 (blue), using density countours. Importantly, this selection represents one possible realization of the sample in G20 extracted from the TNG50 galaxies in the mass and SFR ranges of the ALPINE survey. This means that there could be other selections of simulated galaxies that also fit the same criteria of the sample in G20, but include slightly different individual galaxies. In this figure, the spacing between the blue contours (observations) appears larger than the spacing between the red contours (simulations). This suggests that the observed sample has a more spread-out distribution of galaxies in the joint space of M_* and SFR compared to the simulated sample. Areas where the blue and red contours overlap suggest regions of similar M_* and SFR values between both samples.

Fig. 1. Comparison of the distribution of the SFRs and M* values between the simulated (red) and observed (blue) samples (75 and 50 galaxies, respectively). The joint distribution is visualized using density contours. The contour spacing represents the density of data points, with closer contours indicating regions of higher data density. The blue bars in the bottom-right corner represent the average 1-sigma error for the SFR and M* of the observations. The histograms along the x and y axes provide marginal distributions of the M* and SFR, respectively, for both observed and simulated samples.

3.3. [CII] cooling rate and flux estimation

In order to obtain an estimate of the [CII] emission for each selected halo, we used the analytical [CII] cooling rate (Λ_CII) from Bisbas et al. (2022). In their study, Bisbas et al. (2022) performed [CII] synthethic observations using the GRIFFIN smooth particle hydrodynamics simulations (Lahén et al. 2020). These simulations have a resolution of 4 M_⊙ per gas particle and use a modified version of GADGET-3 (Springel 2005) described in Hu et al. (2014). They assume a nonequilibrium model of cooling and chemistry that traces H₂, H⁺, CO, H, C⁺, O, and free electrons (Nelson & Langer 1997; Glover & Mac Low 2007), only when the gas temperature is < 3 × 10⁴ K. For temperatures above this threshold, equilibrium cooling tables based on 12 different metal species (H, He, N, C, O, Si, Mg, Fe, S, Ca, Ne, and Zn) are used (from Wiersma et al. 2009). Furthermore, the simulations also include routines to account for the interstellar UV radiation field and mechanisms for stellar feedback, such as photoionization, photo-electric heating, and supernovae.

Based on the direct Λ_CII estimations from the simulations, Bisbas et al. (2022) proposed a theoretical approach to get Λ_CII, based on analytical expectations. They demonstrated that the Λ_CII values obtained from the analytical expressions are in agreement with the majority of the Λ_CII values obtained directly from their simulations. This approach assumes that the gas is optically thin and is applicable for T > 150 K. The analytical cooling rate can be expressed as:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{CII}}=\frac{R_{\mathrm{c},\mathrm{dex}}\times 2e^{-91.25/\mathrm{T}}}{R_{\mathrm{s}}+R_{\mathrm{c},\mathrm{dex}}(1+2e^{-91.25/\mathrm{T}})}{R}_{\mathrm{s}}{E}_{\mathrm{CII}}{n}_{{\mathrm{C}}^{+}}(\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}),\end{array}$$(1)

where T is the temperature of each gas cell and the excitation temperature of CII transition is 91.25 K, E_CII is the energy hν_CII = 1.25988 × 10⁻¹⁴ erg, n_C⁺ is the number density of C⁺ particles, R_s is the spontaneous emission equal to the corresponding Einstein A-coefficient R_s = 2.291 × 10⁻⁶ s⁻¹, and R_c, dex is the total de-excitation rate given by:

$$\begin{array}{c}\hfill R_{\mathrm{c},\mathrm{dex}}=R_{\mathrm{c},10}(\mathrm{H})n_{\mathrm{H}}+R_{\mathrm{c},10}({\mathrm{H}}_2)n_{{\mathrm{H}}_2}+R_{\mathrm{c},10}(\mathrm{e}-)n_{\mathrm{e}},\end{array}$$(2)

where R_c, 10(e−), R_c, 10(H₂) and R₁₀(H) are collisional de-excitation rates with e−, H, and H₂ as colliding partners, respectively. We consider the contribution of collisions with H₂ negligible. Indeed, previous works argue that the contribution of molecular gas to the total [CII] luminosity is minor (e.g., Franeck et al. 2018; Tarantino et al. 2021; Bisbas et al. 2022). Therefore, we focus on [CII] emission originating from PDRs where sufficient amounts of HI and e− are present as collision partners (e.g., Olsen et al. 2015). This assumption is also based on the fact that TNG50 does not follow the formation and evolution of molecular gas, consistent with the simulation setup.

Then, R_c, 10(e−) and R_c, 10(H) are (from Leiden Atomic and Molecular Database, Schöier et al. 2005):

$$\begin{array}{cc}\hfill R_{\mathrm{c},10}(\mathrm{e}-)& =2.426206\times {10}^{-7}{\left(\frac{T}{100}\right)}^{-0.345}({\mathrm{cm}}^{-3}{\mathrm{s}}^{-1})\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill R_{\mathrm{c},10}(\mathrm{H})& =3.113619\times {10}^{-10}{\left(\frac{T}{100}\right)}^{0.385}({\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}).\hfill \end{array}$$(4)

To estimate Λ_CII (Eq. 1), all required parameters except n_C⁺ can be estimated or assumed from the outputs of TNG50-1 (i.e., temperatures and densities).

To gain insight into the physical properties of the selected galaxies, Fig. 2 presents a temperature-density (T − n_H) phase diagram within R₂₀₀, of one z = 5 system. The color bar indicates the gas mass distribution, where the yellow (blue) regions correspond to areas of high (low) gas mass concentration within the halo. The region where most of the mass is concentrated (∼T < 10^4.5 K, 10^−2.5 < n_H < 10⁻¹ cm⁻³) is dominated by ISM gas. It is important to note that all dense gas cells with n_H > 0.1 cm⁻³ are by construction star-forming in these simulations. The temperatures of this star-forming gas is determined by the effective equation of state (Springel & Hernquist 2003) and they do not represent the physical temperatures of the gas (Pillepich et al. 2019). Instead, the star-forming gas cells (n_H > 0.13 cm⁻³) are assumed to be at their cold-phase temperature of 1000 K (Springel & Hernquist 2003), as the cold-phase is dominant in mass with respect to the warmer component in this regime. This approach is common in modeling emission lines arising from the cold-phase in TNG simulations (e.g., Nelson et al. 2021).

Fig. 2. Phase diagram for one of the selected TNG50 galaxies at z = 5, weighted and color-coded by the gas mass. The diagram encompasses all gas cells within R200. For densities above the star-formation threshold nH = 0.1, the temperatures are determined by the effective equation of state (Springel & Hernquist 2003).

In order to estimate the values of n_C⁺ we need the C⁺ abundances for each gas cell. As an initial simplified approach, we assume a scenario where all (i.e., the total) neutral carbon is converted into C⁺. Subsequently, we explore other cases in which the C⁺ abundances are obtained using the spectral synthesis code CLOUDY version 17.03 (Ferland et al. 2017). For the CLOUDY calculations, we either assumed the redshift dependent UVB of Khaire & Srianand (2019) as the only photoionization source, or the combination of UVB and radiation from nearby young stars (ages < 10 Myr). For the young stars we used the spectra model of Cerviño et al. (2002), and assumed that only a fixed 5% of the stellar radiation can escape the molecular birth cloud (i.e., escape fraction f_esc = 0.05). The normalization of the stellar radiation field (ϕ) in units of M_⊙ year⁻¹ kpc⁻² at each gas particle position ${\overrightarrow{r}}_j$ was computed in the optical thin approximation (no absorbtion along the line-of-sight beyond the constant at-the-source f_esc):

$$\begin{array}{c}\hfill \varphi ({\overrightarrow{r}}_j)=\frac{f_{\mathrm{esc}}}{{10}^7\mathrm{yr}}{\displaystyle \sum _i}\frac{M_i}{{({\overrightarrow{r}}_i-{\overrightarrow{r}}_j)}^2},\end{array}$$(5)

where M_i and ${\overrightarrow{r}}_i$ are the mass and position of stellar particles with ages less than 10 Myr, following Obreja et al. (2019). To compute the ionization fraction χ = n_C⁺/n_C, we thus run one zone plane-parallel CLOUDY models for fixed hydrogen number density log(n_H), temperature log(T), metallicity log(Z/Z_⊙), normalization of the young star radiation field log(ϕ) and redshift (z = 3, 4, 5 for the UVB), to cover all the physical conditions of gas cells in the TNG50 galaxies. No dust and hence no metal depletion is assumed in accordance with the TNG50 simulation setup. Subsequently, we used 3D or 4D barycentric interpolation in either χ[log(Z), log(n_H), log(T)] or χ[log(ϕ), log(Z), log(n_H), log(T)] CLOUDY tables, to get the C⁺ ionization fraction of each gas cell for the UVB and UVB plus young star photoionization models at fixed redshift. Figure 3 shows examples of input spectra used in the CLOUDY calculations, and Appendix B gives an example of a CLOUDY input file.

Fig. 3. Example of different input spectra used in the CLOUDY calculations. The solid light-gray line represents the UVB SED of Khaire & Srianand (2019), and the dashed gray line the CMB, both at z = 5. The vertical black lines represent the ionization energies of C→C+ (dotted) and C+→C++ (dashed). The colored lines represent the contribution to the SED due to young stars (Cerviño et al. 2002) at different normalizations, plus the contribution of the UVB. The units of the radiation field ϕ are M⊙ year−1 kpc−2.

Once the C⁺ abundances are computed, all the quantities in Eq. (1) are available. We can then obtain the [CII] luminosities per cell as L_[CII] = Λ_CII × V, with V the volume of each cell, and estimate their respective flux densities following Solomon et al. (1992):

$$\begin{array}{c}\hfill I_{[\mathrm{CII}]}=\frac{L_{[\mathrm{CII}]}}{1.04\times {10}^{-3}{\nu }_{\mathrm{obs}}{D}_{\mathrm{L}}^2}(L_{\odot })(\mathrm{Jy}\mathrm{km}{\mathrm{s}}^{-1}),\end{array}$$(6)

where L_[CII] is the [CII] line luminosity in L_⊙ units, ν_obs is the observed frequency in units of GHz, and D_L corresponds to the luminosity distance in Mpc.

3.4. [CII] 2D flux maps and surface brightness profiles

After getting L_[CII] for each gas cell, we constructed [CII] surface brightness maps. For this, we obtained 2D projections by integrating the emission within a slice encompassing a range of [−R₂₀₀, +R₂₀₀] (or ∼[−40, +40] kpc at the median z = 5) from the galaxy’s center⁵. As a conservative approach, we did not apply any cut in velocities, ensuring the inclusion of the central galaxy in our analysis. Beyond the chosen range, structures at larger distances are not expected to greatly contribute to the [CII] signal, which is dominated by the central galaxy. In contrast, the observed maps apply a velocity cut of [−200, +200] km s⁻¹. For reference, these velocities are similar to the virial velocity for the average halo hosting the selected TNG50 galaxies (v_vir ∼ 220 km s⁻¹).

In order to make these 2D maps comparable with the observations, we applied a 2D Gaussian smoothing with a FWHM equal to the stacked synthetized beam obtained in the ALMA observations (0.9 arcsec or ∼6 kpc at z = 5; see Sect. 2). Figure 4 shows an example of [CII] 2D surface brightness maps for one of the selected galaxies at z = 5, at full resolution⁶ (left panel) and with smoothing (right panel). We have not added noise to our modeled smoothed SB maps because our primary focus is to assess whether we can achieve similar flux levels as observed. Future simulations of ALMA observations with SIMALMA⁷ are needed to incorporate noise levels and achieve a more realistic representation.

For each galaxy, we then extract circularly averaged radial profiles from their 2D smoothed SB maps. The profiles were obtained by considering radial bins of 0.15 arcsec, as done in the observations. Subsequently, total averaged profiles were obtained by stacking the individual profiles. To accomplish this, we utilized a total of 20 orientations per galaxy and performed 1000 iterations of the experiment, randomly selecting one orientation per galaxy and calculating the final average based on these realizations.

4. Results

In this section we present the results obtained from our post-processing model, which generates [CII] emission in high-redshift galaxies using the TNG50 simulations. By comparing our model predictions with observational data from the ALPINE survey, we assess the validity and performance of our approach, and try to provide some intuition on the aforementioned [CII] emission scenarios.

Fig. 4. [CII] emission surface-brightness map for one of the selected TNG50 galaxies at z = 5, assuming C+ abundances with UVB and local young stars as ionization sources. Left: the map is presented at the full-resolution of the simulation. Right: the map has been convolved with the average ALMA beam of the G20 observations (bottom left corner). The white contour represents the average 2σ level of the observations (0.08 Jy km s−1 beam−1).

4.1. Relation between [CII] luminosity and star formation rate

Several studies have found that the [CII] luminosity is tightly correlated with the SFR, with a nearly linear relation (e.g., De Looze et al. 2014; Olsen et al. 2017; Lagache et al. 2018; Schaerer et al. 2020). However, at high-z the scatter around the local relation tends to increase (e.g., Casavecchia et al. 2024). This fact has been primarily attributed to variations in galaxy ISM physical properties, such as gas density and metallicity (e.g., Lagache et al. 2018). In this context, we explore the [CII]–SFR relation obtained from our post-processing model, and compare the obtained results with the observational data from the ALPINE survey. Therefore, the subsequent analysis checks the consistency between our model predictions and the observations.

The results are shown in Fig. 5, where both the simulated (left panel) and observational (right panel) data are color-coded by stellar mass. For the simulated galaxies, the x-axis shows the SFR averaged over the last 100 Myr, and the y-axis the [CII] luminosity obtained by assuming UVB and young stars as photoionization fields, and by integrating inside the area within the 2σ level of the observations (stacked flux map, r ∼ 1.3 arcsec). Comparing the two panels, it is evident that our results are in agreement with the observations, when considering the simulated galaxies that can be detected by the observations (L_[CII] > 10^7.8 L_⊙). In addition, for both sets of data we can appreciate a trend in the stellar masses along the relation. However, when comparing the distributions of normalized distances⁸ between M_*, L_[CII] and SFR through a Kolmogorov–Smirnov test, it is found that the two samples are not drawn from the same distribution (p-value = 0.002). A possible explanation for this disparity is the significant scatter in M_* and SFR values in the observed data compared to TNG50 (i.e., Fig. 1). Our results are also roughly in agreement with the expected relationship for local starburst galaxies (blue line, De Looze et al. 2014) and for high-z galaxies (red dashed line, Lagache et al. 2018). We note that three of our simulated galaxies lie below the expected relationships, exhibiting luminosities below the observational limit. Since these galaxies represent non-detections in the observations of G20, we exclude them from any future analysis.

Fig. 5. [CII] luminosity versus SFR for the TNG50 galaxies (left panel) and the ALPINE survey (right panel). The shaded blue area represents the 1σ relationship of De Looze et al. (2014) for starburst local galaxies. The dashed red line shows the relationship of Lagache et al. (2018) for z = 5 star-forming galaxies. The average 1σ error for L[CII] and SFR in the right panel are 0.05 dex and 0.4 dex, respectively.

Overall, the results of this analysis provide validation for the predictions made by our theoretical model on galactic scales. In the next sections we explore what these models predict for the [CII] emission in the CGM.

4.2. [CII] Surface brightness profiles

4.2.1. Averaged profiles of all gas

The study of radial [CII] surface brightness profiles in high-redshift galaxies enables a deeper understanding of the spatial distribution of [CII] emission, the properties of the cold ISM, the interplay between star formation and feedback, and the properties of the CGM (e.g., Fujimoto et al. 2019, 2020; G20). In Fig. 6 we present the resulting [CII] surface brigthness profiles of our simulated galaxies, for the three different cases of C⁺ abundances calculated (represented by different colors; see the legend), as explained in Sect. 3.3. The black squares show the averaged profile from G20, the solid lines show the average (including z = 4, 5 and 6) profiles from TNG50, and the shaded areas encompass 95% of the distribution of the averaged profiles through multiple realizations of the experiment, accounting for different orientations of the galaxies (see Sect. 3.4). The radial distances are shown up to ∼3 arcsec (or ∼19 kpc at z = 5), which corresponds to the maximum radius analyzed in the observational profiles.

Fig. 6. Averaged [CII] SB profiles obtained for our sample of 72 TNG50 galaxies at z = 4, 5, and 6 (solid lines), with shaded areas enclosing a 95% of different results obtained by different galaxy orientations. The different line colors represent the TNG50 profiles obtained assuming C+ abundances equal to the total C abundance (dark pink), C+ abundance obtained considering the UVB as incident field (light blue) and C+ abundance considering UVB and local young stars as ionization sources (orange). As a comparison, the black squares represent the results from the ALPINE survey.

For all cases, the most central part of the profile is dominated by the PSF of the observations (up to a radius of ∼1.5 arcsec), while the flatter part at larger scales correspond to the inner CGM of these galaxies. In the scenario where all carbon (C) is assumed to be converted into [CII] (dark pink line), comparing with the observations we get an overestimation of the emission at central scales (∼8 times higher) and an underestimation at scales larger than ∼1.5 arcsec (∼6 times lower). When considering UVB as ionizing source (light blue line), we can see that our predicted [CII] flux is also higher (∼7 times) than the observations at central scales, and lower (∼7 times) at large scales. In the last case, which assumes UVB and young stars as ionizing sources (orange line), the simulated flux is comparable to the observed one at central scales (in agreement with our analysis in Sect. 4.1), but it is even lower at larger scales (∼10 times lower). We note that the additional ionizing flux removes C⁺ as it undergoes photo-ionization to C⁺⁺, resulting in a reduction of the central flux. Hereafter we focus on this last case which provides more realistic predictions for each galaxy emission (i.e., Fig. 5), and discuss what could cause the discrepancy on CGM scales. Indeed, the results of this model suggest that there is extended [CII] emission within the inner CGM. However, the flux normalization of these extended [CII] halos compared to the observations is not matched. In what follows we analyze in detail the predicted SB profiles and explore the contributions from different gas phases.

4.2.2. Contribution of star-forming gas

By examining the contribution of star-forming gas (i.e., gas with n_H > 0.13 cm⁻³; see Sect. 3.1) to the total averaged [CII] surface brightness profile, we can refine our understanding of the underlying physical mechanisms responsible for the observed emission. Figure 7 (left panel) shows the total contribution of the dense star-forming gas (magenta dot-dashed line) to the total [CII] SB profile (orange line; including UVB and local young stars). The non star-forming gas contribution is shown by the brown dashed line. The magenta and orange lines are indistinguishable, indicating that the total [CII] emission is almost completely dominated by the star-forming gas. Since the [CII] emission in the ISM is primarily emitted from PDRs where ionized carbon atoms are heated by UV radiation from nearby young stars (e.g., Stacey et al. 2010), this result is expected at central galactic scales (up to ∼1.5 arcsec), where most of the star-formation takes place. On the other hand, understanding the origin of star-forming gas dominating in the inner CGM (0.1 ≤ R/R₂₀₀ < 0.6, or ∼1.5 ≤ r < 4 arcsec at z = 5) is more complex.

Fig. 7. Contribution and distribution of the star-forming gas. Left: Averaged contribution of the star-forming gas (dot-dashed magenta line) and non star-forming gas (dashed brown line) to the total [CII] emitting gas (orange line). The shaded areas enclose 95% of different results obtained by different galaxy orientations. Right: SFR density map for a TNG50 galaxy at z = 5. The box goes out to ∼40 kpc or ∼R200 from the central galaxy.

Indeed, our selected TNG50 galaxies present star-forming gas extended up to R₂₀₀. This can be appreciated in Fig. 7 (right panel), that shows a projection of star-forming gas for one galaxy at z = 5 (same galaxy shown in Fig. 4). In Sect. 5.1 we discuss the possible origin of this large-scale star-forming gas.

4.2.3. High and low-SFR subsamples

In order to understand the relationship between the [CII] emission and SFR, G20 divided their sample in two subsamples: high and low-SFR galaxies, using as threshold 25 M_⊙ year⁻¹. Then, to determine the typical size of the stacked [CII] line core (i.e., between [−200, 200] km s⁻¹; see Sect. 2), they computed SB radial profiles for both subsamples (see Fig. 8). They find that the subsample of high-SFR galaxies show a more extended profile (radial scales > 10 kpc at z = 5), while the profile of low-SFR galaxies drops around r ∼ 1.5 arcsec (or r ∼ 10 kpc). They claim that this result supports the scenario where star-formation feedback is responsible for the presence of [CII] halos.

Fig. 8. Averaged [CII] SB profiles (solid lines) assuming C+ abundances from the UBV and local young stars’ case. The blue line is the averaged profile for the high-SFR subsample, and the green one for the low-SFR subsample. The black triangles and squares represent the results for high and low-SFR galaxies of the observations, respectively. The shaded areas encompass the 2.5 and 97.5 percentiles of the distribution of all profiles (for each galaxy and each rotation).

In order to further compare with the observations, we did the same analysis dividing our sample in high and low-SFR using the same threshold of 25 M_⊙ year⁻¹. This comparison is shown in Fig. 8, where the triangles (squares) represent the high (low) SFR subsample of G20, and the blue (green) solid curve represents the high (low) SFR subsample of the TNG50-1 galaxies. Our high and low-SFR subsamples contain 21 and 50 simulated galaxies, respectively. The shaded areas encompass the percentiles 2 and 98 of the whole distribution of averaged profiles per galaxy per rotation. The black symbols are the averaged profiles obtained by the observations. At central galactic scales (r < 1.5 arcsec), the profile of the high-SFR galaxies is higher (∼2 times) than the one of low-SFR galaxies, which is expected (e.g., Fig. 5). Between 1.5 < r < 2 arcsec both profiles are in agreement. Then, at larger scales (i.e., CGM scales), the profile of the high-SFR subsample overcomes the one of low-SFR galaxies. However this is a slight difference compared to what is found in the observations, where the high-SFR galaxies profile is ∼4 times higher than the low-SFR galaxies profile. Hence, according to our model, [CII] halos are expected with similar intensities in both high and low-SFR galaxies. A possible reason for this discrepancy is that the simulations may not accurately capture the complex physical processes responsible for producing [CII] halos in galaxies, such as the effects of feedback from star formation or active galactic nuclei. A detailed discussion is presented in Sect. 5.1.

4.2.4. Contribution of inflows and outflows

Understanding the role of inflows and outflows in galaxy formation and evolution is essential to fully comprehend the physical processes that govern the growth and properties of galaxies. In this context, the identification of the [CII] contribution from each of these components allow us to disentangle the impact of inflows and outflows on the overall [CII] emission properties of the studied galaxies. In particular, cold (T ∼ 10⁴ K) inflows and outflows have been proposed as a possible origin of [CII] halos (Fujimoto et al. 2019). In the case of cold gas inflows, the infalling streams can produce shock heating, resulting in the production of [CII]. Conversely, extended [CII] emission may be generated by ionized carbon, powered by outflows originating from AGN or star formation processes. In this context, we estimated the gas radial velocities “v_rad” with respect to the center of the halos, specifically the position at the minimum gravitational potential. We denote as outflow the gas with v_rad > 0, and inflow as v_rad < 0 km s⁻¹.

The results of this analysis are shown in Fig. 9, where we separate the contribution to the [CII] radial profile from inflowing and outflowing gas. The blue dashed line indicates the contribution of inflowing gas, while the red dashed line indicates outflows. We can appreciate that on galactic scales (up to ∼1.5 arcsec) the contribution of inflows and outflows is balanced. On larger scales the contribution of the infalling gas dominates, increasing with radial distance.

Fig. 9. Averaged contribution of the infalling gas (dashed blue line) and outflowing gas (dashed red line) to the total [CII] emitting gas (orange line). The shaded areas represent the 2σ uncertainty obtained using 1000 samples of galaxies with randomly picked orientations, as explained in Sect. 3.4.

In order to further investigate [CII] as a kinematical tracer, the left panel of Fig. 10 shows the gas radial velocity versus radial distance normalized by R₂₀₀, weighted by [CII] luminosity, for a stacking of all the selected halos at redshift 4, 5, 6. The histograms in the right panels show the total distribution of inflows and outflows at different radial ranges: 0 ≤ R/R₂₀₀ < 0.1, 0.1 ≤ R/R₂₀₀ < 0.6, and 0.6 ≤ R/R₂₀₀ < 1. As a reference, the average R₂₀₀ for the TNG50 selected galaxies at z = 5 is ∼40 kpc (or ∼6 arcsec). For the first radial bin (inner galactic scales), our predictions indicate a balance between inflows and outflows. For the second radial bin (inner CGM), we note that the contribution of inflows exceeds the contribution of outflows by 40%. At larger scales (third radial bin), the contribution of inflows (98%) dominates the [CII] emission.

Fig. 10. Contribution of inflows and outflows. Left: Gas radial velocity versus radial distance (normalized to the virial radius), weighted by [CII] luminosity, for a stacking of all the selected TNG50 galaxies (i.e, at redshifts 4, 5, and 6). Velocities greater than 0 denote outflows, and less than 0 inflows. Right: the histograms show the total distribution of inflowing and outflowing gas, at different radial bins (see legend).

These results suggest that infalling gas may play a more significant role in the [CII] emission at larger scales, potentially indicating the presence of gas accretion processes that contribute to the buildup of gas reservoirs in high-redshift galaxies. These results are intriguing due to the expectation that inflowing gas from the intergalactic medium typically possesses low metallicity (e.g., Pallottini et al. 2014). Consequently, it is anticipated that this gas would have a minimal impact on [CII] halos. Therefore, these results raise interesting implications for our understanding of metal-enriched gas inflows. We further explore these results in Sect. 5.

4.2.5. Contribution of satellites

Another potential origin of [CII] halos is related to satellite galaxies (e.g., Fujimoto et al. 2019). This scenario suggests that the presence of satellite galaxies surrounding central star-forming galaxies gives rise to the observed extended [CII] emission. In this context, we analyzed the contribution of gas originating from sources other than the central galaxy, specifically from satellite galaxies, to the total [CII] emission. Our sample excludes systems that have “observable” companion galaxies, that is, with SFR > 1.5 M_⊙ and at a projected distance < 3 arcsec from the central galaxy, to be in agreement with G20. However, we analysed the contribution of low-SFR satellite galaxies (that would not be detectable by observations) that might still be present at small and large scales.

In this context, Fig. 11 shows the contribution of satellites (purple dashed line) and the contribution of the central galaxies (dashed gray line) to the total [CII] profile (orange line, averaged for our sample of 72 galaxies at z = 4, 5 and 6). As expected, we note that for the inner part of the galaxies, i.e between 0 and ∼2 arcsec, the contribution from the central galaxies dominate. However, there is also a very low contribution from the satellites because of projection effects. We also note that beyond ∼2 arcsec the contribution of satellites starts to dominate. This contribution increases up to ∼100% at R₂₀₀ (∼7 arcsec).

Fig. 11. Averaged contribution of gas within the satellites (dashed purple line) and central galaxies (dashed gray line) to the total [CII] emitting gas (orange line). The shaded areas enclose 95% of different results obtained by different galaxy orientations.

Interestingly, these results suggest that gas that is part of or associated with satellites galaxies may play an important role in observations of [CII] halos, especially when analyzing stacked data, in contrast to what has been argued in observational works (e.g., Fujimoto et al. 2019, 2020).

To further explore the contribution of satellites, Fig. 12 shows the same analysis as Fig. 10, but only for the gas from satellites, neglecting the contribution of the central galaxies. From the 2D histogram (left panel), it is evident that companion galaxies (clumpy structures) are already present on small scales, ∼0.1 R/R₂₀₀ (or 4 kpc/0.7 arcsec at z = 5, with low contribution to the total [CII] emission). From the histograms at the right panel, we note that after 0.6 R/R₂₀₀ almost all (99%) of the gas from satellites is inflowing toward the central halo. These results will be further elaborated and discussed in the next section.

Fig. 12. Same as Fig. 10, but for a stacking of all satellites (rejecting the contribution of the central galaxies) of all the selected TNG50 halos.

5. Discussion

In this section we explore different properties associated with extended [CII] emission based on the results found in this work. Subsequently, we aim to interpret these properties and discuss potential scenarios concerning the origin of [CII] halos.

5.1. Origin of large-scale star-forming gas

The CGM is a complex and dynamic environment influenced by different physical processes operating across various scales. These processes include gas accretion, feedback from star formation and black holes, and interactions with the surrounding environment (see Tumlinson et al. 2017 and references therein). In this context, at high-z the presence of star-forming gas extending beyond the body of the central galaxies is not unexpected in the TNG model (e.g., Donnari et al. 2019). In particular, Pillepich et al. (2019) studied the evolution of stellar and gaseous components in star-forming TNG50 galaxies at 0 < z < 6. They found that, for all masses and redshifts, the star-forming gas (as traced by Hα) settles into disky or elongated morphologies. At high-z, they also find that the star-forming gas often exhibits a complex structure, characterized by significant asymmetries on large scales (beyond the stellar body), which can be attributed to the prevalence of inflows and outflows during this epoch.

In Sect. 4.2.2, we show that the star-forming gas in the TNG50 galaxies is extended up to ∼R₂₀₀, dominating the [CII] SB profiles from inner galactic to CGM scales. One option is that this extended star-forming gas originates from the central galaxy. Galaxies undergoing active star formation can generate strong outflows that expel gas into their surroundings. This gas could then cool down and lead to star formation in the halo. Another option is that this star-forming gas comes from satellites galaxies. Interactions can trigger star formation and redistribute gas within and around galaxies, or gas can be stripped or ejected from satellite galaxies. Furthermore, the extended star-forming gas may be fueled by accretion of cold gas from the intergalactic medium along filamentary structures.

Considering the kinematic analysis performed in Sect. 4.2.4 together with the contribution of satellites in Sect. 4.2.5, we conclude that the large-scale (r > 1.5 arcsec) star-forming gas is mostly infalling and belongs to the satellites galaxies. Specifically, from Fig. 12 we obtain that at R/R₂₀₀ > 0.6 almost all the gas belonging to satellites is infalling (right third histogram). Another illustrative example is presented in Fig. 13, that shows an averaged SFR-weighted radial velocity map of the star-forming gas for a galaxy at z = 5. Inflows are represented in blue and outflows in red. It is evident from the figure that the majority of the large-scale star-forming gas primarily consists of infalling gas.

Fig. 13. Average SFR-weighted radial velocity map for a galaxy at z = 5. Only the star-forming gas is shown. The box goes out to ∼40 kpc or ∼R200 from the central galaxy.

As a consistency check, we also analyzed the metallicity distribution of the star-forming gas. The chemical signatures within the CGM are influenced by different mechanisms such as galactic outflows from supernovae and black holes, and the ejection or stripping of metals from satellite galaxies (e.g., Tegmark et al. 1993; Shen et al. 2013; Nelson et al. 2019b; Torrey et al. 2019; Péroux et al. 2020). Figure 14 shows an averaged SFR-weighted metallicity map (right panel) of the star-forming gas for a galaxy at z = 5. Its corresponding averaged radial profile is shown in the left panel. As expected, the gas that is in the central body of the galaxy exhibits higher metal enrichment (by almost 2 orders of magnitude) compared to the large-scale accreting gas. Additionally, we note that the infalling gas streams in the CGM are already metal polluted at this redshift, with typical values in the range of 10⁻¹–10⁻² Z_⊙. As a reference, the dominant metallicity of outflowing gas for z ∼ 4 galaxies (M_* ∼ 10^9.5 M_⊙) has been found to be ∼10^−0.7 Z_⊙ (Nelson et al. 2019b), while gas from the intergalactic medium is expected to have a low-metallicty of ∼10⁻³ Z_⊙ at z ∼ 4 − 6 (Pallottini et al. 2014). Overall, this result is expected since metal-enriched gas is necessary for the presence of C⁺ ions and consequently [CII] emission.

Fig. 14. Gas metallicity distribution for a galaxy at z = 5. Left: Average radial profile of the star-forming gas metallicity (weighted by SFR). Right: Average SFR-weighted metallicity map. Only the star-forming gas is shown. Its corresponding radial profile is shown in the left panel. The box size radius is ∼40 kpc or ∼R200.

5.2. Role of satellite galaxies in producing [CII] halos

Most observational works studying [CII] halos have focused on isolated galaxies (i.e., with no observable minor or major mergers, e.g., Fujimoto et al. 2019, 2020; G20; Herrera-Camus et al. 2021; Fudamoto et al. 2022). However, contamination from possible faint undetected satellites cannot be discarded. The presence of satellites may affect the gas dynamics, heating processes, or feedback mechanisms, influencing the overall [CII] emission properties. In this context, some attempts from observations have been made to estimate the contribution of such undetected satellites to the observed [CII] halos. For instance, Fujimoto et al. (2019) investigated the possible contribution of satellites to the extended [CII] emission they found in a stacking of data of star-forming galaxies at z = 5 − 7. They estimated radial values of the L_CII/SFR ratio, finding the highest ratios in the outer regions. These values are in disagreement with those expected for local dwarf galaxies. The authors interpret this as an indication that satellites are likely not the origin of the observed [CII] halo. On the other hand, Ginolfi et al. (2020b) studied one merging system at z ∼ 4.6, from the ALPINE survey. They find a [CII] halo extending up to r ∼ 15 kpc, and argue that its origin is mainly interstellar gas stripped by strong gravitational interactions, with a possibly less significant contribution from galactic ouflows and star formation from small faint satellites. Schimek et al. (2024) presented a theoretical effort of simulating [CII] emission in a zoom-in simulation of a system at z = 6.5 undergoing a major merger. Their results indicate that the [CII] emission in the CGM is higher in an accreting filament and tidal features from the merging galaxies. Recently, Di Cesare et al. (2024) studied the [CII] emission in the CGM of six major mergers from the ALPINE survey. They find extended [CII] emission (≳20 kpc) in all the systems. By comparing with simulated galaxies, they infer that the origin of this emission is tidal stripped gas and the presence of unresolved star-forming satellites.

From the previous section and the results presented in Sect. 4.2.5, we find that star-forming gas from satellites or companion galaxies contribute substantially to the extended [CII] emission beyond ∼2 arcsec. This suggests that faint satellites, which may be undetected in observations, and their associated gas could be important contributors to the [CII] halos observed around galaxies. However, as we are not obtaining the full emission level of the observed [CII] profile on CGM scales, we cannot confirm that. A possibility for this missing flux is that there might be other physical processes affecting the extended [CII] emission that are not well reproduced by TNG50 (see the next section). To test what would be the maximum contribution of satellites in our model, we derived the average [CII] SB profile for a larger sample of TNG50 galaxies (169) that covers the physical properties range of in the sample of G20 (see Sect. 3.2), but without restricting to similar distributions in SFR and stellar mass. Importantly, we select these galaxies without removing the systems with mergers. The result is presented in Fig. 15, where the C⁺ abundances were estimated by assuming UVB and young stars as ionizing sources. Interestingly, the level of [CII] emission at large scales is ∼8 times higher than what is found in our cleaned sample (see Fig. 6), being closer to the observations. This supports the scenario in which the presence of satellites and their interactions with the central galaxies could have a significant impact on the [CII] emission in the inner CGM. This effect could be especially important in observational results based on stacking techniques.

Fig. 15. Averaged [CII] SB profiles obtained for a sample of 169 TNG50 galaxies at z = 4, 5, and 6 (orange line) which include close satellites. The shaded area encloses 95% of different results obtained by different galaxy orientations. The C+ abundances have been estimated considering UVB and local young stars as ionization sources.

In summary, our model, which shows that [CII] traces star-forming gas, suggest that satellites (including faint and undetected ones) can play an important role in shaping the extended [CII] emission in the inner CGM of high-z galaxies. However, further investigation is needed to fully understand the specific role of satellites in reproducing the exact levels of the observed extended [CII] emission. Even when including galaxies with close satellites, we are not able to fully match the observations. Simulations capable to resolve less massive galaxies (TNG50 resolves galaxies with M_* ≳ 10⁷ M_⊙; see Sect. 3.1) would be essential. Future theoretical works exploring the detailed properties of satellites, their gas content, and interactions with the central galaxies could provide insights into the mechanisms responsible for the observed [CII] halos of high-z galaxies.

5.3. Role of outflows in [CII] emission of high and low-SFR galaxies

Based on a spectral and cube stacking analysis, G20 argue that [CII] halos trace circumgalactic gas that has been previously enriched by past outflows driven by star formation activity. They found extended [CII] emission only for the subsample of galaxies with high-SFR (SFR > 25 M_⊙ yr⁻¹), implying that galaxies with more active star formation should present more extended [CII] emission. This interpretation is further supported by semi-analytical models (Pizzati et al. 2020, 2023). For instance, Pizzati et al. (2023) implemented a model to predict [CII] halos originating from outflowing gas, which incorporates key factors such as the outflow mass loading factor⁹ (η), the parent galaxy SFR, and the circular velocity of the dark-matter halo. In particular, the value of η is related to the radiative cooling of the wind and its capability to reach low temperatures (T ∼ 10² − 10⁴ K). Conversely, the gravitational potential of the dark-matter halo can substantially decelerate the gas expansion. These authors compare their results with ALPINE individual galaxies (from Fujimoto et al. 2020) finding that detected [CII] haloes are a natural outcome of starburst-driven outflows. Their model also suggests that low-mass systems also present [CII] halos, however their emission levels are too faint to be detected with the current observational capabilities.

In Sect. 4.2.3, we analyzed the [CII] SB profiles of our high-SFR and low-SFR subsamples, finding no major difference between both profiles at CGM scales, in contrast with the observational results of G20. Both subsamples present extended [CII] emission; however, the emission of the high-SFR subsample should be higher (up to ∼20 times) to match the observations. To have more insights regarding the reason of this discrepancy, we constructed the [CII] radial SB profile of the outflowing gas only and look for differences in both SFR subsamples (for completeness, the same is done for the inflowing gas). This exercise is shown in Fig. 16, where the dashed line indicates the average for high-SFR galaxies (SFR > 25 M_⊙ yr⁻¹), and the dotted line for low-SFR galaxies (SFR < 25 M_⊙ yr⁻¹). The radial axis extends to ∼7 arcsec (∼R₂₀₀). We note that for galactic scales (r < 2 arcsec), the SB of the outflowing gas is ∼2 times higher for the high-SFR galaxies. Then, at r ∼ 2 arcsec there is a drop in the profile of the high-SFR subsample, which then raises again, overcoming the outflowing gas profile of low-SFR galaxies at r ∼ 2.8 arcsec, until ∼R₂₀₀ (by a factor of ∼10). This indicates that outflows in the high-SFR galaxies are more energetic and hence more extended, contributing to a larger observable [CII] SB for scales larger than r ∼ 2.8 arcsec. However, this increased outflow contribution is still much smaller than the contribution from infalling star-forming gas in our model. If outflows were the main contributors to the [CII] halos, their phase diagram in the TNG50 simulation are not accurate enough to produce observable levels of emission. A detailed discussion about the implications of the TNG feedback model is presented in Appendix A.

Fig. 16. Averaged [CII] SB profiles for the outflowing (inflowing) gas only, represented by the red (blue) color. Our original TNG50 sample has been divided in two subsamples: high-SFR (dashed lines) and low-SFR (dotted lines). The C+ abundances have been estimated considering UVB and local young stars as ionization sources.

In summary, our predicted [CII] emission in TNG50 galaxies suggests that [CII] halos originate through a combination of contribution of star-forming gas (associated with satellites galaxies on CGM scales) and outflows from galactic scales. The contribution of inflowing gas from satellites dominates beyond ∼2 arcsec. However, compared to the observations we do not match the level of [CII] emission on CGM scales. This discrepancy could be related to the limited resolution of the simulations, issues in the feedback model of TNG, and/or other physical processes not implemented in our model. The next section details the caveats of our work.

5.4. Caveats

There are certain limitations and considerations related to the simulations and model assumptions that need to be taken into account when interpreting the results and conclusions presented in the previous sections:

• Cold dense gas clumps in the CGM: The TNG50 simulations do not fully resolve cold high-density gas within the CGM, potentially resulting in an underestimation of the [CII] emission. It is possible that the presence of cold clumps in the CGM would contribute to higher [CII] emission levels than those captured in this work. These cold clumps are known to endure within outflows, which can be intensified by cosmic rays, a factor not incorporated in IllustrisTNG (e.g., Rathjen et al. 2021). The presence and survival of dense CGM gas is supported by observations (e.g., Cantalupo et al. 2014; Arrigoni Battaia et al. 2015; Hennawi et al. 2015; Vidal-García et al. 2021) and high-resolution “cloud-crushing” simulations (e.g., Gronke & Oh 2018, 2020; Kanjilal et al. 2021). In particular, Nelson et al. (2020) attempted a resolved study of small-scale cold CGM gas in TNG50, finding the presence of cold clouds whose number and mass increase with resolution. Despite these efforts, such processes are still largely unresolved, and cosmological models of galaxy formation are not converged in their CGM properties (e.g., Oppenheimer et al. 2021). In addition, most of the [CII] emission comes from thin PDRs (e.g., Gullberg et al. 2015), which remain unresolved in our model. Given the critical densities associated with C⁺ and, more significantly, the dependence of Λ_CII on particle number density (Λ_CII ∝ n²; see Eqs. 1 and 2), it becomes necessary to incorporate higher densities, potentially through the application of density distributions on sub-grid scales (e.g., Olsen et al. 2021; Schimek et al. 2024; Khatri et al. 2024).

• Coldest temperature: The temperature of the star-forming gas cells is assumed to be fixed at 1000 K (cold mode, Springel & Hernquist 2003). While this specific value is motivated by the subgrid model used in the simulations, in reality star-forming gas is characterized by a distribution of temperatures which spans also lower values (e.g., Kainulainen et al. 2009; Draine 2011; Olsen et al. 2021). Colder dense gas (currently not modeled in the TNG50 suite), down to the [CII] excitation temperature of ∼91 K, could enhance the production of [CII] emission.

• SFRs in large-box cosmological simulations: Reproducing the highest levels of SFR estimated in observations is a well known issue for large-box cosmological simulations (e.g., Furlong et al. 2015; Donnari et al. 2019; Pallero et al. 2019; Aoyama et al. 2019; Katsianis et al. 2020). The low statistics or absence of very high SFR galaxies in the simulations could affect our analysis. In particular, the number of high-SFR galaxies studied here is smaller than the number of low-SFR galaxies. Our high-SFR subsample consists of 21 galaxies (versus 50 for the low-SFR subsample), and are in the range of 25−270 M_⊙ year⁻¹, while in the observations they analyze 20 high-SFR galaxies in the range of 25−623 M_⊙ year⁻¹ (versus 20 low-SFR). The most star-forming systems in the observations, those at ∼600 M_⊙ year⁻¹, are dominating the averaged [CII] halo emission. However, in TNG50 such systems are simply not available.

• Turbulence dissipation in the CGM: The CGM is expected to be a turbulent environment, influenced by different perturbations associated with galaxy interactions, accretion, and outflows (e.g., Vidal-García et al. 2021). Turbulence dissipation affects the dynamics, thermal properties, and chemical enrichment of the CGM (see Faucher-Giguère & Oh 2023 and references therein). In the cold CGM, turbulent widths have typical values of ∼10−30 km s⁻¹, which indicates random motions with velocity dispersions of σ ∼ 1000 km s⁻¹ (e.g., McCourt et al. 2018; Rudie et al. 2019). Currently, these processes are not followed by the TNG suite. Considering the above, the inclusion of sub-grid scale high-velocity turbulences could be essential to enhance the frequency of collisions between C⁺ and its colliding partners, leading to an increased [CII] emission in the CGM (e.g., Godard et al. 2014).

• H₂ as collisional partner: In our post-processing model, we have neglected H₂ as colliding partner because of the reasons explained in Sect. 3.3, and especially because its physics is not modeled by the simulations. [CII] emission due to collisions with H₂ is characterized by the highest critical density when compared to collisions with e− and H (see Sect. 1), and such densities are not present on CGM scales in the TNG50 simulations. Additionally, H₂ abundance is assumed to be dependent on metallicity (e.g., Krumholz et al. 2018); hence, minimal H₂ is expected on halo scales where metallicities are lower than those toward galactic scales. However, there is growing observational evidence of the presence of molecular gas outside the main body of high-z galaxies, i.e reaching the inner CGM (e.g., Ginolfi et al. 2017; Falgarone et al. 2017; Emonts et al. 2019, 2023; Li et al. 2020; Vidal-García et al. 2021). This suggests that H₂ may indeed be an important collisional partner out to the inner CGM and contributor to the [CII] emission in those regions.

• Assumptions in the CLOUDY calculations: To obtain the C⁺ abundances in the simulated systems, we used CLOUDY tables by assuming a fixed escape fraction of ionizing photons of 5%. However, it is important to note that this escape fraction might be variable at different radial scales, as well as the optical thickness of the medium. Some works suggest that f_esc should be very small (< 1%) to produce the observed [CII] halos (e.g., Pizzati et al. 2020, 2023), and it may also vary for different galaxy mass ranges (e.g., Xu et al. 2016). In addition, the CLOUDY calculations and the simulations setup neglects the effects of dust. However, dust is known to be important in galaxy evolution (e.g., Draine 2004). For instance, dust is key for the chemistry of the gas by capturing selected elements onto its grains and by catalyzing formation of the H₂ molecule. Dust grains can contribute to gas heating mainly through the photoelectric effect (Weingartner & Draine 2001), and gas cooling mainly through free-particle capture onto the grain surface and collision of ions with grains (e.g., Burke & Silk 1974; Montier & Giard 2004). Specifically, in PDRs, where [CII] has been shown to be the main coolant, the cold dust emission spatially correlates with [CII] emission, indicating the importance of gas heating through the photoelectric effect (e.g., Habart et al. 2001).

6. Summary

To investigate the nature of [CII] halos detected in observations around galaxies at z > 4, we conducted a post-processing analysis on the TNG50 simulation. We aimed to model [CII] emission in a sample of simulated galaxies at z ∼ 4 − 6. Specifically, by incorporating C⁺ abundances derived by assuming UV background and young stars as radiation sources, we estimated [CII] emissivities and generated mock observations, comparing our results with the ALPINE [CII] emitting galaxies investigated by G20. This approach allowed us to study the possible origin of extended [CII] emission in high-z galaxies and to discuss the contribution of different factors to this emission. The main results of this work are summarized as follows:

(i) Our model reproduces the L_CII versus SFR relationship found by previous works in local and high-z galaxies. This result provides validation for the predictions from our model on galactic scales.

(ii) We derived an averaged [CII] SB profile and compared it to the observed data. Our predicted [CII] emission levels on central scales (r < 1.5 arcsec or ∼10 kpc at z = 5) are in agreement with the observations. However, at CGM scales, our results underestimate the [CII] emission levels by a factor of ∼10.

(iii) The [CII] emission in our systems is dominated by the dense star-forming gas, as expected. This gas extends up to the virial radius.

(iv) From a kinematic analysis, we find that on inner galactic scales (r < 0.7 arcsec or 5 kpc at z = 5), the [CII] emission traces a balance between inflows and outflows (contribution of 51 and 49%, respectively).

(v) The contribution of infalling gas from satellites dominates the [CII] emission in the CGM, specifically after r ∼ 2 arcsec (or ∼13 kpc at z = 5). This suggests that faint and undetected satellites and their associated gas may play a significant role in shaping the observed extended [CII] emission.

(vi) We find further tension with the observations when comparing the SB profiles for the low- and high-SFR galaxies, finding no significant difference in the extended [CII] emission levels (up to r ∼ 3 arcsec or ∼20 kpc at z = 5) between the two subsamples. We propose that this could be an effect of the feedback model and the limited resolution of the TNG simulations.

(vii) According to our results, the origin of [CII] halos is a combination of the contribution of gas from satellites and outflows from the central galaxies. However, the contribution of infalling gas by galaxies completely dominates at 4 < r < 7 arcsec (0.6 < r/R₂₀₀ < 1).

This theoretical effort highlights the complex nature of extended [CII] emission in the inner CGM of high-redshift galaxies. Further investigations based on state-of-the-art simulations and observations are necessary to fully understand the specific role of satellites and other physical processes (e.g., galactic winds and turbulences) in reproducing the observed [CII] halos. In particular, cosmological simulations with a higher resolution in the ISM and CGM are needed. While running cosmological simulations with higher refinements is currently computationally expensive and requires tuning, one could achieve improved results with the use of more detailed subgrid prescriptions based on radiative transfer calculations run on high-resolution turbulent boxes (Buhlmann et al., in prep.). On the other hand, future observations from large submillimeter single-dish telescopes (e.g., AtLAST; Klaassen et al. 2020) in combination with high-resolution interferometric observations (e.g., CRISTAL; Posses et al. 2024) will be able to improve our general knowledge of the CGM and specifically of [CII] halos in high-z objects and verify whether there are any biases in the current small observational samples (e.g., a selection bias favoring more luminous [CII] halos at a given SFR).

Acknowledgments

The authors thank the anonymous referee for their detailed comments that helped improving the manuscript. The authors like to thank Dylan Nelson for insightful comments after the presentation of an earlier version of this work. A.O. has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 443044596.

References

Appendix A: TNG feedback model and [CII] emission

The TNG model assumes that star formation feedback drives galactic outflows, which are launched from star-forming gas with a wind velocity (v_w) that depends on the local dark matter velocity dispersion (σ_DM, e.g., Vogelsberger et al. 2013). The outflow mass loading is determined by the supernova energy and the assumed wind speed, while the metal content is a fraction of ISM metallicity. The model uses a kinetic wind scheme where wind particles are created and decoupled from the gas until they leave the dense ISM. Outside the dense medium, the wind particles recouple with the gas, transferring their mass, momentum, metals, and thermal energy to the surrounding medium. Wind particles recouple to the gas cell when: 1) the density drops below a threshold (n_H = 0.05 × density threshold for star formation = 0.006 cm⁻³), or 2) after a maximum travel time has been reached (0.025 × Hubble time at respective z) (Pillepich et al. 2018).

For the TNG50 galaxies at z = 5 (median of the sample) studied here, the average distance¹⁰ at which the wind particles would recouple with the gas from the condition 1), is ∼5 arcsec (or ∼32 kpc) and ∼4.6 arcsec (or ∼30 kpc) for the high and low-SFR subsamples, respectively. On the other hand, following condition 2), the maximum distance¹¹ after the maximum travel time is ∼2.8 arcsec (or ∼18 kpc) and ∼2.4 arcsec (or ∼15 kpc) for the high and low-SFR subsamples, respectively. These results indicate that for the case of our selected halos, the wind particles recouple with the gas cells according to the second criteria, which is first met. The different radial distances estimated for both subsamples confirm that the outflows of the high-SFR galaxies propagate to larger distances compared to the low-SFR galaxies (Fig. 16). On the other hand, the distance of ∼18 kpc (∼2.8 arcsec) obtained for the high-SFR subsample is comparable to the distance that is reached by the outflows in the observations of G20 (r ∼ 15 kpc). It also coincides with the upturn of the outflowing gas SB profile, ocurring immediately after the drop in this profile (∼2 arcsec; Fig. 16). This suggests that, as the wind particles recouple with the gas and transfer their properties, it can lead to a change in the gas properties and potentially in an increase in the [CII] emission. Therefore, if all of the observed extended [CII] emission is attributed to stellar winds, it would be necessary to incorporate a new subgrid prescription able to enhance the wind contribution (i.e., the outflow mass loading factor) significantly. This increase could be achieved by including a prescription for an entrained clumpy, cold dense phase that is currently missing (e.g., Hamann et al. 2013; Liang et al. 2016; Faucher-Giguère et al. 2016; McCourt et al. 2018; Gronke & Oh 2018, 2020).

Finally, it is important to discuss the effect of AGN feedback. While the observed samples do not report the evidence for AGN activity in the targeted galaxies down to current observational limits, all our selected TNG50 galaxies have AGNs (L_BOL ≲ 10^45.5 erg s⁻¹) affecting the surrounding gas with thermal mode feedback (i.e., high accretion rates, average Ṁ/Ṁ_Edd = 0.5). The average AGN feedback energy is higher for the high star-forming sample compared to the low star-forming sample (ΔE = 10^43.9 erg s⁻¹ and 10^43.3 erg s⁻¹, respectively). This stronger AGN feedback in the most star-forming galaxies may lead to higher gas temperatures, resulting in similar or lower [CII] emission compared to the low star-forming galaxies. Importantly, an AGN would also affect the C⁺ abundances with its hard radiation. This effect is not taken into account in our calculations, but would increase the tension with observations as we expect Carbon to be more ionized when accounting for an AGN contribution.

Appendix B: Example of input CLOUDY file

In the current version of CLOUDY (C17.03), the UVB models of Khaire & Srianand (2019) can be easily accessed with the table KS18 redshift [value], while for our specific choice of young stellar radiation field we pass the spectrum with the table SED "[name_spectral_file]" and provide the normalization as the logarithm of the ionizing photon flux phi(H) [value] corresponding to each radiation field normalization ϕ. Lines starting with "#" in the following input file example denote comments. The output file *.C_ionf contains the ionization fractions for all the levels of the carbon atom.

# set the hydrogen number density
constant density
hden [value] log
# set the temperature
constant temperature, t = [value] K log
# set the metallicity
metals [value] log
# set the incident radiation fields
table KS18 redshift [value]
CMB redshift [value]
table SED "[name_spectral_file]"
phi(H) [value]
# set the stopping criteria
stop zone 1
# set the output
set save prefix "[name_input_file]"
save element carbon last ".C_ionf"
iterate to convergence

For the UVB-only models, we use the same input file, but remove the two lines table SED "[name_spectral_file]" and phi(H) [value].

All Figures

Fig. 1. Comparison of the distribution of the SFRs and M* values between the simulated (red) and observed (blue) samples (75 and 50 galaxies, respectively). The joint distribution is visualized using density contours. The contour spacing represents the density of data points, with closer contours indicating regions of higher data density. The blue bars in the bottom-right corner represent the average 1-sigma error for the SFR and M* of the observations. The histograms along the x and y axes provide marginal distributions of the M* and SFR, respectively, for both observed and simulated samples.

In the text

	Fig. 2. Phase diagram for one of the selected TNG50 galaxies at z = 5, weighted and color-coded by the gas mass. The diagram encompasses all gas cells within R200. For densities above the star-formation threshold nH = 0.1, the temperatures are determined by the effective equation of state (Springel & Hernquist 2003).
In the text

Fig. 3. Example of different input spectra used in the CLOUDY calculations. The solid light-gray line represents the UVB SED of Khaire & Srianand (2019), and the dashed gray line the CMB, both at z = 5. The vertical black lines represent the ionization energies of C→C+ (dotted) and C+→C++ (dashed). The colored lines represent the contribution to the SED due to young stars (Cerviño et al. 2002) at different normalizations, plus the contribution of the UVB. The units of the radiation field ϕ are M⊙ year−1 kpc−2.

In the text

Fig. 4. [CII] emission surface-brightness map for one of the selected TNG50 galaxies at z = 5, assuming C+ abundances with UVB and local young stars as ionization sources. Left: the map is presented at the full-resolution of the simulation. Right: the map has been convolved with the average ALMA beam of the G20 observations (bottom left corner). The white contour represents the average 2σ level of the observations (0.08 Jy km s−1 beam−1).

In the text

Fig. 5. [CII] luminosity versus SFR for the TNG50 galaxies (left panel) and the ALPINE survey (right panel). The shaded blue area represents the 1σ relationship of De Looze et al. (2014) for starburst local galaxies. The dashed red line shows the relationship of Lagache et al. (2018) for z = 5 star-forming galaxies. The average 1σ error for L[CII] and SFR in the right panel are 0.05 dex and 0.4 dex, respectively.

In the text

Fig. 6. Averaged [CII] SB profiles obtained for our sample of 72 TNG50 galaxies at z = 4, 5, and 6 (solid lines), with shaded areas enclosing a 95% of different results obtained by different galaxy orientations. The different line colors represent the TNG50 profiles obtained assuming C+ abundances equal to the total C abundance (dark pink), C+ abundance obtained considering the UVB as incident field (light blue) and C+ abundance considering UVB and local young stars as ionization sources (orange). As a comparison, the black squares represent the results from the ALPINE survey.

In the text

Fig. 7. Contribution and distribution of the star-forming gas. Left: Averaged contribution of the star-forming gas (dot-dashed magenta line) and non star-forming gas (dashed brown line) to the total [CII] emitting gas (orange line). The shaded areas enclose 95% of different results obtained by different galaxy orientations. Right: SFR density map for a TNG50 galaxy at z = 5. The box goes out to ∼40 kpc or ∼R200 from the central galaxy.

In the text

Fig. 8. Averaged [CII] SB profiles (solid lines) assuming C+ abundances from the UBV and local young stars’ case. The blue line is the averaged profile for the high-SFR subsample, and the green one for the low-SFR subsample. The black triangles and squares represent the results for high and low-SFR galaxies of the observations, respectively. The shaded areas encompass the 2.5 and 97.5 percentiles of the distribution of all profiles (for each galaxy and each rotation).

In the text

	Fig. 9. Averaged contribution of the infalling gas (dashed blue line) and outflowing gas (dashed red line) to the total [CII] emitting gas (orange line). The shaded areas represent the 2σ uncertainty obtained using 1000 samples of galaxies with randomly picked orientations, as explained in Sect. 3.4.
In the text

Fig. 10. Contribution of inflows and outflows. Left: Gas radial velocity versus radial distance (normalized to the virial radius), weighted by [CII] luminosity, for a stacking of all the selected TNG50 galaxies (i.e, at redshifts 4, 5, and 6). Velocities greater than 0 denote outflows, and less than 0 inflows. Right: the histograms show the total distribution of inflowing and outflowing gas, at different radial bins (see legend).

In the text

	Fig. 11. Averaged contribution of gas within the satellites (dashed purple line) and central galaxies (dashed gray line) to the total [CII] emitting gas (orange line). The shaded areas enclose 95% of different results obtained by different galaxy orientations.
In the text

	Fig. 12. Same as Fig. 10, but for a stacking of all satellites (rejecting the contribution of the central galaxies) of all the selected TNG50 halos.
In the text

	Fig. 13. Average SFR-weighted radial velocity map for a galaxy at z = 5. Only the star-forming gas is shown. The box goes out to ∼40 kpc or ∼R200 from the central galaxy.
In the text

	Fig. 14. Gas metallicity distribution for a galaxy at z = 5. Left: Average radial profile of the star-forming gas metallicity (weighted by SFR). Right: Average SFR-weighted metallicity map. Only the star-forming gas is shown. Its corresponding radial profile is shown in the left panel. The box size radius is ∼40 kpc or ∼R200.
In the text

	Fig. 15. Averaged [CII] SB profiles obtained for a sample of 169 TNG50 galaxies at z = 4, 5, and 6 (orange line) which include close satellites. The shaded area encloses 95% of different results obtained by different galaxy orientations. The C+ abundances have been estimated considering UVB and local young stars as ionization sources.
In the text

	Fig. 16. Averaged [CII] SB profiles for the outflowing (inflowing) gas only, represented by the red (blue) color. Our original TNG50 sample has been divided in two subsamples: high-SFR (dashed lines) and low-SFR (dotted lines). The C+ abundances have been estimated considering UVB and local young stars as ionization sources.
In the text