© The Authors 2023

1. Introduction

The cosmic microwave background (CMB) photons we observe have been influenced by their interaction with the baryonic matter as they travel through the Universe. These interactions modify the CMB photons and are called secondary anisotropies of the CMB. The study of these secondary anisotropies allows for the extraction of information about the large-scale structure (LSS) of the Universe as well as about the primary anisotropies, the matter content of the Universe. One of the secondary anisotropies, the thermal Sunyaev–Zeldovich (tSZ) effect, is due to the inverse Compton scatting of the CMB photons by the hot electrons in the intracluster medium (ICM) of galaxy clusters. This effect is a crucial tool for cluster detection and was first used by the South Pole Telescope (SPT; Staniszewski et al. 2009) and the Atacama Cosmology Telescope (ACT) Collaborations (Menanteau et al. 2010) and later by the Planck Collaboration (Planck Collaboration XXVII 2016) to construct cluster catalogues. The tSZ effect can be used to measure cluster properties, such as temperature (Erler et al. 2018). It is also crucial for cosmology, as it has a strong dependence on matter clustering σ₈ and matter density (Komatsu & Seljak 2002; Planck Collaboration XXII 2016; Bolliet 2018). By fitting the tSZ power spectrum extracted from data with models of the signal and its contaminants, assuming a Gaussian likelihood function, and applying a Markov chain Monte Carlo (MCMC) analysis with a fixed mass bias b, one can constrain these cosmological parameters (Planck Collaboration XXI 2014). This implies that a small improvement in the tSZ power spectra statistical significance or a reduction of the bias due to residual contamination from other signals will lead to a significant improvement in measurements of σ₈ and Ω_m. Therefore, a clean, uncontaminated tSZ map is essential.

The cosmic infrared background (CIB) is the emission from stars and galaxies absorbed and re-emitted by dust in the infrared domain. It is the main contaminant in tSZ reconstruction and is also a tracer of the matter distribution and the LSS (Hurier 2015). The CIB contamination, or leakage, originates from the spatial correlation between CIB and tSZ and is also due to the component separation techniques, such as internal linear combination (ILC; Bennett et al. 2003; Eriksen et al. 2004; Park et al. 2007; Kim et al. 2008; Delabrouille et al. 2009), not being able to perfectly separate the signal from the many other astrophysical emissions present in the sky. Previous studies on CIB leakage in the Planck tSZ map made by Planck Collaboration XXII (2016) showed that the amplitude of the leakage increases when the CIB originates from higher redshift populations. The dependence of the leakage on the redshift of the CIB population comes from the fact that the separation methods of components primarily focus on reducing the main contaminant, which at high frequencies is Galactic thermal dust, and is thus less efficient at reducing the high-z CIB contaminants. This could lead to a bias in tSZ estimation for high-z clusters, as mentioned in Vikram et al. (2017). Another study by Alonso et al. (2018) showed that the error bars on the tSZ profile are not much improved by the deprojection (i.e. nulling) of the CIB with more constrained component separation methods such as constrained internal linear combination (CILC; Remazeilles et al. 2011). Dust contamination can reduce the number of detected galaxy clusters by up to ∼9% within z ∈ [0.3, 0.8] (Melin et al. 2018). Similarly, Zubeldia et al. (2023) demonstrated that CIB contamination in the tSZ signal can lead to ∼20% fewer clusters being detected and proposed a publicly available code to produce CIB-free cluster catalogues.

In the past three decades, several space missions, such as the COsmic Background Explorer (COBE; Mather et al. 1994), the Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2003), and Planck (Planck Collaboration I 2014), have observed the sky at microwave frequencies and extracted full-sky maps of the CMB, tSZ, CIB, and other astrophysical components. Upcoming ground-based experiments, such as those using the Simons Observatory Large Aperture Telescope (SO LAT; Ade et al. 2019) and the Fred Young Submillimeter Telescope (FYST; CCAT-prime Collaboration 2023), will allow for unprecedented measurements of the tSZ and CIB thanks to their high spatial resolution and high sensitivities.

In this work, our focus is on quantifying the benefits of combining future experiments, such as SO LAT and FYST to reduce the CIB leakage in ILC-reconstructed tSZ maps. A cleaner tSZ map will allow better constraints on the cosmological parameters σ₈ and Ω_m to be drawn, more accurate cluster counts and detection, and more accurate measurements of the properties of clusters, such as their temperature. Without mitigation, CIB leakage along the line of sight (LoS) could systematically bias measurements of the tSZ effect and lead to higher tSZ measurements, or it could lower such measurements by filling the characteristic tSZ decrement at frequencies below 217 GHz (Schaan et al. 2021). Obtaining less contaminated signals by combining microwave experiments is common and will be part of the observation strategies of future experiments. For example, previous studies, such as Ade et al. (2019), have shown that SO LAT will achieve precise measurement of the tSZ power spectrum when combined with Planck data. As Planck observations are not contaminated by the atmosphere, this combination allows SO LAT to mitigate such contamination on large scales. Similarly, Melin et al. (2021) constructed a bigger cluster catalogue by combining Planck and the South Pole Telescope (SPT-SZ) data. Our work also encourages the combination of all available data from previous and upcoming CMB experiments.

This paper is structured as follows: In Sect. 2, we formalise the tSZ and CIB description in our study. In Sect. 3, we describe the modelling of high-resolution maps of the different microwave sky emissions. In particular, we first explain how we generated Galactic foregrounds, the extragalactic components, and their processing. Then, we detail the generation of noises and the handling of beams and masks. This is followed by a conclusion of the output of the sky modelling. In Sect. 4, we summarise the methods used to extract the tSZ effect, starting with the ILC and followed by CILC. Next, we explain beam convolution handling and the treatment of full-sky maps, and we provide an overview of the methods used. We present a detailed description of our results and an overview in Sect. 5. Section 6 contains a summary and discussion of the above-mentioned results.

2. Formalism

2.1. The thermal Sunyaev–Zeldovich effect

The thermal tSZ effect, predicted by Sunyaev & Zeldovich (1970, 1972), is the inverse Compton scattering of CMB photons on the hot electrons present in the ICM. A consequence of this effect is the distortion of the CMB blackbody spectrum, due to the shift of the low frequency (corresponding to low energies) CMB photons to higher frequencies (corresponding to higher energies). The temperature change of the CMB due to the tSZ effect can be expressed as:

(1)

where T_CMB = 2.72548 ± 0.00057 K (Fixsen 2009) is the temperature of the CMB and y is the dimensionless Compton-y parameter measuring the electron pressure integrated over the LoS proper distance l,

(2)

where σ_T is the Thomson cross-section, m_e the electron mass, and P_e(l) = n_e(l)k_BT_e(l) the electron pressure composed of the electron density n_e, the Boltzmann constant k_B, and the electron temperature T_e. In Eq. (1), f(x_ν) is the spectral shape of the tSZ (Sunyaev & Zeldovich 1972, 1980):

(3)

where x_ν = hν/k_BT_CMB, h being the Planck constant and ν the frequency. The frequency spectrum of the tSZ is characterised by a decrement in the intensity change for frequencies lower than ν ∼ 217 GHz and an increment at higher frequencies. To express the frequency variations induced by the tSZ effect in terms of an intensity variation, we use the conversion factor to pass from units of thermodynamic temperature (K_CMB; Planck Collaboration IX 2014) in Eq. (1), to intensity units (Jy sr⁻¹),

(4)

where B_ν(T) is the Planck blackbody spectrum at a fixed frequency ν, I₀ = 2(k_BT_CMB)³/(hc)² ≈ 270 MJy sr⁻¹, c is the speed of light, and

(5)

Therefore, the tSZ intensity change is given by

(6)

More details on the Sunyaev–Zeldovich (SZ) effects can be found in the review by Mroczkowski et al. (2019). In this study, we place ourselves in a simplistic approach and do not take into account the relativistic correction of the tSZ effect (Wright 1979; Itoh et al. 1998; Pointecouteau et al. 1998; Chluba & Sunyaev 2012), even though Planck data have shown that this leads to an underestimation of the amplitude of its power spectrum (Remazeilles et al. 2018).

2.2. The cosmic infrared background

The CIB is the extragalactic integrated emission produced by all stars and galaxies during their formation and throughout their life. The CIB frequency spectrum peaks in the infrared regime because the majority of sources are far and their light is redshifted by the expansion of the Universe but also because some of this light is absorbed by dust and re-emitted into the infrared regime. This existence of the CIB was predicted by Partridge & Peebles (1967) and first detected by Puget et al. (1996), Fixsen et al. (1998), Hauser et al. (1998). Contrary to the tSZ, one frequency spectrum is not sufficient to fully model it. The CIB spectral energy distribution (SED) depends on the redshift of sources. It can be approximated as a modified blackbody spectrum:

(7)

where A_dust is the amplitude of the dust emission, β is the dust emissivity spectral index, and T_dust is the temperature of the dust grains. The CIB is a tracer of the entire history of star formation in the Universe. One characteristic of the CIB of particular importance for our study is its non-zero correlation with the tSZ signal. Indeed, some of the unresolved dusty star-forming galaxies constituting the CIB are found in galaxy clusters. Therefore, the spatial correlation between tSZ and CIB is non-zero, and it was measured in Reichardt et al. (2012, 2021), Dunkley et al. (2013), George et al. (2015), Choi et al. (2020b).

3. Modelling the microwave sky

We set up a pipeline called Skymodel to create mock high-resolution maps of the microwave sky that would be observed in future experiments. Our pipeline is based on Hierarchical Equal Area isoLatitude Pixelisation (HEALPix) template maps of the microwave sky from existing simulations as well as the most recent observations of the microwave sky. The microwave sky components are usually separated into two categories: the emissions coming from the Galaxy, called Galactic foregrounds, and the extragalactic components. Our Skymodel deals differently with these two categories, using an existing Python package for the Galactic foregrounds and HEALPix maps from mock simulations for the extragalactic ones.

3.1. Galactic foregrounds

Galactic foregrounds are the various emissions originating from the Milky Way and are thus in the foreground of extragalactic emissions. We considered the following Galactic foregrounds: Galactic synchrotron; free-free emission; anomalous microwave emission (AME), which is supposedly produced by spinning dust (Dickinson et al. 2018); and thermal emission from Galactic dust grains.

All of the above-mentioned components were generated using the Python Sky Model (PySM; Thorne et al. 2017). The PySM¹ package is based on Galactic foreground template maps obtained from Planck Collaboration XXV (2016) that were extracted using the COMMANDER algorithm (Eriksen et al. 2004, 2008; Planck Collaboration IX 2015; Planck Collaboration IV 2020) on Planck 2015 raw data, data from WMAP, and 408 MHz all-sky survey by Haslam et al. (1981, 1982), reprocessed by Remazeilles et al. (2015). The original data maps were at a HEALPix (Górski et al. 2005) pixel resolution of pix_size = 13.74′ (N_side = 256). The resolution of a HEALPix pixel is defined as the square root of the pixel area. The PySM upgrades them to 6.9′ (N_side = 512) and adds small-scale fluctuations generated randomly by extrapolating the power spectra of those maps (Miville-Deschênes et al. 2007). We generated each foreground component at the desired frequency using the basic spectral model of PySM (i.e. “s1”, “f1”, “d1”, and “a1”). These maps were then upgraded to a higher N_side = 4096 and smoothed with a 10′ full width at half maximum (FWHM) Gaussian to mitigate pixelation artefacts (see Fig. 1).

Fig. 1. Projections of the PySM sky at 353 GHz composed of free-free, synchrotron, AME, and Galactic dust. The upper row shows cutouts around the Galactic south pole, and the lower row shows cutouts around the Galactic centre. Each cutout is composed of 600 × 600 pixels of the size 0.86′, making 8.6° ×8.6° maps. The left column shows cutouts from the original PySM sky, which feature pixelation artefacts due to their pixel size of 6.9′ (Nside = 512). The right column shows modified cutouts obtained by first updating the original map to a higher resolution, 0.86′ (Nside = 4096), by oversampling using the HEALPy function ‘ud_grade()’ and then smoothing with a Gaussian kernel of a width that is close the original pixel size of the PySM maps (FWHM = 10′). This reprocessing of the PySM maps removes the pixelation artefacts while only marginally lowering the spatial resolution.

3.2. The extragalactic components

The extragalactic components include the CMB, tSZ, CIB, and the kinematic Sunyaev–Zeldovich (kSZ) effect. These components were generated at various frequencies using template maps provided by the Websky extragalactic CMB mocks (Stein et al. 2020), which are HEALPix full-sky maps of the extragalactic sky at N_side = 4096 (pix_size = 0.86′). The large-scale structure of the Websky simulations was obtained through a peak-patch and 2nd-order Lagrangian perturbation theory (2LPT) simulation with 12288³ particles and a box size of 15.4 h⁻¹ Gpc. The authors assumed a flat Λ cold dark matter (ΛCDM) cosmology with h = 0.68, Ω_cdm = 0.261, Ω_b = 0.049, Ω_Λ = 0.69, n_s = 0.965, and σ₈ = 0.81. The simulations can be downloaded from their website². An overview of the provided template maps, their frequencies, and units can be seen in Table 1.

To generate the tSZ map, Stein et al. (2020) projected the pressure profiles from hydrodynamical simulations (Battaglia et al. 2012) onto the mass-Peak Patch halo catalogue. They did not consider field contributions, and the minimum halo mass considered to compute their pressure profile was ∼1 × 10¹³ M_⊙. Their tSZ power spectrum is in good agreement with the Planck data (Stein et al. 2020). Websky provides all-sky, dimensionless Compton-y maps (linked to the tSZ effect through Eqs. (1) or (6)), which can be multiplied by the tSZ spectrum in Eqs. (3) and (5) to obtain frequency maps of the tSZ. The tSZ signal was simulated by Websky in its non-relativistic limit.

The CIB was generated using the halo occupation distribution model proposed by Shang et al. (2012) and with the parameters given in Viero et al. (2013) to fit the CIB power spectrum with the Herschel observations. Within this model, the rest-frame SED of a source is a modified blackbody at low frequencies and a power law at high frequencies. The CIB’s SED depends on the source redshift, halo mass, and frequency at which it is observed. The Websky halo catalogue is populated with central and satellite galaxies. The field contribution was not considered. The Websky simulations generated CIB maps at the Planck High-Frequency Instrument (HFI) frequencies (see Table 1), and in particular, the 545 GHz CIB power spectrum was normalised to fit the Planck 545 GHz CIB power spectrum (Planck Collaboration I 2014) at ℓ = 500.

Because of the redshift dependence of the CIB distribution, a decorrelation was expected between the CIB at different frequency channels. The decorrelation was determined in the Websky simulations (Stein et al. 2020) by computing the Pearson product-moment correlation coefficients on the power spectrum of the CIB maps, evaluated at the Planck HFI frequency channels. The correlation matrix between the different CIB maps provided by the Websky simulations is illustrated in Fig. 2. The Websky CIB decorrelation values are in agreement with Planck Collaboration I (2014) and Lenz et al. (2019).

Fig. 2. Correlation matrix between the Websky CIB maps provided at the Planck HFI frequencies. The axis shows the frequencies of the maps in GHz, and the dimensionless colour bar shows the correlation coefficient.

In this work, we extrapolated the CIB model to arbitrary frequencies by fitting each pixel of the existing Websky CIB maps with the following modified blackbody:

(8)

where ν₀ = 353 GHz. This was done using the Scipy “curve_fit()” function to derive three all-sky parameter value maps for the variables A_CIB, T_dust, and β. These full-sky maps were injected in Eq. (8) to produce Websky-like CIB intensity maps at new frequencies.

In the first approximation, we assumed this exact modified black body across frequencies, which is not fully realistic since the CIB is composed of different sources at different redshifts, and their SED is not the same. The Websky simulation includes cross-correlation between the tSZ and CIB, which is in agreement with the Planck data results (see Fig. 11 of Stein et al. 2020).

For the CMB component, the Websky simulations provided us with two files containing random realisations of the primary T/Q/U a_ℓm coefficients, one being unlensed and one being lensed. The lensed file was obtained using the unlensed a_ℓm file and the CMB lensing convergence from 0 < z < 1100 in the Born approximation. From those a_ℓm, we could reconstruct a lensed and an unlensed CMB map in μK_CMB. Using those maps and the conversion factor (Eq. (4)), we could obtain a CMB map at any desired frequency. In the first approximation, we used the unlensed CMB. For the kSZ component, the situation is analogous to the tSZ one. The Websky templates provided frequency-independent maps of the kSZ, which can be seen as a ‘y_kSZ-map’. We multiplied the maps by the conversion factor in Eq. (4) to get a kSZ map at any frequency. The Websky simulations do currently not provide radio point source (RPS) templates, but the Websky team will add them in a future release. This extragalactic component is therefore not accounted for in our study. The different extragalactic full-sky map templates were all converted to the same unit, that is, either temperature units (K_CMB) or intensity units (Jy sr⁻¹) units, using Eq. (4).

3.3. Future experiments

We simulated the microwave sky as will be seen by two future ground-based telescopes. One of these telescopes is the FYST of the Cerro Chajnantor Atacama Telescope prime (CCAT-prime) Collaboration (CCAT-prime Collaboration 2023). It is an upcoming 6 m diameter telescope that will be located at 5600 m on the mountain of Cerro Chajnantor in Chile. Due to its very low precipitable water vapour, this exceptional site will provide one of the best atmospheric transmissions for ground-based microwave observations to date (Bustos et al. 2014). The crossed-Dragone optical design proposed by Niemack (2016) will deliver a high-throughput, wide-FoV telescope capable of rapidly scanning large areas of the sky. The first light is expected in 2024.

The second telescope of interest will be part of the Simons Observatory (SO; Ade et al. 2019). The SO will be composed of three refracting 0.4 m small aperture telescopes (SATs) and one cross-Dragone 6 m large aperture telescope (LAT). The four instruments will be located at 5300 m on Cerro Toco in the Atacama Desert in Chile.

We studied what benefits may result from the combination of FYST with the SO LAT. To do so, we used PySM and the Websky mock sky to simulate the microwave sky at frequencies and sensitivities, meaning the smallest signal a telescope can detect above the random background noise, as well as the beam sizes of future observations conducted with these telescopes, see Table 2.

3.4. Noises, masks, and beams

Our Skymodel pipeline includes two types of noise: white noise, to reproduce the instrumental noise, and atmospheric red noise, which is present for ground-based experiments. The noise implementation is based on the SO Collaboration (Ade et al. 2019) noise curves forecasted and adapted for CCAT-prime by Choi et al. (2020a):

(9)

where ℓ_knee = 1000 and α_knee = −3.5. The values for N_red and N_white for the SO LAT and FYST are included in the Skymodel. For this reason, the atmospheric noise model can only be computed at valid SO LAT frequencies or FYST frequencies (see Table 2) and takes the form of a series of power spectra, as defined by Eq. (9) and shown in Fig. 3. Using those power spectra, the Skymodel generates full-sky map realisations of the atmospheric noise or atmospheric and instrumental noise. This is done through the HEALPy function ‘synfast()’, which generates map realisations from the power spectrum while assuming the field to be Gaussian. This was done for each frequency band. In reality, neighbouring frequency bands will observe similar regions of the sky at the same time. Thus, the atmospheric noise is partially correlated. We adopted a 90% correlation between each pair of neighbouring frequency channels as suggested in Ade et al. (2019). Fixing the ‘synfast()’ seed allowed for the production of a 100% correlated map for two neighbouring channels. Another seed and a correction factor can be applied to the power spectrum to produce a down-weighted correlated map and an uncorrelated map of one of the two channels. This new map and the fully correlated map can be linearly combined to produce a 90% correlated map using

(10)

Fig. 3. SO LAT (blue lines) and FYST (red lines) sum of the atmospheric red and instrumental white noise spectra based on the noise model presented by Ade et al. (2019) and adopted by Choi et al. (2020a). The noise curves have been beam corrected.

where M is the final map at the frequency ν, which will be 90% correlated with its neighbouring channel. The HEALPy synfast() function is represented by h, and is the atmospheric noise curve, which is 100% correlated when produced from the same seed as its neighbouring channel and 0% when produced from another seed.

The frequency bands are broad and extended over some frequency range, and their transmission is not perfect but varies over their extent. In the first approximation, we considered the frequency bands to be delta functions with perfect transmission.

In addition to the maps of astrophysical emission and noise, the Skymodel uses survey masks. The masks are stored at N_side = 256 and were upgraded to N_side = 4096 in order to correspond to the Websky map resolution. In this study, we used a mask that excludes the Galactic plane with a sky fraction of f_sky = 0.6 (see Fig. 4). This mask was derived in Erler et al. (2018) by masking the 40% brightest pixels of a Planck Galactic emission map. This is an initial approximation, as in reality SO LAT and FYST will observe a smaller fraction of the sky.

Fig. 4. Galactic foreground mask derived from Planck data by Erler et al. (2018; fsky ≈ 0.6). The white part is the masked region, and the blue part is what was used in our analysis. Overplotted is the HEALPix tessellation scheme for Nside = 4.

The Skymodel pipeline also simulates the effect of telescope beams using the Gaussian FWHMs given in Table 2. The product of the Skymodel is a full-sky HEALPix map containing all the Galactic foreground components, the extragalactic components, instrumental white noise, and atmospheric red noise at a given frequency and assuming a given beam size (see Appendix D).

4. Method to extract thermal Sunyaev–Zeldovich map

To study the power spectrum measurement of the tSZ effect for the SO LAT and FYST telescopes, a tSZ map needs to be extracted from the simulated microwave sky. The sky is composed of various Galactic and extragalactic emissions that are considered to be contaminants of our signal of interest. Component separation methods permit the recovery of a signal from the multi-component microwave sky. Internal linear combination (ILC) is a multi-frequency component separation technique that was first used by Bennett et al. (2003) to extract CMB maps from WMAP data and later refined by Eriksen et al. (2004). Multiple flavours of ILC exist. In their paper, Tegmark et al. (2003) suggested applying it in Fourier space rather than in map space. Needlet ILC (NILC; Delabrouille et al. 2009) processes the maps in both spatial and harmonic space by proposing to filter each multi-frequency map by a set of ℓ-dependent window functions called needlets in order to apply multiple ILCs on different multipole intervals. This makes the ILC weights ℓ-dependent, allowing them to take advantage of the scale-dependent distribution of astrophysical components to minimise the noise and better extract the signal of interest. The Modified ILC Algorithm (MILCA; Hurier et al. 2013) differentiates itself by allowing multiple constraints on the Astrophysical emissions in order to reduce their contamination. Scale-discretised directional wavelet ILC (SILC; Rogers et al. 2016) filters the map with ℓ windows that are direction dependent.

In this work, we used ILC and CILC methods in a map-based space. It is also possible to use these methods directly in Fourier space, which is computationally faster but does not allow for spatial separation of signals. Fourier space ILC/CILC only has information on the power spectra of the component in ℓ-space and no spatial information on the localisation of the signal over the full-sky map.

4.1. Internal linear combination

Internal linear combination is widely used in tSZ studies (Bennett et al. 2003; Fornazier et al. 2022). Its main advantage is that it is a so-called “blind” method because only the SED of the component of interest must be known, and no information on the contaminants is needed. The ILC supposes that (1) the maps are composed of a linear combination of astrophysical components and noise, and (2) the astrophysical components are uncorrelated.

A set of N_obs observed map m_i(p) at a frequency i and a pixel p can be expressed as

(11)

where a_i is the amplitude of the frequency spectrum of the signal of interest s at each frequency i. The noise at a given frequency i and pixel p is represented as n_i(p).

To retrieve the signal of interest s, an estimator is defined as

(12)

where ω_i is the ILC weight at the frequency i. The ILC weights are defined such that they follow two conditions. One of them is that in order to obtain an unbiased estimate, the ILC weights must have a unit response to the frequency spectrum of the signal of interest s characterised by a_i at each frequency i (see Eq. (13)):

(13)

The other condition is that the minimisation of the presence of the unwanted astrophysical contaminants and the instrumental noise in the reconstructed map of the signal s is achieved through the condition of minimum variance of the reconstructed map, VAR(). As shown in Eriksen et al. (2004), because the components are uncorrelated, we have VAR(, where is the value of the covariance of the observed maps between the frequencies j and k. Thus, the minimum variance condition is given by

(14)

The minimum variance under the constraint of Eq. (14) can be achieved using the Lagrange multiplier method to solve the above equation. The solution gives the following ILC weights:

(15)

where W is the vector containing the ILC weights, and A is the vector containing the frequency spectrum information of the signal of interest and is called ‘the mixing vector’. The covariance matrix of the set of maps is represented by . The map containing our signal of interest is given by

(16)

From that, we can define the residual noise (r_N), that is, the noise left over by the ILC in the reconstructed map. This noise can be defined as the difference between the ILC-recovered map and the true signal:

(17)

The contaminants are not absent from the recovered signal map. Moreover, the low variance contaminants might not be reduced in the resulting map. This becomes a problem when assumption two of the ILC is not fully verified and our target signal is spatially correlated with another emission, such as for the tSZ and CIB. This is the case between CIB and tSZ. The remaining CIB signal in the map can bias the ILC estimate of the tSZ. One way to reduce this effect is through CILC.

4.2. Constrained internal linear combination

The CILC method was first introduced by Remazeilles et al. (2011). It adds one extra constraint to the ILC approach. If the spectral shape of one of the contaminants is known, it can be taken into account in the modelling such that instead of Eq. (11), it is supposed that

(18)

where s(p) is the signal of interest and q(p) is the contaminant we want to suppress. Respectively, a_i and b_i are the intensity of the frequency spectrum evaluated at a frequency i of the signal and contaminant. All the remaining astrophysical contaminants and the instrumental noise at a given frequency i and pixel p are contained by n_i(p). The CILC weights are defined such that

(19)

The CILC weights have a unit response to the SED of the signal of interest s and a null response to the SED of the astrophysical contaminant q we want to suppress so that this contaminant is absent from the reconstructed signal.

The condition of minimum variance (see Eq. (14)) stays the same, and only the expression of the weights is changed to take into account the two conditions imposed by the CILC,

(20)

where W is the vector containing the CILC weights, B is the mixing vector of the contaminant, A is the mixing vector of the signal, and is the covariance matrix of the set of maps. Those weights are then applied to the set of multi-frequency maps using Eq. (16) to recover a map containing the signal of interest free from the nulled contaminant and with minimum contamination from other astrophysical signals and instrumental noise. The CILC method can be generalised to null more than one contaminant if their SED is known. However, these additional constraints reduce the degree of freedom available for the condition of minimum variance to be verified (Hurier et al. 2013), thus generating increased noise for the CILC compared to the ILC. Similar to the ILC, the residual noise can be defined as the difference between the CILC-recovered map and the signal of interest (see Eq. (17)), except that the weights are given by Eq. (20).

4.3. Beam convolution

The set of multi-frequency maps was smoothed with a Gaussian kernel to the lowest frequency channel resolution so that they could be linearly combined by the ILC. No meaningful physical information is available on scales smaller than the resolution of each map. When linearly combining the maps, the upper limit of meaningful physical information was therefore given by the lowest resolution map.

More sophisticated component separation methods, such as NILC, allows the original resolution of each frequency map to be kept. Indeed, the NILC weights are a function of the multipole ℓ for each frequency map and can therefore be set to zero when the scale exceeds the resolution of the map. Setting the weights to 0 above a given ℓ, allows the lower-resolution maps to contribute to the signal reconstruction only up to their beam-limited scale, leaving the higher-resolution maps to be the contributors of signal reconstruction on smaller scales.

We note that when simulating the microwave sky, only the astrophysical components and atmospheric noise are limited by the instrument beam. Instrumental noise is not beam limited. The inclusion of the beam was done using HEALPy to smooth the maps a_ℓ, m in the multipole space and multiply them with a ℓ-dependent Gaussian kernel. Our set of multi-frequency maps can be expressed as

(21)

where are the a_ℓ, m coefficient of the smoothed map at the frequency i. The Gaussian beam in the multipole space at the frequency i is indicated as b_i(l). Next, a_ℓ, m(s) is the a_ℓ, m decomposition of component of interest map, and a_ℓ, m(c_i) represents all astrophysical contaminants and the atmospheric noise at a frequency i. Finally, a_ℓ, m(n_i) is the instrumental noise a_ℓ, m at the frequency i. In general, to apply an ILC, all maps have to be smoothed down to the lowest resolution. This is necessary to avoid adding information from other higher resolution, higher frequencies maps that are smaller than the beam size of some of the maps into the recovered map. The smoothing of the maps to the lowest map resolution was done using HEALPy, and this process can be formalised as

(22)

where are the a_ℓ, m coefficient of the map at a frequency i smoothed down to the lowest resolution. The Gaussian kernel at this desired ‘lowest’ resolution is indicated at b^f. The consequence of the necessary smoothing to the lowest resolution is an increase of the instrumental white noise contribution in each frequency map by a factor b_f(ℓ)/b_i(ℓ) because for all the frequency channels, except the lowest one, we have b_f(ℓ) > b_i(ℓ).

4.4. Full-sky treatment

The microwave full-sky maps are composed of many different astrophysical signals. Some of those signals are spatially localised. The Galactic foregrounds are mainly localised around the Galactic plane. When applying the map-based ILC to a full-sky map, instead of applying it straightforwardly to the multi-frequency set of maps, thus getting only one set of weights for the full sky and equally treating regions that are strongly contaminated by Galactic foregrounds and regions that are less contaminated, we tessellated the sky into different zones and applied an ILC on each of these zones separately. This resulted in a set of different weights for each region so that they are always optimised for the level of contamination present in their zone. This has been done previously in various studies (Tegmark et al. 2003). The most common tessellation approach is to separate the sky into zones of Galactic contamination.

To tessellate the sky, we used the native nested HEALPix tessellation (see Fig. 4). This nested scheme tessellates the sky originally into 12 equal-area fields. The sides of each of those fields can be further divided to create more equal-area cells of smaller sizes. Using HEALPix nested scheme to create sub-fields (Górski et al. 2005), we chose to tessellate the sky using N_side = 4, which gives fields, each covering an area of ≈215 deg². In general, tessellation can create border effects along the limits of the fields. In Eriksen et al. (2004), the border effects were attenuated by smoothing the ILC weights. In this work, we do not correct the border effects.

The Galactic plane was masked (see Fig. 4) with a f_sky ≈ 0.6 Galactic mask map to reduce foreground contamination and better extract the tSZ effect. To compute power spectra on a masked map, we used PyMaster (Alonso et al. 2019), which corrects the power spectra for the masking of part of the sky. In this work, we do not apodise the masks, as it is computationally very intensive at our high resolution and does not make a significant difference in the resulting power spectrum (see Fig. A.2). The full treatment which includes, creation, tessellation, masking and ILC is performed at the map level. Only when a tSZ full-sky map was extracted was the power spectrum of the map computed using spherical harmonics.

4.5. Debiasing

The instrumental white noise biasing of the tSZ power spectrum extraction was estimated for a 4000 h observation time (Choi et al. 2020b). Therefore, the mean of this noise could be computed and removed to debias the tSZ reconstruction. This was done by generating random realisation of the instrumental white noise for each frequency at the sensitivities given in Table 2 and computing the projected noise by multiplying the ILC weights from the tSZ reconstruction with the noise maps and summing them. The power spectrum of this noise was then subtracted from the ILC-extracted Compton-y (y_ILC) power spectrum. All the following power spectrum results were debiased for instrumental white noise.

Similarly, as the CMB power spectrum is well-known, we could generate random realisation from the Websky CMB power spectrum to debias for the CMB residual noise. In the first approximation, this study presents an idealised case where the CMB power spectrum is perfectly known and its residual is removed. With real data, this CMB debiasing will be limited by the uncertainties on the recovered CMB power spectrum.

4.6. Overview of steps

In this section, we describe our tSZ extraction method and illustrate it in Fig. 5.

Fig. 5. Scheme of the method we use in this paper to extract the tSZ signal. The input of the sky model is template maps coming from the Websky simulations and PySM. The sky model processes them to produce maps of the different sky components at a given frequency, sensitivity, and beam. This set of multi-frequency maps {Mi1, …, Min} was smoothed to a common resolution and tessellated, and the Compton-y signal was extracted by the ILC. The PyMaster algorithm was used on the resulting masked map to extract the power spectrum.

First, using the Skymodel based on the Websky simulation and PySM, we generated multi-frequency maps of the microwave sky at the frequencies, sensitivities, and beam sizes of SO LAT and FYST (see Table 2). The set of maps is in the HEALPix format and uses intensity MJy sr⁻¹ units.

Second, the simulated maps were all smoothed down to the resolution of the lowest frequency channel used in our analysis.

Third, a Galactic mask that excludes 40% of the sky (see Fig. 4) was applied to the set of multi-frequency maps.

Fourth, each map was tessellated into smaller equal-area regions using the HEALPix nested format.

Fifth, an ILC was applied to each multi-frequency set of the tessellated area to obtain the Compton-y patches. The patches were reassembled to obtain the full-sky Compton-y map. The covariance matrix was computed only on the unmasked data by our ILC.

Sixth, the power spectrum of the recovered Compton-y map was computed using PyMaster and debiasing, and the residual noises were computed using Eq. (17).

5. Results

To quantify how much CIB could leak into the reconstructed tSZ maps for future experiments and access the detectability of the tSZ power spectrum, we used the microwave sky model detailed in Sect. 3. We simulated a set of multi-frequency maps at the frequencies, sensitivities, and beam sizes for SO LAT and FYST (see Table 2) with atmospheric noise. We applied an ILC or CILC component separation method to extract a tSZ map and looked at the ratios between the different noises and signal power spectra. We restricted ourselves to the scales ℓ ∈ [500, 5000]. We chose this lower limit to avoid the cosmic variance and the upper limit because of the 2.2′ resolution of our HEALPix maps, which allowed us to resolve scales up to ℓ ∼ 4909.

5.1. Internal linear combination results

By construction, the weights of an ILC must verify the first condition, which is to have a unit response to the SED of the tSZ effect (see Eq. (13)). This first condition means that the element-wise product of weights (ω_i) times the mixing vector (a_i) is an indicator of the fractional contribution of each frequency band to the recovery of the tSZ signal. It can be negative since the frequency spectrum of the tSZ effect is negative below 217 GHz. The weights attributed by the ILC can themselves be negative for optimisation purposes, with the only constraint being that the sum over all frequencies must be one. The value of the weights is given by the mean of the 192 weights (one per tessellated field), and the error bars are given by the standard deviation of those weights. The top panel of Fig. 6 shows that for SO LAT alone (orange dots), the 93 GHz channel is the main contributor to the tSZ reconstruction by a factor ω₉₃a₉₃ ≈ 1.4. This channel overestimates the tSZ signal, and the SO LAT 145 GHz channel was used by the ILC to correct this. As a result, the 145 GHz channel contributes roughly by a factor ω₁₄₅a₁₄₅ ≈ −0.4 to the tSZ reconstruction. This result is not surprising, as those frequency bands are located at the tSZ decrement that is characteristic of the signal in the microwave sky and where the values of the mixing vector are negative. When SO LAT and FYST are combined (SO LAT+FYST; blue dots), the SO LAT 93 GHz and 145 GHz channels remain the main contributors to the tSZ reconstruction and have similar values, ω₉₃a₉₃ ≈ 1.45 and ω₁₄₅a₁₄₅ ≈ −0.45. Other channels from SO LAT and FYST do not contribute to the reconstruction. When FYST is alone (green dots), the main contributor to the tSZ reconstruction is the 280 GHz channel, by a factor ω₂₈₀a₂₈₀ ≈ 1.1, and the 350 GHz channel, by ω₃₅₀a₃₅₀ ≈ 0.05. The 220 GHz and 405 GHz channels were used to mitigate the over evaluation of the tSZ signal with the respective factors ω₂₂₀a₂₂₀ ≈ −0.05 and ω₄₀₅a₄₀₅ ≈ −0.1.

Fig. 6. Internal linear combination weights for SO LAT-like (orange) and FYST-like (green) experiments and SO LAT+FYST (blue) when the sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top figure: fractional contribution of each frequency channel to the recovered y-map. Bottom figure: frequency channels that contribute to maximising the tSZ signal and minimising the noise. The error bars show the standard deviation around the mean weight of all of the 192 tessellated fields for each frequency map.

The second condition that has to be verified by the ILC weights is to ensure that the variance of the recovered map has been minimised (see Eq. (14)). The bottom panel of Fig. 6 shows the values of the weights. In combination with the top panel, we concluded that the weights associated with some frequency bands that do not contribute to the tSZ signal reconstruction but have a high absolute value contribute to the noise minimisation. In particular, the 220 GHz channel located close to the tSZ null is often attributed a high weight by the ILC even though it does not contribute to the tSZ reconstruction. This is because it is the only map containing all the noise contaminants and almost no signal. It is thus the perfect map for noise minimisation. But this is not the case in Fig. 6, except for when FYST is alone (green dots), because of the presence of atmospheric noise in our mock sky. The strong atmospheric noise contaminates the 220 GHz band, making it less exploitable by the ILC (see the comparison with Fig. A.1 for a case without atmospheric noise). Figure 6 shows that in the presence of all the microwave emissions and atmospheric noise, the SO LAT frequency channels that probe the tSZ decrement, 93 GHz and 145 GHz, drive both the tSZ reconstruction and the noise minimisation of the tSZ signal.

However, the weights do not tell us anything about the CIB contamination or the detectability of the tSZ signal. For that, we turned to the power spectrum and the residual ILC noise power spectra. The ILC output is a map containing the Compton-y signal and residual noises (Eq. (17)). Using PyMaster, the power spectra of this map and the residual noises can be computed and compared to the input-pure Compton-y map power spectra from the Websky simulations. We defined as the power spectra of the quantity X or of the residual noise r_X of the quantity X. For SO LAT alone (see Fig. 7 top-left panel) the dominant noise residual on small scales is the CIB residual (r_CIB), which starts to dominate around ℓ ∼ 3700. The cumulative foreground residual (r_f) and the kSZ residual (r_kSZ) are both subdominants to the input signal by around one order of magnitude over the scales ℓ ∈ [500, 5000]. The atmospheric noise residual (r_Atmo) dominates the input signal for ℓ < 1400 and ℓ > 4500.

Fig. 7. Power spectra of the ILC residual noises compared to the input Compton-y power spectrum from the Websky simulations (red). The simulated sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top panels: power spectra of various quantities. The grey curves show the ILC noises residual power spectra, denoted by rX (X being the kSZ); the cumulative Galactic foregrounds (f); and the atmospheric noise (Atmo). The black curve shows the power spectrum of the CIB residual noise. The blue curve shows the sum of all grey and black curves, that is, the total ILC noise residual debiased for instrumental white noise and the CMB. The top-left panel shows the results for the simulated SO LAT data, while the top-right panel shows the results for the simulated data for SO LAT and FYST combined. Bottom panel: ratio of the input Compton-y power spectrum (red) and the total noise residual (blue) for SO LAT and SO LAT+FYST. All power spectra have been beam corrected and bin averaged over a window Δℓ = 50.

In the combined case of SO LAT+FYST (see right panel of Fig. 7), the situation is similar to the case where SO LAT is alone except that the atmospheric noise residual (r_Atmo) is higher when FYST is present. This could be due to higher atmospheric noise being present in the high-frequency channels of FYST (see Fig. 3). The CIB noise residual power spectrum is significantly lower for SO LAT+FYST. It starts to dominate above the input signal around ℓ ∼ 4200, while for SO LAT alone it was around ℓ ∼ 3700. We concluded that adding FYST to SO LAT reduces the CIB noise residual.

The lower panel of Fig. 7 shows the extent the input Compton-y power spectrum is dominated by the total ILC noise residual as a function of the scale. The total noise residual (r_tot) power spectrum dominates over the input signal on nearly all scales except around ℓ ∈ [2100, 2700] for SO LAT alone. This is because of our debiasing of the CMB on those critical scales. Overall, SO LAT alone has a better input/noise ratio on scales ℓ < 2500 because of its slightly lower atmospheric noise residual, and SO LAT+FYST does better than SO LAT alone on scales ℓ > 3000 because of its lower CIB noise residual.

5.2. Constrained internal linear combination results

Constrained internal linear combination can be used to ‘null’ the contribution of a component of the microwave sky using its frequency spectrum. Nullifying several components was tested.

Before any debiasing, the CMB is the main contaminant of the recovered Compton-y power spectra on large scales, so it seemed natural to deproject (i.e. ‘null’) it. However, the additional constraints on the ILC weights (see Eq. (18)) led to a smaller solution space for the weights to extract the signal, which resulted in a ‘noise penalty’, that is, higher noise in the recovered signal map. In this particular case of a mock sky containing all Galactic foregrounds and extragalactic emission, excluding RPSs, and applying a CILC to deproject CMB, the ILC weights are so constrained to null CMB that they become less efficient at reducing the second main contaminant, the atmospheric noise contamination, which blows up. As a consequence, the total residual noise dominates the input signal even more, by up to one order of magnitude on small scales (see Appendix B). Deprojecting the CMB is therefore not a recommended option for our study. This effect was already reported in Fig. 36 of Ade et al. (2019), but in their study, the blowing up of the noise was mitigated on large scales by the addition of Planck data.

Because of its high-frequency coverage, we expected FYST to better probe the CIB, which is the dominant contaminant above ν > 300 GHz, and therefore better target CIB contamination and remove it. Deprojecting the CIB using CILC seemed the natural option to take full advantage of the potential constraining power of the FYST high-frequency channels. We note that a perfect deprojection of the CIB is not possible because its signal is not composed of only one SED. We adopted the modified blackbody spectrum given in Eq. (8), with A_CIB = 1, T_dust = 10 K, and β = 1.6. Parameter values were calculated by computing the average for each of the parameter maps derived in Sect. 3.2.

The top panel of Fig. 8 shows that the CIB noise residual (r_CIB) power spectrum is reduced by the CILC deprojection for both SO LAT and SO LAT+FYST, but it is not zero because the CIB is composed of multiples spectra and only one is used for the deprojection. It also shows that the atmospheric noise residual (r_Atmo) power spectrum increases for SO LAT and decreases for the combined SO LAT+FYST when the CIB is deprojected (see Fig. 7 for comparison). This can be explained by the higher degrees of freedom available with the higher number of frequency channels for SO LAT+FYST, making the CILC noise penalty have a less penalising effect in comparison to when SO LAT is alone. The lower panel of the figure shows that SO LAT+FYST offers a large ℓ ∈ [1700, 3500] window over which the Compton-y signal is not noise dominated, while for SO LAT alone, the signal remains noise dominated (exclusively by the atmospheric noise). For SO LAT alone or for SO LAT+FYST combined, compared to the ILC case, the ratio between the input target signal and the resulting noise is not much higher in general (Fig. 7). This was also seen when applying our pipeline to Planck data, as the recovered Compton-y map amplitude is roughly the same as the ILC and CILC. This is not surprising because we also saw that the ILC and CILC weights are very similar. This observation of the CILC not providing a much lower noise bias Compton-y map compared to the ILC was already reported by Alonso et al. (2018).

Fig. 8. Power spectra of the CILC residual noises when deprojecting CIB compared to the input Compton-y power spectrum from the Websky simulations (red). The simulated sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top panels: power spectra of various quantities. The grey curves show the ILC noise residual power spectra, denoted by rX (X being the kSZ); the Galactic foregrounds (f); and the atmospheric noise (Atmo). The black curve shows the power spectrum of the CIB residual noise. The blue curve shows the sum of all grey and black curves, that is, the total ILC noise residual debiased for instrumental white noise and the CMB. The top-left panel shows the results for the simulated SO LAT data, while the top-right panel shows the results for the simulated data for SO LAT and FYST combined. Bottom panel: ratio between the input Compton-y power spectrum (red) and the total noise residual (blue) for SO LAT and SO LAT+FYST. All power spectra have been beam corrected and bin averaged over a window Δℓ = 50.

5.3. Overview

Figure 9 gives an overview of the relative amplitude changes for the power spectrum of residual noises of each astrophysical component when combining SO LAT and FYST and shows a comparison to SO LAT alone in the two previous cases applying an ILC (in blue) and a CILC deprojecting the CIB (in orange). To see how much SO LAT+FYST reduces residual noise power spectra in comparison to SO LAT alone, we took the ratio of the residual power spectra of each component in both cases:

(23)

Fig. 9. Value of the ratios of the power spectra of X for SO LAT+FYST compared to SO LAT(). The ILC-extracted Compton-y signal (yILC) and rX represent the residual noise of the quantity X. The shaded regions are the density distribution of the ratios over the ℓ range. The bar is the mean averaged over the ℓ range of the ratio distribution. In blue is the ILC case (see Sect. 5.1) and in orange is the CILC case when CIB is deprojected (see Sect. 5.2).

where is the power spectrum of the component X or of the residual noise r_X for SO LAT+FYST and is the same quantity but for SO LAT alone. For the extracted Compton-y map (y_ILC), the blue violin plot shows a mean improvement of around 4% for SO LAT+FYST compared to SO LAT alone but with a large scatter. Especially at ℓ < 2800, SO LAT performs better (see bottom panel of Fig. 7). When deprojecting the CIB (orange violin plot), we moved to a mean improvement of around 8% in the ILC-extracted Compton-y map for SO LAT+FYST compared to SO LAT alone. In this case, the scatter is also much smaller than before because of the high number of degrees of freedom offered by the combined SO LAT+FYST frequency bands, and the CILC noise penalty is much lower for SO LAT+FYST compared to SO LAT alone, making the combination strictly superior to extract a less biased tSZ power spectrum in the CILC case.

For the residual CIB noise (r_CIB), the power spectrum is on average around 35% lower when combining SO LAT and FYST in the ILC case. In the CILC case, it drops to ∼7%, which might not be surprising because the CILC focusses on nulling as much CIB as possible no matter the frequencies. The addition of FYST does reduce the CIB residual even more but not as much as without the null constraint. Another reason for the drop in improvement when combining SO LAT and FYST in the CILC case might be the higher atmospheric noise at the high frequencies where the CIB signal dominates, undermining the advantage of FYST to constrain and remove CIB.

In the ILC case, the CMB residual noise (r_CMB) power spectrum is reduced by ∼22% when combining SO LAT with FYST compared to SO LAT alone. Because the ILC cannot separate CMB and kSZ due to their identical spectral frequency distribution, the residual kSZ noise (r_kSZ) is also suppressed when combining SO LAT and FYST, on average by around 18%. When deprojecting CIB with a CILC, adding FYST to SO LAT also reduces the CMB and kSZ residual noises but only by ∼5% over all scales on average.

The sum of all galactic foreground noise residuals is defined as the cumulative foreground residual (r_f). Its noise power spectrum is improved on average by ∼10% for SO LAT+FYST compared to SO LAT in the ILC case. In the CILC case, the cumulative foreground residual power spectrum is reduced by ∼30% when combining FYST with SO LAT. This could be due to the additional constraint on the CIB dust spectrum, nulling part of the Galactic dust contribution, which is more easily nulled with FYST high-frequency coverage.

The FYST high-frequency channels are more contaminated by the atmosphere. For this reason, the atmospheric red noise residual (r_Atmo) power spectrum is ∼30% higher when SO LAT and FYST are combined rather than when SO LAT is alone in the ILC case. When deprojecting the CIB, the atmospheric noise residual becomes slightly lower ∼1% for SO LAT+FYST compared to SO LAT alone. This is due to the noise penalty generated by the additional constraint in the CILC worsening the atmospheric noise for LAT alone compared to SO LAT+FYST. Because SO LAT+FYST benefits from additional frequency bands and therefore more degrees of freedom, it is more resilient to the inevitable noise increase caused by additional spectral constraints.

6. Summary and discussion

Using full-sky HEALPix template maps from the Websky simulations (Stein et al. 2020) and PySM, we created a pipeline that allowed us to simulate the high-resolution maps of the different microwave sky emissions at any beam, frequency, and sensitivity. Based on the Choi et al. (2020a) noise modelling, our pipeline also simulates atmospheric red noise at the locations of the Simons Observatory and FYST.

We simulated the sky that would be seen during experiments at these locations and used a map-based ILC and CILC to extract back the thermal Sunyaev–Zeldovich (tSZ) effect from this multi-frequency, multi-component sky. To minimise Galactic foreground contamination, we used a Galactic mask derived from Planck (see Fig. 4). To optimise our ILC weights in order to take advantage of the spatial distribution of some of the microwave emissions, we tessellated the sky using the HEALPix equal-area nested tessellation scheme. We made use of CILC to deproject the CIB by using a simple modified blackbody spectrum to represent the CIB frequency variation and constrain the ILC weights to nullify it. We then used PyMaster to compute the power spectrum and its associated residual noise power spectra, which are the noise coming from all the other microwave sky emissions left by ILC and CILC in the recovered map.

We note that our study does have certain shortcomings that should be taken into account. For example, the Websky simulations do not yet offer a map of RPSs, and these extragalactic components are therefore absent from our study. The RPSs are an important contaminant of the tSZ signal, as they are located in galaxy clusters. If RPSs are not considered, they can fill up the tSZ decrement and bias an estimation towards a lower estimated tSZ signal. A solution to this problem is to mask the RPSs, but this also lowers the estimated tSZ power spectrum, as it removes parts of the tSZ signal. This effect especially impacts measurements of the tSZ signal below 90 GHz, which we do not study here. A study from Holder (2002) showed that RPS subtraction can lead to an underestimation of the tSZ power spectrum of between ∼30% and ∼50%, depending on the scale. We did not consider the relativistic correction of the tSZ effect, which can lead to a significant underestimation of the signal (Remazeilles et al. 2018). This correction is left to future studies. Regarding the method, our tessellation scheme can lead to border effects that were reported by Eriksen et al. (2004) but not accounted for. Our pipeline was tested on Planck data with and without tessellation and only a marginal (i.e. < 5%) difference could be seen in the power spectrum amplitude. Applying ILC in map-based space also has the disadvantage of the necessary resolution downgrade of all maps to the lower beam channel, which is 2.2′ for SO LAT, thus reducing the benefits of FYST’s higher resolution channels. This also prevented us from using Planck data that are essential to lowering atmospheric and CMB contamination on large scales. This was addressed by performing CMB debiasing of the residual CMB noise in the extracted Compton-y map. A Fourier space or needlet ILC implementation is left for future work.

We find that when using an ILC to extract tSZ, the combination of SO LAT with FYST reduces the CIB residual noise power spectrum by ∼35%, on average, on scales ℓ ∈ [500, 5000] when compared to SO LAT alone. For the CMB and kSZ components, the addition of FYST to SO LAT reduces the residual noise by ∼22% on average, with a strong scale dependence for the CMB, as the addition of FYST mainly helped on large scales. For the kSZ component, the reduction observed when adding FYST to SO LAT is not scale dependent. Adding FYST also helps reduce the Galactic foreground contamination by, on average, ∼10% on scales ℓ ∈ [500, 5000]. However, FYST’s high-frequency channels are more affected by atmospheric noise due to poorer transmission. Therefore, the residual noise power spectrum coming from the atmosphere is on average 30% higher when SO LAT and FYST are combined compared to when SO LAT is alone. This greatly affects the total residual noise budget, lowering the benefits of adding FYST to SO LAT for the tSZ recovery to only an average ∼4%, and on larger scales, SO LAT alone performs better. We find that over ℓ ∈ [2100, 2700], the ILC-extracted tSZ power spectrum is not noise dominated for SO LAT alone. But on smaller scales, it is extremely biased, and its noise budget is dominated by the CIB. When SO LAT is combined with FYST, the CIB noise residual dominates much later ℓ ∼ 4200, but because the atmospheric noise is higher, there is no window over which the ILC-extracted tSZ power spectrum is not noise dominated.

We note that the CIB residual noise reduction that FYST brings to SO LAT is not something that SO LAT could achieve on its own by simply having a longer survey time. This was tested by reducing the instrumental white noise in SO LAT bands by two, and we observed a reduction of the ILC residual instrumental noise in the reconstructed Compton-y map, but the other noise components remained the same. This ‘boosted’ version of SO LAT alone did not have a lower CIB residual power spectrum compared to the planned SO LAT or SO LAT plus FYST combination.

Because tSZ and CIB have a non-zero correlation in the Websky simulations and because FYST’s high frequencies probe the CIB better, we used CILC to deproject the CIB component, modelling it with a simple one-component modified blackbody spectrum. Compared to the ILC case, adding FYST to SO LAT lowers the CIB noise residual by only ∼7%. The CMB and kSZ noise residuals are also reduced by ∼5% when combining SO LAT and FYST when using CILC. However, the cumulative foregrounds are suppressed by ∼30%, and atmospheric noise is a bit lower for the SO LAT and FYST combination. Despite a lower reduction of the CIB residual noise when combining SO LAT and FYST, the CIB deprojection case offers a large ℓ ∈ [1800, 3500] window where the CILC-extracted Compton-y map is not noise dominated after instrumental noise and CMB debiasing.

Acknowledgments

The authors would like to thank Kaustuv Basu for his guidance and advice throughout the project. Srinivasan Raghunathan, Jacques Delabrouille, and Frank Bertoldi for their time, valuable insights, and numerous helpful discussions. The referee for useful comments and numerous improvements to the draft. The editor Laura Pentericci for her advice and help on the submission process. Michael Coronado for his language editing work which significantly helped improve the manuscript. The CCAT-prime Collaboration for their feedback on this work and Colin Hill for his help. Laila Linke for her useful comments on this draft. We furthermore acknowledge support from the International Max-Planck Research School (IMPRS) and the Bonn-Cologne Graduate School of Physics and Astronomy (BCGS). The simulations of the CIB used in this paper were developed by the Websky Extragalactic CMB Mocks team, with the continuous support of the Canadian Institute for Theoretical Astrophysics (CITA), the Canadian Institute for Advanced Research (CIFAR), and the Natural Sciences and Engineering Council of Canada (NSERC), and were generated on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund – Research Excellence, and the University of Toronto.

References

Appendix A: Planck data

A.1. Internal linear combination and constrained internal linear combination weights

We applied our pipeline to Planck High-Frequency instrument (HFI) Data Release 2³ (DR2; Planck Collaboration XXII 2016) in order to see how well we could recover the true y using an ILC. We also considered a CILC to deproject CIB by constraining the weights so that they have a null response to the frequency spectrum of the CIB (see Eq. (20)). However, the CIB deprojection is not perfect, because the CIB frequency spectrum depends on the population of the sources and their redshift.

Looking at the ILC weights gave us information about which channel contributes the most to extracting the tSZ signal and minimising the noise. These were the two conditions that the weights had to satisfy (see Eq. (13) and (14)). The ILC weights ponderated by their respective frequency mixing vector value showed the fractional contribution of each frequency channel to the recovered y. The fractional contribution of the ∼220 GHz band was null because the value of the mixing vector at this frequency is null.

For a map-based ILC applied to Planck data (see blue points in Fig. A.1), we concluded from the weight that the 143 GHz channel is the one that contributes the most, with a factor of ω₁₄₃a₁₄₃ ≈ 0.7, to the reconstructed final map as well as to minimising the noise. The 217 GHz channel does not contribute to the final reconstructed y map because the tSZ signal is null at this frequency. However, its high ω₂₁₇ value makes it the main contributor to noise minimisation. This was expected, as the tSZ signal is absent at this frequency and the map contains all the noise contaminants. Next are the contributions of the 353 GHz channel and the 545 GHz channel. The contribution values of those channels are consistent with previous studies (Erler et al. 2018). In the case where a CILC is applied to Planck data to deproject the CIB (orange points, Fig. A.1), we saw that the weights are almost identical to the ILC ones. This is expected, as Alonso et al. (2018) also found that deprojecting CIB with a CILC does not result in much less CIB contamination in the reconstructed tSZ profile, arguing that this shows the ILC already does a sufficient job of cleaning the CIB.

Fig. A.1. Internal linear combination (blue) and CILC (orange) weights with CIB deprojected on Planck data. Top figure: Fractional contribution of each frequency channel to the recovered y-map. Bottom figure: Frequency channels that contribute to maximising the tSZ signal and minimising the noise. The error bars show the standard deviation around the mean weight of all of the 192 tessellated fields for each frequency map.

A.2. Processing of the maps

We tested the consistency and effect of each of the post-processing steps applied on the map (see Fig. A.2). The blue curve in the figure shows the ILC-extracted power spectrum computed with a Galactic mask and without tessellation (see Fig. 4), while the lighter blue dotted curve shows the same but with apodisation. The black curve shows the case with the tessellation of the sky and the ILC being applied to each rhombus separately. The grey dotted line shows when apodisation was also applied. For our particular case of tSZ extraction, the mask apodisation only lowers the power spectra of the maps by less than 5%. Comparing the power spectrum with tessellation (black) and without (blue) showed that tessellation helps reduce the noise on large scales, up to ℓ ≈ 2000, by a factor of approximately two. At ℓ = 2500, it increases the noise by around 15%. This can be understood as the size of the rhombus being small enough to break down regions highly contaminated by unwanted astrophysical components on large scales and to weigh down such regions in order to better recover the tSZ signal. But the rhombus size might be too big to optimally deal with smaller region noise. Moreover, with 192 total rhombuses over the sky, we might be creating border effects on smaller scales and thus increasing the noise residual on those scales. The scales (ℓ > 2200) are also affected by the ∼9.66’ beam.

Fig. A.2. Power spectra of the Compton-y map using our ILC pipeline on Planck DR2 with different configurations. In blue is the power spectrum of the extracted map with a simple fsky = 0.6 Galactic mask. In black is the same thing but with a tessellation of the sky. The dotted lines (blue and black) indicate the inclusion of apodisation in each case.

Appendix B: Constrained ILC: Deprojecting the cosmic microwave background

When applying an ILC to extract back a Compton-y map, the CMB is the dominant noise residual (r_CMB) on large scales (see Fig. B.1). Using CILC to deproject the CMB, the residual CMB noise is suppressed by three orders of magnitude. However, the additional constraints imposed on the weights to deproject the CMB led to a smaller solution space and therefore increased noise. This led to a one order of magnitude increase of the atmospheric residual noise (r_Atmo), which worsens the estimate of the extracted Compton-y power spectrum, especially on small scales, by more than one order of magnitude. This is commented on extensively in Sec. 5.2 and was tested using our Skymodel pipeline (see Sec. 3).

Fig. B.1. Power spectra of ILC/CILC (full lines) and their noise residuals (dotted lines). Specifically, the orange curve is the power spectrum of the ILC-extracted Compton-y map. The lighter orange lines are the noise residuals. Specifically, they are the ILC residual CMB noise (rCMB) and rAtmo for the atmospheric noise. The blue lines represent the case with CILC where the CMB is deprojected. The solid, darker blue curve is the power spectrum of the CILC-extracted Compton-y map. The lighter blue curves are the noise residuals.

Appendix C: ILC weights

In our study, we tesselated the sky into 192 rhombi of equal area using the HEALPix nested scheme. This sky separation allowed us to better adapt to the spatially dependant contaminant of some Astrophysical emissions, such as the Galactic foregrounds. The ILC we applied on each zone resulted in one weight per frequency and per rhombus. An example of a full sky weights map produced and used by our ILC can be seen in Fig. C.1.

Fig. C.1. Maps of the weights for each patch of the sky where the ILC was applied for the 860 GHz frequency band. The weights were computed from simulations that included all Galactic foregrounds and extragalactic components (see Table 1) at the sensitivities, beam sizes, and frequencies of SO LAT and FYST (see Table 2) and with 90% correlated atmospheric red noise generated using the noise curves given by Choi et al. (2020a), (see Fig. 3).

Appendix D: Simulated maps

In this work, we used the Skymodel pipeline based on PySM and the Websky simulations to generate mock maps of the microwave sky at the frequencies, sensitivities, and instrument beam size of the upcoming SO LAT and FYST (see Table 2). The atmospheric red noise was generated using the noise curves of Choi et al. (2020a) and the red noise of the neighbouring frequency channels of a given telescope were correlated at 90% using Eq. 10. The resulting simulated microwave sky maps used can be seen in Fig. D.1.

Fig. D.1. Mock HEALPy maps of the microwave sky simulated using the Skymodel pipeline (see Sec. 3). All Galactic foregrounds and extragalactic components (see Table 1) at the sensitivities, beam sizes, and frequencies of SO LAT and FYST (see Table 2) and the 90% correlated atmospheric red noise between neighbouring channels generated using the noise curves given by Choi et al. (2020a) were included (see Fig. 3).

All Tables

All Figures

Fig. 1. Projections of the PySM sky at 353 GHz composed of free-free, synchrotron, AME, and Galactic dust. The upper row shows cutouts around the Galactic south pole, and the lower row shows cutouts around the Galactic centre. Each cutout is composed of 600 × 600 pixels of the size 0.86′, making 8.6° ×8.6° maps. The left column shows cutouts from the original PySM sky, which feature pixelation artefacts due to their pixel size of 6.9′ (Nside = 512). The right column shows modified cutouts obtained by first updating the original map to a higher resolution, 0.86′ (Nside = 4096), by oversampling using the HEALPy function ‘ud_grade()’ and then smoothing with a Gaussian kernel of a width that is close the original pixel size of the PySM maps (FWHM = 10′). This reprocessing of the PySM maps removes the pixelation artefacts while only marginally lowering the spatial resolution.

In the text

	Fig. 2. Correlation matrix between the Websky CIB maps provided at the Planck HFI frequencies. The axis shows the frequencies of the maps in GHz, and the dimensionless colour bar shows the correlation coefficient.
In the text

	Fig. 3. SO LAT (blue lines) and FYST (red lines) sum of the atmospheric red and instrumental white noise spectra based on the noise model presented by Ade et al. (2019) and adopted by Choi et al. (2020a). The noise curves have been beam corrected.
In the text

	Fig. 4. Galactic foreground mask derived from Planck data by Erler et al. (2018; fsky ≈ 0.6). The white part is the masked region, and the blue part is what was used in our analysis. Overplotted is the HEALPix tessellation scheme for Nside = 4.
In the text

Fig. 5. Scheme of the method we use in this paper to extract the tSZ signal. The input of the sky model is template maps coming from the Websky simulations and PySM. The sky model processes them to produce maps of the different sky components at a given frequency, sensitivity, and beam. This set of multi-frequency maps {Mi1, …, Min} was smoothed to a common resolution and tessellated, and the Compton-y signal was extracted by the ILC. The PyMaster algorithm was used on the resulting masked map to extract the power spectrum.

In the text

Fig. 6. Internal linear combination weights for SO LAT-like (orange) and FYST-like (green) experiments and SO LAT+FYST (blue) when the sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top figure: fractional contribution of each frequency channel to the recovered y-map. Bottom figure: frequency channels that contribute to maximising the tSZ signal and minimising the noise. The error bars show the standard deviation around the mean weight of all of the 192 tessellated fields for each frequency map.

In the text

Fig. 7. Power spectra of the ILC residual noises compared to the input Compton-y power spectrum from the Websky simulations (red). The simulated sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top panels: power spectra of various quantities. The grey curves show the ILC noises residual power spectra, denoted by rX (X being the kSZ); the cumulative Galactic foregrounds (f); and the atmospheric noise (Atmo). The black curve shows the power spectrum of the CIB residual noise. The blue curve shows the sum of all grey and black curves, that is, the total ILC noise residual debiased for instrumental white noise and the CMB. The top-left panel shows the results for the simulated SO LAT data, while the top-right panel shows the results for the simulated data for SO LAT and FYST combined. Bottom panel: ratio of the input Compton-y power spectrum (red) and the total noise residual (blue) for SO LAT and SO LAT+FYST. All power spectra have been beam corrected and bin averaged over a window Δℓ = 50.

In the text

Fig. 8. Power spectra of the CILC residual noises when deprojecting CIB compared to the input Compton-y power spectrum from the Websky simulations (red). The simulated sky contains all extragalactic components (see Table 1), Galactic foregrounds, and atmospheric noise. Top panels: power spectra of various quantities. The grey curves show the ILC noise residual power spectra, denoted by rX (X being the kSZ); the Galactic foregrounds (f); and the atmospheric noise (Atmo). The black curve shows the power spectrum of the CIB residual noise. The blue curve shows the sum of all grey and black curves, that is, the total ILC noise residual debiased for instrumental white noise and the CMB. The top-left panel shows the results for the simulated SO LAT data, while the top-right panel shows the results for the simulated data for SO LAT and FYST combined. Bottom panel: ratio between the input Compton-y power spectrum (red) and the total noise residual (blue) for SO LAT and SO LAT+FYST. All power spectra have been beam corrected and bin averaged over a window Δℓ = 50.

In the text

Fig. 9. Value of the ratios of the power spectra of X for SO LAT+FYST compared to SO LAT(). The ILC-extracted Compton-y signal (yILC) and rX represent the residual noise of the quantity X. The shaded regions are the density distribution of the ratios over the ℓ range. The bar is the mean averaged over the ℓ range of the ratio distribution. In blue is the ILC case (see Sect. 5.1) and in orange is the CILC case when CIB is deprojected (see Sect. 5.2).

In the text

Fig. A.1. Internal linear combination (blue) and CILC (orange) weights with CIB deprojected on Planck data. Top figure: Fractional contribution of each frequency channel to the recovered y-map. Bottom figure: Frequency channels that contribute to maximising the tSZ signal and minimising the noise. The error bars show the standard deviation around the mean weight of all of the 192 tessellated fields for each frequency map.

In the text

	Fig. A.2. Power spectra of the Compton-y map using our ILC pipeline on Planck DR2 with different configurations. In blue is the power spectrum of the extracted map with a simple fsky = 0.6 Galactic mask. In black is the same thing but with a tessellation of the sky. The dotted lines (blue and black) indicate the inclusion of apodisation in each case.
In the text

Fig. B.1. Power spectra of ILC/CILC (full lines) and their noise residuals (dotted lines). Specifically, the orange curve is the power spectrum of the ILC-extracted Compton-y map. The lighter orange lines are the noise residuals. Specifically, they are the ILC residual CMB noise (rCMB) and rAtmo for the atmospheric noise. The blue lines represent the case with CILC where the CMB is deprojected. The solid, darker blue curve is the power spectrum of the CILC-extracted Compton-y map. The lighter blue curves are the noise residuals.

In the text

Fig. C.1. Maps of the weights for each patch of the sky where the ILC was applied for the 860 GHz frequency band. The weights were computed from simulations that included all Galactic foregrounds and extragalactic components (see Table 1) at the sensitivities, beam sizes, and frequencies of SO LAT and FYST (see Table 2) and with 90% correlated atmospheric red noise generated using the noise curves given by Choi et al. (2020a), (see Fig. 3).

In the text

Fig. D.1. Mock HEALPy maps of the microwave sky simulated using the Skymodel pipeline (see Sec. 3). All Galactic foregrounds and extragalactic components (see Table 1) at the sensitivities, beam sizes, and frequencies of SO LAT and FYST (see Table 2) and the 90% correlated atmospheric red noise between neighbouring channels generated using the noise curves given by Choi et al. (2020a) were included (see Fig. 3).

In the text