The informativeness of [C ii] line-intensity mapping as a probe of the H i content and metallicity of galaxies at the end of reionization Open Access

ABSTRACT

Line-intensity mapping (LIM) experiments coming online now will survey fluctuations in aggregate emission in the [C ii] ionized carbon line from galaxies at the end of reionization. Experimental progress must be matched by theoretical reassessments of approaches to modelling and the information content of the signal. We present a new model for the halo–[C ii] connection, building upon results from the Feedback In Realistic Environments simulations suggesting that gas mass and metallicity most directly determine [C ii] luminosity. Applying our new model to an ensemble of peak-patch halo light-cones, we generate new predictions for the [C ii] LIM signal at z ≳ 6. We expect a baseline 4000-h LIM survey from the Cerro Chajnantor Atacama Telescope facility to have the fundamental sensitivity to detect the [C ii] power spectrum at a significance of 5σ at z ∼ 6, with an extended or successor Stage 2 experiment improving significance to 48σ at z ∼ 6 and achieving 11σ at z ∼ 7.5. Cross-correlation through stacking, simulated against a mock narrow-band Lyman-break galaxy survey, would yield a strong detection of the radial profile of cosmological [C ii] emission surrounding star-forming galaxies. We also analyse the role of a few of our model’s parameters through the pointwise relative entropy (PRE) of the distribution of [C ii] intensities. While the PRE signature of different model parameters can become degenerate or diminished after factoring in observational distortions, various parameters do imprint themselves differently on the one-point statistics of the intrinsic signal. Further work can pave the way to access this information and distinguish different sources of non-Gaussianity in the [C ii] LIM observation.

1 INTRODUCTION

As the studies of galaxy formation and cosmology advance, galaxies at high redshift and in particular at the end of the epoch of reionization (EoR) are as important as ever. From the perspective of cosmology, understanding how the Universe exits the EoR relies on understanding the physical conditions, abundances, and cosmic structure of the first stars and galaxies that carry out the reionization process. From the perspective of galaxy formation, to understand the nature of these first galaxies and proto-galaxies is to understand the starting point of cosmic star formation history, with star formation activity continuing to rise all the way through the ‘cosmic noon’ of z ∼ 2–3 (see e.g. the review of Madau & Dickinson 2014). The EoR thus sits at a key intersection of cosmology and extragalactic astrophysics, between structure formation and galaxy formation. Studying the population of galaxies from the EoR is fundamental to understanding the Universe both as we know it today and as it was almost 14 Gyr ago.

Although targeted or pencil-beam surveys are already detecting individual objects from the EoR, conventional spectroscopic galaxy surveys will struggle to probe EoR galaxies at the efficiency required for this population to act as a true large-scale cosmological probe. Thus, a full understanding of the EoR demands alternative ways to measure the clustering of star formation and gas at high redshift.

Intensity mapping (IM), or line-intensity mapping (LIM), is a technique analysing low- to medium-resolution spatial–spectral observations to probe a three-dimensional (3D) cosmological volume for the signature of large-scale structure as traced by spectral lines associated with specific atomic or molecular species. The strength of such analyses resides in their potential to probe large cosmological volume efficiently while measuring the contribution of all galaxies, including those impractical for standard galaxy surveys to detect individually (see reviews by e.g. Kovetz et al. 2017, 2019; Bernal & Kovetz 2022).

While atomic hydrogen will be most abundant particularly at early cosmic epochs, the past decade has seen studies of LIM concepts leveraging other spectral lines associated with dust-obscured star formation activity to complement IM surveys of the 21-cm atomic hydrogen line (see Fig. 1). This includes the rotational transitions of carbon monoxide (CO) observable at centimetre wavelengths, one of the most abundant luminous tracers of molecular gas and thus the subject of significant work in LIM modelling (Righi, Hernández-Monteagudo & Sunyaev 2008; Visbal & Loeb 2010; Lidz et al. 2011; Pullen et al. 2013; Li et al. 2016; Moradinezhad Dizgah & Keating 2019; Sun et al. 2019; Yang et al. 2021, 2022) and experiment (Keating et al. 2016, 2020; Cleary et al. 2022; Keenan, Keating & Marrone 2022).

At shorter wavelengths, the 157.7-µm [C ii] line from ionized carbon (with rest-frame frequency ν_rest = 1900.5 GHz) becomes one of the most promising tracers of high-redshift galaxies. In our local volume, [C ii] is a key far-infrared (FIR) cooling line that makes up as much as 1 per cent of the bolometric FIR luminosity in local gas-rich and highly luminous galaxies (Stacey et al. 1991). Given mounting observations with the Atacama Large Millimetre/submillimetre Array (ALMA) of individual [C ii] lines at z ∼ 4–8 (see e.g. Knudsen et al. 2016; Pentericci et al. 2016; Smit et al. 2018; Yan et al. 2020; Fujimoto et al. 2021), the existence of an appreciable population of [C ii] emitters in the late EoR appears certain. The brightness of the cosmological [C ii] emission at this relatively metal-poor epoch remains uncertain, however, and has been the subject of a substantive body of prior modelling work (e.g. Gong et al. 2012; Uzgil et al. 2014; Silva et al. 2015; Yue et al. 2015; Serra, Doré & Lagache 2016; Dumitru et al. 2019; Chung et al. 2020; Karoumpis et al. 2022).

At time of writing, a range of EoR [C ii] pathfinders have begun construction or even operation, including the CarbON [C ii] line in post-rEionization and ReionizatTiOn (CONCERTO) spectrometer (Fasano et al. 2022; Monfardini et al. 2022), the Tomographic Ionized-carbon Mapping Experiment (TIME; Crites et al. 2014; Sun et al. 2021; Crites 2022), and the Epoch of Reionization Spectrometer (EoR-Spec) instrument on the Cerro Chajnantor Atacama Telescope (CCAT) facility (Nikola et al. 2022; CCAT-Prime Collaboration 2023). With all these projects in active development and/or analysis phases, and early science results anticipated within this decade, this work aims to push the modelling of the [C ii] LIM signal forward in two important ways: first by reassessing the nature of the underlying halo–[C ii] connection and second by assessing the information content of the one-point statistics of the [C ii] intensity fluctuations to be observed.

The model of this work distinguishes itself from many previous models by mediating the halo–[C ii] connection not solely through star formation rate (SFR), but through the neutral atomic hydrogen (H i) mass and gas-phase metallicity in the interstellar medium (ISM) of the galaxies in dark matter haloes, based on insights from the Feedback In Realistic Environments (FIRE) project examined by Liang et al. (2024). In formulating this halo model, we reproduce certain trends like the evolution of [C ii] luminosity against FIR luminosity, and arrive at a model consistent with current observational constraints on high-redshift [C ii] abundances.

Applying our halo model to an ensemble of approximate cosmological simulations, we aim to answer the following questions:

• What is the detectability of [C ii] emission in near- and far-future mm-wave LIM surveys, in either autocorrelation or cross-correlation?

• How do different model parameters affect the voxel intensity distribution (VID) of the line-intensity fluctuations, as represented by the differential relative entropy?

• How do observational effects like the instrument response and thermal noise affect the observability of these distinct signatures in relative entropy?

The paper is structured as follows. After a summary of context around the halo–[C ii] connection and one-point statistics in Section 2, we present our new [C ii] model and simulation methods in Section 3. We then show our results and forecasts in Section 4, and conclude with forward-looking discussion in Section 5.

Unless otherwise stated, we assume base-10 logarithms, and a Lambda cold dark matter cosmology with parameters Ω_m = 0.286, |$\Omega _\Lambda = 0.714$|⁠, Ω_b = 0.047, H₀ = 100h km s⁻¹ Mpc⁻¹ with h = 0.7, σ₈ = 0.82, and n_s = 0.96, to maintain consistency with previous simulations used by Chung et al. (2023). Distances carry an implicit h⁻¹ dependence throughout, which propagates through masses (all based on virial halo masses, proportional to h⁻¹) and volume densities (∝h³).

2 CONTEXT IN PREVIOUS WORK

2.1 A FIRE-tested formula: moving away from [C ii] as a tracer of star formation

The SFR is a highly attractive lever arm for [C ii] modelling at high redshift. With [C ii] being a dominant FIR cooling line mostly unaffected by dust extinction, [C ii] emission is a natural mechanism to balance heating from newly forming stars. Equally, [C ii] emission is attractive as a proxy for SFR at high redshifts, where a signpost for obscured star formation activity is needed for a full understanding of cosmic star formation history.

However, previous work has already cast significant doubt that the connection between [C ii] luminosity and SFR is reliably universal (Herrera-Camus et al. 2015; Vallini et al. 2015; Croxall et al. 2017), as other factors like metallicity appear to modulate the observed relation. A full review of such work is well beyond the scope of this paper, and we refer the reader to Liang et al. (2024) for further details. The upshot, however, is that such doubts motivate a re-examination of the nature of the halo–[C ii] connection based on a physically motivated picture of [C ii] emission. Such work is particularly important if we are to leverage synergies between mm-wave LIM and 21-cm tomography at high redshifts (Padmanabhan 2023).

To this end, we make use of insights from the work of Liang et al. (2024), who conducted a comprehensive analysis of the relationship between the [C ii] luminosity and SFR of simulated star-forming galaxies at z ∼ 0–8 generated by the FIRE project. The mock galaxy sample includes 507 FIREBox galaxies (including 91 at z ≥ 6), which reside not in zoom-in boxes but in a |$(15h^{-1}\,$| cMpc)³ cosmological box with prescriptions for hydrodynamics and feedback (Feldmann et al. 2023), as well as 10 MassiveFIRE zoom-in simulations (6 of which have snapshots at z ≥ 6; see references in Liang et al. 2024). Through post-processing of these simulations to derive the interstellar radiation field and the ionization and excitation states of gas clouds, Liang et al. (2024) derived [C ii] luminosities for the FIRE galaxies studied, as well as other properties like the SFR and the bolometric infrared (IR) luminosity L_IR.

Not only do the FIRE simulations reproduce the empirically observed overall correlations between the [C ii] luminosity, SFR, and L_IR, but Liang et al. (2024) also reproduce the empirically observed evolution of the [C ii]–SFR or [C ii]–L_IR relation with both redshift and L_IR. Of interest is the finding that dense molecular regions contribute negligibly to [C ii] emission in the FIRE simulations, which would imply that the link between [C ii] luminosity and SFR is truly not through directly tracing the cool star-forming molecular gas (as is the case for CO, for instance). Rather, the [C ii] emission primarily coincides with the neutral H i regions between the ionized and molecular gas, where ionized carbon is appreciably present and the gas density is sufficiently high to encourage cooling through [C ii] line emission.

The ultimate result of Liang et al. (2024) was a physical picture of [C ii] emission where four factors modulate the scaling between [C ii] and SFR: the gas depletion time-scale, the gas-phase metallicity, the gas density, and the fraction of gas in ionized and neutral regions. This work will leverage this proposed formula to move beyond [C ii] as simply a proxy for SFR, and build a new halo model instead relating [C ii] luminosity to the expected metallicity and neutral hydrogen mass in each halo, using this model to build new forecasts for upcoming [C ii] LIM experiments as will take place on the CCAT facility.

2.2 The importance of one-point statistics in summarizing a highly non-Gaussian signal

In traditional cosmological contexts like studies of the cosmic microwave background, we are used to the unreasonable effectiveness of Gaussian random fields as a description of cosmological fields. What reflects this is the common use of the power spectrum as the primary summary statistic of observations of density fluctuations, including LIM.

However, structure formation is very much a non-linear process resulting in non-Gaussian structures. In LIM contexts, Breysse et al. (2017) first explicitly proposed the use of the probability distribution function of observed voxel intensities – i.e. the VID – as a way to capture the highly non-Gaussian information present in line-intensity fluctuations, without relying on higher order statistics whose estimation and covariance pose significant challenges to detection and interpretation. Work continues on understanding the VID and extending one-point statistics in LIM contexts, on both theoretical and numerical fronts (Breysse, Anderson & Berger 2019; Ihle et al. 2019; Sato-Polito & Bernal 2022; Breysse, Chung & Ihle 2023; Chung et al. 2023).

It is, however, the logarithm of the VID, rather than the VID itself, that is most closely related to the information content of the LIM observation. Lee et al. (2024) put forward a relevant insight in this context, using one-point statistics to quantify the non-Gaussian information injected into the cosmic IR background (CIB) by gravitational lensing. In particular, that work put forward quantities related to the relative entropy – or equivalently to the Kullback–Leibler divergence (KLD; Kullback & Leibler 1951) – between the intensity probability distributions of the lensed and unlensed CIB.

It is indeed natural to consider the information content in the VID – i.e. the probability density function (PDF) P(I) – as the relative entropy from a model PDF Q(I) to the actual P(I):

However, Lee et al. (2024) considered it equally interesting to examine the integrand, the pointwise relative entropy (PRE) per intensity bin (dubbed ‘weighted relative entropy’ in that work). We will consider ‘per intensity bin’ (or ‘per log-intensity bin’ in some cases) to be implicit and often refer simply to

[We will often switch between using dS_rel as shorthand for either dS_rel/dI as defined above or dS_rel/d(log I); the usage will be clear based on, e.g., the binning of the discrete histogram approximating the underlying continuous P(I).] Note that many works (e.g. MacKay 2003) define ‘relative entropy’ as the negative of our definition in equation (1), as that is equal to the KLD. However, our definition naturally leads to the PRE from Q to P being the Shannon information for P minus that for the ‘starting’ Q, weighted by P.

Using the PRE – the log-space difference of one distribution against another, weighted by the true distribution – provides a fuller description of the effect of a given variable on the intensity PDF. The integrated quantity of relative entropy will quantify the net statistical distance by which that variable pushes the PDF away from an alternative expectation, but the PRE can provide full templates of how the difference arises at different ranges of intensity. Thus, distinguishing different modes of distributional distortion would disentangle a range of astrophysical and cosmological signatures that can imprint non-Gaussianities in the observable Universe.

Towards the conclusion of this work, rather than explicitly forecasting the detectability of these signatures, we will use the PRE to consider the information content in the VID related to a number of astrophysical model parameters. This in turn will guide future work that will more explicitly probe the detectability of the relevant signatures as well as signatures of primordial non-Gaussianity in the clustering and abundances of [C ii] emitters.

3 METHODS

3.1 Line-luminosity model: a novel halo–[C ii] connection

As we reviewed in Section 2, through the analysis of the relationship between SFR and [C ii] emission in ∼500 simulated galaxies, Liang et al. (2024) derived a linear scaling relationship for the ratio between SFR and the global [C ii] luminosity |$L_\rm{{[}C\,{\small II}{]}{}}$| of a galaxy:

where |$f_\rm{{[}C\,{\small II}{]}{}}$| is the fraction of the total galactic gas mass contained in H i and H ii regions, |$t_{\rm {dep}}\equiv M_{\rm {gas}}/\rm {SFR}$| is the galaxy’s gas depletion time-scale (equal to the total gas mass divided by the SFR), and |$\overline{Z}_{\rm {gas}}$| and |$\overline{n}_{\rm {gas}}$| are the average gas-phase metallicity and gas particle number density, respectively.¹

By the end of this section, we will be able to transform this relation into a model of [C ii] luminosity as a function of halo mass and redshift. We first establish H i mass as a primary intermediary property in Section 3.1.1, before adding dependence on metallicity in Section 3.1.2 and a stochastic step in Section 3.1.3.

3.1.1 L_[C ii] as a function of H i mass

By multiplying both sides of equation (3) by the SFR, we recover a proportionality relationship between |$L_\rm{{[}C\,{\small II}{]}{}}$| and M_gas:

For now, we neglect the dependence on metallicity and gas density, and consider the proportionality between |$L_\rm{{[}C {\small II}{]}{}}$| and |$f_\rm{{[}C {\small II}{]}{}}M_\text{gas}$|⁠. The latter quantity effectively ends up being the gas mass contained in H i and H ii regions. The original motivation for this being the relevant quantity in the work of Liang et al. (2024) is that it is entirely acceptable to ignore the contribution of molecular (H₂) regions to [C ii] emission, as most of the carbon there is in a neutral rather than ionized state. In particular, however, simulations show that it is chiefly the H i regions that dominate the determination of the [C ii] luminosity. This is because, as Liang et al. (2024) discuss, even the [C ii] emission from H ii regions – which often contribute more to the total [C ii] luminosity than H i regions – technically originates from the interface of the H ii gas with the H i-rich regions.

Based on the above considerations, we approximate the H i gas as the main driver of |$f_\rm{{[}C {\small II}{]}{}}M_{\rm {gas}}$|⁠, such that the [C ii] luminosity scales linearly with the H i gas mass M_H i, with some proportionality coefficient |$A_\rm{{[}C {\small II}{]}{}}$|⁠:

Villaescusa-Navarro et al. (2018) examine haloes and galaxies in cosmological simulations to provide a relationship between the total halo mass M_h and M_H i for M_h ∼ 10⁹–|$10^{13}\ {\rm M}_{\odot }$|⁠, featuring an exponential cut-off at lower halo masses and a power law at higher halo masses:

where M_min is the cut-off mass (essentially in units of M_⊙), M_0,H i is the overall normalization, and α_H i is the power-law index at high halo mass. Villaescusa-Navarro et al. (2018) determine values for these parameters at z ∈ {0, 1, 2, 3, 4, 5}. While we are looking into haloes of z ∼ 6–8, the relation evolves little at the highest redshifts where Villaescusa-Navarro et al. (2018) perform their fitting. Based on the z = 5 parameter values, we use α_H i = 0.74, and scale the other parameters down to 70 per cent of their values at z = 5, so that M_0,H i = 1.9 × 10⁹ M_⊙ and M_min = 2.0 × 10¹⁰. This downscaling is an approximate extrapolation of the redshift evolution of these parameters at z ≳ 4.

Recall that our goal is to express the [C ii] luminosity of a halo as a function of its mass M_h. Since we have M_H i(M_h), the missing ingredient is the coefficient |$A_\rm{{[}C {\small II}{]}{}}$| in equation (5). We calibrate this proportionality at an SFR of |$\rm {SFR}^{*}\sim 1\ M_\odot$| yr⁻¹, not uncommon for FIREBox simulated galaxies at z ∼ 6.

Throughout this work, we obtain averages of certain properties as a function of halo mass and redshift based on the model of Behroozi, Wechsler & Conroy (2013a, b), which connects halo mass to stellar mass and SFR through empirical models matching outputs from cosmological simulations to observational constraints. Here, we interpolate across outputs from this model to find the halo mass |$M_{\rm h}^{*}$| corresponding to |$\rm {SFR}^{*}$| at z ∼ 6, and subsequently the corresponding H i mass |$M_{\rm H\, {\small I}}^{*}\equiv M_{\rm H\, {\small I}}(M_{\rm h}^{*})$| based on equation (6). We find |$M_{\rm H\, {\small I}}^{*} \sim 2\times 10^9\ \mathrm{ M}_\odot$|⁠. As per Liang et al. (2024),

meaning that we should have |$L_\rm{{[}C {\small II}{]}{}}(M_{\rm h}^{*})\sim 10^7\ \mathrm{ L}_{\odot }$|⁠. This sets the value of |$A_\rm{{[}C {\small II}{]}{}}=0.005$|⁠, meaning that when considering only the proportionality of [C ii] luminosity with H i mass, the relation between halo mass and [C ii] luminosity is

with M_H i(M_h) given by equation (6).

3.1.2 L_[C ii] as a function of metallicity

The gas-phase metallicity Z determines the amount of carbon in neutral and ionized gas available to emit in [C ii] in the first place, and thus is fundamental to our model. The proportionality that we have derived above effectively assumes a constant global gas-phase metallicity across all haloes, but this is far from realistic as we expect significant evolution of metallicity on average with both halo mass and redshift. We would in fact like to calculate the metallicity Z expected in each halo and use it as an additional variable that modulates the halo–[C ii] connection, so that

It is common to relate the metallicity of galaxies to their stellar masses M_* and SFR, to the point where this relation is considered ‘fundamental’ (at least out to some moderately high redshift) and dubbed the fundamental metallicity relation (FMR) and studied using a range of galaxy selections and parametrizations (Lara-López et al. 2010; Mannucci et al. 2010; Hunt et al. 2012; Cresci, Mannucci & Curti 2019; Curti et al. 2020). We use the model of Behroozi et al. (2013b) again, using their fitting formula for average stellar mass 〈M_*〉(M_h, z) as a function of halo mass and redshift, and interpolating their outputs for 〈SFR〉(M_h, z).

Note in particular that the stellar mass–halo mass (SMHM) relation parametrization is

where M₁ is the characteristic halo mass, ϵ is a fitting factor, and f(x) is defined as

where α is the power-law index at low mass, γ_f is a fitting index, and δ is a fitting factor. Most of the fitting parameters include some redshift dependence, including

with best-fitting values of α₀ = −1.412 and α₁ = 0.731. We refer the reader to Behroozi et al. (2013b) for further details on the other parameters.

To find metallicities, Heintz et al. (2021) use the FMR parametrization of Curti et al. (2020), and we do the same:

where |$M_0(\mathrm{SFR}) = 10^{m_0} \times \mathrm{SFR}^{m_1}$|⁠, with best-fitting parameters |$\tilde{Z}_0$| = 8.779, m₀ = 10.11, m₁ = 0.56, γ = 0.31, and β = 2.1. Here, the metallicity |$\tilde{Z}=12+\log {(\mathrm{O/H})}$| is not in units of solar metallicity. Curti et al. (2020) note that the solar oxygen abundance corresponds to |$\tilde{Z}=8.69$|⁠, meaning that

All that remains is to recalibrate the proportionality between [C ii] luminosity and M_H i, as we did while building the first version of our model. For the halo mass |$M_{\rm h}^{*}$| corresponding to SFR|$^{*}\sim 1\ {\rm M}_{\odot }\,\text{yr}^{-1}$| at z ∼ 6, we should expect a metallicity of Z*/Z_⊙ ∼ 0.3 based on the relations outlined above. One can already see that |$\alpha _\rm{{[}C {\small II}{]}{}}\approx A_\rm{{[}C {\small II}{]}{}}/(Z^{*}/\mathrm{ Z}_\odot)\sim 0.02$|⁠. To be more precise, we have

with M_H i(M_h) still given by equation (6) but the expected metallicity value Z(M_h, z) now given by equations (13) and (14) combined with 〈SFR〉(M_h, z) and 〈M_*〉(M_h, z) based on Behroozi et al. (2013b). Additionally, we assume a minimum halo mass of |$M_{\rm h}=10^{10}\ \mathrm{ M}_\odot$| required for [C ii] emission, with the implication being that we do not expect lower mass haloes to host carbon and/or diffuse gas in sufficient abundance for there to be appreciable cooling through [C ii] emission (cf. e.g. table 2 of Chung et al. 2020 and discussion in their section 3.1).

3.1.3 Adding lognormal scatter to metallicity

A halo model cannot capture all of the baryonic physics interior to each halo driving the [C ii] emission, but we can approximate such physics as a stochastic process affecting the [C ii] luminosity. In this case, we assume that the carbon abundances are most likely to be driven away from the expected value, and randomly scatter the metallicity Z(M_h, z) by a lognormal distribution with standard deviation σ_Z = 0.4 (in units of dex), which in turn will scatter the [C ii] luminosity with the same distribution since |$L_\rm{{[}C {\small II}{]}{}} \propto Z$|⁠. We choose the fiducial value for σ_Z based on previous work by Vizgan et al. (2022), which showed that the log-linear correlation between [C ii] luminosity and H i mass exhibits a log-space root-mean-square scatter of 0.39 dex.

To add scatter to metallicity, we compute a positive scaling factor ζ_Z randomly drawn from a lognormal distribution with the following PDF:

where σ = σ_Z ln(10) and μ = −σ²/2, with μ chosen to preserve the linear mean value of Z for a given halo mass and redshift. Then, from the expected metallicity value Z(M_h, z) given halo mass and redshift, we may obtain a randomized metallicity simply by multiplying by ζ_Z ∼ P(ζ_Z; σ_Z). Therefore, when ‘painting’ [C ii] luminosities on to dark matter haloes in a cosmological simulation, we calculate |$L_\rm{{[}C {\small II}{]}{}}$| from the H i mass and the randomized metallicity ζ_ZZ found for each halo:

3.2 Noise model: fundamental sensitivities of Stage 1 and 2 LIM experiments

Along with the signal, we generate forecasts of observational experiments. To do so, we convolve the pure signal by a Gaussian approximation to the main beam, and add Gaussian noise to the result. This work only aims to examine the fundamental sensitivity that [C ii] LIM will be able to achieve in the near- to medium-term future, and so we omit considerations of either non-Gaussian sky noise (which will be the subject of significant near-future work from the first generation of [C ii] experiments) or interloper emission (rejection of which is already the subject of considerable theoretical efforts – see e.g. Breysse, Kovetz & Kamionkowski 2015; Cheng et al. 2016; Lidz & Taylor 2016; Sun et al. 2018; Cheng, Chang & Bock 2020; Chen & Pullen 2022).

To add Gaussian noise to our signal, we obtain the sensitivity per map voxel. We refer the reader to other works like appendix B of Chung et al. (2020) and references therein for further details on noise-related quantities; but for mm-wave observations, the convention is to express sensitivities in terms of the noise-equivalent intensity (NEI):

where |$\rm {NEFD}$| is the noise-equivalent flux density (the noise level per instrumental beam) and Ω_beam is the solid angle in the sky covered by the beam. Approximating the beam as having a Gaussian profile with full width at half-maximum of θ_FWHM, the beam solid angle is given by

Given the NEI, the actual noise level per voxel is given by

where N_feeds is the number of spectroscopic feeds, and t_vox is the average observing time per volume element (or voxel) of the survey volume. For a spectrometer with simultaneous coverage of its observing band, this would simply be the observing time per sky pixel, but a scanning spectrometer will incur a penalty based on the number of steps required to cover the entire observing band.

We will centre our observational forecasts on the CCAT facility, which will deliver the largest field of view out of all currently funded EoR [C ii] experimental projects. Based on the most recent specifications at ν ∼ 280 GHz (CCAT-Prime Collaboration 2023), the two EoR-Spec instrument modules proposed for Prime-Cam, a seven module camera, will target an angular resolution of θ_FWHM = 48 arcsec, a frequency resolution of δ_ν = 2.8 GHz (or ν/δ_ν ≈ 100), and NEFD values of 64, 81, and 120 mJy s^1/2 per beam in first-, second-, and third-quartile weather (respectively) at a source elevation of 45°.

Like the CCAT-Prime Collaboration (2023) work, we assume that the survey spans two fields of 4 deg², and that EoR-Spec requires 42 steps to cover the full observing band. While the original CCAT-Prime Collaboration (2023) forecasts assume 6912 beams and 80 per cent yield (which would imply ≈5530 active detectors) and a weighted average NEI across the varied weather conditions, we assume 5040 active detectors on the EoR-Spec focal plane and an NEFD per beam of 72.5 mJy s^1/2 (a simple arithmetic mean of the first- and second-quartile NEFD values). These two differences end up mostly cancelling each other, as we assume slightly fewer active detectors but slightly better NEFD per beam, and we end up with sensitivity figures within 10 per cent of those obtained by CCAT-Prime Collaboration (2023).

We consider two experimental scenarios as follows:

• The first scenario is a CCAT-like experiment, specifically one conforming to specifications of the proposed baseline CCAT deep spectroscopic survey (CCAT-DSS) of CCAT-Prime Collaboration (2023) with 2000 h of observing time per 4-deg² field. Assuming a pixelization of the survey volume with a pixel size of |$\theta _\text{FWHM}/\sqrt{8\ln {2}}=20 \ \mathrm{ arcsec}$|⁠, this gives 58 s within each pencil beam of 20 arcsec × 20 arcsec, but t_vox ≈ 1.4 s after accounting for the scanning spectrometer penalty (which is equivalent to assuming t_vox = 58 s but N_feeds = 5040/42 = 120).

• The second scenario is a future Stage 2 successor LIM experiment, with 4.8 × 10⁶ spectrometer-hours per field (versus the 2000 × 120 = 2.4 × 10⁵ spectrometer-hours achieved with the baseline CCAT survey), or t_voxN_feeds = 38 h. Advancements in on-chip spectrometers and readout technologies should relieve some of the efficiency limitations of the current generation of spectrometers and enable 10⁶–10⁷ spectrometer-hour surveys to deploy within the next decade (Karkare et al. 2022). The CCAT facility is fundamentally capable of hosting a Stage 2 mm-wave LIM experiment of this calibre, as Prime-Cam occupies only the central |$\sim \!\! 40{{\ \rm per\, cent}}$| of the total field of view delivered by the telescope, and future CCAT instruments could deploy many more spectrometer modules. We do not necessarily assume that these modules will be scanning spectrometers like the EoR-Spec modules, but assume that the same NEFD is fundamentally achievable since the values we assume are already weather-limited.

3.3 Generating cosmological line-intensity mock data cubes

Having pulled out our halo model for [C ii] from the FIRE simulation analysis of Liang et al. (2024), and defined experimental parameters for near-future [C ii] surveys, it only remains to simulate observations of the [C ii] signal. We use a suite of approximate cosmological simulations to generate catalogues of haloes spanning the survey volume, and then generate 3D cubes of simulated observations by ‘painting’ luminosities on to these haloes.

3.3.1 Cosmology to haloes: the peak-patch method

The peak-patch (or mass peak-patch) method is a model for identifying haloes in the initial dark matter fluctuation field. Moving beyond the original peak theory of Bardeen et al. (1986) and the excursion set method of Bond et al. (1991), the peak-patch method is a convolutional method, identifying mass peaks in the initial (Lagrangian) density field after processing through a hierarchy of smoothing filters at different scales, and finding the final properties of peaks as collapsed haloes by solving for homogeneous ellipsoidal collapse and using second-order Lagrangian perturbation theory to resolve flow dynamics. First expounded and validated in the work of Bond & Myers (1996a, b, c), the peak-patch method has recently returned to work following re-implementation and re-validation by Stein, Alvarez & Bond (2019). The peak-patch method generates large ensembles of independent cosmological realizations with much greater ease than conventional N-body methods, making it suitable for signal forecasts as well as covariance estimation for novel statistics. This is particularly true in the context of LIM (Ihle et al. 2019; Chung et al. 2023), which demands large light-cone runs resolving relatively low-mass haloes, since the peak-patch method can intrinsically evolve the initial mass peaks to a continuous range of final redshifts instead of having to stitch together snapshots at discrete time-steps.

We use a suite of peak-patch simulations generated to probe the possible use of novel statistics to apply to CO LIM at EoR redshifts and previously detailed in a number of papers, including Chung et al. (2023). These simulations were designed to push to lower masses by a factor of several compared to the peak-patch simulations generated at z ∼ 3 for Ihle et al. (2019), and so have a size of (960 Mpc)³ and a resolution of 5640³ cells. In principle, this should correspond to a minimum resolvable halo mass of |$6.4\times 10^9\ \mathrm{ M}_\odot$|⁠; after adjusting masses to match abundances predicted by the halo mass function (HMF) of Tinker et al. (2008) with empirical high-redshift corrections from Behroozi et al. (2013b), we find that our simulations should be complete down to |$M_{\rm h}\approx 1.2\times 10^{10}\ \mathrm{ M}_\odot$|⁠, well matched to our assumed minimum emitter halo mass of |$10^{10}\ \mathrm{ M}_\odot$|⁠. [Note that such adjustments are not strictly necessary particularly at lower z than those considered here; Stein et al. (2019) showed that peak-patch simulations agree out of the box with the Tinker et al. (2008) HMF within 10 per cent for |$M_\text{h}\gtrsim 10^{11}\ \mathrm{ M}_\odot$| and z ≤ 2.]

The simulations are light-cone runs with the centre of the box placed 8668 Mpc away from the observer; although the minimum requirement for CO LIM simulations was to cover z = 5.8–7.9, the light-cones actually cover a slightly wider interval of z = 5.6–8.2. The box size then corresponds to a sky area of 6° × 6° at the highest redshifts accessed. As we only need to simulate patches of 2° × 2°, it is theoretically possible to slice each simulation into 9 semi-independent sub-volumes for a total of 2430 light-cones instead of 270. However, this work does not focus on aspects like covariances of statistics, and 270 light-cones is already a sufficient ensemble to reliably identify average expectations and variances. We instead favour using wholly independent volumes and use the central 4 deg² of each light-cone.

3.3.2 Haloes to maps: mocking [C ii] clustering

To consider how the [C ii] line-intensity field responds to the presence of our simulated haloes, we use modified versions of limlam_mocker² and lim³ codes that interface with each other to paint the halo–[C ii] ‘response function’ of equation (17) on to our peak-patch light-cones, as well as appropriately handle convolution of the beam profile and addition of noise. The code is more generally applicable to other atomic and molecular line species, with the only prerequisite being an appropriate halo model for line emission.

The code uses the peak-patch simulation outputs – which catalogue each halo’s mass, Eulerian positions, comoving distances, and redshift – as well as a function to relate halo properties to line luminosities [in our case, the |$L_\rm{{[}C {\small II}{]}{}}(M_{\rm h},z)$| function as described in Section 3.1] and the observational parameters described in Section 3.2, including the survey extent in angular and frequency coordinates, the NEFD, θ_FWHM, N_feeds, δ_ν, and the total survey time (from which the code derives t_vox).

For all haloes in the simulation that fall within the specified observational volume, we compute the [C ii] luminosity of each halo, but any mass-luminosity function, for any atomic or molecular line, can be prescribed. We then bin all haloes and thus their luminosities into discrete voxels based on the pixelization of the sky and the frequency channelization. The conversion of the total luminosity |$L_{\rm{{[}C {\small II}{]}{}}, \rm {vox}}$| in each voxel to the voxel intensity |$I_\rm{{[}C {\small II}{]}{}}$| is straightforward:

where V_vox is the comoving volume of each voxel, and the prefactor c/[4πν_restH(z)] depends on the speed of light c, the rest-frame line emission frequency ν_rest, and the Hubble parameter H(z) at the emission redshift.

From our peak-patch light-cone simulations, we produce [C ii] intensity cubes of size θ_S = 2° along each angular coordinate – i.e. right ascension (RA) and declination – and ν_obs = 270 ± 20 GHz in frequency space. Based on the CCAT beamwidth of θ_FWHM = 48 arcsec or (1/75)°, we require a pixel size equal to the scale of the corresponding Gaussian profile, |$\theta _\text{FWHM}/\sqrt{8\ln {2}}$|⁠, yielding |$(\theta _\mathrm{ S}/\theta _\text{FWHM})\sqrt{8 \ln {2}}=353$| pixels on each transverse side. Along the line of sight in frequency space, we have a bandwidth of 40 GHz and a spectral resolution ν/δ_ν ≈ 100 or δ_ν = 2.8 GHz, which makes for 15 frequency channels.

The observing frequencies encompass [C ii] emission at z = 5.6–6.6; at the central frequency, we obtain V_vox ≈ 22 Mpc³ and c/[4πν_restH(z)] ≈ 7.2 × 10⁻⁴ Jy sr⁻¹ Mpc³ |$\mathrm{ L}_\odot ^{-1}$|⁠. For limited purposes, we also produce intensity cubes at other redshifts, including at ν_obs = (226 ± 14) GHz using otherwise similar parameters.

We consider the raw [C ii] cube in some cases, but in most cases we will consider the cube after convolution with the angular Gaussian profile approximating the CCAT beam, and with the addition of Gaussian noise as described in Section 3.2.

3.4 Sanity checks: luminosity functions and [C ii] deficits

Before moving on to the primary results of our simulations, we show the results of a couple of sanity checks of the model, first against other z ∼ 6 [C ii] luminosity functions in the literature and then against the so-called [C ii] deficit. The latter is not so prominent in the range of luminosities primarily relevant to our work, but there is none the less weak downward evolution of |$L_\rm{{[}C {\small II}{]}{}}/L_\text{IR}$| with increasing L_IR.

As our peak-patch light-cones have halo masses adjusted to match the Tinker et al. (2008)–Behroozi et al. (2013b) form of the HMF, it is sufficient to use that HMF model to generate randomly drawn masses given some large cosmological volume, convert these masses to [C ii] luminosities (with randomly drawn metallicity values), and obtain the resulting luminosity function. We show this in Fig. 2 alongside selected prior simulations and observations.

Our model output compares favourably to the prediction of Popping et al. (2016), which uses a semi-analytical model to predict line-luminosity functions at z = 0–6; that work noted its difficulty in reproducing bright CO and [C ii] emitters at high redshifts, a difficulty that we do not share to the same extent. We also predict similar abundances of |$L_\rm{{[}C {\small II}{]}{}}\sim 10^8$|–|$10^9\ \mathrm{ L}_\odot$| objects as the model of Lagache et al. (2018), although considerably fewer bright objects and somewhat more faint objects. [The comparison of Béthermin et al. (2022) in turn shows that the model of Chung et al. (2020) predicts a luminosity function largely similar to that of Lagache et al. (2018), as does their own model when using the Lagache et al. (2018) SFR–[C ii] relation.]

We also show a similar level of consistency (i.e. within an order of magnitude) with observational constraints, which is remarkable given that we have not explicitly fitted to them in great detail beyond the specific pivot point of |$L_\rm{{[}C {\small II}{]}{}}\sim 10^7\ \mathrm{ L}_\odot$| (which is in fact quite far from where most observational constraints lie).

• The ALMA Spectroscopic Survey Large Program (ASPECS LP) survey was unable to go beyond upper limits on the luminosity function at z ∼ 6 (Decarli et al. 2020); our simulated source abundances lie safely below the upper limits.

• Our simulated luminosity function lies systematically below the measurements of the ALMA Large Program to INvestigate [C ii] at Early times (ALPINE) survey (Yan et al. 2020), an ALMA survey of ultraviolet (UV)-selected galaxies. The target selection of ALPINE may be a source of bias, but the discrepancy against our model is certainly not as radical as against that of Popping et al. (2016) (which posited a similar explanation for inconsistency with a lower limit at z ∼ 4).

• The ALMA Reionization Era Bright Emission Line Survey (REBELS) has also surveyed [C ii] emitters at z ∼ 7, but has not publicly released explicit luminosity function constraints; in its place, we compare the fitting of Padmanabhan (2023) to the work in preparation against our model, and find consistency within a factor of 2–3 for |$L_\rm{{[}C {\small II}{]}{}}\gtrsim 10^8\ \mathrm{ L}_\odot$|⁠. For comparison, publicly discussed REBELS [C ii] emitters (e.g. in Ferrara et al. 2022) have luminosities of |$10^{8.6}\ \mathrm{ L}_\odot$|⁠.

We also compare against constraints on the cumulative number density of [C ii] emitters above certain detection thresholds.

• An alternative ASPECS LP analysis by Uzgil et al. (2021) derived upper limits on source abundances down to luminosities of |$\sim\!\! 10^8\ \mathrm{ L}_\odot$|⁠; we show their upper limit for sources at z = 6–7 and demonstrate consistency with this limit. Although the original work of Uzgil et al. (2021) shows a number of other limits, the qualitative comparison against their other limits is the same.

• In addition, two surveys of lensing clusters have each yielded a single detection of a [C ii] emitter at z ∼ 6. The ALMA Lensing Cluster Survey (ALCS) (Fujimoto et al. 2021) mapped strongly lensed regions across 33 different clusters, covering a total of 88 arcmin²; the Deep UNCOVER-ALMA Legacy High-z (DUALZ) (Fujimoto et al. 2023) specifically targets Abell 2744, in both a wide 24-arcmin² survey and a deep 4-arcmin² survey.As in the comparison against ALPINE data, our predictions lie systematically on the lower side of these limits. However, we caution for both ALCS and DUALZ that these constraints not only rely on a single detection in small effective survey areas of ∼10⁰–10² arcmin², but also rely strongly on the mass model of the lensing cluster, both in the determination of the intrinsic luminosities that were or should have been detectable, and in the determination of the effective comoving volume surveyed.

For the |$(L_\rm{{[}C {\small II}{]}{}}/L_\text{IR})$|–L_IR relation, we look to the simulated FIRE galaxies as a comparison point, rather than observational trends.⁴ Real-life resolved sources at z > 5 are mostly exceptionally bright sources with |$L_\text{IR}\gtrsim 10^{12}\ \mathrm{ L}_\odot$| (and typically |$L_\rm{{[}C {\small II}{]}{}}\gtrsim 10^9\ \mathrm{ L}_\odot$|⁠), whereas it is almost entirely the complementary population that is most abundant and will dominate the large-scale clustering component of the [C ii] LIM signal.

To obtain the IR luminosity for our haloes, we use the same IR–SFR conversion as Liang et al. (2024):

This is the relation of Kennicutt (1998), recalibrated for a Kroupa (2002) initial mass function.

Recall that we formulated our final model in steps, and initially considered in Section 3.1.1 a model with only explicit dependence on H i mass (implying either fixed metallicity across all haloes or lack of any metallicity dependence). We show the |$(L_\rm{{[}C {\small II}{]}{}}/L_\text{IR})$|–L_IR relation for this model in the upper panel of Fig. 3. At high L_IR, Liang et al. (2024) note that t_dep ∼ M_gas/SFR drives the so-called [C ii] deficit as there is less gas available to cool in [C ii] relative to the amount of star-forming activity. In the regime of |$L_\text{IR}\lesssim 10^{12}\ \mathrm{ L}_\odot$|⁠, however, the predicted evolution in |$L_\rm{{[}C {\small II}{]}{}}/L_\text{IR}$| with L_IR is far too steep when attempting to explain [C ii] emission with gas mass alone, and in particular overpredicts [C ii] emission in low-mass objects relative to the FIRE simulations.

The lower panel of Fig. 3 shows the |$(L_\rm{{[}C {\small II}{]}{}}/L_\text{IR})$|–L_IR relation for our final model, using both H i mass and metallicity to inform the [C ii] luminosity. As Liang et al. (2024) note, metallicity is clearly key in explaining a different kind of [C ii] deficit, which is the decline in |$L_\rm{{[}C {\small II}{]}{}}/L_\text{IR}$| with increasing redshift at fixed L_IR. In general, less IR-luminous galaxies are less metallic and less abundant in carbon, which counteracts the increased relative availability of diffuse gas. The softened evolution of |$L_\rm{{[}C {\small II}{]}{}}/L_\text{IR}$| with L_IR is far more consistent with the simulated FIRE galaxy population. While the deficit is softened, we still predict a strong break in the |$(L_\rm{{[}C {\small II}{]}{}}/L_\text{IR})$|–L_IR relation at |$L_\text{IR}\gtrsim 10^{10}\ \mathrm{ L}_\odot$| away from the naïve expectation of constant |$L_\rm{{[}C {\small II}{]}{}}/L_\text{IR}$|⁠.

It is possible that at the very highest IR luminosities shown in Fig. 3, the dependence on metallicity overcompensates against the declining abundance of gas in general, and leads to an overprediction of [C ii]. Inversely, metallicity dependence suppresses our z ≈ 8 predictions at low and high L_IR relative to FIRE predictions. A crude parametrization of |$f_\rm{{[}C {\small II}{]}{}}$| or a re-evaluation of the FMR at z ≳ 6 would allow the model to capture the behaviour of the brightest and highest redshift sources, but this task is beyond the scope of this work as |$L_\text{IR}\gtrsim 10^{12}\ \mathrm{ L}_\odot$| sources will be too rare to significantly affect the [C ii] LIM signal and as a precise replication of the FIRE outputs is not the primary focus of this work.

4 SIMULATION RESULTS

4.1 Line-intensity maps and power spectra

As previously outlined, we generate 270 light-cones of simulated [C ii] emission at z ≳ 6. For most of this work, we will be focusing in particular on the part of the signal originating at z = 5.6–6.6 [corresponding to ν_obs ∈ (250, 290) GHz]. We plot the intensity map for a wider interval spanning z ∈ (5.8, 7.6) in Fig. 4, which shows the [C ii] emission tracing the large-scale clustering of our simulated dark matter haloes, but with an intensity that falls off rather rapidly beyond z ∼ 7.

No actual experiment will perfectly recover the signal as shown in Fig. 4. As previously stated, the primary focus of LIM forecasts is the power spectrum, encapsulating the variance of the [C ii] line-intensity fluctuations at different scales, which should be detectable above noise. The power spectrum also lets us calculate a Wiener filter given the expected power spectrum and the actual noisy observation, thus constructing a minimum-variance linear estimator of the original [C ii] line-intensity map. We consider implementation details for the power spectrum and Wiener filter calculations in Appendix  A, and focus here on the results.

Fig. 5 shows the average power spectrum across all 270 light-cones, compared to both previous models (Silva et al. 2015; Serra et al. 2016; Dumitru et al. 2019; Chung et al. 2020; Kannan et al. 2022; Karoumpis et al. 2022; Sun et al. 2023) and sensitivities expected from future [C ii] LIM experiments (see Chung et al. 2020 for calculation details). At large scales, our predictions are similar to those of previous models if relatively pessimistic. However, our model is decidedly more pessimistic at small scales, indicating that we tend to predict significantly less shot noise compared to previous models, which is to say we do not predict a substantial population of rare bright [C ii] emitters. In this sense, our predictions tie closest to those of Kannan et al. (2022) and Sun et al. (2023).

Greater stochasticity in [C ii] luminosity than the fiducial expectation could account for some of this difference, but not easily. Fig. 6 shows the [C ii] power spectrum for different amounts of lognormal scatter in metallicity ranging from 0.0 to 1.0 (in units of dex), where the fiducial value is σ_Z = 0.4 dex. Higher scatter in metallicity introduces more random small-scale fluctuations of the signal, and as a result, the clustering component of the power spectrum becomes subdominant to shot noise at lower k for higher σ_Z. For an arbitrarily high level of scatter, the signal becomes entirely dominated by shot noise on all scales. However, the values required are extreme, with σ_Z ≳ 0.7 dex required to achieve a level of shot noise relative to clustering similar to other model predictions. (As for why such values would be extreme relative to the fiducial σ_Z = 0.4 dex, consider that a change in the standard deviation of a Gaussian variable from σ = 0.4 to σ = 0.8 would decrease the correlation coefficient against another variable by a factor of 2, all other covariances being held equal.)

In fact, the low level of shot noise arises in large part as a fundamental feature of our halo model: by its very motivation our model prescribes non-linearly suppressed [C ii] emission from more massive objects relative to their SFR, mimicking the so-called [C ii] deficit. By contrast, the other models use either a linear relation or mildly non-linear power law to obtain [C ii] luminosity from either SFR or IR luminosity. A key exception is the prediction of Sun et al. (2023), which uses the limfast code (Mas-Ribas et al. 2023) to generate a model of [C ii] emission motivated by similar physical variables as in our model, albeit simulated at the molecular cloud level rather than at the dark matter halo level. Another exception exists in one of the models of Karoumpis et al. (2022), which works from the results of Vallini et al. (2015) to account for the effect of metallicity on the [C ii] luminosity. However, metallicity accounts for the redshift evolution of the [C ii]–SFR relation given the metal-poor nature of early galaxies, while the primary driver of the [C ii] deficit at high luminosities is the gas mass. We also see the reluctance of our model to prescribe very bright objects in the luminosity function shown in Fig. 2, especially in comparison to the model of Lagache et al. (2018). All this suggests that the shape of the power spectrum and the relative balance of the clustering and shot noise components may be a key discriminator between competing physical pictures of [C ii] emission. Although such pictures must include a wide range of processes that impact the halo–[C ii] connection at both low and high masses, the shape of the power spectrum may yet test different explanations of the [C ii] deficit.

Prospects for detecting the power spectrum to this end are promising. When summed in quadrature across all |$k\lesssim 1$| Mpc⁻¹, the total signal-to-noise ratio for the z ∼ 6 [C ii] power spectrum is ≈5σ for the baseline CCAT-DSS experiment. Although a detection at z ∼ 7.5 would be out of reach for CCAT-DSS, our Stage 2 LIM experimental parameters would be sufficient to achieve an ≈11σ overall detection of the z ∼ 7.5 [C ii] power spectrum, as well as an improved ≈48σ detection at z ∼ 6. The improvement would allow measurement of the shape of the [C ii] power spectrum beyond a simple initial detection, probing the contribution of populations of faint and bright galaxies to the total signal.

A Wiener-filtered line-intensity map provides a more intuitive illustration of what this detectability of the power spectrum means in terms of sensitivity to [C ii] intensity fluctuations at different scales. We process the raw line intensities from an example realization, as shown in the upper panel of Fig. 7, through a Wiener filter given the Stage 2 forecast, ending up with the result in the lower panel of the same figure demonstrating recovery of prominent large-scale features (again, see Appendix  A for more information on the Wiener filter processing).

4.2 Stacking [C ii] on optical galaxies

We perform a mock stacking analysis to consider recoverability of the average [C ii] signal associated with locations of individual objects. Specifically, we assume that the haloes in our simulations contain not only [C ii]-emitting gas, but also Lyman-break galaxy (LBG) populations that can be detected through photometric drop-out (as the Lyman break in the continuum flux falls between two photometric filters), and simulate stacks on a hypothetical narrow-band LBG survey that can access thousands of objects over the sky area of a single 2-deg² [C ii] LIM field.

We assume drop-out identification of LBGs at z ∈ (5.9, 6.1), and assign LBG counts to each halo in this redshift range. Harikane et al. (2016) consider a halo occupation distribution model for a range of LBG populations detected in Hubble deep imagery and large-area Subaru/Hyper Suprime-Cam (HSC) data, including a z ≈ 5.9 LBG population detected as i-drop-outs with threshold rest-frame UV aperture magnitude of 28.4, or absolute UV magnitude of M_UV < −19.1. Although this limiting magnitude is mostly driven by the Hubble data, which only span hundreds of square arcminutes, we assume that future ultra-deep photometric data will be capable of reaching this limit over areas on the scale of square degrees. This is not unreasonable based on the original HSC Subaru Strategic Program specifications (Aihara et al. 2018), which specified a depth of i ≃ 28 for the HSC-UltraDeep (3.5 deg²) layer; as of the third data release (Aihara et al. 2022), the HSC data had reached a depth of i ≃ 26.9 in the Deep layer, with the UltraDeep layer being ≃0.8 mag deeper.

We briefly recap the halo occupation parametrization of Harikane et al. (2016) parametrization here with minor changes to notation. The average number of LBGs in each halo is given as an assumed star formation duty cycle f_duty = 0.6 times the number of both central and satellite LBGs when active:

The average central LBG count 〈N_cen〉 is given by a sigmoid that reaches 0.5 at some characteristic mass M_min,LBG and has some width σ_log M:

where for the z ∼ 6 LBG population, Harikane et al. (2016) obtain log(M_min/M_⊙) = 11.03 as the best-fitting value. The number of satellite LBGs is then this scaled by a power law:

Harikane et al. (2016) use results of prior empirical studies of halo occupation (Kravtsov et al. 2004; Zheng et al. 2005; Conroy, Wechsler & Kravtsov 2006; Martinez-Manso et al. 2015) to fix a number of parameters. This includes σ_log M = 0.2, α_LBG = 1, and

as well as

After obtaining 〈N〉(M_h), we make a Poissonian draw to determine the ‘actual’ number of LBGs hosted in each halo. This results in a typical count of (4.3 ± 0.4) × 10³ LBGs at z ∈ (5.9, 6.1) in each 2-deg² light-cone. The corresponding comoving number density of ∼5 × 10⁻⁴ Mpc⁻³ is lower than that found by Harikane et al. (2016) in data, but consistent with our own recalculation at z = 5.9 using the HMF and the prescription for 〈N〉(M_h), suggesting that our results are likely consistent with their fitting to the data (modulo cosmology differences).

We generate the actual mock stack on LBGs by cutting out a thumbnail from the [C ii] map centred on the host halo location, 50 pixels on each side (≈17 arcmin × 17 arcmin), and then averaging the thumbnails weighted by the LBG count of all host haloes.

We first show in Fig. 8 the radial profile of the intensity of the stack away from the centre, as recovered from both noiseless simulations and realistic forecasts for [C ii] observations with CCAT-DSS and the Stage 2 concept. The profile appears to have two components, namely a more compact component localized near the stack centre corresponding to the [C ii] emission from the host halo, and a more diffuse component extending out to larger angular distances corresponding to emission from the surrounding large-scale structure. We find a slight upward bias in the CCAT-DSS forecast for the second component at the level of 10 per cent, but otherwise the fundamental sensitivity of CCAT-DSS already allows for strong recovery of the radial profile of [C ii] emission around z ∼ 6 LBGs. With Stage 2 sensitivities, the variance in the recovered profile is comparable to the corresponding variance for the noiseless stack, suggesting that instrumental noise in the [C ii] observation is not the dominant source of uncertainty in the radial profile in that scenario.

We consider a phenomenological model of the radial profile as the sum of Gaussian and power-law profiles, with r here being the radial distance in pixels (recalling that 1 pixel covers 20 arcsec along each angular dimension):

Table 1 shows the recovered values of all five free parameters {I_D, r₀, δ_D, I_G, ϱ} from both the noiseless stack and the average expected stack from the CCAT-DSS and Stage 2 LIM forecasts, given respective errors. Note that while we naïvely expect ϱ = 2 pixel² based on the beam profile being a Gaussian with scale of 1 pixel, the stacking procedure effectively convolves this with the pixel window function during the stack to widen the actual profile. Otherwise, the recovered parameter values are generally consistent with what we see in Fig. 8, with the slight bias in the CCAT forecast for the large-scale component reflected in the lower recovered value of δ_D.

For illustrative purposes, we also show a single realization of the stacking exercise in Fig. 9, including radial profiles recovered from that particular realization re-projected into a thumbnail. Note, however, that the Stage 2 LIM sensitivities allow very strong recovery of the [C ii] signal correlated with LBG locations without the need to isotropize the signal by averaging in radial bins. In fact, although recovery of the radial profile is already very good with CCAT-DSS sensitivities, an isotropized analysis even with Stage 2 sensitivities will miss out on the finer details of clustering of the [C ii] signal, as shown by the structure evident in the difference map between the recovered isotropized thumbnail and the true noiseless stack. Future work should explore the possibility of improving the information content of LIM stacks through anisotropic techniques like oriented stacking (see e.g. Lokken et al. 2022, 2023).

4.3 Relative entropy in one-point statistics

Returning to summary statistics of the [C ii] observation on its own, we now consider the ability of one-point statistics to distinguish the non-Gaussian signature of different parameters.

Recall from Section 2 our definition of PRE per intensity bin:

Consider the ‘true’ P(I) to be our fiducial [C ii] model, and consider the alternate Q(I) to be a version of the model with a change to some parameter λ. The PRE due to finding an observation with the fiducial value λ₀ of the parameter instead of the alternate value of λ₀ + Δλ is

We may consider how this PRE varies with λ (or equivalently its displacement Δλ away from λ₀). Locally, this is simply given by the derivative of the PRE:

or, taking Δλ → 0 while noting that Q(I; λ = λ₀) = P(I),

Therefore, the differential PRE incurred by a change in a given model parameter is simply the derivative of the log probability with respect to that parameter, weighted by the fiducial probability.

We numerically estimate the PRE differential with respect to a number of parameters that we encounter in Section 3.1, namely

• |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠, the overall normalization of |$L_\rm{{[}C {\small II}{]}{}}(M_{\rm h},z)$| in equation (17), with fiducial value of 0.024;

• σ_Z, the metallicity lognormal scatter in units of dex, with fiducial value of 0.4 dex;

• α_H i, the power-law index of the M_H i(M_h) relation, with fiducial value of 0.74;

• log M_min, the exponential cut-off mass scale in the M_H i(M_h) relation, with fiducial value of 10.3;

• α₀, controlling the power-law index of the SMHM relation of Behroozi et al. (2013b), with fiducial value of −1.412; and

• γ, a power-law index in the FMR parametrization of Curti et al. (2020), with fiducial value of 0.31.

First, we consider the PRE in the noiseless signal simulation, before the introduction of either the angular beam or Gaussian noise. The upper panel of Fig. 10 shows the resulting PRE derivatives (per log-intensity bin) with respect to the six chosen parameters. Due to the significantly smaller fiducial value of |$\alpha _\rm{{[}C {\small II}{]}{}}$| relative to those of the other parameters, we appropriately rescale the PRE derivatives against the other parameters for easy comparison.

Since the PRE derivative is of the same sign as the derivative of the log-VID, the effect of increasing a given parameter λ is to displace voxel counts towards ‘higher entropy’ from bins of lower d(dS_rel)/dλ to bins of higher d(dS_rel)/dλ. For instance, by increasing |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠, we uniformly increase the brightness of the entire [C ii] intensity signal, thus pushing voxels from lower to higher intensities. For our fiducial model, this displacement appears to happen across a cross-over point of |$\sim\!\! 10^2$| Jy sr⁻¹.

Unfortunately, the effect of many parameters is largely similar to that of |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠, and their PRE derivatives appear mostly the same as |$\mathrm{ d}(\mathrm{ d}S_\text{rel})/\mathrm{ d}\alpha _\rm{{[}C {\small II}{]}{}}$| after a rescaling and/or sign flip. Decreasing any one of M_min, α₀, or γ increases the brightness of the signal in a way that cannot be distinguished from increasing |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠, at least not through the informational divergence in the one-point statistics of the map. In Fig. 10, we have specifically chosen rescalings that accentuate this degeneracy. (Appendix  B shows the same derivatives with less drastic rescalings and against the natural log of the parameter value; the qualitative conclusions remain the same.)

We do, however, find some parameters whose effect appears distinct from that of |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠. For example, increasing α_H i makes the M_H i(M_h) relation steeper about the pivot point of M_h/M_⊙ = M_min, thus increasing H i mass for haloes with M_h/M_⊙ ≳ M_min but decreasing it for haloes with M_h/M_⊙ ≲ M_min. This results in voxels with middling intensities becoming either brighter or dimmer, depending on the haloes they contain. Similarly, increasing σ_Z can either increase or decrease the metallicity, pushing voxels to both brighter and dimmer intensities.⁵ The characteristic intensity scales and effect sizes associated with these variables are different, meaning that d(dS_rel)/dα_H i and d(dS_rel)/dσ_Z appear quite distinct from each other as well as from |$\mathrm{ d}(\mathrm{ d}S_\text{rel})/\mathrm{ d}\alpha _\rm{{[}C {\small II}{]}{}}$|⁠.

However, much of this becomes insignificant against the effect of introducing the angular beam and Gaussian noise to the observation, which results in the PRE derivatives (per linear intensity bin) shown in the lower panel of Fig. 10. The PRE derivative with respect to α_H i is now indistinguishable in shape from that with respect to |$\alpha _\rm{{[}C {\small II}{]}{}}$|⁠, while the PRE derivative with respect to σ_Z still appears distinct but is significantly reduced in amplitude, suggesting its significantly reduced informational impact on the VID.

4.3.1 Discussion point 1: observational distortions

It is clear that both smoothing of the signal by the beam and contamination from thermal noise cause significant loss and distortions in the informational sensitivity of the VID to our model parameters. To better understand how this occurs, we further vary the beamwidth θ_FWHM and noise scale σ_n, and examine the resulting changes in the PRE derivative.

We summarize our findings in Fig. 11. First, we consider the effect of the beam between the outer left panel and the inner left column of panels. Broadly speaking, the effect of convolving the field with a beam of finite size is a coarse pooling. This pooling dilutes the intensity in the brightest voxels across its neighbours, suppressing the heavy right tail that would originally be present in the VID and compressing the effect of each parameter to a narrower range of intensities. The effect is more complex at lower intensities, and so the distortion of the PRE derivative ultimately depends at some level on the clustering of the signal and the effect thereon of each model parameter.

Next, moving to the inner right column of panels of Fig. 11, we consider the effect of introducing Gaussian noise, considering a few intermediate scenarios between the noiseless case and the Stage 2 forecast with |$\sigma _{\rm n}=3.2\times 10^3\,$| Jy sr⁻¹. Inevitably, the PDF of the total intensity is the convolution of the individual signal and noise intensity PDFs, so introducing Gaussian noise means smoothing out the signal PDF by a Gaussian profile of size σ_n. Given that we were plotting the PRE per log-intensity bin across a range of 1–10⁴ Jy sr⁻¹, the apparent result of convolution with the noise PDF is unsurprising. This is in fact the point where, as we move to considering the PRE derivative per linear intensity bin, the distinct signatures of some parameters homogenize, and overall the PRE derivatives clearly decrease significantly in magnitude as noise spreads the actual information content of the model parameter across more intensity bins. A lower value of σ_n ≈ 600 Jy sr⁻¹ – effectively corresponding to a Stage 4 LIM survey – better preserves the distinct nature of some PRE derivatives (e.g. against σ_Z), but the noise in a Stage 2 LIM experiment ends up having a more drastic effect.

By a combination of effective averaging of the signal and smoothing of the actual VID, the beam and Gaussian noise in our survey simulations tend to Gaussianize the observation and decrease the amount of information-theoretical divergence imparted to the VID by the change in any given model parameter.

4.3.2 Discussion point 2: paths forward for optimal information recovery from the VID

While observational distortions do wipe away a significant amount of non-Gaussianity, it is possible to consider paths to extract the residual non-Gaussian signatures from the distorted VID. Detailed work on this front is beyond the scope of this paper, but we will briefly discuss a few possible avenues.

Appropriate filtering or processing of the observational data could conceivably improve detectability of the distinct signatures of different parameters. For instance, the relative entropy analysis of Lee et al. (2024) used not the full simulated CIB map, but maps filtered in different bands of harmonic multipoles. Similarly, Fourier or wavelet analyses could reduce the impact of noise.

However, it seems difficult to extract an information-theoretical free lunch from the noisy LIM observation, at least on its own. Certainly, the simplest noise-reducing filter, a rolling average, will ultimately tend to Gaussianize the signal even as it reduces the noise, which will be undesirable for certain signatures of skew or kurtosis. More advanced statistical methods could conceivably be better matched to the specific problem and perform better at isolating signal from noise in a non-linear, non-Gaussian fashion, but would likely require large amounts of training and thus be at least implicitly constrained by prior expectations and other external information.

One of the best empirical paths forward likely involves explicitly introducing external information and probing one-point statistics of cross-correlations or cross-correlations of one-point or higher order statistics. Although the isotropic stacks of Section 4.2 do not necessarily reveal the clustering around the external source catalogue very prominently on casual inspection, an oriented stack that emphasizes the filaments of the cosmic web may itself have one-point statistics that reveal non-Gaussian signatures of certain model variations. Similarly, the one-point statistics of a line-intensity map reconstruction through a non-linear extension of conventional techniques like linear covariance-based filtering (Chung 2023) could also access non-Gaussianities.

Another technique being developed for LIM applications is to deconvolve noise from signal through operations on the characteristic functions dual to the VID or counts-in-cells (Breysse et al. 2019, 2023; Chung et al. 2023). For reference, Fig. 12 shows the spectral densities of the PRE derivatives. Purely as a qualitative heuristic, the spectral densities suggest that at smaller scales (i.e. higher values of the Fourier dual of the histogrammed intensity), the PRE does contain some small amount of information content imprinted in different ways by the distinct parameters of our model. Therefore, by leveraging the mutual one-point statistical information between a [C ii] LIM survey and a correlated observation, Fourier-space deconvolved distribution estimation (DDE; Breysse et al. 2023) can potentially access discriminatory information content by rejecting disjoint noise and systematics.

Finally, we note that accessing the information content of line-intensity maps is subject to further observational distortions that we have not accounted for – e.g. scale-dependent sky noise, removal of interloper emission, and other mapmaking pipeline artefacts – many of which will involve non-linear operations. Tailored simulations of specific experimental scenarios using large ensembles of numerical simulations of the kind we have considered here will thus continue to be critical to evaluation of information content, covariance estimation, and tests of inference pipelines.

5 CONCLUSIONS

By now, we have answered all of the questions we put forward in Section 1 of this paper:

• What is the detectability of [C ii] emission in near- and far-future mm-wave LIM surveys, in either autocorrelation or cross-correlation? We expect a confident detection of the [C ii] power spectrum at z ∼ 6 with CCAT, followed by detection at z ∼ 7.5 with a Stage 2 LIM experiment. Stacking against a drop-out LBG survey can yield strong detections of the [C ii] signal near star-forming galaxies, as demonstrated with mock source densities corresponding to survey depths similar to those already achieved in existing surveys over multiple square degrees of sky.

• How do different model parameters affect the VID of the line-intensity fluctuations, as represented by the differential PRE? The PRE derivatives of many parameters become apparently degenerate with that of a simple uniform scaling of the signal brightness, but some parameters like the metallicity scatter σ_Z provide a distinct non-Gaussian signature that redistributes probability from middling intensities to both lower and higher intensities, or vice versa.

• How do observational effects like the instrument response and thermal noise affect the observability of these distinct signatures in relative entropy? By Gaussianizing the signal and smoothing out the VID, resolution effects and thermal noise both blur the signature in the log-VID, suppressing the observable PRE differential significantly and/or rendering all parameters degenerate with the overall brightness of the [C ii] fluctuations.

We only stand at the beginning of the larger endeavour of [C ii] LIM, which will continue to require significant logistical and theoretical effort for optimal analysis and interpretation. Future work should investigate techniques like oriented stacking or the Fourier-space DDE – or indeed entirely novel statistical methods – leveraging external cross-correlation or mutual information to isolate signal from interlopers and noise. As we stated towards the end of the previous section, such methods may even be key to accessing sufficient information content to discriminate between different kinds of astrophysical and even cosmological non-Gaussian signatures. This includes signatures of primordial intermittent non-Gaussianities propagating through large-scale structure and galaxy formation, which we will examine in future work through the lens of [C ii] LIM among other probes (Carlson et al., in preparation). In turn, these techniques and non-Gaussianities will require further study through rigorous simulation of both signal and noise.

In this aspect, a full understanding of observational effects will be important, but equally important is a full understanding of the [C ii] signal, or at least as full an understanding as one can get at this juncture. This work represents the state of the art in terms of a halo model grounded on physical insights from small-scale simulations, but even detailed simulations of individual galaxies are only as good as their inputs. Testing such models with [C ii] LIM will be desirable in future; this may demand the use not of manual tailored analyses of a simulation set (as in Liang et al. 2024) but of outputs of statistical methods like emulation or unsupervised dimensionality reduction in order to generate predictions for cosmological [C ii] fluctuations from a wide range of galaxy-scale models. Doubtless our understanding of the halo–[C ii] connection will continue to evolve with further observations of high-redshift [C ii], though for now constraints on [C ii] at z ≳ 6 rely on serendipitous detections and specific selections. It may happen that the first detections of cosmological [C ii] through LIM, in either autocorrelation or cross-correlation, will arrive at constraints free from such selection effects to inform further future LIM experiments.

ACKNOWLEDGEMENTS

Research in Canada was supported by NSERC and CIFAR. Parts of these calculations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund – Research Excellence; and the University of Toronto. Other parts of these calculations were performed on high-performance computing nodes at CITA; the authors thank John Dubinski and M. J. Huang for tireless maintenance of these resources.

The authors thank Jason Sun for kind and constructive comments on this manuscript during the time that it was in review, and for providing data for Fig. 5.

PH acknowledges support from CITA over the course of this research, including through the CITA Summer Undergraduate Research Fellowship and an NSERC Undergraduate Student Research Award.

The authors also thank Patrick Breysse and Clara Chung, whose work with JRB under the 2020 Summer Undergraduate Research Programme in astronomy and astrophysics at the University of Toronto acts as the foundation of this work.

DTC was supported by a CITA/Dunlap Institute postdoctoral fellowship. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. The University of Toronto operates on the traditional land of the Huron-Wendat, the Seneca, and, most recently, the Mississaugas of the Credit River. DTC also acknowledges support through the Vincent and Beatrice Tremaine Postdoctoral Fellowship at CITA.

LL acknowledges financial support from the Swiss National Science Foundation (grant number P2ZHP2_199729) and the University of Toronto Faculty of Arts and Science.

Finally, we thank the anonymous referee who reviewed the manuscript and provided useful comments.

This research made use of astropy,⁶ a community-developed core python package for astronomy (Astropy Collaboration 2013, 2018, 2022), as well as the numpy (Harris et al. 2020) and scipy (Virtanen et al. 2020) packages. This research also made use of NASA’s Astrophysics Data System Bibliographic Services.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: IMPLEMENTATION DETAILS FOR POWER SPECTRA AND WIENER FILTER CALCULATIONS

In order to carry out operations involving the power spectrum, we project the grid of our mock data cubes from observation space into comoving space. Given our cosmology and central redshift, it is straightforward to convert the x, y (in angular space), and z (in frequency space) coordinates of each voxel in our cube into comoving transverse and longitudinal positions (x_co, y_co, and z_co). We consider distortions from the light-cone geometry to be negligible for our simulated sky area and redshift range. We can compute the total span (Δx_co, Δy_co, and Δz_co) and the comoving distance per voxel (dx_co, dy_co, and dz_co) along each axis. Given the signal |$t(\boldsymbol {r})$| as a function of voxel coordinates |$\boldsymbol {r} = (x_\text{co},y_\text{co},z_\text{co})$|⁠, the power spectrum is then given as a function of the dual wave vector |$\boldsymbol {k}$|⁠:

with |$\operatorname{FFT}(\boldsymbol {k})$| denoting the 3D discrete Fourier transform (DFT) of the signal |$t(\boldsymbol {r})$|⁠.

In terms of positional indices (n_x, n_y, n_z) running from (0, 0, 0) to (N_x − 1, N_y − 1, N_z − 1) and wave vector indices (K_x, K_y, K_z) running through the same range, the DFT under our convention (the convention used by numpy) is

with the triple sum being over all voxel positions |$\boldsymbol {r}$| described by (n_x, n_y, n_z).

We map the wave vector indices (K_x, K_y, K_z) back to the comoving wave vector |$\boldsymbol {k} = (k_x, k_y, k_z)$| dual to |$\boldsymbol {r}$| based on the total width (Δx_co, Δy_co, Δz_co) and sampling (dx_co, dy_co, dz_co) along each dimension. The factor of the voxel volume in equation (A1) applied to |$\operatorname{FFT}(\boldsymbol {k})$| rescales the DFT to be an approximation of the continuous Fourier transform, appropriate for use in the power spectrum calculation.

At each point of our Fourier transform cube, we compute the wavenumber |$k=\sqrt{k_x^2 + k_y^2 + k_z^2}$|⁠. We define a one-dimensional (1D) k-grid between 0 and the maximum k value across the 3D |$\boldsymbol {k}$|-space. We obtain an average of the 3D power spectrum |$P(\boldsymbol {k})$| weighted by the number of Fourier modes N_m(k) belonging to each k-bin, and thus retrieve the spherically averaged power spectrum of our data cube.

Fig. A1 shows the comparison between the 1D power spectrum of the pure signal P_s(k) and the noise power spectrum amplitude |$P_\text{noise}=\sigma _\text{n}^2V_\text{vox}$|⁠, as well as the power spectrum of the forecast observation P_f(k) = P_s(k) + P_noise. The sensitivities shown in Fig. 5 are found effectively by taking |$P_\text{noise}N_\text{m}^{-1/2}(k)$|⁠, while the signal-to-noise ratio at each k is found as |$P_\text{s}(k)N_\text{m}^{1/2}(k)/[P_\text{s}(k)+P_\text{noise}]$|⁠. Again, we refer the interested reader to Chung et al. (2020) for further details.

The Wiener filter is the optimal linear estimator for the signal given known signal and noise power spectra. Each Fourier mode is scaled by the following weighting function:

Given the inverse DFT operation |$\operatorname{IFFT}$|⁠, we can calculate the Wiener-filtered signal |$\widetilde{t}(\boldsymbol {r})$| from the forecast |$T(\boldsymbol {r})$| incorporating both signal |$t(\boldsymbol {r})$| and noise |$\sim\!\! \mathcal {N}(0,\sigma _\text{n})$|⁠:

(Note that our notation is intentionally slightly loose for simplicity, making implicit the mappings between discrete indices and continuous coordinates.) As Fig. A1 shows, this results in a filtered signal with a power spectrum given by

which approaches the original signal for signal-dominated k, and is heavily suppressed for noise-dominated k.

APPENDIX B: PRE DERIVATIVES RELATIVE TO MODEL PARAMETER FIDUCIAL VALUES

Fig. 10 in the main text shows the PRE derivatives with respect to the parameter values, and rescaled to emphasize a key degeneracy in the information content of multiple parameters. Fig. B1 is an alternate visualization, albeit with very much the same effect. Here, instead of d(dS_rel)/dλ for each parameter λ considered, we plot |$\lambda \, \mathrm{ d}(\mathrm{ d}S_\text{rel})/\mathrm{ d}\lambda$|⁠. This is equal to the PRE derivative against the natural log of the parameter, or the PRE differential for a relative unit change in a parameter instead of an absolute unit change. We also only apply sign changes to better illustrate the same degeneracy observed between parameters. Otherwise, we do not apply rescalings as aggressively as in Fig. 10, except for replacing log M_min with log M_min − 10.