# CONCERTO: Extracting the power spectrum of the [C II] emission line

M. van Cuyck, N. Ponthieu, G. Lagache, A. Beelen, M. Béthermin, A.

Gkogkou, M. Aravena, A. Benoit, J. Bounmy, M. Calvo, et al.

► To cite this version:

M. van Cuyck, N. Ponthieu, G. Lagache, A. Beelen, M. Béthermin, et al. CONCERTO: Extracting the power spectrum of the [C II] emission line. Astronomy and Astrophysics - A&A, 2023, 676, pp.A62. 10.1051/0004-6361/202346270. hal-04128387

HAL Id: hal-04128387

https://hal.science/hal-04128387

Submitted on 24 Aug 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# CONCERTO: Extracting the power spectrum of the
[$C_{\text{II}}$] emission line

M. Van Cuyck<fnref n="1" />, N. Ponthieu<fnref n="2" />,<fnref n="1" />, G. Lagache<fnref n="1" />, A. Beelen<fnref n="1" />, M. Béthermin<fnref n="1" />,<fnref n="3" />, A. Gkogkou<fnref n="1" />, M. Aravena<fnref n="4" />,
A. Benoit<fnref n="5" />, J. Bounmy<fnref n="6" />, M. Calvo<fnref n="5" />, A. Catalano<fnref n="6" />, F. X. Désert<fnref n="2" />, F.-X. Dupé<fnref n="1" />, A. Fasano<fnref n="1" />, A. Ferrara<fnref n="7" />, J. Goupy<fnref n="5" />,
C. Hoarau<fnref n="6" />, W. Hu<fnref n="1" />, J.-C Lambert<fnref n="1" />, J. F. Macías-Pérez<fnref n="6" />, J. Marpaud<fnref n="6" />,
G. Mellema<fnref n="8" />, A. Monfardini<fnref n="5" />, and A. Pallottini<fnref n="7" />

Received 28 February 2023 / Accepted 20 May 2023

## ABSTRACT

**Context.** CONCERTO is the first experiment to perform a [CII] line intensity mapping (LIM) survey on the COSMOS field to target $z > 5.2$. Measuring the [CII] angular power spectrum allows us to study the role of dusty star-forming galaxies in the star formation history during the epochs of Reionization and post-Reionization. The main obstacle to this measurement is the contamination by bright foregrounds: the dust continuum emission and atomic and molecular lines from foreground galaxies at $z \lesssim 3$.

**Aims.** We evaluate our ability to retrieve the [CII] signal in mock observations of the sky using the Simulated Infrared Dusty Extra-galactic Sky (SIDES), which covers the mid-infrared to millimetre range. We also measure the impact of field-to-field variance on the residual foreground contamination.

**Methods.** We compared two methods for dealing with the dust continuum emission from galaxies (i.e. the cosmic infrared background fluctuations): the standard principal component analysis (PCA) and the asymmetric re-weighted penalized least-squares (arPLS) method. For line interlopers, the strategy relies on masking low-redshift galaxies using the instrumental beam profile and external catalogues. As we do not have observations of CO or deep-enough classical CO proxies (such as $L_{IR}$), we relied on the COSMOS stellar mass catalogue, which we demonstrate to be a reliable CO proxy for masking. To measure the angular power spectrum of masked data, we adapted the P of K EstimatoR (POKER) from cosmic infrared background studies and discuss its use on LIM data.

**Results.** The arPLS method achieves a reduction in the cosmic infrared background fluctuations to a sub-dominant level of the [CII] power at $z \sim 7$, a factor of >70 below our fiducial [CII] model. When using the standard PCA, this factor is only 0.7 at this redshift. The masking lowers the power amplitude of line contamination down to $2 \times 10^{-2}$ Jy$^2$ sr$^{-1}$. This residual level is dominated by faint undetected sources that are not clustered around the detected (and masked) sources. For our [CII] model, this results in a detection at $z = 5.2$ with a power ratio $[C_{\text{II}}]/(\text{residual interlopers}) = 62 \pm 32$ for a 22% area survey loss. However, at $z = 7$, $[C_{\text{II}}]/(\text{residual interlopers}) = 2.0 \pm 1.4$, due to the weak contrast between $[C_{\text{II}}]$ and the residual line contamination. Thanks to the large area covered by SIDES-Uchuu, we show that the power amplitude of line residuals varies by 12–15% for $z = 5.2–7$, which is less than the field-to-field variance affecting $[C_{\text{II}}]$ power spectra.

**Conclusions.** We present an end-to-end simulation of the extragalactic foreground removal that we ran to detect the [CII] at high redshift via its angular power spectrum. We show that cosmic infrared background fluctuations are not a limiting foreground for $[C_{\text{II}}]$ LIM. On the contrary, the CO and [CII] line contamination severely limits our ability to accurately measure the [CII] angular power spectrum at $z \gtrsim 7$.

**Key words.** galaxies: star formation – galaxies: high-redshift – galaxies: ISM – large-scale structure of Universe – cosmology: observations

## 1. Introduction

The Epoch of Reionization (EoR) started with the birth of the first stars and the formation of the first galaxies, whose ultra-violet light ionized the neutral hydrogen of the intergalactic medium. As time went by, the bubbles of ionized gas increased in number, expanding and merging, eventually reionizing the whole Universe. While we are beginning to understand the

overall chronology of the EoR, the full story of reionization remains to be written. However, recent and planned observational advances promise progress in this respect over the coming decade.

The Planck satellite has provided an accurate measurement of the Thomson scattering optical depth, $\tau$, using the cosmic microwave background (CMB) angular power spectrum (APS; Planck Collaboration Int. XLVII 2016). It shows that the

<footnote n="1">¹ Aix-Marseille Univ., CNRS, CNES, LAM, Marseille, France
e-mail: mathilde.van-cuyck@lam.fr</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="4">⁴ Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Univ. Grenoble Alpes, CNRS, LPSC/IN2P3, 38000 Grenoble, France</footnote>

<footnote n="7">⁷ Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy</footnote>

<footnote n="8">⁸ The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova 10691 Stockholm, Sweden</footnote>
Universe was reionized at 10% at $z_{10\%} = 10.4 \pm 1.8$. On the other end of the reionization epoch, studies of the ionized intergalactic medium through quasar and gamma-ray burst absorption in the Lyman-$\alpha$ and Lyman-$\beta$ hydrogen transitions set the late end at $z \approx 6$ (e.g., Chornock et al. 2013; Eilers et al. 2019; Wang et al. 2020; Bosman et al. 2022). As of today, the total ionizing photon budget required to reionize the Universe is still not probed by the deep surveys that have constrained the faint-end slope of the ultraviolet luminosity function in the early Universe (e.g., Ellis et al. 2013; Oesch et al. 2015; Carilli et al. 2016; Livermore et al. 2017; Drake et al. 2017). This may change in the near future with deep observations from the James Webb Space Telescope (JWST).

The CarbON CII line in post-rEionization and ReionizaTiOn (CONCERTO) project (CONCERTO Collaboration 2020) seeks to constrain the dust-obscured star formation history and its spatial distribution at high redshifts ($z > 5.2$). To that aim, CONCERTO is the first experiment to conduct a line intensity mapping (LIM) survey of the $[C_{II}]$ line at the end of the EoR and during the post-reionization era. A LIM survey consists of mapping the surface brightness in a given field of view as a function of position and observed frequency (see Kovetz et al. 2017 and Bernal & Kovetz 2022 for a review). Targeting an emission line then allows the redshift information to be recovered. Line intensity mapping is naturally sensitive to faint sources, missed by classical galaxy surveys, that potentially contribute the most to the ionizing photon budget (e.g., Finkelstein et al. 2012; Dayal et al. 2013; Atek et al. 2015). The CONCERTO project then studies surface brightness fluctuations in the Fourier space using the APS.

The $[C_{II}]$ 157.7 µm line is a line of choice for LIM at high redshifts. It is one of the brightest far-infrared lines and is unattenuated by the dust content of galaxies. Furthermore, it has been shown by both observations (e.g., Gullberg et al. 2015; Schaerer et al. 2020) and theoretical works (e.g., Pallottini et al. 2017; Lagache et al. 2018) that $[C_{II}]$ luminosity is a reliable tracer of star formation in the redshift range of interest. $[C_{II}]$ emission from high-redshift galaxies is dominated by the emission from photo-dissociation regions, which are regions of cool neutral gas penetrated by ionizing ultraviolet radiations from hot young stars (Pineda et al. 2014) that surround molecular clouds (Vallini et al. 2017; Pallottini et al. 2022). In addition, the transparent millimetre atmospheric window enables access from the ground to the $[C_{II}]$ emission originating at $z \gtrsim 4$. The main challenge for any kind of LIM experiment remains the foreground contamination, which can completely dominate the signal of interest. For the CONCERTO $[C_{II}]$ survey, the brightest foreground comes from fluctuations in the cosmic infrared background (CIB). The CIB is the line-of-sight-integrated emission of all dusty galaxies over the entire history of the Universe, from 8 to 1000 µm (Lagache et al. 2000). Further, there is a second foreground due to interloping emission lines, which are atomic and molecular emission lines redshifted in the same frequency band. For the $[C_{II}]$ survey in the CONCERTO frequency range, the main interlopers are the CO rotational ladder and the two $[C_{I}]$ fine structure lines from foreground galaxies at $z \lesssim 3$. In this paper we present in detail a deconfusion method designed for the $[C_{II}]$ survey conducted with CONCERTO.

The treatment of continuum foregrounds has mainly been studied in the context of 21 cm LIM. Most of the techniques rely on both its spectral smoothness and its strong fluctuations to reconstruct and subtract it (McQuinn et al. 2006; Furlanetto et al. 2006; Morales et al. 2006; Chang et al. 2010; Parsons et al. 2012; Liu & Tegmark 2012; Switzer et al. 2015;

Yue et al. 2015) or used to avoid the contaminated angular scales in the (3D) Fourier space (Chapman et al. 2016). Regarding the interlopers, the cross correlation between the data and an alternative tracer probing the same cosmic volume can be used to disentangle the signal of interest from the interlopers. This option has been used for popular lines in intensity mapping (Lidz et al. 2009; Visbal & Loeb 2010; Masui et al. 2013; Croft et al. 2016; Chung et al. 2019) as well as for the $[C_{II}]$ line (Gong et al. 2012; Switzer et al. 2019; Padmanabhan et al. 2022; Chung 2023; Pullen et al. 2023) below $z < 3$. However, there are not enough ancillary probes tracing the EoR at $z > 6$ to cross correlate with the $[C_{II}]$ LIM data (but see Chang et al. 2015; Comaschi et al. 2016). A variety of deconfusion methods are thus under study. By exploiting the anisotropies in the APS of interlopers projected on the $[C_{II}]$ frame (Cheng et al. 2016; Lidz & Taylor 2016; Gong et al. 2020), it is possible to separate out the interloper contribution in the power spectrum, but a high sensitivity and a large-area survey are required to detect these anisotropy features. Cheng et al. (2020) propose a deconfusion method in the phase space for reconstructing the one-point statistics of CO-contaminated data, which can only be used if the CO spectral line energy distribution (SLED) in galaxies varies by less than 20% and with a reasonable noise level. A deep learning approach for de-blending $H_\alpha$ at $z = 1.3$ from $[O_{III}]$ at $z = 2$ was also investigated in Breysse et al. (2015); it achieved a 91% precision but was not applied to the $[C_{II}]$ case.

Among the interloper separation methods, masking is conceptually the simplest one. The $[C_{II}]$ signal is not correlated with the signal from interlopers as they originate from lower redshifts. Hence, removing interloper-contaminated pixels should give an unbiased measurement of the $[C_{II}]$ power spectrum by lowering the interlopers' amplitude in the total APS. Masking efficiency has been studied through its effect on the clustering and shot noise power spectrum levels (Gong et al. 2014; Sun et al. 2018). However, due to the lack of realistic far-infrared to millimetre sky simulations to apply this method to, many side effects are poorly studied, such as the correlation between faint and high-mass galaxies, the correlation between the signal and the mask, and the mixing of $k$ modes induced by the mask. Yue et al. (2015) applied the masking to mock maps, but they did not investigate the systematics as cosmic variance or residual contamination.

In this paper we study all the steps for an end-to-end component separation method applied to the $[C_{II}]$ LIM survey conducted with CONCERTO in order to extract the $[C_{II}]$ APS. We are interested in the detection significance and accuracy of the $[C_{II}]$ APS reconstruction that can be achieved in a typical CONCERTO field after foreground removal. We also study the effect of the field-to-field (FtF) variance on the masking results for CONCERTO-like field sizes. This variance affects both the interloper and the $[C_{II}]$ signals.

For this, we used Simulated Infrared Dusty Extragalactic Sky (SIDES), a newly developed realistic empirical simulation of the dusty infrared-millimetre sky up to $z = 7$, which is fully described in Béthermin et al. (2017, 2022) and Gkogkou et al. (2023), hereafter B17, B22, and G22, respectively. The component separation methods are applied directly to these mock observations, as they would be applied to the real CONCERTO observations. For the CIB contamination (continuum emission), our component separation is based on two different approaches. The first is principal component analysis (PCA) adapted to LIM data (e.g., Bigot-Sazy et al. 2015; Yohana et al. 2021), and the second is a baseline estimator called asymmetric reweighed penalized least squares (arPLS; Baek et al. 2015).
For the interloper emission lines, as CONCERTO carries out the [CII] survey on the widely studied Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) field, we rely on a custom-made masking strategy. In this observational context, our main focus is on the 2D APS, which is the initial quantity derived from the data. The APS of each species can be directly compared, unlike the spherically averaged power spectra, which require projection onto the [CII] frame. Additionally, as detailed in B22, the lowest radial modes will not be explored, due to CONCERTO's low spectral resolution, and the majority of the captured signal will be predominantly present in the well-sampled transverse modes.

The paper is organized as follows. We first present the CONCERTO mock observation cubes from SIDES (Sect. 2). In Sect. 3 we address the removal of CIB anisotropy contamination. In Sect. 4 we present the masking strategy based on the use of the stellar mass, M*, as a CO proxy and a beam width criterion. In Sect. 5 we present the masking results on the interloper maps alone and the interloper-contaminated [CII] maps. Finally, we summarize our results and conclude in Sect. 6. The validation of our APS estimator is detailed in Appendix A. Throughout the paper we assume a Planck Collaboration XIII (2016) cosmology and a Chabrier (2003) initial mass function, as in SIDES (B17).

## 2. SIDES simulation

The SIDES simulation, which mimics the CONCERTO extra-galactic sky, is a central element of our study. SIDES accurately replicates a series of statistical properties, spanning from mid-infrared to millimetre wavelengths, including the redshift distributions and number counts of galaxies, line luminosity functions, and CIB fluctuations (see B17, B22, and G22 for a complete description of the simulation and for comparisons with observational constraints). It produces realistic LIM mock data constituting our input sky used to evaluate the component separation method efficiency.

In Sect. 2.1, we briefly summarize the important SIDES characteristics for our use. We further review the map and cube-making process in Sect. 2.2.

### 2.1. The SIDES catalogue

SIDES produces a catalogue of clustered galaxies with their sky coordinates and redshift (up to z = 7), together with a large set of physical parameters, as for example their star formation rate (SFR), stellar mass (M*), or flux in several instrument bandpasses. All relations used to get these physical parameters are derived empirically from observational trends and include an appropriate scatter. First, from a dark matter (DM) simulation, SIDES produces a galaxy light cone by linking the stellar mass M* to the DM halos' mass, using abundance matching. Two versions exist: one using the 1.4 × 1.4 deg² DM light cone from the Bolshoi-Planck simulation (B17, B22) and one using the 117 deg² DM light cone from the Uchuu simulation (Ishiyama et al. 2021, G22). Next, galaxies are randomly assigned to star-forming or passive galaxies, according to the observed evolution of the star forming galaxy fraction by Davidzon et al. (2017). Passive-galaxy contributions are then neglected because they are reported to be very faint in the far-infrared/sub-millimetre (e.g., Whitaker et al. 2021). Star-forming galaxies are classified as main sequence (MS) or star-bursting (SB) galaxies according to a redshift-evolving fraction described in Béthermin et al. (2012). The SFR of each galaxy is derived using its redshift, M*, and type (MS or SB).

Then, a SED is attributed to each galaxy based on its type, total infrared luminosity $L_{IR}$ and mean radiation field intensity $\langle U \rangle$. $L_{IR}$ is derived from the direct use of the Kennicutt (1998) conversion factor and SFR. As for the mean radiation field $\langle U \rangle$, it follows a redshift evolution described in Béthermin et al. (2017) and reflect the evolution of the dust temperature with redshift.

For CO, the luminosity of the J = 1 transition is first computed using the $L_{IR} - L_{CO(1-0)}$ correlation from Sargent et al. (2014). For SB galaxies, an offset of -0.46 dex is added to $L_{IR}$. Then, for the rest of the transitions, SLED templates from Bournaud et al. (2015) for MS and from Birkin et al. (2021) for SB are attributed. For MS, contributions to the final SLED of clumpy and diffuse SLEDs are defined by the relation between $\langle U \rangle$ and CO(5-4)/CO(2-1) flux ratio from Daddi et al. (2015).

The two [CII] fine-structure lines are derived from empirical relations between line luminosity ratios. First the relation between [CII](1-0)/$L_{IR}$ and CO(4-3)/$L_{IR}$ is calibrated to derive $L_{[CII](1-0)}$. Then, the relation between [CII](2-1)/[CII](1-0) and CO(7-6)/CO(3-4) is calibrated to derive $L_{[CII](2-1)}$. The calibration of these relations is discussed in B22.

Finally, two models are implemented to compute the [CII] luminosity, following two $L_{[CII]}$-SFR relations. The first one is the empirical relation calibrated in the local universe from De Looze et al. (2014), hereafter DL14. The second one is from Lagache et al. (2018), hereafter L18, obtained via a semi-analytical model, and containing a weak redshift dependence. Although the exact $L_{[CII]}$-SFR relation remains debated at high redshifts, the recent observational results of the ALPINE survey (Schaerer et al. 2020) shows that high-redshift galaxies tend to follow the local $L_{[CII]}$-SFR relation. For this reason, we take [CII] following the DL14 relation as our fiducial model.

### 2.2. Cube making

From the SIDES catalogue, we create spatial-spectral data cubes for each component: continuum (CIB), CO, [CII] and [CII]. The adopted pixel size is 5 arcseconds and the frequency axis ranges from 125 to 305 GHz with an absolute spectral resolution of 1 GHz. Voxels (the spatial-spectral elements of a cube) are populated with the corresponding intensities. Lastly, the beam smearing is applied to each frequency channel by convolution in Fourier space with a Gaussian beam kernel and periodic boundary conditions. The Gaussian beam profile in one frequency channel is

$$G_{\text{beam}} = \frac{1}{2\pi\sigma^2} e^{-\frac{r^2}{2\sigma^2}}, \qquad (1)$$

where $\sigma$ is the standard deviation of the Gaussian beam profile linked to the frequency $v$ of the channel by

$$\sigma = \frac{FWHM}{2\sqrt{2\ln(2)}} = \frac{1}{2\sqrt{2\ln(2)}} \frac{1.22c}{vD}, \qquad (2)$$

with FWHM is the full width at half maximum of the instrument beam, c is the celerity of light and D = 12 m is the diameter of the APEX telescope where CONCERTO is installed. All the components, CIB, CO, [CII], and [CII], are summed together to create a foreground-contaminated cube. This cube represents the mock observation of the sky as observed by CONCERTO. The aim of the deconfusion is to retrieve the [CII] signal embedded in the foreground-contaminated data. We deal with the two types of foreground one by one, starting with the continuum foreground.

## 3. Continuum removal

The CIB consists of a smooth, frequency-coherent baseline in the mock electromagnetic spectra, resulting from the sum of dust
continuum emission from galaxies along the line of sight (los) plus a part of continuum emission from the neighbouring los due to the beam smearing. The CIB has strong angular fluctuations that are the dominant contribution to the measured power spectrum (in particular, see Fig. 15 of B22). In this section we compare the efficiency of two different methods for continuum removal: the PCA and the asymmetric re-weighted penalized least squares (arPLS). The analysis was conducted on frequencies ranging from 125 to 305 GHz. For the [CII] signal, as the SIDES simulation only covers the redshift range $z \le 7$, the analysis was done only for $\nu > 237$ GHz. We first assessed the continuum removal on SIDES. SIDES contains both the CIB (mean continuum level) and its spectral fluctuations, which are about 10–20% of the CIB. But for CONCERTO, even if in theory it would be sensitive to the CIB, the data processing removes the mean levels. We thus also performed the continuum removal after subtracting the CIB at each frequency of the SIDES simulation (i.e. removing the mean level).

### 3.1. Principal component analysis

Principal component analysis is a non-parametric method of dimensionality reduction, based on constructing linear combinations of variables in the dataset. Geometrically, these new variables, the so-called principal components, represent the directions of maximum variance in the dataset. The principal components do not have a direct physical interpretation. They are only linked to the variance, not to an underlying physical origin. Most of the variance in the dataset is described by the first principal components, ranked by the amount of variance they carry. Dimensionality reduction corresponds to discarding principal components of lower significance, when projecting the dataset onto this new basis. As mentioned above, the CIB produces the strongest fluctuations at all frequencies in CONCERTO. It is therefore likely that its contamination is condensed in the first few principal components and can be accurately removed. We use for this work the PCA implementation from the package scikit-learn (Pedregosa et al. 2011). We create a matrix, $X$, that contains the collection of spectra from the foreground-contaminated mock cube of Sect. 2.2. The dimensions of $X$ are $N_v \times N_p$, where $N_v$ is the number of frequency channels and $N_p$ is the number of pixels. For clarity, we use Greek symbols for the frequency indices and a Latin symbol for the spatial indices ($p$ for pixel). The covariance matrix, $C$, of the spectral data reads

$$C_{vv'} = \tilde{X}_{vp} \tilde{X}_{pv'}, \qquad (3)$$

where $\tilde{X}$ is defined as

$$\tilde{X}_{vp} = X_{vp} - \mu_p, \qquad (4)$$

and $\mu_p$ is the average in the pixel $p$ over all frequencies. The covariance matrix, $C$, is square, symmetric, and has dimension $N_v$. Next, we computed the eigenvector matrix, $\tilde{E}$, and the block matrix, $\lambda$, that contains the eigenvalues from the diagonalization of $C$, such that

$$C_{vv'} = E_{\nu\alpha}^{-1} \lambda_{\alpha\beta} E_{\beta v'}. \qquad (5)$$

Eigenvectors and eigenvalues are sorted according to the power enclosed in each principal component. We select the $k$ ($k < N_v$) first principal components and build the block matrix $\lambda'$ with their corresponding eigenvalues. This allows the derivation of what is usually called the projection templates:

$$T_{pv} = \tilde{X}_{vp}^T E_{\nu\alpha}^{-1} \frac{1}{\sqrt{\lambda'_{\alpha\nu}}}, \qquad (6)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Spectrum of CIB (continuum) emission in one pixel as a function of frequency (black line). The reconstructed spectra obtained with PCA using the first (red), the first two (orange), and the first three (blue) principal components are over-plotted.</figcaption></figure>

and the continuum is therefore estimated as

$$\tilde{X}_{vp}^{\text{continuum}} = E_{\alpha\nu}^{T} \sqrt{\lambda_{\alpha\beta}} T_{p\beta}^{T} + \mu_{p}. \qquad (7)$$

Eventually, the line spectrum is given by

$$X_{vp}^{\text{lines}} = X_{vp} - \tilde{X}_{vp}^{\text{continuum}}. \qquad (8)$$

Figure 1 shows a representative example of reconstructed continuum spectra with PCA. When using the first two components, the continuum is not well reconstructed, meaning that the continuum features are not totally contained in the first two principal components. When adding the third principal component, the reconstructed spectrum does get closer to the intrinsic continuum spectrum but line emission starts to show up. Hence, it is not possible to reconstruct more accurately the continuum by adding more principal components.

Figure 2 shows the effect of the reconstructed continuum subtraction on the continuum power spectrum amplitude, using the first two principal components. Around $z = 5.2$, the residual continuum is 29 times lower than our fiducial [CII] model and 7 times lower than the L18 [CII] model. But the higher the redshift, the lower the contrast. At $z = 7$, the residual continuum is at a similar level than the fiducial [CII] model and 10 times higher than the L18 model.

When removing the CIB (i.e. the mean level of each frequency channel of the cube), the lines contaminate the continuum estimate starting from the third principal component, as previously. However, removing the CIB decreases the continuum variance, leading to a larger residual continuum power, which is a factor of 11 above the lines signal for $\nu > 237$ GHz. Hence, PCA does not have the cleaning efficiency required up to redshift 7, and we turned to another method often quoted in the literature that fits the continuum with slowly varying functions.

### 3.2. Advanced baseline estimator: arPLS

We test in this section the advanced baseline estimator named arPLS, introduced by Baek et al. (2015) and based a Whittaker-smoothing-based algorithm. Originally developed to serve in a
<figure><img src="image_3.png" /><figcaption>Fig. 2. APS amplitude as a function of frequency for the continuum (solid black line), the [CII] generated with the DL14 or L18 SFR-L$_{[CII]}$ relations (green and light green lines, respectively), the residual continuum obtained with PCA (Sect. 3.1) using the first two principal components (orange line), and the residual continuum obtained with arPLS (Sect. 3.2, blue line). The power spectrum amplitude is averaged at $k = 0.26 \text{ arcmin}^{-1}$ in frequency channels of 1 GHz width.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Relative error of arPLS to the ‘true’ continuum as a function of frequency. The mean is taken over the spectra of all the los of the mock cube.</figcaption></figure>

context of Raman spectroscopy, this method was applied to time series in an astronomical context by Zeng et al. (2021). Here, we adapted the use of arPLS to spectral vectors. This estimator makes use of the frequency-coherent feature of CIB. It adjusts a slowly varying function to the smooth component of a (electromagnetic) spectrum under the assumption that it varies slowly

with frequency. This smooth component (e.g., the black smooth baseline in the spectrum shown in Fig. 1) is the spectrum of CIB, the sum of all dust continuum emission. Moreover, as the frequency-coherent part of the spectrum is fitted, all other frequency-coherent contamination present as Galactic dust and the CMB are fitted simultaneously. Galaxies are different for each los of the mock cube and can drastically change the slope of the baseline so we need to run the baseline estimator on a pixel-by-pixel basis. The smoothed function is fitted to the data by iteratively minimizing a penalized least squares function:

$$S(b) = (y - b)^T W(y - b) + Rb^T D^T Db, \quad (9)$$

where $y$ is the spectrum to be fitted, $b$ the smooth baseline to be found, $W$ the block weight matrix with the weights $w$ for each data point, and $D$ the difference matrix. It is a quasi blind method as the only parameter to be set is the regularization parameter $R$ controlling the trade-off between the agreement of $b$ to the data and its smoothness. In this work, we set it to 1. To find the line peak regions, weights are asymmetrically given to data points:

$$w_i = \begin{cases} f(y_i - b_i, m, \sigma), & y_i \ge b_i \\ 1, & y_i < b_i, \end{cases} \quad (10)$$

with

$$f(d^-, m, \sigma) = \frac{1}{1 + \exp(2(d - (-m + 2\sigma))/\sigma)}, \quad (11)$$

where $d^- = y^- - b$ represents the distance between the current estimate of $b$ and the data points that are positioned exclusively below it. The $m$ and $\sigma$ are the mean and spread $y^-$ of $d$.

This weighting scheme is designed to be robust to reasonable noise level blurring the baseline. Here, we first want to test the ability of this estimator to retrieve the baseline of the simulated LIM electromagnetic spectra along the CONCERTO frequency band, without any instrument noise.

Figure 3 shows the relative error between the continuum fitted with arPLS and the intrinsic continuum (from the CIB cube). While the baseline is slightly overestimated at all frequencies, the mean error stays below 0.3% for $v \ge 237$ GHz ($z_{[CII]} \le 7$).

Figure 2 shows that the power amplitude of continuum is brought below the fiducial [CII] by a factor of >72 and below the L18 [CII] by a factor >4, for $z = 5. - 72$. This allows a clear detection of the [CII] power even in the fainter [CII] case of L18.

For $v \le 237$ GHz, the mean relative error and its dispersion become larger as there are more and more CO lines. Around 130 GHz, at the end of the frequency band, the accuracy decreases even more since there are fewer neighbouring frequency channels to correlate with. Despite this, the fluctuations of the residual continuum are about 27 times lower than the CO power amplitude and permit a further investigation of the CO power spectrum at those frequencies (below 237 GHz).

Finally, we tested arPLS after CIB removal. In that case, the residual power continuum has the same amplitude as before (blue line in Fig. 2) but it is noisier, with a standard deviation of 10 times larger. In any case, the residual contamination is still lower than the L18 [CII] power spectrum above 237 GHz.

We conclude that this method is appropriate for removing the continuum emission from galaxies for CONCERTO, incidentally better than a standard PCA, and is worth being investigated with noise realizations. In what follows, we detail the interloper component separation method.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Stellar mass distribution of sources as a function of different CO proxies available in SIDES. The stellar mass distribution of SIDES sources is shown in grey scale as a function of different proxies up to $z < 3.5$ (from left to right): $M_* - SFR$, $M_* - L_{IR}$, and $M_* - S_{850}$. The best fits for each relation in each redshift bin (given in Table 1) are over plotted, from pink for the first bin ($0 < z < 0.35$) to cyan for the last bin ($2.75 < z < 3.5$).</figcaption></figure>

# 4. Removing the contamination of line interlopers

Once the continuum power amplitude is successfully lowered, the main source of contamination becomes the line interlopers, that is to say, the emission lines redshifted in the same frequency band as [CII]. Dominant interlopers are the CO rotational ladder and two [C1] fine structure lines from foreground galaxies at $z < 3$. The [O1] 145 µm, [NII] 121.9 µm and NII 205.2 µm lines are also present in the CONCERTO frequency range but their intensities are negligible with respect to the [CII] intensity (Silva et al. 2015). Consequently, only CO and [C1] interlopers are considered.

The LIM survey exploits the confusion regime from the far-IR to (sub-)millimetre, where the beam smearing is significant. Because of the confusion regime and together with the expected instrument noise level for CONCERTO (CONCERTO Collaboration 2020) we do not expect to detect lines from individual foreground galaxies. It becomes thus impossible to blindly mask the brightest pixels using a fix intensity threshold. Therefore, to select interloper-dominated voxels, we need to rely on external catalogues of the COSMOS field (Ilbert et al. 2013; Laigle et al. 2016; Nayyeri et al. 2017; Weaver et al. 2022) providing accurate positions and redshifts of the foreground sources. In this section we present the masking strategy based on the ancillary catalogue of detected sources in the COSMOS field. We first justify our choice of using an external catalogue in stellar mass $M_*$ and its use as a CO proxy to select foreground sources. We then detail the masks we build from this ancillary catalogue.

## 4.1. Stellar mass as a CO power proxy

For typical redshifts of foreground sources in COSMOS ($z \sim 1-2$), the selection wavelengths of the COSMOS catalogues trace the stellar emission, that traces the stellar mass, rather than the ultraviolet, that traces the unobscured SFR. Therefore, the stellar mass is usually more complete in those catalogues than the SFR, from which moreover the obscured contribution is often missing.

Figure 4 shows the distribution in $M_*$ for galaxies $z \le 3.5$ as a function of other quantities that can also serve as CO proxy available in SIDES: SFR, $L_{IR}$ and flux at 850 µm $S_{850}$. From these distributions, we obtain the best fit relation for objects in each redshift bin of the completeness relation, taken from Laigle et al. (2016, hereafter L16) and reported in Table 1. The

completeness relation of the COSMOS catalogues describes, as a function of redshift, the stellar mass down to which 90% of sources are detected. Then from the best fits, we get equivalent depths to the stellar mass depths $M_*^{L16}$. Results are listed in Table 1. One should reach these depths to get a catalogue as complete in these proxies as in stellar mass. In comparison, $L_{IR}$ values inferred from the Multi-Band Imaging Photometer (MIPS)/Spitzer at 24 µm (Le Floc'h et al. 2005; Bavouzet et al. 2008) are $1 \times 10^{10} L_{\odot}$ at $z = 0.3$, $3 \times 10^{10} L_{\odot}$ at $z = 0.6$, and $1 \times 10^{11} L_{\odot}$ at $z = 1$. It is one order of magnitude above the $L_{IR}$ depth needed. About SFR, no survey has such a completeness in both ultraviolet and IR to compute the total SFR. Geach et al. (2016) report the Submillimetre Common-User Bolometer Array 2 (SCUBA-2) 850 µm confusion limit to be 0.8 mJy beam$^{-1}$ with a completeness below 10% at this flux. Completeness is close to 100% above 5 mJy beam$^{-1}$, which is largely above (between 1 and 3 orders of magnitude) the $S_{850}$ requested depth.

To further motivate this masking strategy using SIDES, we look at the contribution of each stellar mass population to the CO power spectrum. Especially, we want to make sure that (1) most of the CO power is produced by sources above the L16 completeness relation, so that they can actually be detected and masked in real data, and (2) that CO power is produced by a reasonable amount of sources, so that masking them would not cost too much in survey area.

Figure 5 shows the cumulative contribution to the CO power amplitude per increasing logarithmic step of stellar mass at $z_{[CII]} = 5.2$. By looking at the cumulative contribution rather than the non-cumulative contribution, we take into account for each stellar mass population (i.e. at each step) both its auto-correlation and its cross-correlation with the less massive populations. The cross-correlation contribution is actually important: the non-cumulative contributions (i.e. the auto-correlations alone) of all stellar mass populations represent 60% of the total CO power. The remaining 40% is due to cross-correlations between the different stellar masses. Figure 5 shows that the most massive galaxies, by their auto and cross correlations, contribute the most to the CO power. This is also the case at higher redshifts: at $z_{[CII]} = 6.5$ and $z_{[CII]} = 7$. The correspondence between [CII] redshifts and interloper redshifts is given in Table 2.

For the three redshifts, galaxies more massive than $\ge 10^{8.1} M_{\odot}$ represent only 22–25% of the total object number but produce most of the CO interloping power (99.7%). On the other hand, 75–77% of objects less massive than $10^{8.1} M_{\odot}$ are responsible for
<table><caption>**Table 1.** Equivalent depth in SFR, $L_{IR}$, and $S_{850\mu m}$ to the $M_*^{L16}$ depth of the COSMOS completeness relation, obtained with the best fits from Fig. 4.</caption><thead><tr><th>Redshift</th><th>log M<sub>*</sub><sup>L16</sup>/M<sub>☉</sub></th><th>SFR [M<sub>☉</sub>.yr<sup>-1</sup>]</th><th>L<sub>IR</sub> [L<sub>☉</sub>]</th><th>S<sub>850μm</sub> [mJy]</th></tr></thead><tbody><tr><td>0.35–0.65</td><td>8.1</td><td>(4.0 ± 1.3) × 10<sup>-2</sup></td><td>(4.0 ± 1.3) × 10<sup>8</sup></td><td>(4.1 ± 1.6) × 10<sup>-3</sup></td></tr><tr><td>0.65–0.95</td><td>9.1</td><td>(8.3 ± 2.7) × 10<sup>-1</sup></td><td>(8.3 ± 2.7) × 10<sup>9</sup></td><td>(1.8 ± 0.6) × 10<sup>-2</sup></td></tr><tr><td>0.95–1.3</td><td>9.3</td><td>1.7 ± 0.6</td><td>(1.7 ± 0.6) × 10<sup>10</sup></td><td>(2.7 ± 0.9) × 10<sup>-2</sup></td></tr><tr><td>1.3–1.75</td><td>9.7</td><td>5.7 ± 1.9</td><td>(5.7 ± 1.9) × 10<sup>10</sup></td><td>(6.9 ± 2.4) × 10<sup>-2</sup></td></tr><tr><td>1.75–2.25</td><td>9.9</td><td>(1.2 ± 0.4) × 10<sup>1</sup></td><td>(1.2 ± 0.4) × 10<sup>11</sup></td><td>(1.2 ± 0.4) × 10<sup>-1</sup></td></tr><tr><td>2.25–2.75</td><td>10.0</td><td>(1.9 ± 0.6) × 10<sup>1</sup></td><td>(1.9 ± 0.6) × 10<sup>11</sup></td><td>(1.7 ± 0.6) × 10<sup>-1</sup></td></tr><tr><td>2.75–3.5</td><td>10.1</td><td>(3.0 ± 1.0) × 10<sup>1</sup></td><td>(3.0 ± 1.0) × 10<sup>11</sup></td><td>(2.3 ± 0.8) × 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Contribution to the total CO APS at $z = 5.2$ of each stellar mass bin and CO transition. The left axis and the dotted black curve show the cumulative number count of CO emitters per bin of stellar mass (in percent) redshifted in the frequency channel 305 ± 0.5 GHz ($z_{[C_{II}]} = 5.2$). The right axis shows the cumulative contribution (in percent) to the total CO power amplitude. Coloured curves correspond to each CO transition, and the red curve corresponds to the sum of all CO transitions. The 90% stellar mass completeness in COSMOS2015 for each CO line is indicated with the same colour code on the x-axis. For comparison, the $[C_{II}]$ level (in percent) is represented with the dashed green line. Power amplitudes are averaged for $0.12 < k < 0.24 \text{ arcmin}^{-1}$.</figcaption></figure>

less than 0.5% of the total CO power. For example, we can highlight that the galaxy population with $10^{10.6} M_{\odot} < M_{*} \le 10^{10.9} M_{\odot}$ is responsible for 24–31% of the CO power while they only represent 0.5–0.8% of all the sources. Hence, the massive sources above $10^{8.1} M_{\odot}$, despite being the less numerous, strongly dominate the contribution to the CO power and belong to the massive class of known objects. This makes masking an efficient method for decreasing the CO contamination by removing a few sources without catastrophically masking the survey area. Hence, the stellar mass $M_{*}$ is a reliable CO proxy for the masking purpose.

To make a realistic masking simulation, we masked only sources in the SIDES catalogue that are above the COSMOS completeness relation. We kept the completeness taken from L16, which is deeper than that from Ilbert et al. (2013) or Nayyeri et al. (2017). Although the completeness reported in

<table><caption>**Table 2.** Emission redshift for the lines at the three observed frequencies.</caption><thead><tr><th>Line</th><th>ν<sub>ref</sub> (GHz)</th><th>z(305 GHz)</th><th>z(253 GHz)</th><th>z(237 GHz)</th></tr></thead><tbody><tr><td>CO(3-2)</td><td>345.81</td><td>0.13</td><td>0.37</td><td>0.46</td></tr><tr><td>CO(4-3)</td><td>461.08</td><td>0.51</td><td>0.82</td><td>0.95</td></tr><tr><td>CO(5-4)</td><td>576.36</td><td>0.89</td><td>1.28</td><td>1.43</td></tr><tr><td>CO(6-5)</td><td>691.63</td><td>1.27</td><td>1.73</td><td>1.92</td></tr><tr><td>CO(7-6)</td><td>806.9</td><td>1.65</td><td>2.19</td><td>2.4</td></tr><tr><td>CO(8-7)</td><td>922.17</td><td>2.02</td><td>2.64</td><td>2.89</td></tr><tr><td>[C<sub>I</sub>](1-0)</td><td>492.16</td><td>0.61</td><td>0.95</td><td>1.08</td></tr><tr><td>[C<sub>I</sub>](2-1)</td><td>809.34</td><td>1.65</td><td>2.20</td><td>2.41</td></tr><tr><td>[C<sub>II</sub>]</td><td>1900.54</td><td>5.23</td><td>6.51</td><td>7.02</td></tr></tbody></table>

<table><caption>**Table 3.** Masked intensity per source and survey area loss at the three observed frequencies, corresponding to $z_{[C_{II}]} = 5.2$, 6.5, and 7 for the masks built with different $r_{\sigma}$.</caption><thead><tr><th>r<sub>σ</sub></th><th>Masked intensity per source</th><th>Masked fraction [%] at 305 GHz</th><th>Masked fraction [%] at 253 GHz</th><th>Masked fraction [%] at 237 GHz</th></tr></thead><tbody><tr><td>1.5</td><td>95%</td><td>7</td><td>11</td><td>14</td></tr><tr><td>2.0</td><td>99%</td><td>11</td><td>15</td><td>22</td></tr><tr><td>2.5</td><td>>99%</td><td>16</td><td>22</td><td>29</td></tr><tr><td>3.0</td><td>>99%</td><td>22</td><td>29</td><td>38</td></tr></tbody></table>

Weaver et al. (2022) is deeper, it only corresponds to a 70% completeness fraction, so the L16 relation represents the best compromise between completeness fraction and depth.

## 4.2. Effective masking

External catalogues provide the position of the sources that we want to mask. We then needed to define the extent of the region to mask around these positions. We based our choice on a beam width criterion. Voxels were masked around the location of the selected sources (and thus around CO or $[C_{I}]$ expected emission) within a radius, $r_{\sigma}$, expressed in units of the Gaussian beam width parameter, $\sigma$ (Eq. (2)). This beam varies with the observed frequency and so does the exact masked area. Table 3 lists the fraction of masked intensity per source and the masked fraction of the maps for masks built with different values of $r_{\sigma}$. Figure 6 displays examples of these masks.

Estimating the APS on masked data is not straightforward because the mask induces aliasing. Several estimators have been developed to overcome this difficulty in the context of the CMB
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Masks obtained following the strategy detailed in Sect. 4 using (from left to right) a masking radius, $r_σ$, of 1.5$σ$, 2.5$σ$, and 3$σ$, respectively. From the top to the bottom, the rows correspond to 305 GHz ($z_{[CII]}$ = 5.2), 253 GHz ($z_{[CII]}$ = 6.5), and 237 GHz ($z_{[CII]}$ = 7), respectively. Masked data appear in black.</figcaption></figure>

measurements. In this work, we used the P of K EstimatorR (POKER; Ponthieu et al. 2011), which was adapted to the case of sky patches of a few square degrees, with arbitrary high angular resolution and complex masks Planck Collaboration XVIII (2011), Planck Collaboration XXX (2014). A detailed discussion on the performance and validation of POKER in the context of this work and the derivation of error bars on APS are presented in Appendix A.

## 5. Interloper separation at $z_{[CII]}$ = 5.2, 6.5, and 7

In this section we assess the efficiency of masking at retrieving the underlying $[C_{II}]_{(DL14)}$ APS in the interloper-contaminated mock LIM data (see Sect. 4.2). As mentioned above, the correction of the mask impact is taken care of by POKER, so we focus here on the residual contribution of CO and $[C_{I}]$ left after masking on the COSMOS field. When using an APS measurement to derive cosmological parameters, one must also account for the FtF variance, namely the overall uncertainty on the relative contribution of CO and $[C_{I}]$ to $[C_{II}]$, in different parts of the sky. This particular point is addressed in Sect. 5.2, and we focus here on the measurement on COSMOS only.

In the following, power amplitudes and power ratios are given for the same $k$ range (as in Table 4), where the clustering

<table><caption>**Table 4.** Power amplitude (in Jy$^2$ sr$^{-1}$) averaged between $k=0.12$ and $0.24$ arcmin$^{-1}$ of the different components, before and after the corresponding component separation method is applied for foreground removal.</caption><thead><tr><th>Frequency</th><th>305 GHz</th><th>253 GHz</th><th>237 GHz</th></tr></thead><tbody><tr><td colspan="4" style="text-align:center; font-weight:bold;">Intrinsic amplitudes</td></tr><tr><td>P(k)<sub>CIB</sub></td><td>1.1 &times; 10<sup>2</sup></td><td>3.6 &times; 10<sup>1</sup></td><td>2.4 &times; 10<sup>1</sup></td></tr><tr><td>P(k)<sub>CO</sub></td><td>6.2 &times; 10<sup>-1</sup></td><td>4.1 &times; 10<sup>-1</sup></td><td>5.1 &times; 10<sup>-1</sup></td></tr><tr><td>P(k)<sub>[CII]</sub></td><td>2.9 &times; 10<sup>-1</sup></td><td>1.1 &times; 10<sup>-1</sup></td><td>1.7 &times; 10<sup>-1</sup></td></tr><tr><td>P(k)<sub>[CIII]</sub></td><td>5.3 &times; 10<sup>-1</sup></td><td>6.0 &times; 10<sup>-2</sup></td><td>3.4 &times; 10<sup>-2</sup></td></tr><tr><td colspan="4" style="text-align:center; font-weight:bold;">CIB residual amplitudes</td></tr><tr><td>PCA</td><td>1.5 &times; 10<sup>-1</sup></td><td>2.8 &times; 10<sup>-2</sup></td><td>2.7 &times; 10<sup>-2</sup></td></tr><tr><td>arPLS</td><td>1.1 &times; 10<sup>-3</sup></td><td>1.5 &times; 10<sup>-4</sup></td><td>1.9 &times; 10<sup>-4</sup></td></tr><tr><td colspan="4" style="text-align:center; font-weight:bold;">r<sub>σ</sub> = 3.0 masked amplitudes</td></tr><tr><td>P(k)<sub>CO+[C<sub>I</sub>]</sub></td><td>(8.7 &plusmn; 2.3) &times; 10<sup>-3</sup></td><td>(1.1 &plusmn; 0.4) &times; 10<sup>-2</sup></td><td>(1.7 &plusmn; 0.6) &times; 10<sup>-2</sup></td></tr><tr><td>P(k)<sub>CO bright</sub></td><td>(9.2 &plusmn; 2.4) &times; 10<sup>-5</sup></td><td>(2.0 &plusmn; 0.7) &times; 10<sup>-5</sup></td><td>(3.3 &plusmn; 1.2) &times; 10<sup>-5</sup></td></tr><tr><td>P(k)<sub>CO faint</sub></td><td>(6.6 &plusmn; 1.7) &times; 10<sup>-3</sup></td><td>(7.7 &plusmn; 2.5) &times; 10<sup>-3</sup></td><td>(1.2 &plusmn; 0.4) &times; 10<sup>-2</sup></td></tr><tr><td>P(k)<sub>[CII] bright</sub></td><td>(5.3 &plusmn; 1.4) &times; 10<sup>-6</sup></td><td>(3.7 &plusmn; 1.2) &times; 10<sup>-6</sup></td><td>(3.1 &plusmn; 1.1) &times; 10<sup>-6</sup></td></tr><tr><td>P(k)<sub>[CII] faint</sub></td><td>(1.9 &plusmn; 0.5) &times; 10<sup>-3</sup></td><td>(2.4 &plusmn; 0.8) &times; 10<sup>-3</sup></td><td>(3.6 &plusmn; 1.3) &times; 10<sup>-3</sup></td></tr><tr><td>P(k)<sub>[CIII]</sub></td><td>(5.3 &plusmn; 1.4) &times; 10<sup>-1</sup></td><td>(7.1 &plusmn; 2.3) &times; 10<sup>-2</sup></td><td>(3.0 &plusmn; 1.0) &times; 10<sup>-2</sup></td></tr></tbody></table>

**Notes.** Error bars on the CO, $[C_{I}]$ and $[C_{II}]$ amplitudes were computed using the dispersion measured on 54 independent fields (with $r_σ = 3.0$). Details are given in Appendix A.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** APS as a function of spatial frequency, *k*. The solid black line shows the APS of the interloper-contaminated map (i.e. CO+[C₁] + [CII]) with no mask. The solid green line shows the intrinsic [CII] APS (from the unmasked contamination-free [CII] map). Coloured crosses show the [CII] APS estimates, obtained from masking the interloper-contaminated map with the different masks from *r*<sub>σ</sub> = 1.5 to 3.0. Coloured points linked by dotted lines show the residual interlopers' APS measured on the masked interloper map (i.e. CO+[C₁] after masking) at (from left to right) *z* = 5.2, 6.5, and 7. The beam APS at the corresponding observed frequency is shown in grey. Vertical lines highlight the range *k* = 0.12–0.24 arcmin<sup>-1</sup>.</figcaption></figure>

dominates. This is because the clustering part of the APS should be easier to interpret than the shot noise, which is dominated by the few strongest sources just below the detection threshold (G22), thus lowering the statistical information.

Results are presented in Fig. 7 at $z_{[C_{II}]}$ = 5.2, 6.5, and 7 for different $r_σ$. With $r_σ$ = 3, the interloper power amplitude is smaller than $2 \times 10^{-2}$ Jy$^2$ sr$^{-1}$ for all three redshifts, that is, it is reduced by a factor of (100 ± 30), (50 ± 16), and (42 ± 15), for each redshift, respectively.

For [CII], our fiducial model that follows the local DL14 L[CII]–SFR relation predicts a drop in power amplitude by a factor of 16 between $z = 5.2$ and 7. In this framework, masking with $r_σ = 2$ at $z_{[CII]} = 5.2$ is enough to obtain a contrast between [CII] and residual interlopers of [CII]/(residual CO+[C₁]) = (30 ± 16) (i.e. 3% of residual contamination). The contrast goes up to (62 ± 32) with $r_σ = 3$. At $z_{[CII]} = 6.5$, the [CII] APS is overestimated. Even masking with $r_σ = 3$ does not achieve a better power ratio than [CII]/(residual CO+[C₁]) = (5.5 ± 3.6) (i.e. 18% of residual contamination after masking). Finally, at $z_{[CII]} = 7$, the contrast between [CII] and residual interlopers is weak due to the low [CII] APS, causing [CII]/(residual CO+[C₁]) = 2.0 ± 1.4 (i.e. 50% of residual contamination after masking) with $r_σ = 3$.

These results raise several questions that we discuss in the next section regarding: how well the [CII] APS is reconstructed, especially at redshift 5.2; the cause of the residual power spectra at $z_{[CII]}$ = 6.5 and $z_{[CII]}$ = 7; and whether increasing $r_σ$ to >3 helps in reducing the residual contamination. To answer these questions, we study separately the masking of [CII], free from contamination, the masking of the bright (i.e. above the mass threshold) sources and the masking of the faint (i.e. below the mass threshold) sources.

## 5.1. Interloper masking

The power spectra from the interloper-contaminated maps measured in the previous section are the sum of the power spectra of the three components ([CII], the bright interlopers and the faint interlopers), all affected by the mask. In Fig. 8, we show what we obtain when using different masks on each component separately, for the three redshifts. Besides, Table 4 compiles the power amplitudes of each component.

We first focus on the APS behaviour when masking the faint and bright sources, as displayed in the second and third rows of

Fig. 8. With $r_σ$ ≤ 1.5, the residual APS shows an overshoot at high *k*. This is due to the residuals on the edges of the bright sources ('wings') that are still dominating the contribution to the APS, when the mask is too narrow compared to the instrument beam. When masking at $r_σ$ ≥ 2, the residual contribution of bright sources drops below $10^{-3}$ Jy$^2$ sr$^{-1}$ and becomes negligible compared to that of the faint sources that stays about $0.8 \times 10^{-2}$ Jy$^2$ sr$^{-1}$. The final power ratio (with $r_σ$ = 3) between faint and bright sources are (faint CO+[C₁])/(bright CO+[C₁]) = (91 ± 47), (470 ± 300), and (480 ± 330) at $z_{[CII]}$ = 5.2, 6.5, and 7, respectively. Interestingly, the third row shows how the APS of faint sources is also affected by the mask, contrary to what was expected in B22. Using SIDES, we indeed see that $r_σ$ = 3 masking lowers the faint power amplitude by a factor of (1.9 ± 0.7)–(2.5 ± 0.6) at $z = 7–5.2$. An explanation for this could be that a fraction of faint sources is clustered around bright sources of the same redshift. Yet, if masking removes bright sources efficiently, it has only a moderate effect on the faint sources. Indeed, increasing $r_σ$ from 1.5 to 3.0 increases the masked area by a factor of 3.1, 2.6, and 2.7 at $z = 5.2$, 6.5, and 7, but reduces the faint signal by only an extra 14%, 10%, and 1%, respectively, while in comparison the bright signal is reduced by over a thousand. As a result, enlarging the masked area around bright sources does not effectively remove more contamination from the faint sources and increases rapidly the loss of survey area.

We also applied the masking procedure to [C₁] and CO maps separately. Sun et al. (2018, hereafter S18) studied the masking of CO emitters having a magnitude $m_K^{AB} \le 22$, which translates to a stellar mass cut evolving with redshift (see Fig. 7 in Sun et al. 2018). This lowers their CO power amplitude by a factor $\ge 100$ at $z = 6.5$ and $k = (3–6) \times 10^{-2}$ h Mpc$^{-1}$, for a 8% loss of survey area. With SIDES, at comparable redshift and spatial frequencies (corresponding to $k = 0.12–0.24$ arcmin$^{-1}$ at $z = 6.5$), masking with this prescription lowers the CO power amplitude by only a factor of 20 for a 22% loss of survey area. This indicates that the difference in the CO power model between S18 and SIDES increases the power faint sources in our mock data. Faint sources left after masking contributes at least 5 times more in SIDES than in S18. Moreover, taking into account the beam smearing when masking increases significantly the masked fraction of the surveyed area. With our
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** APS (in Jy$^2$ sr$^{-1}$) as a function of spatial frequency *k* (in arcmin$^{-1}$) for the three observed frequencies (from left to right): 305, 253, and 237 GHz. This is shown (from top to bottom) for each component separately: all the interlopers (CO+[C>I]), the bright interlopers (i.e. above the mass threshold), the faint interlopers, and [C<sub>II</sub>]. The intrinsic power spectra of each (unmasked) component are shown with solid black lines, and with a solid green line for [C<sub>II</sub>]. The power spectra obtained from the maps that are masked using $r_σ = 1.5–3.0$ are shown with coloured points (from blue to yellow, respectively).</figcaption></figure>

strategy of masking all detected sources that have a magnitude $m_K^{AB} < 24$ in Laigle et al. (2016), the CO power amplitude is lowered by a factor of 50, for a 29% loss of survey area (with $r_σ = 3$). Lastly, CO is still the dominant interloper after masking: the interlopers’ power ratio is (residual CO/residual [C>I]) = (3.2±2.1), on average for the three redshifts studied ($z_{[C_{II}]}$ = 5.2, 6.5 and 7).

With the understanding of how masking affects the interlopers, we show in the last row of Fig. 8 the [C<sub>II</sub>] APS reconstructed from the masked [C<sub>II</sub>] maps of SIDES-Bolshoi. The reconstructed [C<sub>II</sub>] power amplitude is consistent with the intrin-

sic one within $1σ$. However, Appendix A details how the measured [C<sub>II</sub>] APS (using the mask $r_σ = 3$) are biased low, on average by 1%, 7%, and 8% at $z_{[C_{II}]}$ = 5.2, 6.5, and 7, respectively, for 54 Uchuu fields.

In conclusion, the masking strategy with $r_σ = 3$ is able to nicely reconstruct the [C<sub>II</sub>] power spectrum from the SIDES-Bolshoi realization at $z_{[C_{II}]}$ = 5.2, reaching a power ratio [C<sub>II</sub>]/(residual CO) = (80 ± 42) and [C<sub>II</sub>]/(residual [C>I]) = 280 ± 150, with a bias on the [C<sub>II</sub>] APS of 1% due to the mask. At $z_{[C_{II}]}$ = 6.5, the power ratios are [C<sub>II</sub>]/(residual CO) = (7.8 ± 5.1) and [C<sub>II</sub>]/(residual [C>I]) = (25 ± 16), with a 7% bias on
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Mean APS (solid lines) and the $1\sigma$ FtF variance interval (coloured area) for a 2 deg$^2$ field, for [CII], faint CO, faint [C$_I$], and the sum of the last two (green, blue, red, and purple, respectively). This is shown, from left to right, for 300–305 GHz centred on $z_{[CII]} = 5.2$, for 260–265 GHz centred on $z_{[CII]} = 6.5$, and for 235–240 GHz centred on $z_{[CII]} = 7$.</figcaption></figure>

the [CII] APS due to the mask. At $z_{[CII]} = 7$, power ratios are [CII]/(residual CO) = (2.7 ± 1.9) and [CII]/(residual [C$_I$]) = (9.4 ± 6.6), with an 8% bias on the [CII] APS. For all redshifts, residual contamination becomes dominated by the faint sources for $r_{\sigma} \ge 2$ since such sources are moderately affected by the mask.

## 5.2. Field-to-field variance effects on residuals

The detection of the [CII] power spectrum depends on its contrast compared to that of the (masked) faint sources, the latter dominating the residual interloper contamination as shown in the previous section. In addition, these measurements will also be affected by the FtF variance. Indeed, G22 shows that FtF variance cannot be neglected when making a forecast for a LIM experiment, especially for the [CII] LIM survey conducted by CONCERTO on a field the size of COSMOS. We explore this here in more details. We use the Uchuu simulation and the variance model from G22 (their Eq. (8)) to determine the faint CO and faint [C$_I$] APS variance. This equation is used to compute a factor *f* to rescale the variance *σ* for a $\Omega = 1$ deg$^2$ field to the $\Omega = 2$ deg$^2$ field of the SIDES-Bolshoi simulation, assuming that the mean power spectrum $\mu$ does not change between the two field sizes:

$$f = \frac{\sigma(\Omega = 2 \text{ deg}^2)}{\sigma(\Omega = 1 \text{ deg}^2)} \quad (12)$$

The values of *f* for the clustering term and the shot noise term are similar at the 1% level.

Figure 9 shows the mean APS and its associated $1\sigma$ FtF variance interval in a 2 deg$^2$ field, for [CII] and the faint interlopers at $z_{[CII]} = 5.2, 6.5$, and 7. The [CII] power amplitude varies by 40% at $z_{[CII]} = 5.2$ and by 60% at $z_{[CII]} = 6.5$ and 7. The faint (CO + [C$_I$]) APS varies by only 12–15% from $z_{[CII]} = 5.2$ to 7. Thus, the variance of faint sources does not have a strong impact on the contrast between residuals and [CII] in comparison to the [CII] variance. Table 5 lists the level of faint interloper contamination and its ratio to the $1\sigma$ value of the [CII] FtF variance. As a result, residual contamination is smaller than 20% of the [CII] FtF variance up to 287 GHz ($z_{[CII]} = 5.62$) and smaller than 50% of the FtF variance up to 253 GHz ($z_{[CII]} = 6.5$).

<table><caption>**Table 5.** Ratio (in percent) of faint interloper contamination (i.e. 'unmasked faint CO+[C$_I$]') to the [CII] power amplitude when the latter is evaluated at the mean and at the edges of the [CII] $\pm 1\sigma$ FtF variance interval.</caption><thead><tr><th>Frequency band</th><th>At +1σ</th><th>At mean</th><th>At -1σ</th><th>Faint interlopers /1σ</th></tr></thead><tbody><tr><td>305 GHz</td><td>2%</td><td>3%</td><td>4%</td><td>7%</td></tr><tr><td>287 GHz</td><td>6%</td><td>8%</td><td>12%</td><td>21%</td></tr><tr><td>253 GHz</td><td>17%</td><td>27%</td><td>66%</td><td>45%</td></tr><tr><td>237 GHz</td><td>87%</td><td>139%</td><td>348%</td><td>232%</td></tr></tbody></table>

**Notes.** The last column gives (also in percent) the ratio between the faint interloper contamination and the $1\sigma$ value of the [CII] FtF variance.

# 6. Summary and conclusion

We tested a complete foreground deconfusion process for the [CII] LIM survey conducted with CONCERTO using realistic mock data of the sub-millimetre sky provided by SIDES-Bolshoi (B17, B22). To build the SIDES mock data, a galaxy catalogue complete down to $10^7 M_{\odot}$ and $z = 7$ was first obtained using abundance matching between stellar mass and DM halo mass in a simulated DM light cone. Galaxy properties were then derived from their type (MS or SB), redshift, and stellar mass using empirically calibrated relations with appropriate scatters. Finally, mock cubes for CIB, CO, [C$_I$], and [CII] were produced. A foreground-contaminating mock cube was created by adding the separate foreground cubes to the [CII] cube.

We first assessed the contamination by CIB (continuum) fluctuations. We tested the ability of PCA and arPLS to subtract the continuum emission for each *los*, by exploiting its (electromagnetic) frequency-coherent distinctive feature. Principal component analysis is not able to separate the continuum from molecular and atomic lines up to $z_{[CII]} = 7$. The arPLS method fits all the frequency-coherent emissions with a precision of 0.3% in the mock spectra and sufficiently removes the continuum contamination for the two [CII] models provided by SIDES up to $z = 7$. For $5.2 \le z \le 7$, the residual CIB power amplitudes are lower than that of our fiducial DL14 (L18) [CII] model by a factor of 72 (4).

For the deconfusion with the line interlopers, we relied on a masking approach that uses an external COSMOS stellar mass
catalogue, which provides accurate sky positions and redshifts for foreground sources to mask interloper-contaminated voxels. We show that $M_*$ is the deepest CO and $[C_1]$ proxy available in the COSMOS field. We also investigated its use as a CO proxy. The most contaminating sources ($M_* > 10^{8.1} M_\odot$) are the least numerous and are above the L16 completeness relation, allowing us to efficiently remove the interloper contamination with a limited loss of the surveyed area. The spatial extent of the masks is defined with a (frequency-evolving) beam width criterion. To test the masking on our mock data in a realistic way, we masked only sources above the COSMOS (deep) completeness relation described in L16.

For the SIDES-Bolshoi simulation, the power amplitude from bright CO and $[C_1]$ sources was lowered by four orders of magnitudes. Masking is also more efficient than expected in B22 as it also removes a fraction of the faint source signal at the same time, thanks to the clustering. The residual power amplitude left after masking is lower than $2 \times 10^{-2}$ Jy$^2$ sr$^{-1}$, rather independent of frequency, and is dominated by the faint sources that are not masked (especially CO sources). The CO power amplitude is lowered by a factor of 50 at $z = 6.5$ with our masking, while it is lowered by a factor greater than 100 with a smaller masking depth in S18, indicating that the signal of sources below the mass threshold is stronger in SIDES than in S18.

In the framework of our model, which contains a strong decrease in the $[C_{II}]$ APS with redshift, we obtained the following results for the $[C_{II}]$ APS measurement.

At $z_{[C_{II}]} = 5.2$, the power ratio is $[C_{II}]/(\text{residual CO}+[C_1]) = 62 \pm 32$ (2.5% residual contamination) for a 22% surveyed area loss and 1% of bias due to the mask. At this given redshift, the $[C_{II}]$ power spectrum varies by about 40% due to the FtF variance. So, FtF variance is larger than any effect introduced by the masking or by residual contamination (which is $\le 4\%$ of FtF variance).

At $z_{[C_{II}]} = 6.5$, the residual contamination after masking is a factor of $[C_{II}]/(\text{residual CO}+[C_1]) = 5.5 \pm 3.6$ below the $[C_{II}]$ power spectrum amplitude. This residual contamination and the reconstruction by POKER (see Appendix A) led to an overestimation of the $[C_{II}]$ power spectrum by a factor of 1.4 for a 29% loss of survey area. The variance estimate from SIDES-Uchuu gives a residual contamination by the faint sources that should not exceed half of the $[C_{II}]$ FtF variance at this redshift. At $z_{[C_{II}]} = 7$, SIDES predicts a drop in $[C_{II}]$ power spectrum amplitude by a factor of 16 with respect to $z_{[C_{II}]} = 5.2$ and an important FtF variance for a COSMOS-like field. Indeed, the measured residual contamination after masking is only a factor of $2.0 \pm 1.4$ below the $[C_{II}]$ power spectrum amplitude.

Hence, residual contamination is still important above $z_{[C_{II}]} > 6.5$, despite the fact that most (>98%) of the CO power has been masked. Masking even more sources to bring residues to a fainter level is tricky because the $[C_{II}]$ APS is already slightly affected by the mask (i.e. underestimated by 7–8% between $z = 6.5$ and 7; see Appendix A) and the cost per masked source in a surveyed area is lower frequencies (i.e. higher redshifts). To reach the $[C_{II}]$ signal at these redshifts, other or complementary methods for masking will be needed.

Our analysis makes use of the SIDES simulations (Béthermin et al. 2022), in which the $[C_{II}]$ power spectrum is predicted with a much lower amplitude than in earlier models (e.g., Serra et al. 2016). For the first generation of LIM experiments, such as CONCERTO, and given such a low-amplitude $[C_{II}]$ power spectrum, we expect the instrument noise to dominate. Van Cuyck et al. (in prep.) will introduce the noise generator of SIDES and investigate the effects of instrument white

noise on different LIM observables, such as auto-power spectra, cross-power spectra between lines and galaxies, and cross-power spectra between different lines within the LIM dataset. We have shown that our method for removing the continuum emission does not create any systematic residuals, but this will have to be tested with instrument noise. On the other hand, we expect masking to be robust to noise contamination because the interlopers' signal is removed regardless of the noise level. The robustness of cross-correlation measurements between different lines within the LIM dataset in the presence of instrument noise is also worth investigating in the context of component separation. Such a cross correlation could allow for a direct estimate of interlopers' residual power amplitude after masking. The last effect that would need further checks for the masking technique is the inclusion of the catalogues' redshift uncertainties, reported to be $\delta z/(1+z) = 0.025$ below $z \le 3$ for the COSMOS field (Weaver et al. 2022).

Finally, the CO-line emission, while being a strong foreground for the $[C_{II}]$ survey, is of great importance. The CO part of the CONCERTO survey can be used to probe the CO SLED and the gas content of galaxies at $z \le 3$, including at cosmic noon ($2 \le z \le 3$), when star formation in galaxies statistically reaches its peak. The component separation methods and tools developed in this work for $[C_{II}]$ can be adapted and used to exploit CO surveys.

**Acknowledgements.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 788212) and from the Excellence Initiative of Aix-Marseille University-A\*Midex, a French "Investissements d'Avenir" programme. M.A. acknowledges support from FONDECYT grant 1211951, CONICYT + PCI + Max Planck Institute for Astronomy MPG 190030, CONICYT + PCI + REDES 190194 and ANID BASAL project FB210003.

## References

Atek, H., Richard, J., Jauzac, M., et al. 2015, ApJ, 814, 69

Baek, S.-J., Park, A., Ahn, Y.-J., & Choo, J. 2015, Analyst, 140, 250

Bavouzet, N., Dole, H., Le Floc'h, E., et al. 2008, A&A, 479, 83

Bernal, J. L., & Kovetz, E. D. 2022, A&ARv, 30, 5

Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Béthermin, M., Gkogkou, A., Van Cuyck, M., et al. 2022, A&A, 667, A156

Bigot-Sazy, M. A., Dickinson, C., Battye, R. A., et al. 2015, MNRAS, 454, 3240

Birkin, J. E., Weiss, A., Wardlow, J. L., et al. 2021, MNRAS, 501, 3926

Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55

Bournaud, F., Daddi, E., Weiß, A., et al. 2015, A&A, 575, A56

Breysse, P. C., Kovetz, E. D., & Kamionkowski, M. 2015, MNRAS, 452, 3408

Carilli, C. L., Chluba, J., Decarli, R., et al. 2016, ApJ, 833, 73

Chabrier, G. 2003, PASP, 115, 763

Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463

Chang, T. C., Gong, Y., Santos, M., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 4

Chapman, E., Zaroubi, S., Abdalla, F. B., et al. 2016, MNRAS, 458, 2928

Chang, Y.-T., Chang, T.-C., Bock, J., Bradford, C. M., & Cooray, A. 2016, ApJ, 832, 165

Cheng, Y.-T., Chang, T.-C., & Bock, J. J. 2020, ApJ, 901, 142

Chornock, R., Berger, E., Fox, D. B., et al. 2013, ApJ, 774, 26

Chung, D. T. 2023, Phys. Rev. D, 107, 023509

Chung, D. T., Viero, M. P., Church, S. E., et al. 2019, ApJ, 872, 186

Comaschi, P., Yue, B., & Ferrara, A. 2016, MNRAS, 463, 3193

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

Croft, R. A. C., Miralda-Escudé, J., Zheng, Z., et al. 2016, MNRAS, 457, 3541

Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Dayal, P., Dunlop, J. S., Maio, U., & Ciardi, B. 2013, MNRAS, 434, 1486

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

Drake, A. B., Garel, T., Wisotzki, L., et al. 2017, A&A, 608, A6

Eilers, A.-C., Hennawi, J. F., Davies, F. B., & Oñorbe, J. 2019, ApJ, 881, 23

Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7

Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012, ApJ, 758, 93
Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2016, MNRAS, 465, 1789

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

Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49

Gong, Y., Silva, M., Cooray, A., & Santos, M. G. 2014, ApJ, 785, 72

Gong, Y., Chen, X., & Cooray, A. 2020, ApJ, 894, 152

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

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210

Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, ArXiv e-prints
[arXiv:1709.09066]

Lagache, G., Puget, J. L., Abergel, A., et al. 2000, in ISO Survey of a Dusty
Universe, eds. D. Lemke, M. Stickel, & K. Wilke (Springer), 548, 81

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

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Lidz, A., & Taylor, J. 2016, ApJ, 825, 143

Lidz, A., Zahn, O., Furlanetto, S. R., et al. 2009, ApJ, 690, 252

Liu, A., & Tegmark, M. 2012, MNRAS, 419, 3491

Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, Am. Astron. Soc.
Meeting Abst., 229, 141.09

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJ, 763, L20

McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815

Morales, M. F., Bowman, J. D., & Hewitt, J. N. 2006, ApJ, 648, 767

Nayyeri, H., Hemmati, S., Mobasher, B., et al. 2017, ApJS, 228, 7

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104

Padmanabhan, H., Breysse, P., Lidz, A., & Switzer, E. R. 2022, MNRAS, 515, 5813

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

Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621

Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, ApJ, 756, 165

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825

Pineda, J. L., Langer, W. D., & Goldsmith, P. F. 2014, A&A, 570, A121

Planck Collaboration XVIII. 2011, A&A, 536, A18

Planck Collaboration XXX. 2014, A&A, 571, A30

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration Int. XLVII. 2016, A&A, 596, A108

Ponthieu, N., Grain, J., & Lagache, G. 2011, A&A, 535, A90

Pullen, A. R., Breysse, P. C., Oxholm, T., et al. 2023, MNRAS, 521, 6124

Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19

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

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Serra, P., Doré, O., & Lagache, G. 2016, ApJ, 833, 153

Silva, M., Santos, M. G., Cooray, A., & Gong, Y. 2015, ApJ, 806, 209

Sun, G., Moncelsi, L., Viero, M. P., et al. 2018, ApJ, 856, 107

Switzer, E. R., Chang, T. C., Masui, K. W., Pen, U. L., & Voytek, T. C. 2015, ApJ, 815, 51

Switzer, E. R., Anderson, C. J., Pullen, A. R., & Yang, S. 2019, ApJ, 872, 82

Vallini, L., Ferrara, A., Pallottini, A., & Gallerani, S. 2017, MNRAS, 467, 1300

Visbal, E., & Loeb, A. 2010, JCAP, 2010, 016

Wang, F., Davies, F. B., Yang, J., et al. 2020, ApJ, 896, 23

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021, Nature, 597, 485

Yohana, E., Ma, Y.-Z., Li, D., Chen, X., & Dai, W.-M. 2021, MNRAS, 504, 5231

Yue, B., Ferrara, A., Pallottini, A., Gallerani, S., & Vallini, L. 2015, MNRAS, 450, 3829

Zeng, Q., Chen, X., Li, X., et al. 2021, MNRAS, 500, 2969
# Appendix A: Measuring the APS of a masked distribution of point sources with POKER

In this paper we make extensive use of POKER, an algorithm developed to compute the APS of a masked 2D signal on a few square degrees at arbitrary high angular resolution (Ponthieu et al. 2011). POKER was developed to measure the power spectrum of the CIB in the context of experiments such as the infrared astronomical satellite (IRAS) and Planck (Planck Collaboration XVIII 2011) for which the CIB could be seen as an ubiquitous diffuse signal and treated as a random Gaussian field. From the direct Fourier transform of the masked data (pseudo-power spectrum), POKER inverts the mask induced aliasing and provides an estimate of the underlying APS that is unbiased on average when the signal is independent of the mask. When the signal of interest is highly non Gaussian such as in this work, with point sources and voids in between, the outcome of the power spectrum derived on masked area must however be reconsidered, and this, actually regardless of the APS estimator.

For pedagogical purpose, we can start with the simplest example: a signal whose APS is constant for all angular modes. This APS can be that of a single point source or that of uniform white noise. In the case of uniform white noise, even a mask as complex as those shown in this work will leave some signal on the map from which the pseudo-power spectrum can be derived. In the case of a point source, if the source is masked, the pseudo-power spectrum is then zero and no measurement is possible. In this work, the signal can be approximated as resulting from a collection of $N_s$ point sources. If the sources had all the same flux, the intrinsic power spectrum would simply be $N_s$ times the uniform APS of single source of flux $\phi$. The mask covers a fraction $f_{mask}$ of the field but suppresses a fraction of sources $f_{sources}$ that is different because sources are not uniformly distributed, and loosely speaking, several sources can fall in the same hole. In this case, the pseudo-power spectrum is underestimated by $f_{sources}$ and corrected up by $f_{mask}$ by POKER. Figure A.1 illustrates two different examples where the mask and the source locations are taken directly from a 114 arcmin² section of the SIDES-Bolschoi simulation and their flux is fixed to 1: depending on the effective (a priori unknown) clustering of the sources compared to the mask, the reconstructed APS is larger or smaller than the intrinsic $[C_{II}]$ APS. In this case, it is therefore impossible by construction to measure an unbiased APS in the presence of a mask, and the bias will depend on the clustering properties of the signal compared to the mask. We can extend this reasoning to the sum of several distributions of point sources, each distribution with a different flux. Each one of them will lead to an APS estimate biased by a different $f_{sources}/f_{mask}$ so that the final APS is some average of these ratios, weighted by the source flux. Last, the convolution by the point spread function adds another feature. Indeed, it gives a spatial extension to the initial point sources, so that their signal can be only partially masked and more degeneracies show up when computing the absolute level of the pseudo-power spectrum. It is therefore essential to test the APS derivations on realistic source and mask simulations to assess the magnitude of these potential biases.

Figure A.2 presents the results of this test, from which we built the error bars given in Sect. 5.1. For each of the 54 fields of $2 \text{ deg}^2$ cut in the Uchuu simulation, we built a mask with $r_\sigma = 3$ following Sect. 4 and computed the ratio between the out-

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Illustration of the POKER correction for two different 114 arcmin² fields. Top row: CO+[CⅠ] masks corresponding to the two fields, as well as the [CⅡ] sources (red points). Bottom row: APS of the masked [CⅡ] source maps, with or without correction of mask aliasing (labelled ‘Poker’ and ‘Masked’, respectively). On the left-hand side, the mask covers 44% of the sky patch (dashed blue line) but masks only 29% of the [CⅡ] sources (dashed green line). POKER returns an output power spectrum that has the correct shape (flat, up to the sample variance) but that is overestimated by a factor of $(1-0.29)/(1-0.44)=1.27$. This is the dashed red line, and it matches the average level of the output APS. On the right-hand side, the mask covers only 36% of the sky patch but 58% of the sources. Hence, the output APS returned by POKER is underestimated by $(1-0.58)/(1-0.36)=0.65$.</figcaption></figure>

put (masked) and input $[C_{II}]$ APS. We then computed the average of this ratio over all the 54 fields and show that it is not biased low by more than 1%, 7%, and 8% on average for the three redshifts ($z_{[C_{II}]} = 5.2, 6.5$, and 7), up to an angular scale of $k_{lim} = 2.14, 1.94, 1.33 \text{ arcmin}^{-1}$, respectively. The bias increases with redshift because the masked area is larger, as it depends on the beam size. Error bars for an APS obtained using a mask at a given frequency are the dispersion of this ratio at this frequency. The obtained errors bars represent 25%, 33%, and 34% of the masked APS at $z_{[C_{II}]} = 5.2, 6.5$, and 7, respectively, and are also used for smaller $r_\sigma$ masks. Error bars displayed in Fig. A.2 are normalized by the number of realizations. They are constant for all angular scales because they result from an effective overall scaling by the mask to flux ratio (see e.g., Fig. A.1) that varies from one field to another. For reference, the error bars that would come from sampling variance of a diffuse signal are displayed in blue and their amplitude decreases (as expected) with the increasing number of available angular modes at higher $k$.

At smaller angular scales, there is a systematic small-scale bias. In Fig. A.3 we check via simulations that this bias either damps or boosts the estimate depending on the relative size of the mask and the FWHM. Table A.1 gives the spatial frequency at which this bias between the measured and the expected APS becomes larger than 20% and 50%, in the interloper contaminated maps masked at $r_\sigma = 3.0$. These values are well above the beam cutoff, and for consistency, in all this work, we restrict to angular scales $k < k_{20\%}$.
<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Average relative error (black crosses) between the [CII] APS measured in the 54 Uchuu sub-fields masked with $r_{\sigma} = 3$ as a function of spatial frequency, $k$, at $z = 5.2$ (top), $z = 6.5$ (middle), and $z = 7$ (bottom). The mean errors over the $k$ modes are <1%, 7%, and 8% for each redshift, respectively. The dispersion (red error bars) normalized by the square root of the number of realizations is also represented. The mean sizes of these normalized error bars over the $k$ modes are 3%, 4%, and 5% for each redshift, respectively. For reference, the equivalent normalized error bars for a perfect Gaussian diffuse signal from the same number of Monte Carlo realizations are over-plotted in blue.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Output (masked) over input (intrinsic) [CII] APS ratio for various beam sizes but with the same mask. The smaller the beam, the larger the $k$ angular mode up to which the output APS equals the input one. The analytical beam is displayed in solid lines for reference and shows that the effect only appears at high angular scales, where the beam has mostly smeared out any signal.</figcaption></figure>

<table><caption>**Table A.1.** Spatial frequency (in arcmin⁻¹) at which the pseudo power spectrum becomes biased by 20% and 50%, and the corresponding angular size, $r$, as a fraction of the FWHM.</caption><thead><tr><th>frequency</th><th>k<sub>20%</sub></th><th>r<sub>20%</sub> [FWHM]</th><th>k<sub>50%</sub></th><th>r<sub>50%</sub> [FWHM]</th></tr></thead><tbody><tr><td>305 GHz</td><td>2.14</td><td>1.36</td><td>2.35</td><td>1.24</td></tr><tr><td>253 GHz</td><td>1.94</td><td>1.24</td><td>2.14</td><td>1.36</td></tr><tr><td>237 GHz</td><td>1.33</td><td>1.71</td><td>1.46</td><td>1.55</td></tr></tbody></table>