# A Cross-Internal Linear Combination Approach to Probe the Secondary CMB
Anisotropies: Kinematic Sunyaev-Zel’dovich Effect and CMB Lensing

Srinivasan Raghunathan Center for AstroPhysical Surveys, National Center for
Supercomputing Applications, Urbana, IL 61801, USA  Yuuki Omori Department of
Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA
Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA Srinivasan Raghunathan<EMAIL_ADDRESS>Yuuki Omori
<EMAIL_ADDRESS>
###### Abstract

We propose a cross-internal linear combination (cross-ILC) approach to measure
the small-scale cosmic microwave background (CMB) anisotropies robustly
against the contamination from astrophysical signals. In particular, we focus
on the mitigation of systematics from cosmic infrared background (CIB) and
thermal Sunyaev-Zel’dovich (tSZ) signals in kinematic Sunyaev-Zel’dovich (kSZ)
power spectrum and CMB lensing. We show the cross-spectrum measurement between
two CMB maps created by nulling the contributions from CIB (CIB-free map) and
tSZ (tSZ-free map) to be robust for kSZ as the approach significantly
suppresses the total contribution of CIB and tSZ signals. Similarly, for CMB
lensing, we use the approach introduced by Madhavacheril & Hill (2018) but
with a slight modification by using the tSZ-free and CIB-free maps in the two
legs of the quadratic estimator. By cross-correlating the CMB lensing map
created using this technique with galaxy surveys, we show that the biases from
both CIB/tSZ are negligible. We also compute the impact of unmodeled CIB/tSZ
residuals on kSZ and cosmological parameters finding that the kSZ measured
using the standard ILC to be significantly biased. The kSZ estimate from the
cross-ILC remains less affected by CIB/tSZ making it crucial for current and
upcoming CMB surveys such as the South Pole Telescope (SPT), Simons
Observatory (SO) and CMB-S4. With the cross-ILC approach, we find the total
kSZ power spectrum can be measured at very high significance in the coming
years: $35\sigma$ by SPT, $22\sigma$ by SO, and $80\sigma$ by CMB-S4. Finally,
we forecast constraints on the epoch of reionization using the kSZ power
spectrum and find that the duration of reionization, currently unconstrained
by Planck, can be constrained to $\sigma(\Delta z_{\rm re})$ = 1.5 (or) 0.5
depending on the choice of the prior on the optical depth to reionization. The
data products and the associated codes can be downloaded from this
link${}^{\text{\faGithub}}$.

##  1 Introduction

Observations of the millimeter sky offer a wealth of information about the
origin, contents, and evolution of our Universe as well as astrophysical
effects. In this work, we focus on two late-time effects: gravitational
lensing and the Sunyaev Zel’dovich (SZ) effect. The effect of gravitational
lensing – the bending of the path of the CMB photons due to the presence of
matter in the Universe – is an excellent probe of structure growth. The SZ can
be of two kinds: kinematic (kSZ, Sunyaev & Zeldovich 1980) and thermal (tSZ,
Sunyaev & Zel’dovich 1970, 1972) effects. The tSZ arises from the inverse
Compton scattering of CMB photons off free electron in hot intracluster medium
(ICM). As a result, it is a strong probe of both structure formation and the
complex gas physics of the ICM. The kSZ signal originates when moving
electrons with bulk motion scatter the CMB photons introducing a Doppler
shift. This can arise due to spatial differences in ionization fractions in
the high redshift ($z\geq 6$) universe during the epoch of reionization (EoR,
Paul et al. 2021) and also due to bulk motion of haloes containing free
electrons in the local Universe. From now on, we will refer to these two as
patchy-kSZ and homogeneous-kSZ signals. These two kSZ signals can provide key
information on the physics of the EoR (see Choudhury, 2022, for a review) and
the growth of structure, respectively.

Both CMB lensing (Smith et al., 2007; Das et al., 2011; van Engelen et al.,
2012; Planck Collaboration et al., 2014a; Story et al., 2015; Omori et al.,
2017; Wu et al., 2019; Madhavacheril et al., 2023) and tSZ (Planck
Collaboration et al., 2016a; Madhavacheril et al., 2020; Bleem et al., 2021)
have been detected to very high significance by many experiments. They have
also been used to constrain cosmology and astrophysics (Planck Collaboration
et al., 2016a; Bianchini et al., 2020; Douspis et al., 2021; Sánchez et al.,
2022; Omori et al., 2023; Qu et al., 2023). In addition, the tSZ signal has
also been used to produce catalogs of tens of thousands of galaxy clusters
(Bleem et al., 2015; Planck Collaboration et al., 2016b; Hilton et al., 2020;
Huang et al., 2020) which have also been used as cosmological probes (Planck
Collaboration et al., 2016a; Bocquet et al., 2019, for example). Although the
kSZ signal has also been detected (Hand et al., 2012; Soergel et al., 2016;
Hill et al., 2016; Reichardt et al., 2021; Gorce et al., 2022), the
significance of these detections has not been high enough for them to be used
as cosmological or EoR probes (Chen et al., 2023). Of these, Reichardt et al.
(2021) and Gorce et al. (2022) measured the total kSZ power spectrum which
receives contribution from both the patchy and homogeneous terms while the
other works focussed only on the homogeneous term due to the motion of haloes.
The improvement in the quality of CMB data from current and future surveys
(Benson et al., 2014; Bender et al., 2018; The Simons Observatory
Collaboration et al., 2018; CMB-S4 Collaboration, 2019; Sehgal et al., 2019)
significantly improves the expected detection significance of the lensing and
SZ signals. However, with the improvement in the statistical errors, the
systematic errors in the measurements also become important and non-
negligible.

In this work, we focus on systematic biases in the kSZ power spectrum and CMB
lensing reconstruction due to the presence of unwanted astrophysical
foreground signals in the CMB temperature maps. For kSZ, these contaminants
include tSZ and signals from dusty star-forming galaxies (DSFGs, the
integrated light from which makes up the cosmic infrared background or CIB)
and radio galaxies. For CMB lensing, all the above signals including the kSZ
act as sources of bias.

In the past, the standard approach to measure the kSZ power spectrum has been
to use templates, that were predicted from simulations (Shaw et al., 2010,
2012; Battaglia et al., 2013b, a), for jointly fitting and mitigating the
undesired signals (Choi et al., 2020; Reichardt et al., 2021; Gorce et al.,
2022). The use of templates can be harmful as a mis-estimation of the template
can bias the kSZ results significantly as demonstrated in Appendix B. Gorce et
al. (2022) modified the use of kSZ/tSZ templates and built an emulator to
jointly predict the kSZ/tSZ signals using a random forest algorithm trained
numerically for a range of parameters governing cosmology, cluster
astrophysics, and reionization. This helped them to break the degeneracy
between kSZ and tSZ (See Fig. 3 of Gorce et al., 2022) which led to an almost
$\times 2$ improvement in the kSZ signal-to-noise (S/N). However, Gorce et al.
(2022) used a template for the CIB signal.

In the context of CMB lensing, several methods have been proposed to address
the systematic errors due to foreground signals at the expense of the S/N, for
example, by adopting an aggressive masking strategy (van Engelen et al.,
2014), by using shear-only reconstruction (Schaan & Ferraro, 2019), and by
implementing a foreground bias-hardening technique (Namikawa et al., 2013;
Osborne et al., 2014; Darwish et al., 2021; MacCrann et al., 2023).
Madhavacheril & Hill (2018, hereafter MH18) introduced a modified quadratic
estimator (QE) (Hu & Okamoto, 2002) to mitigate tSZ-induced biases by using a
tSZ-free map in one leg of the QE. While this lensing estimator works well to
mitigate lensing bias from tSZ, the tSZ-free maps generally have an enhanced
level of CIB which might be problematic for future surveys, particularly for
cross-correlation studies (Schmittfull & Seljak, 2018; Abbott et al., 2023;
Chang et al., 2023; Omori et al., 2023) that rely on CMB temperature-based
lensing reconstruction.

In this study, we use a different approach and propose to use the cross-power
spectra between two different CMB maps constructed using different internal
linear combination (ILC, Cardoso et al. 2008; Planck Collaboration et al.
2014b) of the different frequency bands. The two CMB maps are different
because they have been derived in a way to null the response of different
foreground signals in them using the constrained-ILC (cILC) approach
(Remazeilles et al., 2011). We refer to this cross-power spectrum approach as
the “cross-ILC technique”. Given that CIB and tSZ form the most important
sources of biases for kSZ, the cross-ILC technique uses the tSZ-free and CIB-
free111CIB signals cannot be fully removed from the map using cILC as they are
made up of multiple populations. As a result, a CIB-free map is generally a
CIB-minimized map and will never be truly devoid of CIB. CMB maps for the
power spectrum estimation. Along these lines, we note here that Kusiak et al.
(2023) recently explored nulling multiple foreground components using large-
scale structure tracers from galaxy surveys. Likewise, for CMB lensing, we use
the MH18 estimator but with a slight alteration of using tSZ-free and CIB-free
CMB maps in the two legs of the QE.

The paper is structured as follows: In §2, we describe the experimental setups
used for forecasting followed by the cross-ILC method and the simulations used
for estimating the systematic errors. We present and discuss the results in §3
for both kSZ and CMB lensing which include the importance of systematic
errors, mitigation strategy using the cross-ILC technique, expected kSZ S/N
and constraints on the EoR. Finally, we discuss the potential applications of
the cross-ILC technique and conclude in §4. In appendixes, we present the
bandpower errors for kSZ; compute biases in kSZ and cosmological parameters
due to residual foregrounds; discuss alternate foreground models; and
demonstrate biases in the cross-correlation of CMB lensing and galaxy surveys
for multiple redshift bins.

The underlying cosmology used in this work was set to Planck 2018 measurements
(TT, TE, EE + lowE + lensing in Table 2 of Planck Collaboration et al. 2020).
For forecasting the S/N, we assume a kSZ power spectrum level of
$D_{\ell}=3\mu{\rm K}^{2}$ where $D_{\ell}=\ell(\ell+1)C_{\ell}/2\pi$.

##  2 Method

###  2.1 Experimental Specifications

Table 1: Band-dependent beam and white noise levels for different experiments
considered in this work. Note that we have made a conservative choice and
ignored frequency channels below 90 GHz in this work.

Experiment | Beam $\theta_{\rm FWHM}$ in arcminutes ($\Delta_{T}$ in $\mu{\rm K{\text{-}}arcmin}$)
---|---
90 GHz | 150 GHz | 220 GHz | 285 GHz | 345 GHz | 410 GHz | 850 GHz
SPT-3G | 1.7 (3) | 1.2 (2.2) | 1 (8.8) | - | - | -
\+ SPT-3G+ | 1 (2.3aaInverse variance weighted noise estimate of SPT-3G 220 GHz and SPT-3G+ 225 GHz channels.) | 0.55 (5.6) | 0.45 (40.2)
SO-Baseline | 2.2 (8) | 1.4 (10) | 1 (22) | 0.9 (54) | -
\+ FYST | 1 (12.3bbInverse variance weighted noise estimate the respective experiment and FYST 285 GHz channels.) | 0.8 (24.8ccInverse variance weighted noise estimate the respective experiment and FYST 345 GHz channels.) | 0.6 (107) | 0.5 (407) | 0.3 ($6.8\times 10^{5}$)
SO-Goal | 2.2 (5.8) | 1.4 (6.3) | 1 (15) | 0.9 (37) | -
\+ FYST | 1 (10.6bbInverse variance weighted noise estimate the respective experiment and FYST 285 GHz channels.) | 0.8 (22.3ccInverse variance weighted noise estimate the respective experiment and FYST 345 GHz channels.) | 0.6 (107) | 0.5 (407) | 0.3 ($6.8\times 10^{5}$)
S4-Wide | 2.2 (1.9) | 1.4 (2.09) | 1 (6.9) | 0.9 (16.88) | -
S4-Ultra Deep$\dagger$$\dagger$footnotemark:  | 3.0 (0.45) | 1.92 (0.41) | 1.32 (1.3) | 1.2 (3.1) | -

††footnotetext: Five meter Three-mirror anastigmat (TMA) telescope design
(CMB-S4 Collaboration, 2019).

In this study, we consider the current survey being conducted on the South
Pole Telescope (SPT) with the SPT-3G camera (Benson et al., 2014; Bender et
al., 2018), and the upcoming CMB experiments Simons Observatory (SO, The
Simons Observatory Collaboration et al. 2018) and CMB-S4 (CMB-S4
Collaboration, 2019) as our baseline surveys. We use both the wide (S4-Wide)
and the deep (S4-Ultra Deep) CMB-S4 surveys (CMB-S4 Collaboration, 2019). For
SO, we consider both the nominal SO-Baseline and also its alternative SO-Goal
configurations (The Simons Observatory Collaboration et al., 2018). Since
mitigating the contamination from CIB is crucial for kSZ, we also include the
high frequency (HF) channels from SPT-3G+ (Anderson et al., 2022) for SPT-3G
and Fred Young Submillimeter Telescope (FYST, CCAT-Prime Collaboration et al.
2023) for SO. SPT-3G+ is a proposed successor of SPT-3G which will observe the
same 1500 ${\rm deg}^{2}$ footprint of SPT-3G but with three HF bands: 220,
285, and 345 GHz (Anderson et al., 2022).

The experimental configurations (beams and noise levels for each band) of the
surveys are given in Table 1. As indicated in the table, we make inverse
variance weighted combinations of the overlapping bands when the HF channels
are included. The specifications for modeling the atmospheric noise (Dibert et
al., 2022) for SPT and CMB-S4 can be found in Table 2. For SO, we use the
values quoted in The Simons Observatory Collaboration et al. (2018) similar to
the setup described in Raghunathan (2022). Throughout this study, we ignore
data from channels below 90 GHz as they are primarily used to reduce
contamination from the galactic synchrotron signals on large-scales which have
negligible impact on both kSZ and CMB lensing.

Table 2: Atmospheric $1/f$ noise specifications ($\ell_{\rm knee}$, $\alpha_{\rm knee}$) for SPT and CMB-S4 experiments. Band [GHz] | SPT | S4-Ultra Deep | S4-Wide
---|---|---|---
90 | 1200, 3.0 | 1200, 4.2 | 1932, 3.5
150 | 2200, 4.0 | 1900, 4.1 | 3917, 3.5
220 | 2100, 3.9 | 2100, 4.1 | 6740, 3.5
285 | 2100, 3.9 | 2100, 3.9 | 6792, 3.5
345 | 2600, 3.9 | - | -

###  2.2 Simulations

Thanks to the advancement in the field of computational astrophysics, we have
several correlated multi-component simulations (Sehgal et al., 2010; Stein et
al., 2020; Omori, 2022) that are crucial to test and understand the biases in
the measurements from future CMB surveys. We primarily use the Agora
simulation released recently by Omori (2022) to quantify the systematics
introduced by CIB, tSZ, and radio signals in both kSZ and lensing
measurements. We do not note significant differences when replacing the
frequency dependent ILC weights computed using Agora with a foreground model
based on SPT measurements (R21). More details about this validation can be
found in Appendix C.1. As a further check, we replace the Agora with Websky
(Stein et al., 2020) simulations in Appendix C.2. We limit these checks to the
S4-Wide survey.

####  2.2.1 Agora

We provide a brief description of Agora simulations here and refer the
interested readers to the original paper (Omori, 2022) for more details. Agora
is a set of correlated extragalactic maps generated based on the MultiDark
Planck 2 (MDPL2; Klypin et al. 2016) simulation using its dark matter
particles and halo catalog. The simulation consists of high resolution maps of
lensed CMB, tSZ, kSZ, CIB and radio sources as well as lensing maps. The tSZ
map is generated by pasting halo profiles that are fit to the BAHAMAS
hydrodynamical simulation (McCarthy et al., 2017; Mead et al., 2020) onto the
MDPL2 haloes.222Halo catalogs are publicly available
http://halos.as.arizona.edu/simulations/MDPL2/hlists/ The CIB and radio
catalogs are based on the MDPL2 UniverseMachine catalogs (Behroozi et al.,
2019).333Available at
http://halos.as.arizona.edu/UniverseMachine/DR1/MDPL2_SFR/.

Agora has been verified to not only reproduce auto-/cross- frequency spectra
from data but also component separated maps such as Compton-$y$ maps (Omori,
2022). In this study, we run measurements on the uncalibrated444Uncalibrated
here means that the simulated maps have been convolved with the experiment-
dependent frequency band pass functions, but have not been adjusted to match
with the measured power spectra, which is typically a 5-10% offset. frequency
maps.

####  2.2.2 Sky components and masking

For the kSZ study, we use the CIB, radio and tSZ signals from both Agora and
Websky simulations. We identify the locations of dusty and radio point sources
with flux $S_{\rm 150}\geq 3\ {\rm mJy}$ which roughly corresponds to $\geq
10\sigma$ detection limit from the surveys considered in this work. Note that
is a conservative choice and we can reduce the masking threshold down to
$5\sigma$ which corresponds to $S_{\rm 150}\sim 1-2\ {\rm mJy}$ (CMB-S4
Collaboration, 2019; The Simons Observatory Collaboration et al., 2018). We
only mask a single pixel for the point sources. We also mask the location of
clusters with mass $M_{\rm 500c}\geq 2\times 10^{14}\ \mbox{M}_{\odot}$ which
roughly corresponds to the $10\sigma$ detection limit (Raghunathan et al.,
2022; Raghunathan, 2022). For the cluster mask, we do not modify the masking
threshold as a function of redshift. However, the size of the mask changes
based on the redshift of the halo and we define the masking radius to be
$2\theta_{500c}$.

For CMB lensing, besides the astrophysical CIB, radio, and SZ signals, we also
use the lensed CMB maps with the associated convergence field. The point
source masking threshold used for CMB lensing is higher ($S_{\rm 150}\geq 6\
{\rm mJy}$) than the one adopted for kSZ. This choice is primarily to reduce
the mean-field lensing bias in the reconstructed lensing map. The cluster
masking threshold is $10\sigma$.

###  2.3 Internal Linear Combination

Figure 1: The total ILC residuals in the MV CMB map. The dash-dotted gray line
shows the primary CMB TT power spectrum while the solid gray line corresponds
to a rough kSZ level of $D_{\ell}=3\mu{\rm K}^{2}$. CMB dominates on large
scales ($\ell\lesssim 4000$) while the residual noise and foregrounds are
higher than kSZ on small scales. We note that the kSZ S/N for individual modes
is always lower than one for all experiments. Figure 2: The residual CIB
(left) and tSZ (right) signals compared to the expected kSZ level of
$D_{\ell}=3\mu{\rm K}^{2}$ in the MV CMB/kSZ maps from different experiments
considered in this work.

We combine data from different frequency channels ${\rm N}_{\rm ch}$ using the
ILC technique in the harmonic space as

$\displaystyle S_{\ell}=\sum_{i=1}^{\rm N_{\rm ch}}w_{\ell}^{i}M_{\ell}^{i},$
(1)

where $S$ corresponds to the desired sky signal, which in this case is the sum
of CMB and kSZ signals. The multipole dependent weights $w_{\ell}$ for each
frequency channel $i$ given in Eq.(2) are computed to minimize the total
variance from experimental noise and foregrounds to obtain the minimum-
variance (MV) signal map. The weights can also be tuned to produce a minimum-
variance map along with an additional constraint of nulling the contribution
to a specific frequency response. This is the cILC approach. The weights are
constructed as

$\displaystyle w_{\ell}^{\rm cILC}$ $\displaystyle=$ $\displaystyle{\bf
C}_{\rm\ell}^{-1}\mathcal{F}\left(\mathcal{F}^{\dagger}{\bf
C}_{\rm\ell}^{-1}\mathcal{F}\right)^{-1}N,$ (2)

where the matrix ${\bf C}_{\ell}$ has a dimension $\rm N_{ch}\times\rm N_{ch}$
and contains the covariance between simulated maps in multiple frequencies at
a given multipole $\ell$; $\mathcal{F}=[A_{S}\ B_{S}\ C_{S}\ ...\ Z_{S}]_{{\rm
N}_{\rm ch}\times S}$ contains the frequency response vector of the signal
component $S$, either the desired signal $A_{S}$ with dimension $\rm N_{\rm
ch}\times 1$ or the undesired sky components $[B_{S}\ C_{S}\ ...\ Z_{S}]$ that
are being nulled using the cILC (Remazeilles et al., 2011) technique which is
specified using $N=[1\ 0\ 0\ ..\ 0]_{S\times 1}$.

For a standard MV ILC, Eq.(2) simplifies to $w_{\ell}^{\rm MV}=\dfrac{{\bf
C}_{\rm\ell}^{-1}A_{S}}{A_{S}^{\dagger}{\bf C}_{\rm\ell}^{-1}A_{S}}$. In Fig.
1, we show the total (noise and foregrounds) ILC residuals for the standard MV
ILC estimator for different experiments considered in this work. The
foreground signals are from Agora simulations. The primary CMB (gray dash-
dotted curve) is much larger than the kSZ (gray solid curve) on large scales
($\ell\lesssim 4000$) while the small scales are dominated by instrumental
noise and foregrounds.

The instrument noise is generally easy to model. The contribution from radio
point sources in desired frequency band $\nu$ below the masking threshold can
also be modeled relatively easily assuming an underlying source distribution
$dN/dS$ as we show later. On the other hand, the residual CIB and tSZ signals
are harder to model. The contributions from CIB and tSZ in the MV CMB map
relative to the kSZ signal $D_{\ell}=3\mu{\rm K}^{2}$ are presented in Fig. 2
in the left and right panels respectively. As it is evident from the figure,
the residual CIB and tSZ signals are non negligible for SPT-3G (green),
S4-Wide (black), SO-Baseline (blue), and SO-Goal (orange). The CIB residuals
are $\sim\times 2$ lower than kSZ for SPT-3G+ (yellow) and S4-Ultra Deep (red)
although the tSZ residuals become higher than kSZ for $\ell\geq 4000$. This
suggests that using the MV ILC map for kSZ measurements could lead to
significant biases depending on the uncertainties in modeling the CIB and tSZ
signals (see Appendix B).

###  2.4 Cross-ILC

Figure 3: The total CIB and tSZ residuals expected in the S4-Wide survey for
different ILC combinations: Think lines are for auto-spectra and thick lines
are for cross-ILC. We note that the residuals in cross-ILC combination (tSZ-
free $\times$ CIB-free in black) is much lower than kSZ ($D_{\ell}=3\mu{\rm
K}^{2}$) and hence can be a compelling means of measuring the kSZ power
spectrum with future CMB surveys.

As discussed above, it is possible to reduce the contamination from CIB/tSZ
using the cILC method in Eq.(2). This can be achieved in multiple ways. For
example,

  * (A)

Creating a single ILC map by jointly nulling tSZ and CIB using $B_{S}$ and
$C_{S}$ in Eq.(2) where $B_{S}$ and $C_{S}$ are the frequency response of tSZ
and CIB SED being nulled. In this case, we will take the auto-spectrum of this
cILC map.

  * (B)

Creating two cILC maps. The first map is a tSZ-nulled map such that $B_{S}$ is
the frequency response of the tSZ signal. The second map is a CIB-nulled map
and here $B_{S},C_{S}$ correspond to the CIB SED(s) being nulled. In this
case, we will then take the cross-spectrum of the two cILC maps.

  * (C)

This is similar to the second method but the CIB SED(s) for nulling are chosen
such that the total CIB+tSZ residuals are lower in the cross-spectrum. _This
is the approach we adopt in this work._

Approach (A): If CIB signal is made up of a single population (i.e:) if the
frequency response can be described by a single SED, then the first approach
is the optimal for mitigating CIB. However, the CIB signal is made up of
multiple galaxy populations at different redshifts and nulling the tSZ
significantly enhances the contamination from CIB populations that are not
included in Eq.(2). Moreover, this approach increases the noise significantly
by $\times 2-\times 3$ compared to the MV ILC.

Approach (B): The noise in approach (B) is lower than (A). However, it does
not take the cross-correlation CIB$\times$tSZ into account and hence not
ideal.

Approach (C): The third approach (C) is similar to (B), but we optimise the
weights for CIB nulling in a manner than reduces the total CIB+tSZ residuals
in the CMB map. Thus, it also takes the cross-correlation CIB$\times$tSZ into
account. We describe the steps to obtain weights for tSZ and CIB nulling
below.

In Fig. 3, we compare the total CIB+tSZ residuals expected in S4-Wide for
different ILC techniques. Thin lines are the auto-spectra of the ILC maps:
blue is for the MV, red is for tSZ-free, and orange is for a multi-component
CIB-nulled ILC map. Thick lines are for the cross-ILC technique obtained as
the cross-spectrum of two different ILC maps. From the figure, we note that
the cross-ILC spectrum obtained from tSZ-free and CIB-free maps shown in black
returns the best total CIB+tSZ residual and is almost an order of magnitude
lower than kSZ for $\ell\in[3500,5000]$. As a result, we use this cross-ILC
combination as the baseline in the rest of the paper. We provide further
details in §3.1.1.

_Weights for tSZ and CIB nulling:_ The weights for tSZ-free map are obtained
by setting $B_{s}=[-4.36,-2.61,-0.1,2.94,5.72,8.79,29.9]$, which correspond to
the frequency response of the tSZ signal in $\nu=[90,150,220,285,345,410,850]$
GHz bands, in Eq. 2. The CIB signal is generally modeled as a modified
blackbody as $\eta_{\nu}=\nu^{\beta_{\rm CIB}}B_{\nu}(T_{\rm CIB})$ where
$T_{\rm CIB}$ is the temperature, $\beta_{\rm CIB}$ is the emissivity index,
and $B_{\nu}(T)$ is the Planck function. Getting the frequency response of the
CIB signal is slightly tricky as the CIB is made up of multiple populations of
DSFGs with different values of $T_{\rm CIB}$ and $\beta_{\rm CIB}$. As a
result, we do not adopt the same values of $T_{\rm CIB}$ and $\beta_{\rm CIB}$
for all experiments. Instead, we perform a blind four parameter ($T^{1}_{\rm
CIB},\beta^{1}_{\rm CIB},T^{2}_{\rm CIB},\beta^{2}_{\rm CIB}$) grid search to
obtain a two component SED for CIB-nulling and the frequency response of these
two SEDs are plugged into $B_{s}$ and $C_{s}$ of Eq. 2. The SEDs that return
the lowest level of CIB+tSZ residual in the multipole window
$\ell\in[3000,5000]$ (gray band in Fig. 3 and Fig. 4) are chosen. We use this
range for kSZ as the sample variance from CMB fully swamps the kSZ detection
at $\ell\leq 3000$. For lensing, we use a slightly different range
$\ell\in[2000,5000]$. We also checked the results with a single SED or a
three-component SED model and find that the two-component SED works the best
in terms of both $S/N$ and foreground residuals. Note that the above weights
are calculated using Agora foreground model but, as mentioned earlier, we
validate this assumption using a foreground model based on SPT measurements
(R21) and Websky simulations (Stein et al., 2020) in Appendix C.

####  2.4.1 Fisher formalism

Table 3: Fiducial values of the parameters and priors used in this work. All
the applied priors are Gaussian with widths given below in the table.

Parameter | Fiducial | Prior
---|---|---
Cosmological parameters:
Amplitude of scalar fluctuations ${\rm ln}(10^{10}A_{s})$ | 3.044 | -
Dark matter density $\Omega_{c}h^{2}$ | 0.1200
Baryon density $\Omega_{b}h^{2}$ | 0.02237
Scalar spectral index $n_{s}$ | 0.9649
Angular size of sound horizon | 1.04092
at recombination $100\theta_{\ast}$
Reionization optical depth $\tau_{\rm re}$ | 0.0544 | 0.007${}^{a}$${}^{a}$footnotemark:
0.002${}^{b}$${}^{b}$footnotemark:
Foreground parameters:
Residual CIB+tSZ $A_{\rm CIB+tSZ}$ | 1 | 0.1
Mean spectral index of sources $\alpha_{\rm rad}$ | -0.76 | 0.1
Scatter in source spectral indices $\sigma(\alpha_{\rm rad})$ | 0.2 | 0.3
Total kSZ power spectrum:
Amplitude $A_{\rm kSZ}$ | 1 | -
Homogeneous kSZ power spectrum:
Amplitude $A_{\rm kSZ}^{h}$ | 1 | 0.1
Spectral tilt $\alpha_{\rm kSZ}^{h}$ | 0 | 0.1
Reionization kSZ:
Mid-point of reionization $z_{\rm re}^{\rm mid}$ | 7.69 | 1.43${}^{a\dagger}$${}^{a\dagger}$footnotemark:
0.41${}^{b\dagger}$${}^{b\dagger}$footnotemark:
Duration of reionization $\Delta z_{\rm re}$ | 4 | -

aafootnotetext: Planck-like $\tau_{\rm re}$ prior.

bbfootnotetext: LiteBIRD-like $\tau_{\rm re}$ prior.

††footnotetext: Planck-like or LiteBIRD-like $\tau_{\rm re}$ prior translated
to prior on $z_{\rm re}^{\rm mid}$.

We use a Fisher formalism to forecast the expected S/N for the total kSZ power
spectrum. For this purpose, we use the cross-ILC combination for all
experiments and also compare it with the S/N obtained from the MV ILC.
Although MV ILC is expected to return the best S/N, we note from Fig. 2 and
Fig. 3, the total ILC residuals receive significant contribution from CIB and
tSZ; and hence will be prone to foreground-induced biases.

For kSZ forecasts, we consider two approaches. In the first case, we compute
the total kSZ S/N and for this we use a ten parameter model of which six are
the standard $\Lambda{\rm CDM}$ parameters $\theta\in[A_{s},\ h,\ n_{s},\
\Omega_{c}h^{2},\ \Omega_{b}h^{2},\ \tau_{\rm re}]$ and the additional
parameters are: $A_{\rm kSZ}$ which quantifies the amplitude of the kSZ signal
$D_{\ell}=A_{\rm kSZ}\times 3\mu{\rm K}^{2}$; $A_{\rm CIB+tSZ}$ which
represents the amplitude of the total CIB and tSZ residuals; $\alpha_{\rm
rad}$ and $\sigma(\alpha_{\rm rad})$ which represent the mean spectral index
and the scatter when modeling the radio residuals.

In the second approach, we replace the total kSZ amplitude $A_{\rm kSZ}$ with
parameters that govern the physics of reionization namely: $z_{\rm re}^{\rm
mid}$ which is the mid-point of reionization and $\Delta z_{\rm re}$ which is
the duration of reionization. Since, it is hard to disentangle between the two
kSZ components (homogeneous and patchy kSZ) from the total kSZ power spectrum,
we include two more parameters along with suitable priors $A_{\rm kSZ}^{h}$
and $\alpha_{\rm kSZ}^{h}$ that represent the amplitude and the slope of the
homogeneous kSZ power spectrum. Thus in the second case, we fit for 13
parameters.

We use CMB TT/EE/TE power spectra since EE/TE will dominate the $\Lambda{\rm
CDM}$ parameter constraints from future CMB surveys (Calabrese et al., 2014;
Galli et al., 2014). For EE/TE, we use data in the multipole range
$\ell\in[30,5000]$. For TT, we use a slightly different strategy. To better
predict the radio residuals, we use $\ell_{\rm max}^{\rm TT}=5000$. However,
following other similar works in the literature (Calabrese et al., 2014; The
Simons Observatory Collaboration et al., 2018; Alvarez et al., 2021), we set
different values of $\ell_{\rm max}^{\rm TT,kSZ}\in[3000,3500,4000,4500]$ and
do not consider modes above $\ell_{\rm max}^{\rm TT,kSZ}$ for kSZ to account
for the small-scale foreground residuals. We do this by setting the
derivatives of the kSZ power spectrum $\dfrac{\partial kSZ}{\partial A_{\rm
kSZ}}=0$ above $\ell_{\rm max}^{\rm TT,kSZ}$. The fiducial values of the
parameters and the priors are listed in Table 3.

##  3 Results and Discussion

###  3.1 Kinematic SZ

Figure 4: The ratio of the total CIB+tSZ residual compared to the expected
level of kSZ ($D_{\ell}=3\mu{\rm K}^{2}$) for the cross-ILC combination tSZ-
free $\times$ CIB-free from different CMB surveys. In the left, we show the
residuals for the baseline experiments while the right panel is after the
inclusion of HF bands from SPT-3G+ or FYST.

In this section, we first show the improvement in CIB and tSZ residual levels
using the cross-ILC technique for all the experiments. Besides the CIB+tSZ
residuals, the CMB maps also contain the residual signals from instrumental
noise and radio galaxies. We handle these next. This is followed by the
calculation of the expected kSZ SNR using the Fisher formalism.

####  3.1.1 Total CIB+tSZ residuals

In Fig. 4, we present the ratio of the total CIB+tSZ residuals over the
expected kSZ signal ($D_{\ell}=3\mu{\rm K}^{2}$) for the cross-ILC combination
tSZ-free $\times$ CIB-free for different experiments. The left panel contains
the residuals for the baseline experiments. As it is evident, the total
CIB+tSZ residuals are always lower than kSZ for the cross-ILC combination. To
be specific, the residuals are roughly $\times 2$ ($\times 3$) lower than kSZ
in $\ell\in[3000,5000]$ ($\ell\in[3500,5000]$) for all experiments. For
CMB-S4, the residuals are $\leq\times 10$ lower than kSZ in the range
$\ell\in[3500,5000]$. By comparing Fig. 4 to Fig. 2, we can deduce that the
kSZ power spectrum can be measured by all the experiments using the cross-ILC
combination much more robustly compared to the MV ILC estimator.

_Inclusion of HF bands:_ In the right panel of Fig. 4, we assess the
performance after the inclusion of HF bands to the baseline experiments:
SPT-3G+ to SPT-3G and FYST to SO. Since the residuals are already small for
S4-Wide and S4-Ultra Deep, we do not include FYST for CMB-S4. For SPT, we find
a significant improvement in the reduction of CIB+tSZ residuals after
including information from SPT-3G+. The noise, although not shown here, also
reduces for SPT-3G. For SO, however, including FYST does not improve the
residuals. In fact, from the right panel of Fig. 4, we note that the residuals
have slightly increased. On the other hand, addition of FYST does reduce the
overall noise in the ILC maps by $\times 1.5$ to $\times 2$ compared to SO-
only at $\ell\in[4000,5000]$. This behavior is not surprising and it is
because of the higher noise levels of SO and FYST compared to SPT-3G and
SPT-3G+. In the case of the former, the ILC algorithm optimises the weights to
reduce the overall noise primarily while for the latter, the weights are
optimised to reduce the foreground residuals. The noise reduction ultimately
helps in improving the final kSZ S/N for SO as shown in Table 4.

We note here that the total CIB+tSZ residuals can be lowered further for SO +
FYST at the expense of a higher noise by introducing a scaling term for noise
or CIB in the covariance matrix ${\bf C}_{\rm\ell}$ in Eq.(2) used for
computing the ILC weights as prescribed by Bleem et al. (2021). While this
approach is not optimal for S/N, it does help in reducing the residuals.
Following this idea, when we scaling the noise levels of FYST’s HF channels
($\nu\in[345,410,857]$) by $\times 5$, we find the CIB+tSZ residuals at
$\ell\in[4000,5000]$ to go down by $\times 1.5$ for SO-Goal.

####  3.1.2 Handling radio residuals

Given that the power spectrum of the radio galaxies follow a Poisson
distribution, they can modeled easily compared to other foreground signals. In
this work, we use a point source masking threshold of $S_{150}^{\rm max}$ = 3
mJy. This is a conservative choice and the masking threshold can be lowered
further for both current and future experiments. We model the power from
sources below the masking threshold by assuming an underlying source
distribution $dN/dS$ along with a mean spectral index $\alpha_{\rm rad}$ and a
scatter $\sigma(\alpha_{\rm rad})$ as

$\displaystyle C_{\ell_{{\nu_{1}\nu_{2}}}}^{\rm radio}=\bigintsss_{\alpha_{\rm
rad}^{\rm min}}^{\alpha_{\rm rad}^{\rm
max}}d\alpha\bigintsss_{0}^{S_{150}^{\rm max}}dS_{150}\ S_{150}^{2}$ (3)
$\displaystyle\frac{dN}{dS_{150}}\
\left(\frac{\nu_{0}^{2}}{\nu_{1}\nu_{2}}\right)^{\alpha}\mathcal{N}\left[\alpha_{\rm
rad}|\bar{\alpha}_{\rm rad},\sigma(\alpha_{\rm rad})\right],$

where the normal distribution $\mathcal{N}\left(\mu,\sigma\right)$ is used to
parameterize the probability density function of source spectral indices.

Figure 5: Radio residuals using Agora expected from S4-Wide for the masking
threshold $S_{150}=3\ {\rm mJy}$ (yellow curve) used in this work. They are
negative for the cross-ILC tSZ-free $\times$ CIB-free combination. We also
show the analytical modeling of the radio residuals using Lagache et al.
(2020) in the yellow band. The band quantifies the residual power assuming a
scatter in the radio spectral index $\sigma(\alpha_{\rm rad})\in[0,0.35]$. For
reference, we also show the residuals and the models for multiple masking
thresholds: $S_{150}=1\ {\rm mJy}$ in blue, $S_{150}=2\ {\rm mJy}$ in green,
and $S_{150}=4\ {\rm mJy}$ in red. The yellow dotted curve near the bottom
right is the residual expected when a radio spectrum $\alpha_{\rm rad}=-0.76$
is nulled along with tSZ.

The residual radio signal in the ILC combination can now be calculated
analytically using Eq.(4) as

$\displaystyle C_{\ell_{\rm ILC}}^{\rm radio}$ $\displaystyle=$ $\displaystyle
w_{\ell_{A}}{\bf C}_{\rm\ell}^{\rm radio}\ w_{\ell_{B}}^{\dagger}$ (4)

where ${\bf C}_{\rm\ell}^{\rm radio}$ is the ${\rm N}_{\rm ch}\times{\rm
N}_{\rm ch}$ covariance matrix containing the auto- and cross-power spectra
between different frequency bands; and $w_{\ell_{A}}$ and $w_{\ell_{B}}$ are
the frequency dependent weights for the two ILC maps $A$ and $B$.

In Fig. 5 we show the radio residual (yellow curve) for the cross-ILC tSZ-free
$\times$ CIB-free combination expected in the S4-Wide survey after masking
source above $S_{150}=3\ {\rm mJy}$, our baseline masking threshold. The radio
residuals are negative for this cross-ILC combination. For the analytical
prediction of the radio residuals, we use Lagache et al. (2020) radio source
distribution model $dN/dS$ with $\alpha_{\rm rad}=-0.76$ (R21). The yellow
band in the figure corresponds to radio residuals assuming a scatter in the
spectral index $\sigma(\alpha_{\rm rad})\in[0,0.35]$. As can be seen from Fig.
5, the width of the band is: $D_{\ell_{\rm radio}}^{\rm sys}=[0.4,0.8,1.0]$
for $\ell\in[3000,4000,5000]$ and these are smaller than the expected kSZ
level of $D_{\ell}=3\mu{\rm K}^{2}$. For reference, we also show the radio
residuals and the predictions for other masking threshold: $S_{150}=1\ {\rm
mJy}$ in blue, $S_{150}=2\ {\rm mJy}$ in green, and $S_{150}=4\ {\rm mJy}$ in
red. When computing the kSZ SNR in § 3.1.4, we account for the radio residuals
using two parameters $\alpha_{\rm rad}$ and $\sigma(\alpha_{\rm rad})$.

_Nulling radio along with tSZ:_ Another approach to mitigate the radio
residuals is to null the response to radio sources along with the tSZ removal
in the first leg of the cross-ILC. This can be achieved by setting the
frequency responses $B_{s},C_{s}$ in Eq.(2): $B_{s}=[-4.36,-2.61,-0.1,2.94]$
is the frequency response of tSZ and $C_{s}=[2.89,1.0,0.62,0]$ in the
frequency response of radio signals for $\alpha=-0.76$ in the bands
$\nu=[90,150,220,285]$.

We show the reduction in radio power using this method in Fig. 5 for our
baseline masking threshold $S_{150}=3\ {\rm mJy}$ as the yellow dotted curve
(near the bottom right corner). Since not all the radio sources have the same
spectral index $\alpha_{\rm rad}=-0.76$, we can expect to have some amount of
residual signal as can be seen in Fig. 5. As it is evident from the yellow
dotted curve, this residual is an order of magnitude smaller than kSZ and
hence negligible. This curve also includes a band around it with
$\sigma(\alpha_{\rm rad})=0.4$ although the width of the band is small and
hence not visible. Note that including an additional constraint for nulling
radio along with tSZ in cILC will increase the noise level of the resultant
map compared to the tSZ-free map. Hence, we present this method only as a
proof of concept and do not pursue it further in this work.

Table 4: Forecasted S/N of the total kSZ power spectrum $D_{\ell}=3\ \mu
K^{2}$ for different experiments using the
Fisher\text{\faGithub}\text{\faGithub}Jupyter notebooks
s1_get_fisher_matrix.ipynb and s2_analyse_fisher_matrix.ipynb used for
forecasting are publicly available. formalism. These are obtained after
marginalizing over the $\Lambda{\rm CDM}$ parameters and $A_{\rm CIB+tSZ}$.
Planck-like prior on $\sigma(\tau_{\rm re})=0.007$ and a 10% prior on
$\sigma(A_{\rm CIB+tSZ})=0.1$ have been assumed. The radio residual modelling
parameters $\alpha_{\rm rad}$ and $\sigma(\alpha_{\rm rad})$ are fixed for
this table. The associated biases in $A_{\rm kSZ}$ due to unmodeled CIB and
tSZ residuals are presented in Fig. 9 of Appendix B.

Experiment | $f_{\rm sky}$ | Total kSZ SNR = $A_{\rm kSZ}/\sigma(A_{\rm kSZ})$
---|---|---
$\ell_{\rm max}^{\rm TT,kSZ}=3000$ | $\ell_{\rm max}^{\rm TT,kSZ}=3500$ | $\ell_{\rm max}^{\rm TT,kSZ}=4000$ | $\ell_{\rm max}^{\rm TT,kSZ}=4500$
MV | Cross-ILC | MV | Cross-ILC | MV | Cross-ILC | MV | Cross-ILC
SPT-3G | 3% | 10.33 | 8.39 | 19.73 | 12.94 | 27.68 | 16.55 | 27.09 | 19.11
\+ SPT-3G+ | 11.51 | 10.90 | 23.61 | 20.23 | 35.76 | 29.20 | 37.20 | 35.77
SO-Baseline | 40% | 15.26 | 5.13 | 24.64 | 7.86 | 30.06 | 10.71 | 25.80 | 13.13
+FYST | 15.36 | 7.30 | 24.90 | 11.44 | 30.42 | 15.82 | 26.18 | 19.60
SO-Goal | 40% | 16.12 | 5.45 | 27.05 | 8.51 | 35.54 | 11.88 | 32.81 | 14.90
+FYST | 16.21 | 7.54 | 27.33 | 12.04 | 35.98 | 16.98 | 33.29 | 21.37
S4-Wide | 50%$\dagger$$\dagger$footnotemark:  | 41.68 | 32.03 | 78.45 | 48.55 | 109.33 | 67.48 | 107.79 | 83.68
S4-Ultra Deep | 3% | 12.87 | 12.37 | 28.45 | 25.77 | 45.41 | 41.50 | 48.48 | 54.11

††footnotetext: Although S4-Wide will map 67% of the sky, we exclude the
regions that are significantly contaminated by galactic foregrounds. See Fig.
1 of Raghunathan et al. (2022).

####  3.1.3 Bandpower errors

Besides the tSZ+CIB residuals, the CMB maps contain residual noise which also
contributes to the overall variance. The noise can be modeled easily and the
noise residuals $N_{\ell_{\rm ILC}}$ in the ILC maps can be computed by
replacing ${\bf C}_{\rm\ell}^{\rm radio}$ in Eq.(4) with ${\bf N}_{\ell}$,
which is a ${\rm N}_{\rm ch}\times{\rm N}_{\rm ch}$ covariance matrix
containing the auto- and cross-noise spectra between different frequency
bands. ${\bf N}_{\ell}$ is generally diagonal but given that current and
future CMB surveys use multi-chroic pixels, the atmospheric noise will be
correlated between bands that share the same pixel on the focal plane. In this
work, we assume that the atmospheric noise is correlated at 90% level between
the adjacent bands (The Simons Observatory Collaboration et al., 2018; CMB-S4
Collaboration, 2019). Tweaking the correlation level does not affect any of
our results. For the case of white noise, ${\bf N}_{\ell}$ is diagonal.

With the foreground and noise residuals in hand, we can now calculate the
bandpower errors as Eq.(5 Knox 1995)

$\displaystyle\Delta\hat{C}_{\ell}$ $\displaystyle=$
$\displaystyle\sqrt{\dfrac{2}{(2\ell+1)f_{\rm
sky}\Delta_{\ell}}}\hat{C}_{\ell}$ (5)

where

$\displaystyle\hat{C}_{\ell}$ $\displaystyle=$
$\displaystyle\sqrt{\dfrac{\hat{C}_{\ell_{AB}}^{2}+\hat{C}_{\ell_{AA}}\hat{C}_{\ell_{BB}}}{2}}$
(6)

and $\hat{C}_{\ell}\equiv C_{\ell_{\rm ILC}}+N_{\ell_{\rm
ILC}}=C_{\ell}+N_{\ell}$ is the sum of CMB, kSZ, residual foregrounds and
noise in the two ILC maps $A$ and $B$.

It is to be noted that although the cross-ILC combination returns a lower
level of CIB+tSZ residuals compared to MV-ILC, modifying the ILC weights to
suppress the foregrounds introduces a noise penalty. We demonstrate this in
Fig. 8 in Appendix A. As a result the total kSZ SNR in the cross-ILC
combination will be lower than the MV ILC. Despite this noise penalty, the low
CIB+tSZ residuals in the cross-ILC combination will allow us to measure the
kSZ signal robustly compared to the MV as discussed in Appendix B.

####  3.1.4 Fisher forecasts

In Table 4 we present the Fisher forecasts for the total kSZ power spectrum
$D_{\ell}=A_{\rm kSZ}\times 3\ \mu{\rm K}^{2}$ for both the MV- and the cross-
ILC combinations. The S/N is computed as $\sigma(A_{\rm kSZ})/A_{\rm kSZ}$ and
we marginalize over the 6 $\Lambda{\rm CDM}$ parameters and $A_{\rm CIB+tSZ}$.
Parameters governing the radio residuals, $\alpha_{\rm rad}$ and
$\sigma(\alpha_{\rm rad})$, are fixed in this case but we also discuss below
the impact on kSZ S/N when they are left free. We apply a Planck-like prior
$\sigma(\tau_{\rm re})=0.007$ and $\sigma(A_{\rm CIB+tSZ})=0.1$ in this case.
Like mentioned in §2.4.1, we use $\ell\in[30,5000]$ for T and P but modify
$\ell_{\rm max}^{\rm TT,kSZ}$ as shown in the table.

As expected, the kSZ S/N increases when we include information from small-
scales. This is because of the sample variance of the CMB which decreases
exponentially when moving towards smaller scale due to diffusion damping.
While the MV-ILC generally returns a better S/N compared to cross-ILC, having
a precise knowledge of CIB and tSZ residuals is important for the MV-ILC (see
Fig. 2 and Fig. 3 to compare the residual levels to kSZ). We demonstrate this
in Appendix B by calculating the biases in $A_{\rm kSZ}$ estimation due to
unmodeled CIB and tSZ signals. We can also note from the table, that the S/N
saturates or even decreases for the MV-ILC when moving from $\ell_{\rm
max}^{\rm TT,kSZ}=4000$ to $\ell_{\rm max}^{\rm TT,kSZ}=4500$. This is because
of the degeneracy between $A_{\rm kSZ}$ vs $A_{\rm CIB+tSZ}$ and the increase
in the small-scale CIB+tSZ residuals. The case is different for cross-ILC and
we note that the S/N constantly increases with $\ell_{\rm max}^{\rm TT,kSZ}$.
The choice of $\sigma(A_{\rm CIB+tSZ})=0.1$ prior does not affect the cross-
ILC but has non-negligible impact on the results from MV-ILC.

In the rest of the paper, we only focus on the constraints using the cross-ILC
technique. For $\ell_{\rm max}^{\rm TT,kSZ}=4000\ (4500)$, the kSZ power
spectrum can be detected at $\sim 19\sigma$ by SPT-3G. Similarly, SO can
measure the kSZ power spectrum at $13-15\sigma$. The advantage of including HF
bands for SPT-3G and SO is also evident from the table. Including SPT-3G+ will
improve the S/N by $\sim\times 2$ and this could be achieved during the later
part of this decade. For SO, adding FYST improves the S/N by $\sim\times 1.5$.
For S4-Wide, the kSZ SNR is $\sim 70\ (80)$ for $\ell_{\rm max}^{\rm
TT,kSZ}=4000\ (4500)$.

_Fitting for radio:_ When we include $\alpha_{\rm rad}$ and
$\sigma(\alpha_{\rm rad})$, we see roughly 10% reduction in kSZ S/N for
$\ell_{\rm max}^{\rm TT,kSZ}=4000$. For, $\ell_{\rm max}^{\rm TT,kSZ}=4500$,
the reduction in kSZ S/N is higher (30-50%) compared to the baseline case when
radio residual parameters are fixed. This is because of the degeneracy between
$A_{\rm kSZ}$ vs radio residuals when restricting the fitting to the same
multipole ranges. If we include information from even smaller scales, the kSZ
signal and radio residuals have significantly different shapes (see Fig. 5)
which breaks the degeneracy between the parameters. For example, when we set
$\ell_{\rm max}=6000$, the reduction in kSZ S/N when fitting for radio
residuals is lower only by 20% compared to fixing the radio residual model
parameters. The other way of mitigating this is by reducing the masking
threshold or by nulling radio along with tSZ in one of the legs of the cross-
ILC as described in § 3.1.2 and demonstrated in Fig. 5. The choice of the
priors used for radio residuals (see Table 3) does not have a significant
impact on the kSZ S/N. For example, increasing the prior widths by $\times 2$
changes the kSZ S/N only marginally.

_Impact due to marginalization of $\Lambda{\rm CDM}$ parameters:_ Assuming
Planck priors on $\Lambda{\rm CDM}$ parameters has no impact on the kSZ S/N
for S4-Wide. If all the $\Lambda{\rm CDM}$ parameters are fixed, the kSZ S/N
increases kSZ by $20\%$. This is because the cosmological constraints are
primarily driven by TE/EE for the future CMB surveys (Galli et al., 2014)
while kSZ constraints are TT-only. Alternatively, if we calculate the
constraints with TT-only, then we find that applying Planck priors on
$\Lambda{\rm CDM}$ parameters improves the kSZ S/N by 20% while fixing them
enhances the S/N by $\times 2$.

###  3.2 Constraints on the epoch of reionization

Figure 6: Marginalized constraints on reionization parameters $z_{\rm re}^{\rm
mid}$ and $\Delta z_{\rm re}$ from kSZ power spectrum. The left and right
panels correspond to $\sigma(\tau_{\rm re})=0.007$ (Planck-like) and
$\sigma(\tau_{\rm re})=0.002$ (LiteBIRD-like) but translated to priors on the
mid-point of reionization $z_{\rm re}^{\rm mid}$. kSZ power spectrum does not
constrain $z_{\rm re}^{\rm mid}$ and hence it is dominated by the prior for
all experiments. The future surveys can reduce the uncertainty on the duration
of reionization $\sigma(\Delta z_{\rm re})\sim 1.5$ (left panel) which is
currently unconstrained by Planck. The constraints on $\sigma(\Delta z_{\rm
re})$ improves by a $\times 2.5-\times 3$ (right panel) when applying the
LiteBIRD-like $\tau_{\rm re}$ prior.

In this section, we compute the constraints on reionization parameters from
the kSZ power spectrum using the cross-ILC technique. Like mentioned in §
2.4.1, we now replace the total kSZ amplitude $A_{\rm kSZ}$ parameter with
four other parameters namely: $z_{\rm re}^{\rm mid},\Delta z_{\rm re},A_{\rm
kSZ}^{h},\alpha_{\rm kSZ}^{h}$ (see Table 3). We use the standard priors
listed in Table 3 for the foreground parameters: $A_{\rm CIB+tSZ},\alpha_{\rm
rad},\sigma(\alpha_{\rm rad})$.

To this end, we use Abundance Matching Box for the Epoch of Reionization
(AMBER) simulations (Trac et al., 2022; Chen et al., 2023) to predict the
patchy kSZ signal using two parameters: the mid-point $z_{\rm re}^{\rm mid}$
and duration $\Delta z_{\rm re}$ of reionization. The fiducial value for
$z_{\rm re}^{\rm mid}=7.69$, corresponding to $\tau_{\rm re}=0.0544$ (Planck
Collaboration et al., 2020), and we set $\Delta z_{\rm re}=4$ (Chen et al.,
2023). For the local homogeneous kSZ signal, we follow (Alvarez et al., 2021)
and model it as $D_{\ell}^{\rm h-kSZ}=A_{\rm
kSZ}^{h}\left(\dfrac{\ell}{\ell_{\rm\ast}}\right)^{\alpha_{\rm
kSZ}^{h}}\mu{\rm K}^{2}$ with $A_{\rm kSZ}^{h}=1.5$, $\alpha_{\rm kSZ}^{h}=0$,
and $\ell_{\rm\ast}=3000$.

Since it is hard to differentiate between the homogeneous and patchy kSZ
signals from the total kSZ power spectrum, we set priors on them:
$\sigma(A_{\rm kSZ}^{h})=0.1$ and $\sigma(\alpha_{\rm kSZ}^{h})=0.1$ (The
Simons Observatory Collaboration et al., 2018). There are a couple of other
techniques of separating the kSZ contributions from the homogeneous and patchy
components. First is the “dkSZ-ing” (Foreman et al., 2022) approach of
predicting and subtracting the homogeneous kSZ using large-scale structure
(LSS) surveys. This is possible since the homogeneous kSZ is sourced by haloes
in the local Universe which will he highly correlated with the LSS tracers.
While the authors primarily focused on improving the constraints on
cosmological parameters using dKSZ-ing, this technique should also be useful
to extract the reionization kSZ signal, as alluded by authors in the paper.
The second approach is to use kSZ 4-pt information (Smith & Ferraro, 2017;
Alvarez et al., 2021). These are, however, outside the scope of this work and
we leave these methods to be explored in a future work.

In Fig. 6, we present the constraints on $z_{\rm re}^{\rm mid}$ and $\Delta
z_{\rm re}$ after marginalizing over 6 $\Lambda{\rm CDM}$, 3 foreground, and 2
homogeneous kSZ parameters. We use Planck-like $\sigma(\tau_{\rm re})=0.007$
prior in the left panel and LiteBIRD-like $\sigma(\tau_{\rm re})=0.002$ prior
in the right panel. The $\sigma(\tau_{\rm re})$ priors are translated to
priors on $\sigma(z_{\rm re}^{\rm mid})=1.43\ (0.41)$ in the left (right)
panel. Since the kSZ power spectrum does not constrain $\tau_{\rm re}$ or
alternatively $z_{\rm re}^{\rm mid}$, they are dominated by the prior for all
experiments. However, we note in the left panel that SPT-3G and SO can achieve
$\sigma(\Delta z_{\rm re})=2$. Including HF information from SPT-3G+ and FYST,
can reduce the uncertainty to $\sigma(\Delta z_{\rm re})=1.6$ while CMB-S4 can
reduce it further to $\sigma(\Delta z_{\rm re})=1.4$. Choosing the LiteBIRD-
like prior can improve these by $\times 2.5-\times 3$. Widening the width of
the priors $\sigma(A_{\rm kSZ}^{h})$ and $\sigma(\alpha_{\rm kSZ}^{h})$ by
$\times 2$ ($\times 5$) increases $\sigma(\Delta z_{\rm re})$ by 10% (30%).

###  3.3 Temperature-based CMB lensing

Figure 7:  The observed foreground-induced lensing biases due to tSZ, kSZ, and
CIB when cross-correlating CMB lensing convergence $\kappa$ map from S4-Wide
with galaxy overdensity $\delta_{g}$ and shear $\gamma$ from LSST-Y1. The
fiducial maximum CMB multipole is set to $\ell_{\rm max}^{T}=3000$ except for
red the curves where we set $\ell_{\rm max}^{T}=5000$. The standard QE (Hu &
Okamoto, 2002) (blue) is significantly affected by tSZ (left panel) while the
other estimators are immune to tSZ. The bias form kSZ (middle panel) is the
same for all estimators with $\ell_{\rm max}^{T}=3000$ and small at $\lesssim
0.1\sigma$ ($\sim 1-3\%$) level. The bias from CIB follows this order: MH18 QE
(green, MH18) > standard QE (blue, Hu & Okamoto 2002) > cross-ILC QE (yellow),
but note that the CIB biases are boosted by $\times 10$ for clarity. For all
foregrounds, the bias from the cross-ILC is the lowest and hence we also
quantify the bias for $\ell_{\rm max}=5000$ as well (red curve). As it is
evident from the red curve, the foreground-induced biases are at a negligible
level for cross-ILC QE even for $\ell_{\rm max}^{T}=5000$.

MH18 demonstrated that it is possible to produce temperature-based lensing
maps that are immune to biases from tSZ with minimal noise penalty by
utilizing two different maps: one high-resolution low-noise map (either the
lowest noise frequency channel or the minimum-variance map) and a tSZ-free
map. One concern with such an approach is the boosted amplitude of CIB in the
tSZ-free map (see red curve in Fig. 3), which could be picked up strongly if
the other input temperature map also contains CIB residuals.

In this study, we build on the general concept of MH18 and replace the low-
noise/minimum variance temperature map with a CIB-free temperature map. This
allows us to reduce biases from both tSZ and CIB simultaneously, allowing us
to produce a clean CMB lensing map with minimal contamination from
extragalactic foregrounds.

We follow the approach used in the original Agora simulations work (Omori,
2022) and estimate the residual bias from various secondary components by
performing lensing reconstruction on single component maps (i.e., individual
temperature maps of tSZ/kSZ/CIB) processed in the same manner as the full
lensing reconstruction (assuming the presence of all the foregrounds). Rather
than showing the auto-spectrum of the reconstructed single-component lensing
maps, we cross-correlate the reconstructed lensing maps with LSS tracers
(galaxy overdensity and galaxy shear) since foreground biases get picked up
more strongly in such cross-correlation measurements (Fabbian et al., 2019).
Specifically, we use Rubin Observatory Legacy Survey of Space and Time (LSST-
like, LSST Science Collaboration et al. 2009) sample for this test.

The results are shown in Fig. 7 where we compare biases from standard QE
(blue), MH18 QE (green), and the cross-ILC QE (yellow and red). We show biases
relative to the error $\dfrac{(\hat{C}_{\ell}^{AB}-C_{\ell}^{AB})}{\Delta
C_{\ell}^{AB}}$ and also the differential bias
$\dfrac{\hat{C}_{\ell}^{AB}}{C_{\ell}^{AB}}-1$. The errors are calculated
analytically using Eq. 5 where $A$ corresponds to CMB lensing $\kappa$ map
from S4-Wide and $B$ is either the galaxy overdensity $\delta_{g}$ or shear
$\gamma$ from LSST year 1 (LSST-Y1) sample. We use CMB $\ell_{\rm
max}^{T}=3000$ in all cases except for the red curve for which we include
information up to $\ell_{\rm max}^{T}=5000$. We perform this cross-correlation
measurement in multiple redshift bins for different experiments although in
Fig. 7 we only show the results for the redshift bin $z\in[0.56,0.82]$ for
S4-Wide for brevity. The results for other experiments and all the redshift
bins are presented in Fig. 12 in Appendix D.

From Fig. 7, we note that the standard QE (blue) shows significant biases due
to tSZ (left panel) as expected. The other curves are immune to tSZ bias. Note
that the tSZ biases are non-zero due to the finite band passes of the
experiments. See Appendix D for more details. The bias from kSZ (boosted by
$\times 2$ for clarity) is also small $\lesssim 0.1\sigma$ ($\sim 1-3\%$) for
the lensing estimators. The bias from CIB is higher for MH18 QE (green) due to
the enhanced level of CIB in the tSZ-free map. However, note that the CIB bias
is boosted by $\times 10$. For cross-ILC QE, the bias from CIB is also
negligible.

##  4 Conclusion

We presented a cross-ILC approach to robustly extract the kSZ and CMB lensing
signals from the current (SPT-3G) and future (SO and CMB-S4) CMB surveys. The
approach uses the cross-spectrum between CIB-free and tSZ-free cILC maps for
kSZ measurement. For CMB lensing, we pass the CIB-free and tSZ-free cILC maps
in the two legs of the QE.

We showed that the residual bias from CIB and tSZ is minimal for this approach
and $\times 3$ to $\times 5$ lower compared to the expected kSZ signal level
of $D_{\ell}=3\ \mu{\rm K}^{2}$. The CIB and tSZ residuals in this cross-ILC
combination are also significantly lower than the measurements using auto-
spectrum of either the MV or the cILC map. We also quantified the residual
radio signals and demonstrated that they can be modeled or mitigated easily.
We have ignored galactic foregrounds in this work as they are negligible on
small scales compared to the extragalactic foreground signals. Similarly we
have also ignored CO in this work as they are also much smaller compared to
kSZ and other foreground signals (Maniyar et al., 2023). For CMB lensing,
using cross-correlation with galaxy overdensity and shear fields from LSST-Y1
sample, we showed that the modified cross-ILC QE has negligible level of CIB
and tSZ-induced biases.

Using Fisher formalism, we forecasted the expected kSZ S/N for CMB surveys.
For the cross-ILC technique, we showed that the total kSZ power spectrum can
be detected at $\sim 19\sigma$ by SPT-3G and the inclusion of data from its
successor SPT-3G+, improves the S/N by $\times 2$. The expected kSZ S/N for SO
is $\sim 13-15\sigma$ with roughly $\times 1.5$ expected after the inclusion
of information from FYST. For CMB-S4, the expected S/N of the total kSZ power
spectrum is extremely high $70-80\sigma$.

We also estimated the biases due to unmodeled CIB/tSZ residuals and found that
the kSZ measurements from the MV-ILC can be significantly biased. The biases
on the cosmological parameters are negligible, though, since the constraints
are dominated by EE/TE. For the cross-ILC measurement, biases on both kSZ and
cosmological parameters are negligible.

We also forecasted the constraints on the epoch of reionization, $z_{\rm
re}^{\rm mid}$ and $\Delta z_{\rm re}$ from kSZ power spectrum using the
cross-ILC approach. While the error on $z_{\rm re}^{\rm mid}$ is dominated by
the prior on the optical depth, the current and upcoming surveys can achieve
$\sigma(\Delta z_{\rm re})=1.4-2$. This improves by $\times 2.5-\times 3$ when
we replace the prior on optical depth from Planck-like to LiteBIRD-like.
Including the kSZ 4-pt information can improve these constraints further
(Smith & Ferraro, 2017).

Applications: The technique has several potential applications. For example,
it can also be used to mitigate the effects of dust and synchrotron signals,
both in temperature and polarization, using dust-/synchrotron-free maps in the
two legs of the cross-ILC estimator facilitating robust measurements of large-
scale CMB for inflationary B-modes. It can also be used for other higher order
statistics, for example, kSZ 4-pt (Smith & Ferraro, 2017; Ferraro & Smith,
2018; Alvarez et al., 2021), projected field kSZ (Hill et al., 2016; Kusiak et
al., 2021), tSZ bispectrum analysis (Crawford et al., 2014), by plugging in
different foreground-removed ILC maps in each of the legs. We have not shown
these explicitly here and leave them to be explored in a future work.

Data/code availability: All the data products produced in this work and the
associated codes are publicly available and can be downloaded from this
link${}^{\text{\faGithub}}$.

## Acknowledgments

We thank Gil Holder for several insightful discussions throughout the course
of this work. We thank Tom Crawford for important discussions, constructive
comments, and helpful suggestions that had significantly improved this
manuscript. We also thank Alex Van Engelen, Colin Hill, Nicholas Huang,
Abhishek Maniyar, Christian Reichardt, Nathan Whitehorn and Kimmy Wu for
discussions and feedback on the manuscript.

This work was partially supported by the Center for AstroPhysical Surveys
(CAPS) at the National Center for Supercomputing Applications (NCSA),
University of Illinois Urbana-Champaign. This work made use of the following
computing resources: Illinois Campus Cluster, a computing resource that is
operated by the Illinois Campus Cluster Program (ICCP) in conjunction with the
National Center for Supercomputing Applications (NCSA) and which is supported
by funds from the University of Illinois at Urbana-Champaign; the
computational and storage services associated with the Hoffman2 Shared Cluster
provided by UCLA Institute for Digital Research and Education’s Research
Technology Group; and the computing resources provided on Crossover, a high-
performance computing cluster operated by the Laboratory Computing Resource
Center at Argonne National Laboratory.

##  Appendix A Bandpower error comparisons

Figure 8: Bandpower errors for S4-Wide survey from different ILC techniques is
shown in the left panel. The data points include contribution from CMB, kSZ,
residual foregrounds and noise. As expected, MV returns the least error while
tSZ-free maps has the largest error. The cross-ILC combination tSZ-free
$\times$ CIB-free is $\sim\times 2$ worse (better) than MV (tSZ-free). The
right panel shows the bandpower errors expected from the cross-ILC technique
for different experiments considered in this work. We assume
$\Delta_{\ell}=250$ for this figure. From the right panel, it is evident that
the cumulative kSZ S/N for all experiments is high.

We show the bandpower errors from different ILC techniques for S4-Wide in the
left panel of Fig. 8. The cross-ILC technique is worse (better) than MV (tSZ-
free) by $\sim\times 2$ across the entire multipole range shown in the figure.
In the right panel, we show the errors for different experiments considered in
this work from the cross-ILC (tSZ-free $\times$ CIB-free) technique. We assume
$\Delta_{\ell}=250$ to compute the bandpower errors using Eq. 5. We note that
the cumulative kSZ S/N expected from all experiments is high with the cross-
ILC technique. As shown previously in §3.1.1 and Fig. 4, the foreground
residuals are lower for the cross-ILC approach making it optimal for
extracting the kSZ signal from current and future experiments.

##  Appendix B Bias due to residual CIB and tSZ signals

Figure 9: Systematic biases in kSZ power spectrum estimation $\sigma_{\rm
sys}(A_{\rm kSZ})$ due to 5% unmodeled CIB and tSZ residuals in the CMB maps.
The bias in the cross-ILC technique is low compared to MV, which is
significantly biased. As expected, the bias in MV increases when extending
$\ell_{\rm max}$ in the analysis.  Figure 10: Biases in cosmological
parameters due to residual CIB and tSZ signals for S4-Wide survey. The biases
are negligible for both MV and cross-ILC since the constraints are dominated
by EE/TE rather than TT.

In this section we calculate the biases in $A_{\rm kSZ}$ due to 5% unmodeled
residual CIB and tSZ signals in the CMB maps. We use a modified Fisher
formalism as given by Eq.(B2) based on Amara & Réfrégier (2008) for estimating
the biases. The bias $b$ on the parameter $i$ is calculated as

$\displaystyle b_{i}=F^{-1}_{ij}B_{j}$ (B1) $\displaystyle{\rm with\
}B_{j}=\sum_{\ell}C_{\ell}^{\rm sys}\Delta C_{\ell}^{-2}\frac{\partial
C_{\ell}^{XY}}{\partial{\theta_{j}}}.$ (B2)

Here $\Delta C_{\ell}$ is calculated with Eq.(5) but after including the
systematic signal as well with $\hat{C}_{\ell}=C_{\ell}+N_{\ell}+C_{\ell}^{\rm
sys}$. The results are shown in Fig. 9 for both MV (hatched bars) and cross-
ILC techniques (filled bars) as a function of $\ell_{\rm max}$ used for kSZ
extraction. We note from the figure that the MV ILC becomes significantly
biased due to the residual CIB and tSZ signals while the biases in the cross-
ILC technique are low in all cases. As expected, biases in MV ILC increases
when the $\ell_{\rm max}$ is extended since the level of CIB+tSZ residuals
also increase with $\ell_{\rm max}$ as can be seen from the blue curve of Fig.
3.

We also calculate the biases in cosmological parameters due to the CIB/tSZ
residuals and do not find significant biases. This lower level of biases is
because the cosmological constraints are dominated by EE/TE while CIB/tSZ
residuals are expected to be largely unpolarized and hence only injected for
TT only. The biases in cosmological parameters expected due to residual
foregrounds for S4-Wide are shown in Fig. 10.

##  Appendix C Validating Agora foreground model

###  C.1 SPT foreground model for ILC weights

We have used Agora simulations as our baseline foreground model. In this
section we test the robustness of this assumption. To this end, we replace the
$\rm N_{ch}\times\rm N_{ch}$ covariance matrix ${\bf C}_{\ell}$ in Eq.(2)
containing the auto- and cross-power spectra between Agora simulations in
different bands by the foreground model from SPT measurements (R21). Modifying
the covariance matrix alters the ILC weights which will indeed modify the
final residuals. We pass the Agora simulations through these modified weights
and compare the resultant residuals with our baseline approach.

Fig. 11 presents this comparison for different ILC combinations: MV (blue),
tSZ-free (red), CIB-free (orange) and the cross-ILC tSZ-free $\times$ CIB-free
(green). The solid lines correspond to ILC weights from Agora foreground model
(baseline case) and the dashed lines are for SPT foreground model. All curves
show the total CIB+tSZ foreground residuals.

Figure 11: Total CIB and tSZ residuals from Agora simulations for two
different ILC weights. Solid lines correspond to the case when Agora
foreground model is used in the covariance matrix to calculate the weights
while the dashed lines use a foreground model based on SPT measurements. The
colors correspond to different ILC combinations. The close match between the
solid and dashed curves asserts that our baseline approach of using Agora
foreground model is a valid assumption.

It is important to point out that there are some differences between SPT and
Agora foreground model used in this work. Firstly, the SPT foreground model
computed using R21 measurements does not have cluster masked and the masking
threshold for dusty point sources in $S_{150}=6.4\ {\rm mJy}$ while for Agora
model, we have masked clusters with mass $M_{\rm 500c}\geq 2\times 10^{14}\
\mbox{M}_{\odot}$ and use a reduced masking threshold for dusty sources
$S_{150}=3\ {\rm mJy}$. For radio point sources, both the models assume the
same $S_{150}=3\ {\rm mJy}$ threshold. Next, the cross-correlation spectra
between CIB and tSZ is simply modeled using a fixed cross-correlation
coefficient of $\rho_{\rm CIB\times tSZ}=0.078$ (R21) in the SPT model.
Finally, the SPT CIB model is an extrapolation of the 150 GHz power spectrum
assuming two SEDs (one for the Poisson component and the other for the
clustering component of the CIB) and does not account for the CIB
decorrelation between different frequency bands. As a result, we expect some
differences between the solid and dashed curves.

In Fig. 11, the tSZ-free (red) curve is sensitive to the differences in the
CIB and radio models between Agora and SPT while the CIB-free (orange) curve
is sensitive to the differences in the tSZ and radio models. The cross-ILC
tSZ-free $\times$ CIB-free (green) and the MV (blue) curves are affected by
all the three (CIB, radio, and tSZ) foreground modeling differences. Despite
the above mentioned differences between the two foreground models, the
agreement between the solid and dashed curves is excellent in all cases.
Specifically, the difference is negligible for cross-ILC tSZ-free $\times$
CIB-free (green). Based on this result, we claim that our baseline approach of
using Agora simulations to estimate the ILC covariance matrix is reasonable.
We also note that for future CMB surveys, the ILC covariance matrix can use
the bandpowers measured by the respective experiment to derive the ILC weights
rather than Agora.

###  C.2 Websky foreground simulations

Now we test the cross-ILC approach by replacing Agora using the multi-
component Websky555Websky simulations are downloaded from
https://mocks.cita.utoronto.ca/index.php/WebSky_Extragalactic_CMB_Mocks
simulations (Stein et al., 2020). Similar to § C.1, we limit this test to the
S4-Wide survey. We approximate $\nu\in[90,150,220,285]$ GHz with the publicly
available [93, 145, 217, 278] GHz Websky simulations and do not attempt to
interpolate the simulations to match the bands used in this work. We also note
that the Websky simulations assume a delta function for the bandpasses centred
on each frequency band (Stein et al., 2020). Even though the CIB power in
Websky simulations are shown to match Planck 545 GHz data well (Stein et al.,
2020), the CIB power at $\ell=3000$ in the 150 GHz band is $\sim\times 2$
higher (Raghunathan et al., 2022) than the SPT results reported by (George et
al., 2015; Reichardt et al., 2021). Subsequently, we multiply the Websky CIB
maps by 0.75. The applied correction is the same for all bands. We note this
simple scaling alone does not make Websky CIB to match SPT results in other
bands but it it makes the agreement better. This simple $\ell$-independent
scaling also does not work for the large-scale clustered part of the CIB
signal. No scaling is applied for Websky tSZ maps.

With the Wesbky CIB/tSZ maps and the power spectra in hand, we run the blind
search algorithm described in § 2.4 to derive the optimal weights for the
cross-ILC estimator. We note that the SEDs obtained from the blind search does
not match the SEDs obtained for Agora simulations. This is because of the
differences between Websky and Agora simulations, although we find Agora to
match data better as demonstrated in Fig. 11. Nevertheless, adopting the SEDs
from the blind search does reduce the total CIB+tSZ residuals. For example,
without the Websky CIB scaling mentioned above, the total CIB+tSZ residuals
are $\times 1.5$ lower than kSZ between $\ell\in[3000,5000]$. However, after
introducing the CIB scaling, we find the CIB+tSZ residuals to go down by
$\times 2$ at $\ell=3000$ and by $\times 4$ $\in[4000,5000]$ compared to the
expected level of the kSZ signal $D_{\ell}=3\mu{\rm K}^{2}$.

##  Appendix D Lensing biases

In Fig. 12, we show the biases in CMB lensing cross-correlations for the
cross-ILC QE. This is similar to Fig. 7 but here, along with S4-Wide, we also
include results for SPT-3G, SO-Baseline, and SO-Goal. In this figure, we only
show the bias over the error. We also present the results for five redshift
bins: Bin 1 = $z\in[0.2,0.42]$, Bin 2 = $z\in[0.38,0.6]$; Bin 3 =
$z\in[0.56,82]$; Bin 4 = $z\in[0.78,1.12]$; and Bin 5 = $z\in[1.11,1.66]$.
Note that the redshift bins are not independent and have some overlap due to
photo-$z$ errors (The LSST Dark Energy Science Collaboration et al., 2018)
although that is irrelevant in this context since we are mainly focused on the
foreground-induced biases in CMB lensing cross correlations. The biases are
all calculated using Agora simulations (Omori, 2022).

The bias due tSZ (left panel) is small $\sim 0.05\sigma$ but non-zero in some
cases. This may be a little counter intuitive since the cross-ILC lensing QE
uses the tSZ-free map in one of the legs and we should expect the bias from
tSZ to be zero. However, note that the finite band passes of the experiments
can lead to small levels of residual tSZ and that can correlate with the non-
zero tSZ signal in the CIB-free map. Nevertheless, the level of the tSZ bias
is negligible. The bias from kSZ is also at a similar level for all
experiments. The bias due to CIB is at sub-percent level and hence negligible.
Note that we have boosted the CIB biases by $\times 10$ for clarity.

Motivated by these small level of biases, in the bottom panel we present the
results for S4-Wide but after increasing the CMB $\ell_{\rm max}^{T}=5000$
during lensing reconstruction. As discussed earlier in §3.3, including
information from smaller scales increases the lensing SNR significantly. For
example, the lensing SNR with the maximum CMB multipole $\ell_{\rm
max}^{T}=5000$ is roughly $\times 1.7$ better than using $\ell_{\rm
max}^{T}=3000$ for temperature-based CMB lensing. This can also lead to
enhanced level of biases which we test here. However, as can be seen from the
bottom panel, the biases are small $\lesssim 0.1\sigma$ even in this case.
While we do not show it explicitly, we also checked the biases due to radio
sources for the fiducial masking threshold of $S_{150}$ = 6 mJy and find them
also to be at sub-percent level. Thus, we conclude that the cross-ILC lensing
QE is robust against the biases induced by CIB, kSZ, radio and tSZ signals.

Figure 12: Biases due to tSZ (left panel), kSZ (middle panel), and CIB (right
panel) in CMB lensing cross-correlations similar to Fig. 7 but extended to
other experiments and more redshift bins. See text for more details.

## References

  * Abbott et al. (2023) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2023, Phys. Rev. D, 107, 023531, doi: 10.1103/PhysRevD.107.023531
  * Alvarez et al. (2021) Alvarez, M. A., Ferraro, S., Hill, J. C., et al. 2021, Phys. Rev. D, 103, 063518, doi: 10.1103/PhysRevD.103.063518
  * Amara & Réfrégier (2008) Amara, A., & Réfrégier, A. 2008, MNRAS, 391, 228, doi: 10.1111/j.1365-2966.2008.13880.x
  * Anderson et al. (2022) Anderson, A. J., Barry, P., Bender, A. N., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12190, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy XI, ed. J. Zmuidzinas & J.-R. Gao, 1219003, doi: 10.1117/12.2629755
  * Battaglia et al. (2013a) Battaglia, N., Natarajan, A., Trac, H., Cen, R., & Loeb, A. 2013a, ApJ, 776, 83, doi: 10.1088/0004-637X/776/2/83
  * Battaglia et al. (2013b) Battaglia, N., Trac, H., Cen, R., & Loeb, A. 2013b, ApJ, 776, 81, doi: 10.1088/0004-637X/776/2/81
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Bender et al. (2018) Bender, A. N., Ade, P. A. R., Ahmed, Z., et al. 2018, in SPIE Conference Series, Vol. 10708, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX, 1070803, doi: 10.1117/12.2312426
  * Benson et al. (2014) Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in SPIE Conference Series, Vol. 9153, SPIE Conference Series, 1, doi: 10.1117/12.2057305
  * Bianchini et al. (2020) Bianchini, F., Wu, W. L. K., Ade, P. A. R., et al. 2020, ApJ, 888, 119, doi: 10.3847/1538-4357/ab6082
  * Bleem et al. (2015) Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, The Astrophysical Journal Supplement Series, 216, 27, doi: 10.1088/0067-0049/216/2/27
  * Bleem et al. (2021) Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2021, arXiv e-prints.  https://arxiv.org/abs/2102.05033
  * Bocquet et al. (2019) Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55, doi: 10.3847/1538-4357/ab1f10
  * Calabrese et al. (2014) Calabrese, E., Hložek, R., Battaglia, N., et al. 2014, J. Cosmology Astropart. Phys, 2014, 010, doi: 10.1088/1475-7516/2014/08/010
  * Cardoso et al. (2008) Cardoso, J.-F., Le Jeune, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 735, doi: 10.1109/JSTSP.2008.2005346
  * CCAT-Prime Collaboration et al. (2023) CCAT-Prime Collaboration, Aravena, M., Austermann, J. E., et al. 2023, ApJS, 264, 7, doi: 10.3847/1538-4365/ac9838
  * Chang et al. (2023) Chang, C., Omori, Y., Baxter, E. J., et al. 2023, Phys. Rev. D, 107, 023530, doi: 10.1103/PhysRevD.107.023530
  * Chen et al. (2023) Chen, N., Trac, H., Mukherjee, S., & Cen, R. 2023, ApJ, 943, 138, doi: 10.3847/1538-4357/ac8481
  * Choi et al. (2020) Choi, S. K., Hasselfield, M., Ho, S.-P. P., et al. 2020, J. Cosmology Astropart. Phys, 2020, 045, doi: 10.1088/1475-7516/2020/12/045
  * Choudhury (2022) Choudhury, T. R. 2022, General Relativity and Gravitation, 54, 102, doi: 10.1007/s10714-022-02987-4
  * CMB-S4 Collaboration (2019) CMB-S4 Collaboration. 2019, arXiv e-prints.  https://arxiv.org/abs/1907.04473
  * Crawford et al. (2014) Crawford, T. M., Schaffer, K. K., Bhattacharya, S., et al. 2014, ApJ, 784, 143, doi: 10.1088/0004-637X/784/2/143
  * Darwish et al. (2021) Darwish, O., Sherwin, B. D., Sailer, N., Schaan, E., & Ferraro, S. 2021, arXiv e-prints, arXiv:2111.00462, doi: 10.48550/arXiv.2111.00462
  * Das et al. (2011) Das, S., Sherwin, B. D., Aguirre, P., et al. 2011, Phys. Rev. Lett., 107, 021301, doi: 10.1103/PhysRevLett.107.021301
  * Dibert et al. (2022) Dibert, K. R., Anderson, A. J., Bender, A. N., et al. 2022, Phys. Rev. D, 106, 063502, doi: 10.1103/PhysRevD.106.063502
  * Douspis et al. (2021) Douspis, M., Salvati, L., Gorce, A., & Aghanim, N. 2021, arXiv e-prints, arXiv:2111.01639.  https://arxiv.org/abs/2111.01639
  * Fabbian et al. (2019) Fabbian, G., Lewis, A., & Beck, D. 2019, J. Cosmology Astropart. Phys, 2019, 057, doi: 10.1088/1475-7516/2019/10/057
  * Ferraro & Smith (2018) Ferraro, S., & Smith, K. M. 2018, Phys. Rev. D, 98, 123519, doi: 10.1103/PhysRevD.98.123519
  * Foreman et al. (2022) Foreman, S., Hotinli, S. C., Madhavacheril, M. S., van Engelen, A., & Kreisch, C. D. 2022, arXiv e-prints, arXiv:2209.03973, doi: 10.48550/arXiv.2209.03973
  * Galli et al. (2014) Galli, S., Benabed, K., Bouchet, F., et al. 2014, Phys. Rev. D, 90, 063504, doi: 10.1103/PhysRevD.90.063504
  * George et al. (2015) George, E. M., Reichardt, C. L., Aird, K. A., et al. 2015, ApJ, 799, 177, doi: 10.1088/0004-637X/799/2/177
  * Gorce et al. (2022) Gorce, A., Douspis, M., & Salvati, L. 2022, A&A, 662, A122, doi: 10.1051/0004-6361/202243351
  * Hand et al. (2012) Hand, N., Addison, G. E., Aubourg, E., et al. 2012, Physical Review Letters, 109, 041101, doi: 10.1103/PhysRevLett.109.041101
  * Hill et al. (2016) Hill, J. C., Ferraro, S., Battaglia, N., Liu, J., & Spergel, D. N. 2016, Physical Review Letters, 117, 051301, doi: 10.1103/PhysRevLett.117.051301
  * Hilton et al. (2020) Hilton, M., Sifón, C., Naess, S., et al. 2020, arXiv e-prints.  https://arxiv.org/abs/2009.11043
  * Hu & Okamoto (2002) Hu, W., & Okamoto, T. 2002, ApJ, 574, 566, doi: 10.1086/341110
  * Huang et al. (2020) Huang, N., Bleem, L. E., Stalder, B., et al. 2020, AJ, 159, 110, doi: 10.3847/1538-3881/ab6a96
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340, doi: 10.1093/mnras/stw248
  * Knox (1995) Knox, L. 1995, Phys. Rev. D, 52, 4307, doi: 10.1103/PhysRevD.52.4307
  * Kusiak et al. (2021) Kusiak, A., Bolliet, B., Ferraro, S., Hill, J. C., & Krolewski, A. 2021, Phys. Rev. D, 104, 043518, doi: 10.1103/PhysRevD.104.043518
  * Kusiak et al. (2023) Kusiak, A., Surrao, K. M., & Hill, J. C. 2023, arXiv e-prints, arXiv:2303.08121, doi: 10.48550/arXiv.2303.08121
  * Lagache et al. (2020) Lagache, G., Béthermin, M., Montier, L., Serra, P., & Tucci, M. 2020, A&A, 642, A232, doi: 10.1051/0004-6361/201937147
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, ArXiv e-prints.  https://arxiv.org/abs/0912.0201
  * MacCrann et al. (2023) MacCrann, N., Sherwin, B. D., Qu, F. J., et al. 2023, arXiv e-prints, arXiv:2304.05196.  https://arxiv.org/abs/2304.05196
  * Madhavacheril & Hill (2018) Madhavacheril, M. S., & Hill, J. C. 2018, Phys. Rev. D, 023534, doi: 10.1103/PhysRevD.98.023534
  * Madhavacheril et al. (2020) Madhavacheril, M. S., Hill, J. C., Næss, S., et al. 2020, Phys. Rev. D, 102, 023534, doi: 10.1103/PhysRevD.102.023534
  * Madhavacheril et al. (2023) Madhavacheril, M. S., Qu, F. J., Sherwin, B. D., et al. 2023, arXiv e-prints, arXiv:2304.05203.  https://arxiv.org/abs/2304.05203
  * Maniyar et al. (2023) Maniyar, A. S., Gkogkou, A., Coulton, W. R., et al. 2023, arXiv e-prints, arXiv:2301.10764, doi: 10.48550/arXiv.2301.10764
  * McCarthy et al. (2017) McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936, doi: 10.1093/mnras/stw2792
  * Mead et al. (2020) Mead, A. J., Tröster, T., Heymans, C., Van Waerbeke, L., & McCarthy, I. G. 2020, A&A, 641, A130, doi: 10.1051/0004-6361/202038308
  * Namikawa et al. (2013) Namikawa, T., Hanson, D., & Takahashi, R. 2013, MNRAS, 431, 609, doi: 10.1093/mnras/stt195
  * Omori (2022) Omori, Y. 2022, arXiv e-prints, arXiv:2212.07420, doi: 10.48550/arXiv.2212.07420
  * Omori et al. (2017) Omori, Y., Chown, R., Simard, G., et al. 2017, ApJ, 849, 124, doi: 10.3847/1538-4357/aa8d1d
  * Omori et al. (2023) Omori, Y., Baxter, E. J., Chang, C., et al. 2023, Phys. Rev. D, 107, 023529, doi: 10.1103/PhysRevD.107.023529
  * Osborne et al. (2014) Osborne, S. J., Hanson, D., & Doré, O. 2014, J. Cosmology Astropart. Phys, 2014, 024, doi: 10.1088/1475-7516/2014/03/024
  * Paul et al. (2021) Paul, S., Mukherjee, S., & Choudhury, T. R. 2021, MNRAS, 500, 232, doi: 10.1093/mnras/staa3221
  * Planck Collaboration et al. (2014a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014a, A&A, 571, A17, doi: 10.1051/0004-6361/201321543
  * Planck Collaboration et al. (2014b) —. 2014b, A&A, 571, A12, doi: 10.1051/0004-6361/201321580
  * Planck Collaboration et al. (2016a) —. 2016a, A&A, 594, A24, doi: 10.1051/0004-6361/201525833
  * Planck Collaboration et al. (2016b) —. 2016b, A&A, 594, A27, doi: 10.1051/0004-6361/201525823
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Qu et al. (2023) Qu, F. J., Sherwin, B. D., Madhavacheril, M. S., et al. 2023, arXiv e-prints, arXiv:2304.05202.  https://arxiv.org/abs/2304.05202
  * Raghunathan (2022) Raghunathan, S. 2022, ApJ, 928, 16, doi: 10.3847/1538-4357/ac510f
  * Raghunathan et al. (2022) Raghunathan, S., Whitehorn, N., Alvarez, M. A., et al. 2022, ApJ, 926, 172, doi: 10.3847/1538-4357/ac4712
  * Reichardt et al. (2021) Reichardt, C. L., Patil, S., Ade, P. A. R., et al. 2021, ApJ, 908, 199, doi: 10.3847/1538-4357/abd407
  * Remazeilles et al. (2011) Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 410, 2481, doi: 10.1111/j.1365-2966.2010.17624.x
  * Sánchez et al. (2022) Sánchez, J., Omori, Y., Chang, C., et al. 2022, arXiv e-prints, arXiv:2210.08633, doi: 10.48550/arXiv.2210.08633
  * Schaan & Ferraro (2019) Schaan, E., & Ferraro, S. 2019, Phys. Rev. Lett., 122, 181301, doi: 10.1103/PhysRevLett.122.181301
  * Schmittfull & Seljak (2018) Schmittfull, M., & Seljak, U. 2018, Phys. Rev. D, 97, 123540, doi: 10.1103/PhysRevD.97.123540
  * Sehgal et al. (2010) Sehgal, N., Bode, P., Das, S., et al. 2010, ApJ, 709, 920, doi: 10.1088/0004-637X/709/2/920
  * Sehgal et al. (2019) Sehgal, N., Aiola, S., Akrami, Y., et al. 2019, in BAAS, Vol. 51, 6.  https://arxiv.org/abs/1906.10134
  * Shaw et al. (2010) Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452, doi: 10.1088/0004-637X/725/2/1452
  * Shaw et al. (2012) Shaw, L. D., Rudd, D. H., & Nagai, D. 2012, ApJ, 756, 15, doi: 10.1088/0004-637X/756/1/15
  * Smith & Ferraro (2017) Smith, K. M., & Ferraro, S. 2017, Phys. Rev. Lett., 119, 021301, doi: 10.1103/PhysRevLett.119.021301
  * Smith et al. (2007) Smith, K. M., Zahn, O., & Doré, O. 2007, Phys. Rev. D, 76, 043510, doi: 10.1103/PhysRevD.76.043510
  * Soergel et al. (2016) Soergel, B., Flender, S., Story, K. T., et al. 2016, MNRAS, 461, 3172, doi: 10.1093/mnras/stw1455
  * Stein et al. (2020) Stein, G., Alvarez, M. A., Bond, J. R., van Engelen, A., & Battaglia, N. 2020, J. Cosmology Astropart. Phys, 2020, 012, doi: 10.1088/1475-7516/2020/10/012
  * Story et al. (2015) Story, K. T., Hanson, D., Ade, P. A. R., et al. 2015, ApJ, 810, 50, doi: 10.1088/0004-637X/810/1/50
  * Sunyaev & Zel’dovich (1970) Sunyaev, R. A., & Zel’dovich, Y. B. 1970, Comments on Astrophysics and Space Physics, 2, 66
  * Sunyaev & Zel’dovich (1972) —. 1972, Comments on Astrophysics and Space Physics, 4, 173
  * Sunyaev & Zeldovich (1980) Sunyaev, R. A., & Zeldovich, Y. B. 1980, MNRAS, 190, 413, doi: 10.1093/mnras/190.3.413
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration, Mandelbaum, R., Eifler, T., et al. 2018, arXiv e-prints, arXiv:1809.01669, doi: 10.48550/arXiv.1809.01669
  * The Simons Observatory Collaboration et al. (2018) The Simons Observatory Collaboration, Ade, P., Aguirre, J., et al. 2018, ArXiv e-prints.  https://arxiv.org/abs/1808.07445
  * Trac et al. (2022) Trac, H., Chen, N., Holst, I., Alvarez, M. A., & Cen, R. 2022, ApJ, 927, 186, doi: 10.3847/1538-4357/ac5116
  * van Engelen et al. (2014) van Engelen, A., Bhattacharya, S., Sehgal, N., et al. 2014, ApJ, 786, 13, doi: 10.1088/0004-637X/786/1/13
  * van Engelen et al. (2012) van Engelen, A., Keisler, R., Zahn, O., et al. 2012, ApJ, 756, 142, doi: 10.1088/0004-637X/756/2/142
  * Wu et al. (2019) Wu, W. L. K., Mocanu, L. M., Ade, P. A. R., et al. 2019, ApJ, 884, 70, doi: 10.3847/1538-4357/ab4186