© The Authors 2023

1. Introduction

The Epoch of Reionization (EoR) started with the birth of the first stars and the formation of the first galaxies, whose ultraviolet 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, τ, using the cosmic microwave background (CMB) angular power spectrum (APS; Planck Collaboration Int. XLVII 2016). It shows that the Universe was reionized at 10% at z_10% = 10.4 ± 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-α and Lyman-β hydrogen transitions set the late end at z ≈ 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≳ 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 ≲ 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_α 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 [C_II] 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 [C_II] 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 [C_II] 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 extragalactic 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 starbursting (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 ⟨U⟩. L_IR is derived from the direct use of the Kennicutt (1998) conversion factor and SFR. As for the mean radiation field ⟨U⟩, 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 ⟨U⟩ and CO(5-4) / CO(2-1) flux ratio from Daddi et al. (2015).

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

Finally, two models are implemented to compute the [C_II] 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 [C_II] 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, [C_I] and [C_II]. 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

$$\begin{array}{c}\hfill G_{\mathrm{beam}}=\frac{1}{2\pi {\sigma }^2}{\mathrm{e}}^{\frac{-r^2}{2{\sigma }^2}},\end{array}$$(1)

where σ is the standard deviation of the Gaussian beam profile linked to the frequency ν of the channel by

$$\begin{array}{c}\hfill \sigma =\frac{\mathit{FWHM}}{2\sqrt{2ln(2)}}=\frac{1}{2\sqrt{2ln(2)}}\frac{1.22c}{\nu D},\end{array}$$(2)

with FWHM is the full with 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, [C_I], and [C_II],, 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 [C_II] 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 [C_II] signal, as the SIDES simulation only covers the redshift range z ≤ 7, the analysis was done only for ν > 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_ν × N_p, where N_ν 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

$$\begin{array}{c}\hfill C_{\nu {\nu }^{\prime }}={\stackrel{\sim }{X}}_{\nu \mathrm{p}}{\stackrel{\sim }{X}}_{\mathrm{p}{\nu }^{\prime }},\end{array}$$(3)

where $\stackrel{\sim }{X}$ is defined as

$$\begin{array}{c}\hfill {\stackrel{\sim }{X}}_{\nu \mathrm{p}}=X_{\nu \mathrm{p}}-{\mu }_{\mathrm{p}},\end{array}$$(4)

and μ_p is the average in the pixel p over all frequencies. The covariance matrix, C, is square, symmetric, and has dimension N_ν. Next, we computed the eigenvector matrix, E, and the block matrix, λ, that contains the eigenvalues from the diagonalization of C, such that

$$\begin{array}{c}\hfill C_{\nu {\nu }^{\prime }}=E_{\nu \alpha }^{-1}{\lambda }_{\alpha \beta }{E}_{\beta {\nu }^{\prime }}.\end{array}$$(5)

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

$$\begin{array}{c}\hfill T_{\mathrm{p}\nu }={\stackrel{\sim }{X}}_{{\nu }^{\prime }\mathrm{p}}^T{E}_{{\nu }^{\prime }\alpha }^{-1}\frac{1}{\sqrt{{{\lambda }^{\prime }}_{\alpha \nu }}},\end{array}$$(6)

and the continuum is therefore estimated as

$$\begin{array}{c}\hfill {\stackrel{\sim }{X}}_{\nu \mathrm{p}}^{\mathrm{continuum}}=E_{\alpha \nu }^T{\sqrt{\lambda }}_{\alpha \beta }{T}_{\mathrm{p}\beta }^T+{\mu }_{\mathrm{p}}.\end{array}$$(7)

Eventually, the line spectrum is given by

$$\begin{array}{c}\hfill X_{\nu \mathrm{p}}^{\mathrm{lines}}=X_{\nu \mathrm{p}}-{\stackrel{\sim }{X}}_{\nu \mathrm{p}}^{\mathrm{continuum}}.\end{array}$$(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.

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.

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 [C_II] model and 7 times lower than the L18 [C_II] model. But the higher the redshift, the lower the contrast. At z = 7, the residual continuum is at a similar level than the fiducial [C_II] model and 10 times higher than the L18 model.

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 arcmin−1 in frequency channels of 1 GHz width.

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 ν > 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 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:

$$\begin{array}{c}\hfill S(b)={(y-b)}^{T}W(y-b)+R{b}^T{D}^{T}Db,\end{array}$$(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:

$$\begin{array}{c}\hfill w_i=\{\begin{array}{c}f(y_i-b_i,m,\sigma ),y_i\ge b_i\hfill \\ 1,y_i<b_i,\hfill \end{array}\end{array}$$(10)

with

$$\begin{array}{c}\hfill f(d^{-},m,\sigma )=\frac{1}{1+\exp (2(d-(-m+2\sigma ))/\sigma )},\end{array}$$(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 σ 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 ν ≥ 237 GHz (z_[CII] ≤ 7).

Fig. 3. Relative error of arPLS to the ‘true’ continuuum as a function of frequency. The mean is taken over the spectra of all the los of the mock cube.

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

For ν ≤ 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 [C_II] 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.

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 [C_II]. Dominant interlopers are the CO rotational ladder and two [C_I] fine structure lines from foreground galaxies at z < 3. The [O_I] 145 μm, [N_II] 121.9 μm and N_II 205.2 μm lines are also present in the CONCERTO frequency range but their intensities are negligible with respect to the [C_II] intensity (Silva et al. 2015). Consequently, only CO and [C_I] 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 ∼ 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 ≤ 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₈₅₀. 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_{\ast }^{\mathrm{L}16}$. 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 × 10¹⁰ L_⊙ at z = 0.3, 3 × 10¹⁰ L_⊙ at z = 0.6, and 1 × 10¹¹ L_⊙ 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⁻¹ with a completeness below 10% at this flux. Completeness is close to 100% above 5 mJy beam⁻¹, which is largely above (between 1 and 3 orders of magnitude) the S₈₅₀ requested depth.

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* − LIR, and M* − S850. 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).

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 [C_II] redshifts and interloper redshifts is given in Table 2.

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[CII] = 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 [CII] level (in percent) is represented with the dashed green line. Power amplitudes are averaged for 0.12 < k < 0.24 arcmin−1.

For the three redshifts, galaxies more massive than ≥10^8.1 M_⊙ 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_⊙ are responsible for less than 0.5% of the total CO power. For example, we can highlight that the galaxy population with 10^10.6 M_⊙ < M_* ≤ 10^10.9 M_⊙ 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_⊙, 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 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_σ, expressed in units of the Gaussian beam width parameter, σ (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_σ. Figure 6 displays examples of these masks.

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.

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 measurements. In this work, we used the P of K EstimatoR (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 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_[CII] = 5.2, 6.5, and 7 for different r_σ. With r_σ = 3, the interloper power amplitude is smaller than 2 × 10⁻² Jy² sr⁻¹ for all three redshifts, that is, it is reduced by a factor of (100 ± 30), (50 ± 16), and (42 ± 15), for each redshift, respectively.

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+[CI] + [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σ = 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+[CI] 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−1.

For [C_II], 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 [C_II] and residual interlopers of [C_II]/(residual CO+[C_I]) = (30 ± 16) (i.e. 3% of residual contamination). The contrast goes up to (62 ± 32) with r_σ = 3. At z_[CII] = 6.5, the [C_II] APS is overestimated. Even masking with r_σ = 3 does not achieved a better power ratio than [C_II]/(residual CO+[C_I]) = (5.5 ± 3.6) (i.e. 18% of residual contamination after masking). Finally, at z_[CII] = 7, the contrast between [C_II] and residual interlopers is weak due to the low [C_II] APS, causing [C_II]/(residual CO+[C_I]) = 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 [C_II] 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 [C_II], 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 ([C_II], 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.

Fig. 8. APS (in Jy2 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+[CI]), the bright interlopers (i.e. above the mass threshold), the faint interlopers, and [CII]. The intrinsic power spectra of each (unmasked) component are shown with solid black lines, and with a solid green line for [CII]. 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).

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⁻³ Jy² sr⁻¹ and becomes negligible compared to that of the faint sources that stays about 0.8 × 10⁻² Jy² sr⁻¹. The final power ratio (with r_σ = 3) between faint and bright sources are (faint CO+[C_I] )/(bright CO+[C_I]) = (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_I] and CO maps separately. Sun et al. (2018, hereafter S18) studied the masking of CO emitters having a magnitude $m_K^{\text{AB}}\lesssim 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 ≥100 at z = 6.5 and k = (3 − 6)  ×  10⁻² h Mpc⁻¹, for a 8 % loss of survey area. With SIDES, at comparable redshift and spatial frequencies (corresponding to k = 0.12 − 0.24 arcmin⁻¹ 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 strategy of masking all detected sources that have a magnitude $m_K^{\text{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_[CII] = 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_II] APS reconstructed from the masked [C_II] maps of SIDES-Bolshoi. The reconstructed [C_II] power amplitude is consistent with the intrinsic one within 1σ. However, Appendix A details how the measured [C_II] APS (using the mask r_σ = 3) are biased low, on average by 1%, 7%, and 8% at z_[CII] = 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_II] power spectrum from the SIDES-Bolshoi realization at z_[CII] = 5.2, reaching a power ratio [C_II]/(residual CO) = (80 ± 42) and [C_II]/(residual [C_I]) = 280 ± 150, with a bias on the [C_II] APS of 1% due to the mask. At z_[CII] = 6.5, the power ratios are [C_II]/(residual CO) = (7.8 ± 5.1) and [C_II]/(residual [C_I]) = (25 ± 16), with a 7 % bias on the [C_II] APS due to the mask. At z_[CII] = 7, power ratios are [C_II]/(residual CO) = (2.7 ± 1.9) and [C_II]/(residual [C_I]) = (9.4 ± 6.6), with an 8 % bias on the [C_II] APS. For all redshifts, residual contamination becomes dominated by the faint sources for r_σ ≥ 2 since such sources are moderately affected by the mask.

5.2. Field-to-field variance effects on residuals

The detection of the [C_II] 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 [C_II] 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 Ω = 1 deg² field to the Ω = 2 deg² field of the SIDES-Bolshoi simulation, assuming that the mean power spectrum μ does not change between the two field sizes:

$$\begin{array}{c}\hfill f=\frac{\sigma (\mathrm{\Omega }=2{\mathrm{deg}}^2)}{\sigma (\mathrm{\Omega }=1{\mathrm{deg}}^2)}.\end{array}$$(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σ FtF variance interval in a 2 deg² field, for [C_II] and the faint interlopers at z_[CII] = 5.2, 6.5, and 7. The [C_II] 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 [C_II] in comparison to the [C_II] variance. Table 5 lists the level of faint interloper contamination and its ratio to the 1σ value of the [C_II] FtF variance. As a result, residual contamination is smaller than 20% of the [C_II] 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).

Fig. 9. Mean APS (solid lines) and the 1σ FtF variance interval (coloured area) for a 2 deg2 field, for [CII], faint CO, faint [CI], 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.

6. Summary and conclusion

We tested a complete foreground deconfusion process for the [C_II] 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⁷ M_⊙ 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 [C_II] were produced. A foreground-contaminating mock cube was created by adding the separate foreground cubes to the [C_II] 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 [C_II] models provided by SIDES up to z = 7. For 5.2 ≤ z ≤ 7, the residual CIB power amplitudes are lower than that of our fiducial DL14 (L18) [C_II] 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_I] proxy available in the COSMOS field. We also investigated its use as a CO proxy. The most contaminating sources (M_* > 10^8.1 M_⊙) 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_I] 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 × 10⁻² Jy² sr⁻¹, 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_[CII] = 5.2, the power ratio is [C_II]/(residual CO+[C_I]) = 62 ± 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 ≤ 4% of FtF variance).

At z_[CII] = 6.5, the residual contamination after masking is a factor of [C_II]/(residual CO+[C_I]) = 5.5 ± 3.6 below the [C_II] power spectrum amplitude. This residual contamination and the reconstruction by POKER (see Appendix A) led to an overestimate 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_[CII] = 7, SIDES predicts a drop in [C_II] power spectrum amplitude by a factor of 16 with respect to z_[CII] = 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 ± 1.4 below the [C_II] power spectrum amplitude.

Hence, residual contamination is still important above z_[CII] > 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 δz/(1 + z) = 0.025 below z ≤ 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 ≤ 3, including at cosmic noon (2  ≤  z  ≤  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.

Acknowledgments

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

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 ϕ. 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.

Fig. A.1. Illustration of the POKER correction for two different 114 arcmin2 fields. Top row: CO+[CI] masks corresponding to the two fields, as well as the [CII] sources (red points). Bottom row: APS of the masked [CII] 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 [CII] 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.

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 deg² cut in the Uchuu simulation, we built a mask with r_σ = 3 following Sect. 4 and computed the ratio between the output (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_[CII] = 5.2, 6.5, and 7), up to an angular scale of k_lim = 2.14, 1.94, 1.33 arcmin⁻¹, 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_[CII] = 5.2, 6.5, and 7, respectively, and are also used for smaller r_σ 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.

Fig. A.2. Average relative error (black crosses) between the [CII] APS measured in the 54 Uchuu sub-fields masked with rσ = 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.

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_σ = 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%.

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.

All Tables

All Figures

	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.
In the text

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 arcmin−1 in frequency channels of 1 GHz width.

In the text

	Fig. 3. Relative error of arPLS to the ‘true’ continuuum as a function of frequency. The mean is taken over the spectra of all the los of the mock cube.
In the text

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* − LIR, and M* − S850. 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).

In the text

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[CII] = 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 [CII] level (in percent) is represented with the dashed green line. Power amplitudes are averaged for 0.12 < k < 0.24 arcmin−1.

In the text

	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.
In the text

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+[CI] + [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σ = 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+[CI] 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−1.

In the text

Fig. 8. APS (in Jy2 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+[CI]), the bright interlopers (i.e. above the mass threshold), the faint interlopers, and [CII]. The intrinsic power spectra of each (unmasked) component are shown with solid black lines, and with a solid green line for [CII]. 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).

In the text

	Fig. 9. Mean APS (solid lines) and the 1σ FtF variance interval (coloured area) for a 2 deg2 field, for [CII], faint CO, faint [CI], 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.
In the text

Fig. A.1. Illustration of the POKER correction for two different 114 arcmin2 fields. Top row: CO+[CI] masks corresponding to the two fields, as well as the [CII] sources (red points). Bottom row: APS of the masked [CII] 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 [CII] 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.

In the text

Fig. A.2. Average relative error (black crosses) between the [CII] APS measured in the 54 Uchuu sub-fields masked with rσ = 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.

In the text

	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.
In the text