© ESO 2024

1. Introduction

In its first billion years, the Universe underwent numerous changes. Following the Big Bang, the hot and ionized plasma cooled down, leading to the formation of neutral hydrogen atoms (e.g., Peebles 1968; Zel’dovich et al. 1969; Seager et al. 2000). This resulted in the Universe being in a state of darkness, known as the “Dark Ages”. The advent of the first stars and galaxies, roughly 300 million years after the Big Bang, re-illuminated the Universe (e.g., Harikane et al. 2022). What makes this period in the history of the Universe even more special is that the appearance of the first light sources concurs with a phase transition in the intergalactic medium, from neutral to ionized (e.g., Becker et al. 2001, 2015; Barkana & Loeb 2001; Fan et al. 2006; Venemans et al. 2013), a period known as the epoch of reionization (EoR).

The connection between the reionization of the Universe and the emergence of the first light sources is well established in current theoretical models (e.g., Zaroubi 2013). Early massive stars located in the first incipient galaxies were the primary drivers of reionization, powering it through their ultraviolet (UV) emission. Despite this understanding, the intricacies of the interaction between early galaxy formation, evolution, and the reionization process remain elusive. Key uncertainties include the energy emitted by the early stars, the fraction of this energy that escaped the early galaxies (e.g., Bouwens et al. 2016), as well as the distribution of those ionizing sources in space and time. Unfortunately, these galaxies are too faint to be detected in large samples with current telescopes, including the state-of-the-art James Webb Space Telescope (JWST, Boylan-Kolchin et al. 2015), which creates challenges for accurately estimating the properties of their population and does not allow for tracing the large-scale structure in which they reside.

A more effective approach than traditional galaxy surveys for tracing both star formation in early galaxies and the large-scale structure in which they reside is line intensity mapping (LIM; see, e.g., Kovetz et al. 2019; Bernal & Kovetz 2022). This method measures the integrated emission of spectral lines from galaxies within a given spatial-spectral resolution element (voxel) using a spectro-imager. By focusing on the singly ionized carbon 158 μm fine-structure line, [CII], which is the brightest line in typical star-forming galaxies (e.g., Stacey et al. 1991) and an indicator of star formation rate (SFR; e.g., De Looze et al. 2014; Le Fèvre et al. 2020; Béthermin et al. 2020; Schaerer et al. 2020; Bouwens et al. 2022), a three-dimensional (3D) tomography of galaxies during the early universe can be obtained (e.g., Karoumpis et al. 2022). Several LIM experiments, such as the EXperiment for Cryogenic Large-Aperture Intensity Mapping (EXCLAIM, Pullen et al. 2023), the Terahertz Intensity Mapper (TIM, Vieira et al. 2020), the CarbON CII line in post-rEionization and ReionizaTiOn (CONCERTO) project (CONCERTO Collaboration 2020), the Tomographic Ionized-carbon Mapping Experiment (TIME, Sun et al. 2021), and the Fred Young Submillimeter Telescope (FYST, CCAT-Prime Collaboration 2023), are underway to target [CII] emission coming from EoR (6 <  z <  9) and post-EoR (3 <  z <  6) redshifts.

Many of these observatories will prioritize measuring the variance in the Fourier mode of these 3D tomographies through their spherically averaged power spectrum (PS). This statistical metric, although relatively simple, offers powerful limits on the underlying line luminosity function, given that the PS correlates with both the first moment (for scales > 10 Mpc) and the second moment (at scales < 10 Mpc) of the line luminosity function (Karoumpis et al. 2022). Such robust constrains of the luminosity functions of [CII] will provide rigorous statistical constraints on the SFR of (post-)EoR galaxies–pivotal, yet elusive parameter in reionization models (e.g., CCAT-Prime Collaboration 2023; Karoumpis et al. 2022).

In this context, where precise measurements of the [CII] PS at different scales during and after the EoR will be possible, it is crucial to study how various cosmic galaxy evolution models influence this PS and can thus be tested by comparison with future LIM measurements. To contribute to this effort, in the first article of this series (Karoumpis et al. 2022), we generated a span of predictions by post-processing the dark matter halo catalog from the IllustrisTNG300 hydrodynamical simulation (TNG300; Pillepich et al. 2018). We then used our predictions to assess the feasibility of detecting the PS of the [CII] line from galaxies at redshifts between 3 and 8 using the spectro-imager of FYST. Our results, which are consistent with empirically motivated predictions (Clarke et al. 2024), demonstrated promising potential for detecting the [CII] PS at the critical comoving length scale of 10 Mpc in four selected FYST bands centered at redshifts 3.7, 4.3, 5.8, and tentatively at redshift 7 (i.e., 410, 350, 280, and 225 GHz). However, those predictions did not account for a significant challenge faced by the LIM technique: the difficulty in distinguishing the contribution of (post-)EoR [CII] galaxies from the infrared (IR) continuum and carbon monoxide (CO) rotational line emission (2 <  J_up <  12) of galaxies located at the same solid angle but at lower redshifts (z = 0 − 5.7; see Fig. 1 and Table 1). While the IR continuum emission is not a major concern for the PS measurement as it is spectrally smooth and its frequency-coherent spectrum can be fitted and removed (Van Cuyck et al. 2023), the contribution of CO emitters presents a real challenge since the variance of their emission is of the same order or even exceeds that of the [CII] galaxies (e.g., Béthermin et al. 2022; Roy et al. 2023).

Fig. 1. Spectral lines observable within four selected EoR-Spec bands (illustrated by red-shaded regions), and originated from galaxies situated at different redshifts. The dashed numbered lines represent transitions within the CO rotational ladder, with the corresponding numbers indicating the Jup for each specific transition. These CO lines act as the foreground contaminant for the [CII] (red thin dashed line) LIM survey.

In this paper, we utilized the halo catalog from our previous work (Karoumpis et al. 2022), and by implementing empirical relations between SFRs and CO line luminosities of galaxies, we modeled the foreground line emission corresponding to the same frequencies as the (post-)EoR [CII] emission. Drawing from these predictions, we evaluated a foreground removal method where voxels containing luminous CO emitters are identified and masked using mock external catalogs with realistic characteristics. To this end, we assumed that the FYST LIM survey will cover a region of the sky benefiting from an external catalog with similar stellar mass completeness limits as the Euclid deep fields (Euclid Collaboration 2022). This external catalog provides the angular position, redshift, stellar mass and SFR of the galaxies and it is used to pinpoint and mask voxels with a high likelihood of containing bright CO emitters. By carefully simulating this technique, we not only assess its effectiveness in mitigating CO contamination, but also examine its influence on the recovered [CII] PS. This includes an analysis of heightened measurement uncertainty resulting from diminished survey volume and the effects of a convolved PS attributable to the window function of the mask.

The paper is organized as follows. Section 2 introduces the formalism adopted for the modeling of the [CII] and CO emission. Section 3 describes the construction of the mock external catalog, and provides an overview of the masking technique. Section 4 examines the statistical properties of the individual CO lines, the impact of masking on the CO PS, the targeted [CII] PS, and the total line CO+[CII] PS. Finally, in Section 5, we critically assess the limitations of our work and suggest potential improvements and additions to our models. Throughout the paper, we adopt a ΛCDM cosmology with the same parameters used in the TNG300 simulation (Pillepich et al. 2018): h = 0.71, Ω_b = 0.046, Ω_m = 0.281, σ₈ = 0.8, Ω_Λ = 0.719, and n_s = 0.963.

2. Generating the mock FYST LIM survey

The Fred Young Submillimeter Telescope (FYST, CCAT-Prime Collaboration 2023) is set to be a state-of-the-art, 6-meter diameter telescope designed for submillimeter to millimeter observations. Placed at 5600 meters on Cerro Chajnantor, FYST will utilize a novel crossed-Dragone optical design (Dragone 1978; Parshley et al. 2018) for fast and efficient wide-field-of-view sky mapping. One of its key instruments, the Prime-Cam receiver, will offer impressive spectroscopic and broadband measurement capabilities, enabling a mapping speed over ten times faster than existing facilities (Vavagiakis et al. 2018). Prime-Cam features seven instrument modules tailored to specific scientific programs. Two of the modules, named Epoch of Reionization Spectrometer (EoR-Spec, Huber et al. 2022), integrate a Fabry-Perot interferometer (Perot & Fabry 1899; Zou et al. 2022) to observe the [CII] line emission at z = 3.5 − 8.05 (210 to 420 GHz) with a resolving power of approximately Δλ/λ = 100. EoR-Spec will target two 5 deg² regions in a survey dedicated to measuring the [CII] PS from (post-)EoR galaxies (CCAT-Prime Collaboration 2023). The chosen survey fields require sensitive auxiliary observations with extensive wavelength coverage to help remove foreground sources (see Sect. 3). Consequently, the Extended-COSMOS (E-COSMOS, Aihara et al. 2018), and Extended-Chandra Deep Field South (E-CDFS, Lehmer et al. 2005) contained in the Euclid Deep Field Fornax (EDF-F, Euclid Collaboration 2022) are the selected fields for the FYST LIM survey. Existing deep coverage in various bands for these fields will soon be augmented by forthcoming imaging and spectroscopy data from present and future state-of-the-art observatories like JWST, Roman, Euclid, and the Large Millimeter Telescope (LMT). For a comprehensive overview of both current and upcoming datasets, see CCAT-Prime Collaboration (2023).

In this context, realistic simulations of the FYST LIM survey are crucial. They allow for the optimization of the survey strategy, the examination of its ability to detect the [CII] LIM signal during the (post-)EoR, and the testing of various foreground mitigation strategies. To this end, we utilize the TNG300 simulation to create a 4° ×4° cone of mock galaxies spanning from z ≈ 0 − 9 (see Sect. 2.1). The [CII] line emission from these galaxies is then incorporated, focusing on the z ≈ 3.4 − 8.3 range to which the FYST observed frequencies are sensitive (Karoumpis et al. 2022, hereafter Paper I; see Sect. 2.2 for a summary). Subsequently, we introduce the CO (J_up = 3 − 12) foregrounds which, in the context of FYST, are emitted by galaxies in the z ≈ 0 − 5.7 range (see Sect. 2.3). The constructed observable cone is then transformed into mock FYST LIM tomographies (see Sect. 2.4), containing both the foreground CO and background [CII] emission, that serves as the basis for evaluating our masking technique.

2.1. The cone of mock galaxies

In this study, we utilized the dark matter (DM) halo catalog presented in Paper I, which contains all the essential results of the TNG300 simulation needed to estimate the line luminosity of the mock galaxies (i.e., angular position, redshift, M_*, SFR) which are defined as the gravitationally bound substructures hosted in the DM halos. This catalog corresponds to an observational cone encompassing a sky area of 4° ×4° and a redshift range of z = 0 − 9, and it was constructed using the TNG300 halo catalogs as building blocks and corrected for resolution effects with the TNG100 DM halo catalogs as a reference (Pillepich et al. 2018). However, for the redshifts of the cone where galaxies emit the brightest CO lines in the frequency range of FYST (i.e., z = 0 − 5.7), we recalculated the resolution corrections, this time using the newly available TNG50 simulation (Pillepich et al. 2019). The reason for this update was that the CO emission originates from redshifts where the TNG100 simulation does not accurately reproduce the observed evolution of the cosmic SFR density (Madau & Dickinson 2014), with the peak epoch appearing at earlier cosmic times in the simulation (z ≈ 3) than the observations (z ≈ 2). Owing to its 17 times better mass resolution, TNG50 closely follows the cosmic SFR density evolution of Madau & Dickinson (2014). The same applies for our cone after the new mass resolution corrections.

To quantify the expected sample variance, we employed an analytical relation derived by Gkogkou et al. (2023) to directly calculate the coefficient of variation of the PS (CV_PS) of the LIM tomography. In their study, Gkogkou et al. (2023) investigated the impact of field-to-field variance (which we refer to as sample variance) on CO and [CII] PS predictions for future LIM experiments. This was achieved by combining a 117 deg² dark matter lightcone with a model of the infrared sky. They developed an analytical formula that estimates the expected CV_PS of the [CII] and CO LIM tomographies, given their dimensions. Upon inputting the survey parameters of a 40 GHz frequency coverage and a 4° ×4° field into their formula, we find that the CV_PS for the clustering (shot noise) component of the CO PS peaks at 410 GHz, with a value of about 7% (4%), while for the [CII] PS, the CV_PS peaks at 225 GHz with a value of about 8% (7%).

Our predictions of line intensities are based on a single realization of this 4° ×4° field and therefore may be subject to sample variance, due to the small number of galaxies traced. This is particularly true at low redshift (near z = 0) where the comoving volume covered by our cone is small, and at high redshift (near z = 7) where star formation occurs predominantly in overdense regions. To quantify the expected sample variance, we employed an analytical relation derived by Gkogkou et al. (2023) to directly calculate the coefficient of variation of the PS (CV_PS) of the LIM tomography. In their study, Gkogkou et al. (2023) investigated the impact of field-to-field variance (which we refer to as sample variance) on CO and [CII] PS predictions for future LIM experiments. This was achieved by combining a 117 deg² dark matter lightcone with a model of the infrared sky. They developed an analytical formula that estimates the expected CV_PS of the [CII] and CO LIM tomographies, given their dimensions. Upon inputting the survey parameters of a 40 GHz frequency coverage and a 4° ×4° field into their formula, we find that the CV_PS for the clustering (shot noise) component of the CO PS peaks at 410 GHz, with a value of about 7% (4%), while for the [CII] PS, the CV_PS peaks at 225 GHz with a value of about 8% (7%). Nevertheless, in Sects. 2.2 and 2.3 we demonstrate that uncertainties on the modeling of the CO and [CII] emission of galaxies introduced much larger uncertainties on the CO and [CII] PS. For example, in the cases of the CO PS at 350 GHz and the [CII] PS at 225 GHz, the relative variability of the models in relation to their mean, (P_max − P_min)/(P_max + P_min), is ≈85% and ≈90% respectively. Therefore, while sample variance affects LIM survey covering only 4° ×4° field, it does not significantly affect the evaluation of the masking technique performed in this paper.

2.2. The [CII] emission

The mock galaxy catalog detailed in Sect. 2.1 underwent various post-processing methods to associate [CII] emission to its galaxies. We restricted these calculations to the redshift range of z ≈ 3.4 − 8.3, where the [CII] emission of galaxies is redshifted into the observing bands of EoR-Spec. All the related methods and calculations are described in detail in Paper I and here we only summarize the most important steps. Firstly, SFRs were attributed to mock galaxies in two ways. One approach used their intrinsic TNG300 SFR, corrected to account for the mass resolution limitation of TNG300 (see Sect. 2.1). The other approach matched the mock galaxy abundance with the observed, dust-corrected ultraviolet luminosity function of high-redshift galaxies (Bouwens et al. 2015), assigning luminosities to the mock galaxies which are subsequently converted into SFRs following Kennicutt (1998). Secondly, the [CII] luminosities of the mock galaxies were estimated from the SFRs using three different SFR-to-L_[CII] relations: from a semi-analytic model of galaxy formation (Lagache et al. 2018), from a hydrodynamical simulation of a high-redshift galaxy (Vallini et al. 2015), and from a high-redshift [CII] galaxy survey (Schaerer et al. 2020).

The various galaxy-to-SFR and SFR-to-L_[CII] relations translate into large variations of the amplitude of the expected [CII] PS (Paper I). Such variations render the assessment of the detectability of the [CII] PS more problematic but at the same time demonstrate the utility of this signal in discriminating between all these models. In the following, our fiducial [CII] model refers to the combination of the recalibrated TNG300 SFR with the SFR-to-L_[CII] relation from Vallini et al. (2015). It was chosen as it produces results that closely align with the median value of our predictions (Paper I).

2.3. The CO emission

The CO molecule is the second most abundant molecule in the Universe, surpassed only by molecular hydrogen (H₂). In the typical conditions of giant molecular clouds, the rotational transitions of the CO that emit line emission are easily excited, making it the most widely used tracer of molecular gas in both the local (e.g., Wang & Hwang 2020; den Brok et al. 2021) and high-redshift universe (e.g., Solomon & Vanden Bout 2005; Daddi et al. 2015; Aravena et al. 2016; Riechers et al. 2019). Its luminosity is proportional to its mass, assuming the number of clouds is small enough for them not to overshadow each other and that they are close to virial equilibrium (e.g., Carilli & Walter 2013). The strong correlation between the CO line luminosities of galaxies and their molecular gas content enables us to forecast these foregrounds using our cone of mock galaxies.

2.3.1. The CO (1–0) luminosity

We estimated the molecular gas content of a galaxy, relying on the well-established relationship between molecular gas and star formation rate. Several studies have shown that the depletion time (t_dep = M_gas/SFR) does not depend significantly on the stellar mass of galaxies and only slightly increases from z ≈ 4 to z ≈ 0 (e.g., Scoville et al. 2017; Tacconi et al. 2018; Kaasinen et al. 2019; Magnelli et al. 2020; Wang et al. 2022). Therefore, changes in the SFR of galaxies with stellar mass and redshift are primarily driven by changes in their gas reservoir. This allowed us to make the simplifying assumption that all galaxies convert their molecular gas into stars at a constant universal rate. To estimate the molecular gas masses of our mock galaxies, we multiplied their SFRs (the corrected intrinsic TNG300 SFRs presented in Sect. 2.1) by t_dep = 570 Myr, corresponding to the depletion time of z ≈ 2 galaxies with a stellar mass of M_* ≈ 10¹⁰ M_⊙ (Tacconi et al. 2018). We verified that even with this simple approach, the cosmic molecular gas mass density inferred from our cone of mock galaxies was compatible with the observations of Riechers et al. (2019), Magnelli et al. (2020), Decarli et al. (2020), and Boogaard et al. (2023) (Fig. 2). This agreement is particularly true around the peak of the molecular gas mass density at z ≈ 2, from where the bulk of the CO foregrounds originate. A more complex model that accounts for changes in the depletion time with redshift and distance from the main sequence (MS; for a definition, see Sect. 2.3.2) produces a similar gas mass density evolution, but to avoid introducing additional parameters, we adopted the simpler model of a constant depletion time.

Fig. 2. Redshift evolution of the cosmic molecular gas mass density of galaxies. The blue horizontal line represents predictions from our cone of mock galaxies, while observational constraints are presented by rectangles. The gray rectangles are for the dust-based estimates of Magnelli et al. (2020), the light gray rectangles are for the COLDz CO survey (Riechers et al. 2019), the light blue rectangle for the CO-based estimates of Boogaard et al. (2023), the red-shaded rectangles are for the ASPECS CO pilot survey (Decarli et al. 2020), while the pink rectangles are for the ASPECS CO LP survey (Decarli et al. 2020). The width of the rectangles indicates the redshift bin size and the height the 1-σ confidence region.

The calculated molecular gas masses were then used to predict the CO line emission of each mock galaxy. We adopted a constant light-to-mass ratio between the specific luminosity of the CO(1−0) line and the molecular gas mass,

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{CO}}=\frac{M_{\mathrm{gas}}}{L{\prime }_{\mathrm{CO}(1-0)}}=3.6M_{\odot }/\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^2,\end{array}$$(1)

as its value for the MS galaxies does not seem to change significantly up to z ≈ 3 (Decarli et al. 2016), and is adopted by the majority of the molecular gas measurements (Daddi et al. 2015; Riechers et al. 2019; Decarli et al. 2020; Lenkić et al. 2020; Chung et al. 2022).

There is observational evidence that starburst galaxies might require adjustments to their α_CO and t_dep. Specifically, a lower t_dep (≈150 Myr as suggested in Tacconi et al. 2018) and a reduced α_CO value of 0.8 (as suggested in Downes & Solomon 1998). Interestingly, the ratio of the proposed depletion times (570/150) is roughly equivalent to the ratio of the α_CO values (3.6/0.8), meaning these changes may roughly offset each other. Therefore, while incorporating these adjustments would not significantly alter the L_{CO(1 − 0)} of our starbursts, it complicates the model further. For this reason, we have chosen to use the simpler model of a constant depletion time and α_CO. Nevertheless, the effect of starburst galaxies is further examined in Sect. 2.3.2.

More generally, although we may underestimate the stochasticity in our scaling relations for the CO(1−0) luminosities due to the simple assumptions of constant t_dep and αCO, the significant uncertainties introduced in Sect. 2.3.2 make any additional scatter in the CO(1−0) luminosities unnecessary for the scope of this study. There, by making some assumptions on the spectral line energy distribution (SLED) of our mock galaxies, we use these CO(1−0) line luminosities to predict the CO line luminosities in transitions that pose the most significant foreground challenges for our LIM survey (from J_up = 3 to J_up = 12).

2.3.2. The CO SLED

The CO SLED of each galaxy is largely defined by the physical conditions (gas density and temperature) of its interstellar medium (ISM; e.g., Narayanan & Krumholz 2014). It is established that due to the diverse range of conditions in the ISM, the CO SLED of a galaxy often requires more than one temperature component to accurately describe it (Valentino et al. 2020). However, simulating the complex interactions of gas physics is beyond the scope of this paper.

Fig. 3. CO luminosity function for various transitions and redshifts. The blue line corresponds to our fiducial model, while the shaded region indicates the range resulting from our optimistic and pessimistic models. Rectangles represent observational data from the ASPECS LP survey (Decarli et al. 2020), with the width indicating the bin size and the height indicating the 1-σ confidence region.

Here we employed a model with a small number of free parameters, yet sufficiently advanced to replicate the CO luminosity function in various transitions observed by the Atacama Large Millimeter/submillimeter Array (ALMA) Spectroscopic Survey (ASPECS, Decarli et al. 2020; Boogaard et al. 2020; Riechers et al. 2020) at frequencies akin to our 225 GHz band (see Fig. 3 and Table 1). To achieve this, we adopted an empirical approach based on the galaxy position in the M_* − SFR diagram, which is a valuable tool for evaluating the properties of a galaxy. The majority of the galaxies on this diagram follows a so-called MS (Elbaz et al. 2007), wherein more massive galaxies tend to have higher SFRs. However, there are distinct exceptions: starburst galaxies, for instance, exhibit exceptionally high SFRs for their stellar masses; or passive galaxies, that have low SFR for their stellar masses. Here, we defined our MSs by fitting a log-linear relation to the M_* − SFR sequence followed by our mock galaxies at different redshifts. Then, we parameterized the position of a galaxy in this diagram using their distance from the MS, ΔMS = log(SFR/SFR_MS), where SFR_MS is the MS SFR at the redshift and stellar mass of the galaxy. Galaxies close to and below the MS have relatively cold ISM (e.g., Carilli & Walter 2013; Magdis et al. 2021), while those above the MS, such as starburst galaxies, have a warmer, denser ISM (e.g., Silich et al. 2009; Carilli & Walter 2013; Kamenetzky et al. 2017). Based on this assumption, we modeled the CO SLEDs of our mock galaxies by linearly combining the observed and well sampled CO SLEDs of a MS-like galaxy and a starburst-like galaxy. We used the CO SLED of the Milky Way (MW), as a characteristic star-forming galaxy with cold ISM conditions, slightly below the MS. On the other hand, we used the CO SLED of NGC 253, a characteristic starburst galaxy with warm and dense ISM condition. This linear combination allowed us to create a continuum of conditions that encompassed a broad range of galaxies from MS to starburst. This approach is supported by observations, which report a correlation between the CO SLED shape of a galaxy and its distance from the MS (e.g., Valentino et al. 2020; Cassata et al. 2020).

The CO SLED of the MW and NGC 253 used in our study are displayed in Fig. 4. The MW observations up to J_up = 5 were taken from Carilli & Walter (2013), and the SLED was extended to J_up = 6 − 10 using the CO(J_up − J_low)/CO(5 − 4) measurement from Wilson et al. (2017). The CO SLED of NGC 253 was taken from Mashian et al. (2015). The difference between the MW and NGC 253 CO SLEDs, which is greater at J_up = 5, is due to the contrasting star formation environments in these galaxies. The MW features a cold, more diffuse molecular gas component, while the intense star formation of NGC 253 creates a warm and dense gas environment, leading to higher excitation of CO molecules and an enhancement mid- to high-J (J_up = 4 − 8) transition emission.

Fig. 4. CO SLED of the Milky Way (black circles) and NGC 253 (black squares). The SLEDs resulting from the linear combination of the two, following Eq. (3), are shown by solid lines, with μ varying from 0 to 1 in steps of 0.1. The red line represents the combination (with μ = 0.4) that best fit the BzK galaxies (red points) and is used here as reference for high-redshift MS galaxies.

We assigned each mock galaxy a SLED resulting from its distance from the MS and the linear combination of these two extreme CO SLEDs. The mixing parameter (μ) was calibrated so that ΔMS = 0 yields a SLED consistent with typical MS galaxies at z = 1.5 (BzK galaxies; Daddi et al. 2015), and ΔMS = 1.5 produces the SLED of NGC 253. In turn, this yields,

$$\begin{array}{c}\hfill \mu =\frac{\mathrm{\Delta }\mathrm{MS}+1}{2.5}·\end{array}$$(2)

The CO SLED resulting from the linear mixing is given by the equation,

$$\begin{array}{cc}\hfill log{L}_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{\prime \mathrm{galaxy}}=& (1-\mu )log{L}_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{\prime \mathrm{MW}}\hfill \\ \hfill & +\mu log{L}_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{{}^{\prime }\mathrm{NGC}253},\hfill \end{array}$$(3)

where $L_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{\prime \mathrm{galaxy}}$, $L_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{\prime \mathrm{MW}}$, and $L_{\mathrm{CO}(J_{\mathrm{up}}-J_{\mathrm{low}})}^{\prime \mathrm{NGC}253}$ are the CO(J_up − J_low) specific line luminosities of our mock galaxy, the MW, and NGC 253, respectively. Figure 4 displays an example of the individual CO SLEDs resulting from Eq. (3) (for μ = 0 − 1 with steps of δμ = 0.1). The linear combination highlighted in red corresponds to ΔMS = 0 (i.e. μ = 0.4), and aligns by design closely with the observed CO SLEDs from BzK galaxies.

This approach where each galaxy has its individual CO SLED constitutes our fiducial model. However, we also considered two extreme models in addition to our fiducial, creating a larger range of possible values for the CO intensity tomographies. In the “low-contamination” model, all galaxies have the SLED of a typical MS galaxy (i.e., μ = 0.4; BzK galaxies; Daddi et al. 2015), while in the “high-contamination” model, all galaxies have the SLED of a typical starburst (i.e., μ = 1; NGC 253).

To verify the accuracy of this approximation, we compared the CO luminosity function inferred from our fiducial, “low-” and “high-contamination” models to that observed by ASPECS (Fig. 3). Our model predictions cover a good range of the 1-σ interval compatible with the ASPECS measurements, with, however, a slight over-prediction of the faint end and under-prediction of the bright end. This former discrepancy could be due to the fact that ASPECS is a line-flux-limited CO survey, which naturally becomes less sensitive and more incomplete at the faint end of the CO luminosity function. The difference at the bright end between our models and the survey has been reported previously in the semi-analytical approach of Popping et al. (2019), Béthermin et al. (2022), Gkogkou et al. (2023). As suggested by Davé et al. (2020), this could be due to the assumption of a constant α_CO for all of our galaxies. In any case, the galaxies of the brightest bin are few in number and therefore constitute a small percentage of the overall number of interlopers that should be masked. This slight underestimation should not therefore significantly impact our results.

2.4. The [CII] and CO mock intensity tomographies

In order to translate the line luminosities of our mock galaxies into mock line intensity tomographies, we assumed that the intensity of a line (either CO J_up = 3 − 12 or [CII]) in a given voxel, positioned at a specific RA, Dec, and ν, is expressed as:

$$\begin{array}{c}\hfill I_{\mathrm{line}}=\frac{1}{\mathrm{\Delta }{\nu }_0{(\mathrm{\Delta }{\theta }_{\mathrm{b}})}^2}\left(\sum _{j\in \mathrm{voxel}}\frac{L_{\mathrm{line}}^j}{4\pi r_j^2{(1+z_j)}^2}\right)\ast G(\mathrm{RA},\mathrm{Dec},\nu ),\end{array}$$(4)

where r_j is the comoving distance of the j-th galaxy, which resides at a redshift z_j and has a line luminosity $L_{\mathrm{line}}^j$; is the symbol of the convolution; and G(RA,  Dec, ν) is the 3D Gaussian function representing the angular and spectral resolution element of the EoR-Spec. To adhere to the Nyquist-Shannon sampling theorem (Shannon 1949), we selected the grid spacing of the mock tomographies to be three times smaller than the full width at half maximum (FWHM) of this 3D Gaussian function. Since the summation of all galaxy luminosities within a voxel essentially represents the integration of the intensity field over the volume of the voxel, we divide twice by Δθ_b = FWHM_RA/3 = FWHM_Dec/3 = FWHM_ang/3 and once by Δν₀ = FWHM_ν/3 to account for the spatial and frequency dimensions of the voxel. This ensures that the calculated intensity is properly normalized by the volume of the voxel.

We assumed that the FYST LIM survey will be divided into four tomographies for each of the four selected frequency bands of the EoR-Spec. These bands are centered at 225 ± 20, 280 ± 20, 350 ± 20, and 410 ± 20 GHz. The average FWHM of the spectral channels for these bands are 2.1, 2.7, 3.6, and 4.4 GHz, respectively. Additionally, the angular beam FWHM for these bands are 0.88, 0.77, 0.65, and 0.62 arcmin, respectively. The intended survey will cover an area of 5 deg². However, for our analysis, we focused on a 16 deg² field, as it offers a testbed with reduced sample variance. It is important to note that when calculating uncertainties (Sect. 4), we considered the sample variance and the white noise of the expected (i.e. two fields of 5 deg²) size and observational time (2000 hours for each field) of the FYST LIM survey. This way we ensure that while the statistical parameters of the foregrounds and the targeted signal approximate the global values, the impact of any low number statistics is effectively captured within the uncertainty estimates.

3. The masking technique

While various alternative methods for addressing submillimeter line foreground contamination have been proposed in the literature, including spectral template fitting (Kogut et al. 2015; Cheng et al. 2020), blind masking (Breysse et al. 2015), and machine learning techniques (Moriwaki et al. 2020; Moriwaki & Yoshida 2021; Zhou et al. 2023), we focus on what appears to be the most effective approach–targeted masking (e.g., Yue et al. 2015; Silva et al. 2015; Sun et al. 2018). This technique employs external catalogs to identify the locations of the brightest foreground sources and mitigates their influence on the targeted PS by excluding the implicated voxels from further analysis.

In our study, the masks were created based on a mock external stellar mass-selected catalog (Sect. 3.1). In order to explore several depths of masking, for each CO line transition in each FYST frequency band, we created two masks of different depth: the “bright” and the “complete” masking arrays. To generate the “complete” masking arrays, we masked every galaxy in our stellar mass-complete external catalog with redshift z_Jup = (ν_Jup/ν_obs ± 20 GHz)−1, where J_up refers to the upper level of the CO transition, and ν_obs is the central observed frequency of one of the four selected EoR-Spec bands.

For the “bright” masking arrays, we employed the same method but limited it to galaxies above the MS, (i.e., ΔMS >  0). This way to identify bright CO contaminants solely based on ΔMS stemmed from the observation that the cumulative CO line contamination from galaxies correlates with the MS boundary (see e.g., Fig. 5). By exclusively masking galaxies above the MS, we effectively limited our masks to those galaxies that significantly contaminate the targeted [CII] signal. This approach optimizes the survey volume available for our PS analysis. The masking arrays, constructed using these two subsamples, are presented in Sect. 3.2, and the methodology behind this masking strategy is discussed in further detail in Sect. 4.2.

Fig. 5. SFR as a function of stellar mass with the bins color-coded to reflect the cumulative CO (6−5) luminosity of the 2D bin. The dashed line indicates the thresholds that define our “bright” subsample with ΔMS >  0 (i.e., galaxies on and above the MS).

3.1. The mock external catalog

The fields that will be targeted by the FYST LIM survey (CCAT-Prime Collaboration 2023) are the E-COSMOS (Aihara et al. 2018) and the E-CDFS (Lehmer et al. 2005), with the latter potentially being extended to cover the whole EDF-F. To mimic the way these surveys would observe our cone of mock galaxies, we simply applied to our mock galaxy catalog redshift-dependent stellar mass limits akin to those affecting an Euclid-deep-like survey (Euclid Collaboration 2020). According to the latest forecasts, such a survey will have multi-wavelength photometric sensitivities that will be equivalent to those obtained today in the COSMOS field. By taking the stellar mass completeness limits from the COSMOS2015 (Laigle et al. 2016) and COSMOS2020 (Weaver et al. 2022) catalogs, we should thus obtain realistic redshift-dependent stellar mass limits to apply to our mock galaxy catalog. These redshift-dependent limits are shown in Fig. 6, starting at ≈0.5 × 10⁹ M_⊙ at low redshifts (z ≈ 0 − 2) and rising to ≈10⁹ M_⊙ at high-redshifts (z ≈ 2 − 4).

Fig. 6. Comparison of the stellar mass completeness limits applied to our mock external catalog (red line) with those of the COSMOS2015 (light gray line Laigle et al. 2016) and COSMOS2020 (dark gray line Weaver et al. 2022) catalogs. The stellar mass completeness limits of the deep Euclid catalogs is expected to lie between these two lines.

In addition to the angular position, redshift information, and stellar mass, we also included in our mock external catalog the SFR for each galaxy. Indeed, we anticipated that the Euclid Deep Fields will provide SFR estimates for the majority of detected galaxies through the spectral energy distribution (SED) fit of all the available optical to near-infrared bands.

We note that the limited angular resolution of the EoR-Spec renders insignificant any potential offset between the CO position and the optically-based position of galaxies in this external catalog, or any potential astrometry inaccuracies in the external catalog. Uncertainties on the angular position of our mock galaxies were thus not included in our mock external catalog. Moreover, given that the Δz associated with our frequency channels largely exceeds the anticipated photometric redshift accuracy of 0.002 × (1 + z) from Euclid (Euclid Collaboration 2021), uncertainties on the line-of-sight distance were also not included in our mock external catalog.

3.2. The masking array

The process of masking voxels containing bright CO contaminant is represented by a 3D array of zeros and ones, with zeros indicating voxels likely containing a bright CO emitter and thus requiring masking. This array is referred to as the binary masking array, B. Here, through empirical testing, we found that the most effective voxel size for masking is that of (FWHM_ν/3)×(FWHM_ang/3)². It is important to note that observed LIM tomographies undergo convolution with the beam of the telescope. This results in intensity from galaxies within a voxel spilling over to neighboring voxels. Additionally, sharp masks, like our binary masking arrays, produce high-frequency components in their Fourier transform that are not well captured when sampling this space discretely. This causes high-frequency details to alias into lower frequencies, distorting the Fourier-transformed data.

A remedy for both issues is to smooth the binary mask using a 3D Gaussian normalized such that its volume integral is equal to 1, a 3D equivalent of the 2D method used in CMB studies (Ponthieu et al. 2011; Kim 2011). While Ponthieu et al. (2011) recommends a Gaussian size double the map resolution, we have opted for a size consistent with the 3D beam of the EoR-Spec. Indeed, due to our broader binning in Fourier space, we do not face the same level of aliasing as CMB studies. This choice aligns with the inherent resolution of the instrument and addresses potential aliasing without significantly reducing the survey volume. However, as highlighted by Kim (2011), Gaussian smoothing can risk foreground contamination, that is, pixels initially zero in the binary mask might become non-zero after smoothing. To counteract this, we multiplied the smoothed mask with the original binary mask.

With the above considerations, we generated our mask as follows: we subtracted B from a unit array, then convolved the outcome with the normalized 3D beam of the EoR-Spec, G(RA,  Dec, ν). By subtracting the convolution result from another unit array and multiplying it by the original binary array, we obtained the “applied masking array”, M. Mathematically, this is expressed as:

$$\begin{array}{c}\hfill M=B\times [\mathbf{1}-(\mathbf{1}-B)\ast G(\mathrm{RA},\mathrm{Dec},\nu )],\end{array}$$(5)

where 1 is an array of ones of identical dimensions to B. Through this approach, the impact M has on the PS of a tomography mirrors the impact B has on the PS of a perfectly deconvolved map. This method bypasses deconvolution, a procedure that is not only computationally intensive but that can also amplify the noise of the observed data.

4. Results

4.1. Line intensities

To inform our masking strategy, we begin by analyzing the predicted intensities of the CO lines at the frequencies of the selected EoR-Spec bands (225 ± 20 GHz, 280 ± 20 GHz, 350 ± 20 GHz, and 410 ± 20 GHz), that is, comparing the mean values (Fig. 7) and CV (i.e., σ_I/μ_I; Fig. 8) of their intensities to those of the targeted line, [CII].

Fig. 7. Mean [CII] and CO line intensities as a function of the observed frequency (or emitted redshift for [CII]). The red dashed line depicts the [CII] fiducial model from Karoumpis et al. (2022). Continuous gray lines represent the prediction for the CO emission, with line shade decreasing as Jup (gray numbers) increases. The thicker black line indicates the mean total CO emission, while the Béthermin et al. (2022) prediction for total CO mean intensity (green line) is included for comparison.

Variations in the mean intensities of spectral lines at specific observed frequencies are primarily due to two factors: the different line luminosities of the galaxy populations emitting each line, and the different distances to each of these populations. In Fig. 7, we see a clear trend in the evolution of line mean intensity with observed frequency. As we shift towards lines with higher rest-frame frequency which therefore originate from higher redshifts, the gradient of their mean intensity as a function of frequency becomes markedly steeper. This steepening is a consequence of more substantial changes between the luminosity distances of galaxies of subsequent frequency channels when dealing with higher redshifts.

Notably the mean intensities of CO(3−2) and CO(4−3) exhibit significant fluctuations at frequencies above 300 GHz. This variability can be attributed to their origin from redshifts z <  0.4, where the volume encompassed by the observational cone is relatively small. Consequently, there is a limited number of emitting galaxies within this volume, leading to intensity fluctuations that are strongly driven by small number statistics. Furthermore, CO(3−2) cannot be observed at frequencies above 345.8 GHz, its rest frequency.

While the mean of the total CO intensity remains relatively constant at all frequencies (≈10³ Jy/sr) and is always dominated by the low-J CO transition from low redshift galaxies, the [CII] intensity shows a notable increase from about ≈10 Jy/sr at 225 GHz to ≈2000 Jy/sr at 410 GHz. We estimate that [CII] becomes the dominant line at > 364 GHz, a prediction that matches well with the results of Béthermin et al. (2022) which quotes a crossover at 358 GHz.

The variability in intensity distribution, as indicated by the CV, is also critical for the efficiency of our masking method. A high CV implies that the line intensity is concentrated in fewer voxels, while a lower CV indicates a more uniform voxel-to-voxel distribution. This aspect becomes increasingly relevant as [CII] signal recovery becomes more challenging at lower frequencies (i.e. higher redshifts).

The differences in CV arise from variations in the size of the galaxy populations contributing to the intensity of each line. This size is contingent on the survey volume, which expands with increasing redshift, and the number density of star-forming galaxies, which peaks around z ≈ 3 in the case of our mock observational cone (in agreement with observations like Conselice et al. 2016). Consequently, the CV typically decreases from z ≈ 0 to z ≈ 4 and then increases for higher redshifts, as shown in Fig. 8 in the case of the [CII] line. For the bright CO lines (i.e., J_up = 4 − 9, see Fig. 7), which are emitted by galaxies at z <  4, only the decreasing trend appears. On the other hand, for the CO lines with J_up >  10, CV is constant with redshift at > 280 GHz since they originate from z ≈ 2 − 3, where the survey volume expands slowly. Additionally, for these J_up, an increasing trend appears at lower frequencies (< 280 GHz) where these lines are emitted by galaxies at z >  4 where the number density of galaxies decreases rapidly.

Fig. 8. Coefficient of variation of the [CII] (red) and CO lines (gray) intensities for the four FYST frequency bands. The gray numbers of decreasing opacity denote the Jup of the CO transition. The secondary x-axis indicating redshift applies solely to [CII] emitters.

Just by comparing the mean and the CV of the line intensities and before even applying any mask, we observed a clear dichotomy. On the one hand, for the recovery of the [CII] emission originating from z <  5 galaxies, the situation is encouraging as the foregrounds have comparable intensity with the targeted signal and the dominating low-J line emission is scarcely distributed in the surveyed volume (i.e., their CV is high). On the other hand, for the recovery of the [CII] emission from z >  5 galaxies, the situation appears less favorable as the foregrounds are on average more than an order of magnitude brighter than the targeted signal and the dominating low-J line emission is homogeneously distributed in the surveyed volume (i.e., CV is low). In the following section, we explore and quantify how these differences impact the effectiveness of the masking technique in removing the CO signal. There, it becomes clear that the differences in the CV of the CO intensity result in high masking efficiency (significant CO reduction with a small fraction of the survey volume masked) for our low-redshift (high-frequency) tomographies, but low masking efficiency (poor CO reduction with a large fraction of the survey volume masked) for our high-redshift (low-frequency) tomographies.

4.2. Evaluating the masking technique efficiency

The effectiveness of a foreground removal technique can only be evaluated by its impact on the statistic of interest. Our study utilizes the spherically averaged PS, computed within the comoving volume of the [CII] reference frame, thus standardizing the analysis frame to avoid constant conversions between different lines. Contrasting with the approach in Paper I, here we compute the PS for tomographies convolved with the 3D beam of the EoR-Spec, rather than incorporating the beam attenuation into the uncertainty of the PS. This adjustment allows us to simulate the effect that the interplay between the transfer functions of the beam and the mask has on the measured PS. We are focusing on the PS at k = 0.02 − 0.32 Mpc⁻¹, a scale that is less affected by sample variance and masking bias. However, our findings hold true at higher k-bins (k = 0.32 − 0.62 Mpc⁻¹ and k = 0.62 − 0.92 Mpc⁻¹), suggesting that the effectiveness of the masking strategy is relatively constant across a broad range of scales.

Mitigating CO contamination while preserving maximal survey volume is far from straightforward due to the variety of possible CO masks and their nuanced effects on both CO and [CII] PS. To unravel these complex behaviors, we have adopted a three-step approach: initially focusing on CO alone to deduce the optimal masking sequence and strategy (Sect. 4.2.1), then examining the impact of masking on the [CII] signal in isolation (Sect. 4.2.2), and finally analyzing the combined effects to fully understand the interplay between masking, CO contamination, and the [CII] signal (Sect. 4.2.3). This comprehensive approach will not only enable us to identify the optimal masking depth, but also to assess the level of contamination at this depth.

4.2.1. Impact of masking on the CO PS

As already discussed in Sect. 3, for each frequency band and for each CO line we generated two distinct masking arrays: a “complete” masking array that masks all the galaxies in our external catalog with z_Jup = (ν_Jup/ν_obs ± 20 GHz)−1 and a “bright” masking array that among the galaxies of the “complete” subsample masks only those with ΔMS >  0.

Given the multitude of available masks, each with a unique impact on the observed data, identifying the optimal masking sequence becomes a critical challenge in our analysis. Therefore, it is important to establish a reliable criterion for determining the most effective sequence of masks. The effectiveness, E(J_up, bright or complete) of masking arrays, M(J_up, bright or complete), is quantified as,

$$\begin{array}{c}\hfill E=\frac{{\mathrm{Standard}\mathrm{deviation}}_{\mathrm{unmasked}}-{\mathrm{Standard}\mathrm{deviation}}_{\mathrm{masked}}}{\mathrm{Number}\mathrm{of}\mathrm{masked}\mathrm{voxels}}·\end{array}$$(6)

The reason for this definition of E is that we need a scale independent value that correlates with the amplitude of the PS per number of masked voxels, like the standard deviation of the intensity normalized by the number of voxels does.

The masking array with the highest E is applied iteratively to the CO-only tomography, continuing until all arrays have been utilized. With each iteration, we recalculate the effectiveness of subsequent masks on the updated data. This iterative process, performed on our CO-only simulations, identifies a sequence of masking arrays for each LIM tomography. It is important to note that this optimal sequence of CO masks remains unchanged at 350 ± 20 GHz and 410 ± 20 GHz if evaluated on a tomography including CO, [CII], and realistic white noise for the FYST LIM survey. In such tomographies, the masking efficiency (E) for each CO line is indeed still dominated by the CO signal. However, at 225 ± 20 GHz and 280 ± 20 GHz, where the [CII] signal and white noise become increasingly important, it is not practical to evaluate the optimal CO masking sequence directly from a realistic tomography including CO, [CII] and instrumental white noise. However, for all our tomographies, we observed that the optimal CO masking sequence does not vary much across the range of CO foreground models implemented in our work. This implies that the optimal CO masking sequences inferred here on our CO-only tomographies are applicable to first order to real observations. Any deviation from the “real” optimal sequence is unlikely to affect the result significantly.

The optimal CO masking sequence as well as their impact on the CO PS as a function of masking depth is depicted in Fig. 9 at k = 0.02 − 0.32 Mpc⁻¹. By construction of our optimal CO masking sequence, the slopes of the lines connecting the CO PS as a function of masking depth starts by having negative values that get progressively less negative as the masking depth increase and so does the masking effectiveness, E. This trend of progressively less negative slope is, however, not always observed in this initial phase due to fact that the definition of E is scale independent and thus not tailored to the specific k ≈ 0.02 − 0.32 Mpc⁻¹ explored in Fig. 9. Overall, this figure demonstrates that in this initial phase, the decrease of the CO PS with masking depth is mostly driven at the high frequency bands by the masking of the CO J_up = 3 (for 350 ± 20 GHz) and J_up = 4 line (for 350 ± 20 GHz and 410 ± 20 GHz) which results in a 50% reduction of the CO PS. At the low frequency bands (225 ± 20 GHz and 280 ± 20 GHz) the decrease is more gradual, yielding to a reduction of one order of magnitude by masking all J_up <  8 lines.

Fig. 9. Effect of masking on the CO-only PS. The range of PS between the “low-contamination” and “high-contamination” models (Sect. 2.3.2) at k = 0.02 − 0.32 Mpc−1 for CO lines (blue) is shown in relation to the percentage of voxels masked within the LIM tomographies for the frequency ranges of 225 ± 40 GHz (associated with z[CII] = 7.4), 280 ± 40 GHz (z[CII] = 5.8), 350 ± 40 GHz (z[CII] = 4.3), and 410 ± 40 GHz (z[CII] = 3.7). The lines represent the optimal CO masking sequence for our fiducial model, while each point and label indicate the Jup and subsample (“b” for “bright” and “c” for “complete”) associated with the masking array applied. The same masking sequence is applied to the “low-contamination” and “high-contamination” models (Sect. 4.2.1). As an illustration and assuming at this stage that the mask does not affect the [CII] PS (but see Sect. 4.2.2), we show the fiducial and range of [CII] PS predictions in red.

As the masking process intensifies, a critical threshold is reached where the amplitude of the CO PS reaches a minimum value (hereafter called the optimal masking depth) and after which it rises despite the increasing masking depth. This occurs when the CO line being masked contributes less to the total intensity of the tomography than the residual emissions from the brighter CO lines that have been masked previously. At this point, the mask begins to function akin to a random mask on the LIM tomography, as it is no longer targeting the brightest voxels. Moreover, since the intensity of the voxels follow a sparse, log-normal distribution (Sect. 4.2.3), this random masking results in the removal of more lower-intensity than high-intensity voxels. However, the correction we apply to the PS for the survey volume lost due to masking assumes that the intensity is uniformly distributed throughout the survey volume. This assumption of uniformity, although simplistic, is necessary because the non-homogeneity of the maps is model-dependent and cannot be accurately predicted without specific models. Nevertheless, it results in an increase in the amplitude of the PS when the mask does not follow the dominant galaxy population.

From Fig. 9, we note that the optimal masking depth as we move to higher frequency tomographies is reached at progressively lower masking depths (60%, 40%, 40% and 30% for the 225 GHz, 280 GHz, 350 GHz and 410 GHz accordingly). Comparing the amplitude of the CO PS at these optimal masking depths with that of the CII PS – assuming, at this stage, that it remains unaffected by the mask (but see Sect. 4.2.2 for further discussion) – allows us to preliminarily evaluate the feasibility of using this targeted masking approach to extract the CII PS signal. We find that in all frequency bands except for the 225 ± 40 GHz, the minimum fiducial CO PS is lower than the corresponding [CII] PS.

4.2.2. Impact of masking on the [CII] PS

In the analysis carried out so far, we have assumed that masking does not affect the PS of the [CII] signal. In this section, we assess the validity of this assumption across different scales.

By examining the influence of the masking on the [CII] PS we deal with an issue already discussed in the context of the CO residuals (Sect. 4.2.1). Due to the sparsity of the signal, the masking is artificially boosting the [CII] PS. Unlike the case of the CO emission, where this effect becomes prominent only after a certain threshold, boosting of the [CII] PS commences from the onset of the masking procedure as the [CII] emission is invariably not correlated with the mask. This boosting of the [CII] PS can henceforth be seen as a bias introduced by the masking process, that contributes additional uncertainty to our analysis. We defined this masking bias, b_mask, as the relative difference at a given scale between the PS amplitudes of the masked and unmasked [CII] tomographies, that is,

$$\begin{array}{c}\hfill b_{\mathrm{mask}}=\left|\frac{P_{\mathrm{masked}}-P_{\mathrm{unmasked}}}{P_{\mathrm{unmasked}}}\right|.\end{array}$$(7)

We evaluated this scale-dependent bias using two extreme models for the [CII] tomography, that is, our model with the highest (hereafter “Luminous”) and our model with the lowest [CII] PS signal (“Faint”; Fig. 10). Then, we evaluated these biases for three masking scenarios, the “light” (7%), “moderate” (21%), and “deep” (41%). In Fig. 10, we present these results for the 350 ± 20 GHz tomography.

Fig. 10. Masking bias representing the relative difference in the PS between masked and unmasked [CII] tomographies for the 350 ± 20 GHz frequency band. Comparisons are made for the most pessimistic (“Faint”, red lines) and optimistic (“Luminous”, black lines) [CII] PS predictions for the masking depths of 7%, 21%, and 41%.

For the “light”, “moderate”, and “deep” scenarios, we find biases of ≈3%, ≈10%, and ≈22% at k = 0.02 − 0.32 Mpc⁻¹ which increase to ≈6%, ≈15%, and ≈28% at k = 0.62 − 0.92 Mpc⁻¹, respectively. The scale dependence of the masking bias depends thus on the level of masking, and it intensifies as the masking becomes more aggressive. At the limiting case of of a masking depth of 41%, b_mask remains below 25% at large scales, but it can escalate to 30% at small scales. A similar scale-dependent bias is observed by Van Cuyck et al. (2023) in their angular [CII] PS predictions. This scale-dependent bias stems from the disparate impact of bright sources on the [CII] PS amplitude at different scales. On larger scales, the PS is predominantly influenced by the clustering signal of numerous star forming galaxies, which itself is proportional to the first moment of the luminosity function (Paper I). However, as we move to smaller scales, the PS is dominated by shot noise, which arises from the discrete nature of the luminous sources and is proportional to the second moment of the [CII] luminosity function (Paper I). At this scale the PS is thus strongly affected by the few, more luminous sources. The low probability of masking these rare bright sources leads to a significant masking bias due to volume corrections calculated assuming an homogeneous source distribution.

Inverting the masking process to correct for masking bias, that is, deconvolving the masking array with the masked intensity maps, would require uncertain modeling of clustering components in order to extrapolate the [CII] emission of the masked voxels. Furthermore, this approach could introduce artifacts when applied to a map with white noise, further complicating the PS analysis. However, it should be noted that this scale-dependent, and sometimes significant, masking bias does not vary greatly (< 5%) depending on the exact [CII] model used. This suggests that reliable scale-dependent corrections for this masking bias can be efficiently calculated using the specific masking strategies employed and [CII] models that are correct to first order. Of course, forward model fits to the observed [CII] PS will be able to account for this bias by construction.

4.2.3. Impact of masking on the CO-Contaminated [CII] PS

We now turn our attention to assessing the impact of masking on the CO-contaminated [CII] PS, P_{(CO + [CII])}, specifically within the k = 0.02 − 0.32 Mpc⁻¹ range (Fig. 11). The masks are applied with the optimal sequence inferred using the CO-only tomographies shown in Fig. 9 (Sect. 4.2.1). Building upon the analysis presented in Sect. 4.2.2, we include the bias introduced by masking on the [CII] signal, which translates into a gradual increase in the [CII] PS as masking depth increases. Finally, as a last piece of information, the uncertainty associated with the [CII] PS for a 2000-hour FYST survey, denoted as σ_[CII], is plotted in Fig. 11 as well. This uncertainty, encompassing both instrumental white noise and sample variance, is defined in Paper I by the equation:

$$\begin{array}{c}\hfill {\sigma }_{[\mathrm{CII}]}=\frac{P_{[\mathrm{CII}]}(k,z)+P_{\mathrm{N}}}{\sqrt{N_{\mathrm{m}}(k)}}·\end{array}$$(8)

Fig. 11. Same plot as Fig. 9 but with the masking bias of [CII] taken into account. We also show the [CII] uncertainty (σ[CII], white line) and the total line (CO+[CII]) PS for the fiducial models (black line).

Assuming that the distribution of masked voxels is uniform, this uncertainty increases with the masking depth proportionally to the decrease in the number of k-modes across all k-bins. A Monte Carlo algorithm validated this uniformity assumption, establishing a consistent relationship across all scales we examined:

$$\begin{array}{c}\hfill \frac{\sigma [\mathrm{CII}]}{\sigma [\mathrm{CII}],\mathrm{masked}}=\sqrt{\frac{N_{k,\mathrm{not}\mathrm{masked}}}{N_k}}=\sqrt{\frac{N_{\mathrm{voxels},\mathrm{not}\mathrm{masked}}}{N_{\mathrm{voxels}}}}·\end{array}$$(9)

By integrating all the elements presented in Fig. 11, we can assess whether the objective of masking has been achieved, that is, to reduce the amplitude of the CO PS below the measurement uncertainty of the [CII] PS (σ_[CII]), while ensuring that this uncertainty remains below the amplitude of the [CII] PS. Two distinct cases are evident in Fig. 11. Firstly, in the higher frequency bands (350 ± 20 and 410 ± 20 GHz), P_{(CO + [CII])} initially undergoes a slight decrease at masking depths of < 10%, attributed to the masking of the rare but bright J_up = 3 and J_up = 4 sources, immediately followed by a gradual increase, mainly due to the [CII] PS masking bias (the CO PS masking bias only playing a role at masking depths > 20%). As already mentioned, since the [CII] emission is invariably uncorrelated with the mask, such an increase in PS from the onset of the masking procedure is characteristic of the dominance of the [CII] signal in the unmasked tomography. At these high-frequency and at the optimal masking depth of < 10%, the [CII] signal dominates and, more importantly, it is well above σ_CII. This guarantees accurate detection of the [CII] PS by the FYST LIM survey in these bands. Secondly, in the lower frequency bands (225 ± 20 and 280 ± 20 GHz), P_{(CO + CII)} decreases progressively with masking depth, reaches a minimum at the optimal masking depth of ≈40 − 60% and finally increases due to a combination of CO and [CII] masking biases. At 280 ± 20 GHz, the decrease in P_{(CO + CII)} is less pronounced than at 225 ± 20 GHz, and in the former case, the masked tomography effectively switches from CO-dominant to [CII]-dominant, whereas in the latter case, the masked tomography remains CO-dominant. At 280 ± 20 GHz, and at the optimal masking depth of 40%, the [CII] signal dominates and lies above σ_CII. FYST LIM will thus detect, albeit marginally, the [CII] PS in this band (signal-to-noise ratio of two and CO contamination of around 15%). In contrast, at 225 ± 20 GHz and at optimal masking depth of 60%, the CO signal still dominates. In this band, an alternative masking approach will be required.

From the observational point of view, the results of Fig. 11 make clear that for a real dataset, the evolution of the PS as function of the masking depth and the location of the optimal masking depth provide key insights into the initial contamination levels of the tomography. A very low value for the optimal masking depth (≲10%), associated with an amplified PS as we surpass this threshold, unambiguously indicates the predominance of [CII] in the initial tomography. In contrast, a value of the optimal masking depth greater than ≳10% indicates the dominance by CO foregrounds in the initial tomography. In this case, accurately assessing the extent of CO contamination at the optimal masking depth is challenging in a real dataset as no clear signatures are visible in the PS and as the [CII] and CO emitters share similar spatial frequency characteristics, complicating their differentiation in Fourier space.

Separation between the [CII] and CO emitters becomes clearer when analyzing how these populations contribute to the average total power of the signal. Specifically, this can be achieved by examining their contributions to the zero displacement of the PS from constructing a histogram that represents the distribution of squared voxel intensities, each weighted by its own squared intensity (Fig. 12). It is crucial to note, however, that in order to create this histogram, knowledge of the absolute intensity of the tomography is needed, while this information is partially lost by our atmospheric and astrophysical continuum foregrounds removal. To estimate a posteriori the magnitude of the mean intensity, we used the PS value at large scales, since we demonstrated in Paper I that k³P_clustering(k ≈ 10⁻² Mpc⁻¹)≈I_mean².

Fig. 12. Pixel distribution histograms weighted with the voxel intensity squared, illustrating the distinction between successful (left; fiducial CO) and unsuccessful (right; pessimistic CO) cases of CO masking of the FYST 280 ± 20 GHz tomography. The top panels show the weighted voxel histograms of these 280 ± 20 GHz tomography masked to the optimal depth and fitted double lognormal distributions. The lower panels detail the progressive masking depth, each curve representing the lognormal fits of the weighted histogram at increasing depth. This highlights the emergence of a bimodal distribution in the successful case versus a single peak distribution in the unsuccessful case. The lognormal fit to the weighted histogram of a white noise tomography equivalent to the upcoming 2000 hours FYST LIM (black dashed line) and a fit to the white noise tomography of a survey 5 times deeper (gray dashed line) are superimposed.

Equipped with this proxy for the mean intensity, we can generate a weighted histogram of the tomography, masked to the optimal depth. When this histogram displays two peaks instead of one, it signifies a transition from a CO-dominated to a [CII]-dominated tomography, achieved by eliminating the brighter CO sources. Conversely, the presence of a single peak signifies that the tomography remains CO-dominated, primarily by the less intense CO sources that were not masked in earlier iterations. Fig. 12 illustrates these two cases for the 280 ± 20 GHz tomography. For the first case we considered our fiducial [CII] and CO PS whereas for the second we combined our fiducial [CII] with our most pessimistic CO PS predictions. For the first case of the fiducial CO PS, as we increase the masking depth, a secondary peak emerges at the brighter end of the histogram, corresponding to the population of [CII] emitters. Notably, at the optimal masking depth, these two peaks have similar heights, marking the transition from CO dominance to [CII] dominance in the PS of this masked tomography. This emergence and predominance of this secondary peak thus signifies the success of the masking process at the optimal masking depth. On the contrary, in the second case, even at the optimal masking depth, the [CII] peak remains hidden in the bright tail of the CO peak. This single-peak histogram, which moves only towards lower voxel intensities as the masking depth increases, signifies the failure of the masking process at the optimal masking depth.

Unfortunately, this rather simple method to distinguish CO-dominated from CII-dominated tomography at the optimal masking depth is strongly affected by the presence of white noise. For example, considering the white noise of the plane 2000 hours of the FYST LIM survey, the peak(s) of the weighted histogram is entirely hidden by it (Fig. 12). Nevertheless, with five times deeper observations, the white noise distribution moves enough towards fainter voxels ($I_{\mathrm{vox}}^2<{10}^4$ (Jy/sr)²) allowing us to identify the emergence of the [CII]-related secondary peak at $I_{\mathrm{vox}}^2>{10}^4$ (Jy/sr)². This white noise distribution effectively acts as a threshold on the x-axis, below which observations become indiscernible.

In essence, our technique unfolds in two steps. We begin by incrementally masking the LIM tomography until we identify the point where the PS at k = 0.02 − 0.32 Mpc¹ is minimized, revealing the optimal depth for masking. Following this, we examine the weighted voxel histogram of the tomography masked to this depth. In this histogram –which is calibrated using the PS measurement and weighted by the square of the voxel intensities– we look for a bimodal distribution as confirmation of successful [CII] signal isolation.

5. Discussion

This study aimed to quantify the contamination of the (post-) EoR [CII] LIM signal by foreground CO emissions for the upcoming FYST LIM survey, utilizing post-processing of the TNG300 hydrodynamical simulation. The strength of our approach lies in its capacity to generate a range of predictions that account for the inherent uncertainties in CO luminosity functions (Riechers et al. 2019; Decarli et al. 2020), and for variations in the CO SLED of galaxies. Based on these predictions, we assessed the effectiveness of a masking technique that leverages an external catalog to identify and subsequently mask bright CO emitters. This led to the identification of a method to determine the optimal masking depth and a novel indicator for successful CO foreground removal.

Our analysis revealed significant variations in CO contamination across different frequency bands of the EoR-Spec, with notably lower contamination at higher frequencies. This variation significantly impacts the efficiency of our masking technique, which shows greater effectiveness at frequencies where initial contamination is lower. Specifically, we found that the targeted [CII] signal could be detected with minimal masking and thus minimal masking bias at 410 GHz and 350 GHz (Fig. 13). In contrast, at 280 GHz and particularly at 225 GHz, our simple CO masking technique faced challenges due to the high level of foreground CO emissions, underscoring the need for alternative strategies at these frequencies.

Fig. 13. Comparison of the [CII] PS for [CII] observations in the redshift range z[CII] = 4.1 − 4.8. The squares denote the combined CO+[CII] PS after masking CO bright voxels, resulting in a 7% reduction in survey volume. The associated error bars on the squares capture the uncertainty in [CII] measurements due to white instrumental noise and sample variance. The blue crosses depict the intrinsic [CII] PS results, binned at the same Δk. The continuous line represents the intrinsic [CII] PS results but with a finer binning.

Our work aligns with the efforts of Van Cuyck et al. (2023) performed in the framework of the CONCERTO experiment, which similarly aimed to test via simulation the extraction of the [CII] signal from the (post-)EoR universe using a comparable masking strategy. Despite the challenges in making direct comparisons due to differences in angular and spectral resolution, both studies highlight the limitations of masking strategies at high redshifts, especially beyond z ≥ 6.5. This consensus emphasizes the complexities of masking strategies and their variable effectiveness across redshift. Diverging from Van Cuyck et al. (2023), we chose to exclude the [CI] line from our simulations due to its highly uncertain intensity values, aligning with methodologies like those suggested by Yue & Ferrara (2019). The absence of [CI] contamination in our simulations is probably largely compensated by the very wide range of CO models, including very pessimistic ones, taken into account in our study. The precise impact of [CI] lines, along with other foregrounds, is expected to be clarified with the advent of first light science observations.

A notable limitation of our study, also highlighted by Gkogkou et al. (2023) and Van Cuyck et al. (2023), concerns the impact of sample variance, of both the CO and [CII] emitters. Our reliance on the single realization of the TNG300 simulation introduces a potential limitation in terms of the accuracy of our results. Nonetheless, given the wide range of models considered in our simulation, we assume that the influence of sample variance on the effectiveness of our masking technique is relatively minor. Nevertheless, recognizing that sample variance could pose a significant challenge in the context of future forward modeling of the [CII] LIM observations, we will explore the development of multiple realizations of wide observational cones by employing machine learning techniques (Garcia et al. 2023).

Our decontamination method itself comes with its limitations. Despite its reliability, informed masking is not the most efficient technique, particularly in terms of survey volume loss, at low frequencies. The challenge remains to develop more refined techniques, that de-confuse the line PS without the cost of losing survey volume. There are several mathematically elegant techniques that could achieve this goal. One such approach would be utilizing the wide frequency coverage of the EoR-Spec to measure the CO PS indirectly from the measurements of the cross-correlation of three lines coming from the same galaxies, leveraging the statistical independence of the rest of the line-components. Another approach involves exploiting the anisotropic nature of the 2D PS of interloper lines (Cheng et al. 2016). Upon projection on the frame of reference of the targeted line, the PS of these interloper lines exhibits distinct anisotropies that can be used to differentiate between the target signal and the foreground contamination. Both approaches are promising, yet they share the limitation of failing to recover the valuable Fourier phase information, which is essential for analyzing histograms related to voxel intensity distribution. A third approach that recaptures the Fourier phase information is using conditional generative adversarial networks to de-confuse the signals at the map level rather than at the PS level. Moriwaki et al. (2020) focused on de-confusing two lines (Hα and OIII), using a large set of simulated data for training their models. Although this approach is expected to be very effective, it is challenging to estimate the uncertainties involved in the de-confusion method of the measurements. Finally, it is worth noting the existence of a fourth approach that is more robust against the model uncertainty and still recaptures the phase information and that relies on the identification and extraction of foreground sources that are characterized by multiple interloping emission lines. This can be done by fitting multiple SLED templates to each line-of-sight of the tomography using the sparse approximation (Cheng et al. 2020). Given the wide frequency range of the EoR-Spec, we consider this method to be the most promising alternative approach beyond targeted masking. This approach will be tested and evaluated in the context of the FYST LIM survey on the third part of our paper series.

While line foregrounds present significant challenges, the importance of atmospheric and astrophysical continuum foregrounds should not be overlooked. Although Van Cuyck et al. (2023) successfully applied the asymmetric re-weighted penalized least-squares technique to remove contamination from the cosmic infrared background, reinforcing our focus on line interlopers as the primary challenge at low frequencies (225 ± 20 GHz, 280 ± 20 GHz), the issue of atmospheric foregrounds remains unresolved and may only be addressed through the use of mock detector data streams. Unfortunately, it is at the high frequencies (350 ± 20 GHz, 410 ± 20 GHz), where CO contamination is not problematic that we expect atmospheric contamination to play a critical role. This juxtaposition makes atmospheric interference potentially the primary barrier to detection at the two higher frequency bands of FYST LIM survey.

6. Conclusion

Our study evaluated the level of contamination from CO rotational line emission when observing the (post-)EoR [CII] PS using the upcoming FYST LIM survey and presented a method for mitigating it. We created mock tomographies of the CO line emitters by post-processing the TNG300 simulation, accounting for the inherent uncertainties in the observed CO luminosity functions and the observed CO SLEDs. The resulting range of CO predictions is in good agreement with the ASPECS observations.

Our analysis revealed that the total mean intensity of CO emissions remains relatively constant, at approximately 1000 Jy/sr, across the observed frequency spectrum. This emission predominantly originates from low-J CO transitions (J_up <  8) in galaxies at lower redshifts (z <  2). In contrast, the [CII] mean intensity exhibits a significant upward trend, increasing from about 10 Jy/sr at 225 GHz to 2000 Jy/sr at 410 GHz, making the [CII] the dominant emission line at frequencies above 364 GHz.

Based on these predictions, we assessed a foreground removal strategy that involves identifying and masking voxels containing bright CO emitters, utilizing mock external catalogs that exhibit realistic features. A key aspect of our analysis was the determination of mask effectiveness, which we defined as the ability of a mask to reduce the standard deviation of the intensity of the signal, normalized by the number of voxels masked. Through iterative application and reassessment of these masks on CO-only tomography, we were able to refine our masking sequence to achieve optimal foreground removal. The optimal sequence for CO masks found here is broadly applicable to real observations, as we see very little variation into this sequence across the broad range of CO model considered here.

Our analysis also highlighted the frequency dependency of the optimal masking depth, which aligns with the point where the PS of the masked tomography reaches its minimum. While less than 10% masking depth is required at higher frequencies (350 ± 20 GHz and 410 ± 20 GHz), a deeper masking of more than 40% is necessary at lower frequencies (280 ± 20 GHz and 225 ± 20 GHz).

At these masking depths, according to our fiducial models, [CII] emission dominates and significantly exceeds the typical instrumental white noise levels at higher frequencies, ensuring a robust detection of the [CII] PS at 350 ± 20 GHz and 410 ± 20 GHz. A tentative detection is also expected at 280 ± 20 GHz, where [CII] marginally prevails over CO contamination. However, at 225 ± 20 GHz, despite deep masking, CO contamination continues to overshadow [CII] emission, highlighting the urgent need for advanced decontamination techniques.

While estimating the optimal masking depth for decontaminating a real LIM tomography is relatively straightforward, discerning the levels of decontamination achieved at this depth presents a more complex challenge. This complexity arises because CO and [CII] emitters share similar spatial frequency characteristics, making it difficult to distinguish between the two solely based on their distribution in Fourier space. A practical solution involves creating a weighted voxel histogram, where each voxel is weighted by the square of its intensity and the histogram is calibrated using the PS measurement of the tomography masked at the optimal depth. Detailed examination of this histogram for signs of bimodality can reveal a successful transition from a CO-dominated to a [CII]-dominated state. This method will, however, only be applicable to data with very low instrumental white noise, about five times lower than that anticipated for 2000 hours of observation with the EoR-Spec.

Our study lays the groundwork for future LIM studies into the high-redshift universe. The development of decontamination techniques capable of extracting the [CII] high-redshift signal as well as the detailed modeling of the astrophysical and atmospheric foregrounds suitable for forward modeling are the necessary next steps.

Acknowledgments

We thank Dongwoo Chung, Patrick Breysse, Yoko Okada, Anirban Roy, Steve Choi, Eiichiro Komatsu, Thoma Nikola, Reinhold Schaaf, Ralf Antonius Timmermann, Kaustuv moni Basu, Vyoma Muralidhara, Maude Charmetant, Nicholas Battaglia, Gordon Stacey and other members of the CCAT science working group for detailed discussions about EoR-Spec; We are grateful to Matthieu Bethermin, Athanasia Gkogkou, Mathilde Van Cuyck, Desika Narayanan, Anthony Pullen, Sylvia Adscheid, Elena Marcuzzo, Prachi Khatri, Benedetta Spina, and Cristiano Porciani for the helpful discussions. This research was carried out within the Collaborative Research Center 956 (project ID 184018867), sub-projects A1 and C4 and within the Collaborative Research Center 1601 (project ID 500700252), sub-projects C3 and C6, funded by the Deutsche Forschungsgemeinschaft (DFG). This research made use of NASA’s Astrophysics Data System Bibliographic Services; Matplotlib (Hunter 2007); Astropy, a community-developed core Python package for astronomy (Astropy Collaboration 2013); PlotDigitizer (http://plotdigitizer.sourceforge.net/).

References

Appendix A: Impact of masking on the PS of individual CO lines

This appendix presents the PS of individual CO (J_up = 3 − 8) spectral lines in the frame of reference of the [CII] line at scales of k = 0.02 − 0.32 Mpc⁻¹ as a function of the percentage of the survey volume masked (Fig. A.1 and A.2) following the optimal masking sequence described in Sect. 4.2. For each CO line PS, two distinct behaviors are observed during the masking process. When the masking step specifically targets the spectral line in question, there is a sharp decrease in the PS. Conversely, when the masking step targets a different spectral line, a slow increase in the PS is observed due to the masking bias (see Sect. 4.2.2).

The behavior of the PS during masking is influenced by how the galaxies emitting the lines are distributed within the 3D tomography. Lines from galaxies that are highly localized within the survey volume, such as CO(3-2) at 350 GHz, exhibit a sharp decrease in their PS when masking targets them. Additionally, there is almost no increase in PS when masking targets other lines due to the minimal residuals left after masking. In contrast, lines from galaxies more evenly distributed throughout the survey volume, such as CO(4-3) at 225 GHz, show a less pronounced decrease in PS when masked, but a stronger increase when other lines are masked, because there remains a substantial amount of residual emission that gets amplified. A key indicator of how uniformly the intensity of each line is distributed in the tomography is the CV of its intensity, as shown in Fig. 8 and described in detail in Sect. 4.1.

Fig. A.1. Effect of masking on the PS of the brightest individual CO lines (Jup = 3 − 8). The lines represent our fiducial model, with the same optimal masking sequence applied as in Fig. 9 and Fig. 11. The results are plotted for the frequency ranges 225 ± 40 GHz (corresponding to z[CII] = 7.4) and 280 ± 40 GHz (z[CII] = 5.8).

Fig. A.2. Same as Fig. A.1, but for the frequency ranges 350 ± 40 GHz (z[CII] = 4.3) and 410 ± 40 GHz (z[CII] = 3.7).

All Tables

All Figures

Fig. 1. Spectral lines observable within four selected EoR-Spec bands (illustrated by red-shaded regions), and originated from galaxies situated at different redshifts. The dashed numbered lines represent transitions within the CO rotational ladder, with the corresponding numbers indicating the Jup for each specific transition. These CO lines act as the foreground contaminant for the [CII] (red thin dashed line) LIM survey.

In the text

Fig. 2. Redshift evolution of the cosmic molecular gas mass density of galaxies. The blue horizontal line represents predictions from our cone of mock galaxies, while observational constraints are presented by rectangles. The gray rectangles are for the dust-based estimates of Magnelli et al. (2020), the light gray rectangles are for the COLDz CO survey (Riechers et al. 2019), the light blue rectangle for the CO-based estimates of Boogaard et al. (2023), the red-shaded rectangles are for the ASPECS CO pilot survey (Decarli et al. 2020), while the pink rectangles are for the ASPECS CO LP survey (Decarli et al. 2020). The width of the rectangles indicates the redshift bin size and the height the 1-σ confidence region.

In the text

	Fig. 3. CO luminosity function for various transitions and redshifts. The blue line corresponds to our fiducial model, while the shaded region indicates the range resulting from our optimistic and pessimistic models. Rectangles represent observational data from the ASPECS LP survey (Decarli et al. 2020), with the width indicating the bin size and the height indicating the 1-σ confidence region.
In the text

	Fig. 4. CO SLED of the Milky Way (black circles) and NGC 253 (black squares). The SLEDs resulting from the linear combination of the two, following Eq. (3), are shown by solid lines, with μ varying from 0 to 1 in steps of 0.1. The red line represents the combination (with μ = 0.4) that best fit the BzK galaxies (red points) and is used here as reference for high-redshift MS galaxies.
In the text

	Fig. 5. SFR as a function of stellar mass with the bins color-coded to reflect the cumulative CO (6−5) luminosity of the 2D bin. The dashed line indicates the thresholds that define our “bright” subsample with ΔMS >  0 (i.e., galaxies on and above the MS).
In the text

	Fig. 6. Comparison of the stellar mass completeness limits applied to our mock external catalog (red line) with those of the COSMOS2015 (light gray line Laigle et al. 2016) and COSMOS2020 (dark gray line Weaver et al. 2022) catalogs. The stellar mass completeness limits of the deep Euclid catalogs is expected to lie between these two lines.
In the text

Fig. 7. Mean [CII] and CO line intensities as a function of the observed frequency (or emitted redshift for [CII]). The red dashed line depicts the [CII] fiducial model from Karoumpis et al. (2022). Continuous gray lines represent the prediction for the CO emission, with line shade decreasing as Jup (gray numbers) increases. The thicker black line indicates the mean total CO emission, while the Béthermin et al. (2022) prediction for total CO mean intensity (green line) is included for comparison.

In the text

	Fig. 8. Coefficient of variation of the [CII] (red) and CO lines (gray) intensities for the four FYST frequency bands. The gray numbers of decreasing opacity denote the Jup of the CO transition. The secondary x-axis indicating redshift applies solely to [CII] emitters.
In the text

Fig. 9. Effect of masking on the CO-only PS. The range of PS between the “low-contamination” and “high-contamination” models (Sect. 2.3.2) at k = 0.02 − 0.32 Mpc−1 for CO lines (blue) is shown in relation to the percentage of voxels masked within the LIM tomographies for the frequency ranges of 225 ± 40 GHz (associated with z[CII] = 7.4), 280 ± 40 GHz (z[CII] = 5.8), 350 ± 40 GHz (z[CII] = 4.3), and 410 ± 40 GHz (z[CII] = 3.7). The lines represent the optimal CO masking sequence for our fiducial model, while each point and label indicate the Jup and subsample (“b” for “bright” and “c” for “complete”) associated with the masking array applied. The same masking sequence is applied to the “low-contamination” and “high-contamination” models (Sect. 4.2.1). As an illustration and assuming at this stage that the mask does not affect the [CII] PS (but see Sect. 4.2.2), we show the fiducial and range of [CII] PS predictions in red.

In the text

	Fig. 10. Masking bias representing the relative difference in the PS between masked and unmasked [CII] tomographies for the 350 ± 20 GHz frequency band. Comparisons are made for the most pessimistic (“Faint”, red lines) and optimistic (“Luminous”, black lines) [CII] PS predictions for the masking depths of 7%, 21%, and 41%.
In the text

	Fig. 11. Same plot as Fig. 9 but with the masking bias of [CII] taken into account. We also show the [CII] uncertainty (σ[CII], white line) and the total line (CO+[CII]) PS for the fiducial models (black line).
In the text

Fig. 12. Pixel distribution histograms weighted with the voxel intensity squared, illustrating the distinction between successful (left; fiducial CO) and unsuccessful (right; pessimistic CO) cases of CO masking of the FYST 280 ± 20 GHz tomography. The top panels show the weighted voxel histograms of these 280 ± 20 GHz tomography masked to the optimal depth and fitted double lognormal distributions. The lower panels detail the progressive masking depth, each curve representing the lognormal fits of the weighted histogram at increasing depth. This highlights the emergence of a bimodal distribution in the successful case versus a single peak distribution in the unsuccessful case. The lognormal fit to the weighted histogram of a white noise tomography equivalent to the upcoming 2000 hours FYST LIM (black dashed line) and a fit to the white noise tomography of a survey 5 times deeper (gray dashed line) are superimposed.

In the text

Fig. 13. Comparison of the [CII] PS for [CII] observations in the redshift range z[CII] = 4.1 − 4.8. The squares denote the combined CO+[CII] PS after masking CO bright voxels, resulting in a 7% reduction in survey volume. The associated error bars on the squares capture the uncertainty in [CII] measurements due to white instrumental noise and sample variance. The blue crosses depict the intrinsic [CII] PS results, binned at the same Δk. The continuous line represents the intrinsic [CII] PS results but with a finer binning.

In the text

	Fig. A.1. Effect of masking on the PS of the brightest individual CO lines (Jup = 3 − 8). The lines represent our fiducial model, with the same optimal masking sequence applied as in Fig. 9 and Fig. 11. The results are plotted for the frequency ranges 225 ± 40 GHz (corresponding to z[CII] = 7.4) and 280 ± 40 GHz (z[CII] = 5.8).
In the text

	Fig. A.2. Same as Fig. A.1, but for the frequency ranges 350 ± 40 GHz (z[CII] = 4.3) and 410 ± 40 GHz (z[CII] = 3.7).
In the text