11institutetext: Université Paris-Saclay, CEA, Département de Physique des
Particules, 91191, Gif-sur-Yvette, France.
11email: alexandre.huchet<EMAIL_ADDRESS>
# A measurement of cluster masses using
Planck and SPT-SZ CMB lensing

Alexandre Huchet  Jean-Baptiste Melin

(Received ; accepted )

We used an unbiased CMB lensing mass estimator on 468 SPT-SZ clusters from the
SPT-SZ and the Planck public data, the first such estimation using combined
ground- and space-based data. We measured the average ratio between CMB
lensing and SZ mass to be $M_{\rm CMBlens}/M_{\rm SZ}=0.98\pm
0.19\rm{\leavevmode\nobreak\ (stat.)}\pm 0.03\rm{\leavevmode\nobreak\
(syst.)}$. The average CMB lensing mass from the combination of the two data
sets is measured at 4.8$\sigma$, which is a significant gain with respect to
the measurement performed on the SPT-SZ only (3.9$\sigma$) or the Planck only
(3.7$\sigma$) data set. We showed that the combination not only takes
advantage of the two different ranges of spatial scales (i.e. Fourier modes)
observed but also exploits the lensing induced correlation between scales
observed by one experiment and the other. This result demonstrates the
importance of measuring a large range of spatial scales for CMB lensing mass
estimation, from arcmin to degrees. This large range of scales will most
probably be provided by the combination of various data sets, such as from the
large and small aperture telescopes of the upcoming Simons Observatory and
future CMB-S4 experiment, and Planck. In this context, the Planck data will
remain a key element for CMB lensing cluster studies in the years to come.

###### Key Words.:

galaxies: clusters: general – cosmic background radiation – Gravitational
lensing: weak – methods: statistical

##  1 Introduction

Clusters of galaxies, located at the nodes of the cosmic web, are very useful
objects for studying both astrophysics and cosmology. They have played an
important role in shaping and reinforcing the $\Lambda$CDM concordance model
and have been used in many works to derive constraints on cosmological
parameters (e.g., White et al. 1993; Carlberg et al. 1996; Bahcall et al.
1997; Allen et al. 2002). Among these approaches, the most promising way to
constrain cosmological parameters with clusters is through cluster counts as a
function of redshift and mass, which are very sensitive to cosmological
parameters, particularly $\Omega_{\rm m}$ and $\sigma_{8}$ (Albrecht et al.
2006; Vikhlinin et al. 2009; Allen et al. 2011). Clusters thus have the
statistical power to be as efficient as other major cosmological probes but
their constraints unfortunately suffer from uncertainty in the relation
between cluster masses and directly observable quantities. This uncertainty on
the cluster mass scale dominates current analyses (e.g., Hasselfield et al.
2013; Planck Collaboration et al. 2016c; Bocquet et al. 2019; Abbott et al.
2020). However, the situation will change in the next few years with the
advent of the large optical facilities _Euclid_
111https://www.esa.int/Science_Exploration/Space_Science/Euclid &
https://www.euclid-ec.org and the Vera Rubin
Observatory222https://rubinobservatory.org & https://lsstdesc.org. These
facilities will improve the accuracy and number of measurements of cluster
masses thanks to weak gravitational lensing, up to redshifts of the order of
one (Laureijs et al. 2011; LSST Science Collaboration et al. 2009). The major
uncertainty for cosmology with clusters at $z<1$ will, therefore, be overcome
in the next decade. Extending cosmological studies with clusters at redshift
$z>1$ will not be possible with the weak gravitational lensing on galaxies due
to the lack of background sources at these redshifts. The cosmic microwave
background (CMB) is hoped to be a new background source for studying the weak
lensing of clusters at higher redshift. Located at $z\sim 1100$ and having
precisely known statistical properties, it makes it possible to measure masses
of galaxy clusters by weak gravitational lensing in the redshift range
$0<z<3$. (Zaldarriaga & Seljak 1999; Seljak & Zaldarriaga 2000; Holder &
Kosowsky 2004; Dodelson 2004; Vale et al. 2004; Lewis & King 2006; Lewis &
Challinor 2006).

The first tools were proposed in the mid-2000s to detect this effect (Maturi
et al. 2005; Hu et al. 2007; Yoo & Zaldarriaga 2008; Yoo et al. 2010), but the
first data sets for which the lensing signal became detectable only appeared
in the early 2010s (Ruhl et al. 2004; Swetz et al. 2011; Planck Collaboration
et al. 2020a). The signal-to-noise for CMB lensing mass measurement is very
low in these data sets. Thus, it is not possible to make individual
measurements and it is necessary to average the signal over several hundred
clusters to pull the signal out of the noise. The first measurements were
carried out almost jointly by the ACT, SPT and Planck collaborations shortly
before the mid-2010s (Madhavacheril et al. 2015; Baxter et al. 2015; Planck
Collaboration et al. 2016c). These detections were made possible thanks to the
pioneering work cited above and improved tools (see Melin & Bartlett 2015).
The improvement of existing tools and the development of new tools are
currently active research fields (Raghunathan et al. 2017; Madhavacheril &
Hill 2018; Raghunathan et al. 2019a; Horowitz et al. 2019; Patil et al. 2020;
Gupta & Reichardt 2021; Levy et al. 2023; Chan et al. 2023). The motivation is
to be ready to analyze data from the new generation of CMB instruments which
will provide data sets that will enable individual mass measurements for the
first time, the Simons Observatory (SO, Ade et al. 2019) and CMB-S4 (Abazajian
et al. 2019). Following its first detection in the mid-2010s, the weak
gravitational lensing on the CMB has been used to measure masses of clusters
selected in the optical and the infrared (Geach & Peacock 2017; Baxter et al.
2018; Raghunathan et al. 2019b; Madhavacheril et al. 2020) as well as to
measure the mass of halos hosted by galaxies or quasars (Raghunathan et al.
2018; Geach et al. 2019). These measurements are performed using the intensity
(i.e. temperature) maps of current data sets, but the polarized signal is
expected to improve significantly the precision of the measurement in future
data sets. In particular, Raghunathan et al. (2019) recently succeeded in
making the first detection of the CMB cluster lensing in the SPTpol data using
only the polarized signal. In parallel to these efforts, first cluster
cosmological analyses using CMB cluster lensing to anchor the mass scale were
performed (Planck Collaboration et al. 2016c; Zubeldia & Challinor 2019),
including a careful assessment of the statistics of the recovered CMB lensing
masses (Zubeldia & Challinor 2020). With the upcoming experiments (SO and
CMB-S4), it will be possible to extract, from the same data set, galaxy
clusters and their masses to high accuracy thanks to CMB cluster lensing,
leading to competitive cosmological constraints (Louis & Alonso 2017;
Madhavacheril et al. 2017; Raghunathan et al. 2022; Chaubal et al. 2022).

In this work, we used an unbiased CMB lensing mass estimator on 468 clusters
from the SPT-SZ cluster catalogue (Bocquet et al. 2019). We worked on flat
tangential maps cut from the all-sky Planck and SPT-SZ maps, and combined
them. We took advantage of the different ranges of spatial scales, or Fourier
modes of the maps, that the two data sets provide. The combined signal-to-
noise ratio for individual mass measurement is comprised between 0.05 and 1.
We thus had to average the 468 mass estimations to obtain the average mass of
the whole sample.

We first introduce the two data sets used in our analysis in Sect. 2. We then
explain our methodology in Sect. 3 and present our simulations in Sect. 4. We
show our results in Sect. 5. We discuss uncertainties in Sect. 6. Finally, we
summarize and conclude in Sect. 7.

Throughout the article, we work in the flat $\Lambda$CDM model with parameters
$H_{0}=70$ km/s/Mpc, $\Omega_{\rm m}=0.3$.

##  2 Data sets

We used two sets of sky maps in our analysis, SPT-SZ and Planck. Our analysis
is performed on the SPT-SZ cluster catalogue only.

###  2.1 SPT-SZ sky maps

The South Pole Telescope (SPT, Ruhl et al. 2004) is a 10 m diameter mirror
facility operating at the Amundsen–Scott South Pole Station in Antartica. It
observed a $\sim$2500 deg2 area of the southern sky between 2008 and 2011, in
three frequency bands centered around 95, 150 and 220 GHz, with a beam size of
1.6, 1.1 and 1.0 arcmin respectively (full width at half maximum, FWHM).This
survey is referred as the SPT-SZ survey (Story et al. 2013). SZ stands for the
Sunyaev-Zel’dovich effect, used to detect the clusters in this survey. The SZ
effect is actually a contaminant in our study, we present it and explain how
we handle it in Sec. 3.4.

We used the ”SPT Only Data Maps” published by Chown et al. (2018). The maps
include all three SPT-SZ frequencies, with a degraded resolution compared to
the original data, that is 1.75 arcmin (FWHM) for the three maps. The SPT-SZ
data contain a point source mask and the transfer function of the instrument
for the three frequencies. All maps are provided in equatorial coordinates,
and in HEALPix format with the resolution parameter $N_{side}=8192$,
corresponding to a pixel size of about 0.43 arcmin. We also extracted the
frequency response of each of the SPT-SZ band from Fig. 10 of Chown et al.
(2018), which are useful to compute the frequency dependance of the thermal SZ
effect in each frequency band $j_{\nu_{i}}$ (see Sect. 3.1).

###  2.2 Planck sky maps

The Planck space mission (Planck Collaboration et al. 2020b) is equipped with
a 1.5 m diameter telescope and two instruments, the Low Frequency Instrument
(LFI) and the High Frequency Instrument (HFI). The latter observed the whole
sky in six frequency bands centered at 100, 143, 217, 353, 545 and 857 GHz at
a resolution of 9.659, 7.220, 4.900, 4.916, 4.675, 4.216 arcmin (FWHM) (Planck
Collaboration et al. 2016a).

Planck HFI maps are full sky but we restrained our analysis to the SPT-SZ
field in order to combine the maps.

The Planck maps are provided in Galactic coordinates, at a HEALPix resolution
$N_{side}=2048$, corresponding to a pixel size of 1.72 arcmin. To facilitate
combining these maps with the SPT-SZ maps, we changed their coordinate system
to equatorial, and upgraded them to $N_{side}=8192$ by zero-padding: we
harmonic transform the maps, add zeros at high multipoles in the harmonic
space, and then perform the inverse harmonic transform. We used the maps from
the Planck Release 2 provided in the Planck Legacy Archive.

###  2.3 Complementarity of the two data sets

The two data sets are complementary. On the one hand, the Planck telescope is
located in space, and its HFI mapped the sky at six different frequencies, but
Planck has somewhat a small mirror for cluster science. On the other hand, the
South Pole Telescope is ground-based, has a large mirror, a large number of
detectors, but the large scales in the maps ($>0.5{\rm deg}$) are filtered out
due to the telescope scanning strategy.

As a result, SPT-SZ has a very good resolution and sensitivity whereas Planck
has an excellent frequency coverage and can map the large angular scales of
the sky. This complementarity should bring information on correlations between
large (Planck) and small (SPT-SZ) scales, leading to a significant improvement
on the final error of our combined CMB lensing analysis compared to the
analyses led on the two independent data sets.

###  2.4 SPT-SZ cluster catalogue

We used the catalogue from Bleem et al. (2015), with masses and redshifts
updated by Bocquet et al. (2019). The catalogue contains 677 objects. We
restricted our analysis to the 516 clusters with measured redshifts, and
subsequently removed 48 clusters located close to the edge of the SPT-SZ field
(see Sect. 3.2). Our final catalogue contains 468 clusters. The SPT-SZ masses
$M_{500}$ (mass enclosed in a sphere of radius $R_{500}$ within which the
average density is 500 times the critical density at the cluster’s redshift)
are estimated using the SZ signal-to-noise ratio versus mass relation for a
flat $\Lambda{\rm CDM}$ model ($h=H_{0}/(100\,{\rm km/s/Mpc})=0.7$,
$\Omega_{\rm m}$=0.3, $\sigma_{8}=0.8$), which is calibrated from external
X-ray and weak lensing data sets. The SPT-SZ masses are corrected from
selection effects. We multiplied the SPT mass by a factor 0.8 to match the
Arnaud et al. (2010) mass definition as was done by Tarrío et al. (2019),
Melin et al. (2021) and Melin & Pratt (2023). We denote by $M_{\rm SZ}$ the
SPT mass rescaled by this factor 0.8. We discuss this choice in the light of
our results in Sect. 7. The distribution of our sample in the mass versus
redshift plane is shown in Fig. 1. Clusters span a large redshift range
$z\in[0,1.7]$ above a mass threshold that varies slowly with redshift $M_{\rm
SZ}\gtrsim 2\times 10^{14}M_{\sun}$.

Figure 1: Mass-redshift distribution of the 468 SPT-SZ clusters used in our
analysis. We display in blue the error bars for the clusters with $M_{\rm
SZ}>7\times 10^{14}M_{\sun}$ and of one out of 20 under this threshold for
clarity.

##  3 Methodology

We use some of the notations from Melin & Bartlett (2015) in the following
subsections.

###  3.1 Concise overview

Our goal is to use an unbiased estimator to measure the mean $M_{\rm
CMBlens}/M_{\rm SZ}$ mass ratio. With this in mind, we did not stack the maps
of all the clusters but instead built a CMB map for each cluster field, from
the nine frequency maps (three SPT and six Planck).

We first cut tangential maps out of the nine sky maps centered at the location
of each cluster. We then used an Internal Linear Combination (ILC, Sect. 3.4)
to build the best lensed CMB map around the cluster. Gravitational lensing of
the CMB induces correlations across its spatial scales. We used these
correlations to build a lensing potential with the minimum variance estimator
from Hu & Okamoto (2002). This lensing potential profile was then compared to
the lensing potential of a Navarro-Frenk-White (NFW, Navarro et al. 1996)
profile using a matched filter (Melin & Bartlett 2015). We can then derive the
individual mass from the NFW profile. Individual errors are large (typically
one to ten times larger than the signal) so we computed a weighted mean mass
of the sample as our measurement.

All those steps were applied on eleven different sets of 468 locations. The
first set is the one with the CMB cluster lensing signal. We called it the
_on_ set. It is centered at the locations of the clusters. The ten following
sets are ten draws at 468 random locations in the SPT-SZ footprint. We called
them the _off_ sets. We ran the same analysis on the _on_ and _off_ sets of
maps, and subtracted at the end the average signal of the ten _off_ sets from
the _on_ set: $\emph{on}-\langle\emph{off}\rangle_{10}$. This method allowed
us to remove spurious signal from instrumental correlated noises, foreground
or background astrophysical sources and any additive systematic bias.
Averaging ten random sets reduces the impact on the final error, which takes
the form
$\sigma=\sqrt{\sigma_{\emph{on}}^{2}+\sigma_{\emph{off}}^{2}}=\sqrt{\sigma_{\emph{on}}^{2}(1+\frac{1}{10})}$.

We assumed that each field contains the primary CMB (lensed!), the thermal SZ
(tSZ) signal, and noise. The noise can be of astrophysical or instrumental
origin. We thus write the maps of one field in Fourier space as

$\left\\{\begin{array}[]{@{}l@{}}m_{\nu_{0}}(\bm{k})=\alpha_{\nu_{0}}(\bm{k})S(\bm{k})+\beta_{\nu_{0}}(\bm{k})y(\bm{k})+n_{\nu_{0}}(\bm{k})\\\
m_{\nu_{1}}(\bm{k})=\alpha_{\nu_{1}}(\bm{k})S(\bm{k})+\beta_{\nu_{1}}(\bm{k})y(\bm{k})+n_{\nu_{1}}(\bm{k})\\\
...\end{array}\right.,$ (1)

or in a more compact way

$\bm{m}(\bm{k})=\bm{\alpha}(\bm{k})S(\bm{k})+\bm{\beta}(\bm{k})y(\bm{k})+\bm{n}(\bm{k}).$
(2)

with $S$ the primary CMB, $y$ the tSZ Compton map of the cluster, and
$n_{\nu_{i}}$ the instrumental and astrophysical noises, that is all the other
components in the frequency band $\nu_{i}$. We also have

$\left\\{\begin{array}[]{@{}l@{}}\alpha_{\nu_{i}}(\bm{k})=b_{\nu_{i}}(k)\times
t_{\nu_{i}}(\bm{k})\\\ \beta_{\nu_{i}}(\bm{k})=\alpha_{\nu_{i}}(\bm{k})\times
j_{\nu_{i}}\end{array}\right.\,,$ (3)

where $b_{\nu_{i}}$ and $t_{\nu_{i}}$ are respectively the isotropic beam and
the transfer function of the instrument in the frequency band $\nu_{i}$. This
transfer function is taken as $t_{\nu_{i}}=1$ for all Planck frequencies but
is anisotropic for SPT-SZ frequencies. Therefore, $\alpha_{\nu_{i}}$ and
$\beta_{\nu_{i}}$ depend on the vector $\bm{k}$ and not only on its modulus
for SPT-SZ data, making the analysis more complex. $j_{\nu_{i}}$ is the
frequency dependence of the tSZ effect, integrated in the $\nu_{i}$ band.

The kinetic SZ (kSZ) signal has the same frequency dependence as the primary
CMB and cannot be separated from it. We estimated its impact on our result in
Sect. 6.1 and corrected from it in our final measurement.

###  3.2 Tangential maps

We performed our analysis in Fourier space from small flat tangential maps. We
hence had to cut them out of the spherical HEALPix maps.

For each cluster, we cut a 10x10 deg2 map centered on it, when possible. This
size matches the size used for the cluster extraction in the Planck official
analysis (Planck Collaboration et al. 2016a). It is large enough to allow for
a good estimation of the covariance matrix of the maps. If the cluster was
close to the border of the field of the SPT-SZ public data, that is less than
17 arcmin ($\sim$40 pixels), we shifted the cluster as close as possible to
the center without including in the tangential maps more than 1% of bad
pixels, that is empty pixels out of the SPT-SZ field. Some clusters being too
close to the border, we had to remove them from our sample. We removed 48
clusters, reducing the SPT-SZ sample from 516 clusters with redshift to 468
clusters with redshift and well within the SPT-SZ public field. The SPT
transfer function is not azimuthally symmetric. We estimated it as in Melin et
al. (2021). We computed it, for each frequency, at the location of the 468 SPT
clusters under the form of an HEALPix map. We then cut 10x10 deg2 maps around
the cluster location, averaged them and symmetrised over the horizontal and
vertical directions to obtain a single transfer function per frequency.
Throughout the analysis, we apodised the maps in real space to avoid artifacts
when using the (2-dimensional) Fourier transform.

###  3.3 Point sources handling

The tangential maps contained point sources that we needed to remove to avoid
contamination of our lensing estimator. To do so, we replaced a circular zone
slightly larger than the FWHM of the beam centered at the point source
location by a constrained Gaussian realisation. This filling had to be
continuous with the surrounding data in order to reduce the impact on our
lensing estimator.

For each cluster, we had 9 maps – 6 Planck and 3 SPT-SZ – that we masked
separately. We first computed the signal-to-noise ratio of the point sources
in each map using a single frequency matched filter. We then iterated on the
sources verifying $S/N>6$ and masked them with a disk of radius $r_{i}$ of 3
(resp. 10) arcmin for SPT-SZ (resp. Planck).

To mask the point sources, we could not simply set the masked regions to zero,
no matter how smooth their edges were. We later measured correlations between
different k-modes (spatial scales) of the map, and these masked regions would
only have created fake correlations if filled with zeros. To ensure that we
were not creating artificial correlations, we used the Hoffman & Ribak (1991)
constrained Gaussian field method (see also Benoit-Lévy et al. 2013).

We first computed a random realisation of a Gaussian temperature field
$\tilde{m}$ from the 2-dimensional power spectrum $P(\bm{k})$ of the map to be
masked $m$. This random realisation was expected to fill the hole and,
although it already had the same variance as $m$, it needed to be continuous
with the neighbouring original map pixels. We thus had to constrain it, using
an annulus encircling the masked zone named the calibrating zone. We chose it
with an outer radius $r_{o}$ of twice the inner radius $r_{i}$, the radius of
the masked zone. We illustrate the disk and annulus in Fig. 2.

Figure 2: Diagram of the two zones used for in-painting. Zone 1 is the in-
painted area (the masked zone), and zone 2 is the area constraining our
Gaussian field to ensure continuity (the constraining zone). The inner radius
$r_{i}$ is 3 (resp. 10) arcmin for SPT (resp. Planck) and the outer radius
$r_{o}$ is twice the inner one.

We now had a zone to mask with a field $m_{1}$ and a constraining zone with a
field $m_{2}$ in the original map, and the same for the Gaussian realisation
($\tilde{m}_{1}$ and $\tilde{m}_{2}$).

The second step involved computing the 2-dimensional correlation function of
$m$, that we noted $\xi(\bm{r})$ with $\bm{r}$ the distance between two
pixels. $\xi(\bm{r})$ is the inverse Fourier transform of $P(\bm{k})$, and
describes the most likely value of a pixel given the values of other pixels at
different distances. From this function, we can obtain the covariance matrix
for chosen sets of pixels. We created two of them, one being the auto-
covariance matrix of the constraining zone $\Sigma_{22}$ and the other one the
covariance matrix between the masked zone and the constraining one
$\Sigma_{12}$. We understand now why we used the correlation function on the
whole map to compute those, as number of pixels in the constraining zone might
not be enough to compute $\Sigma_{22}$ accurately and, since the values of
$m_{1}$ are the ones we have to get rid of, we cannot use only them to compute
$\Sigma_{12}$.

The conditional probability distribution function of a field $m^{\prime}_{1}$
constrained by $m^{\prime}_{2}$ is

$\mathcal{P}(m^{\prime}_{1}|m^{\prime}_{2})=\frac{\mathcal{P}(m^{\prime}_{2}|m^{\prime}_{1})\mathcal{P}(m^{\prime}_{1})}{\mathcal{P}(m^{\prime}_{2})},$
(4)

and it is a Gaussian centered on

$\bar{m}^{\prime}_{1}=\langle\mathcal{P}(m^{\prime}_{1}|m^{\prime}_{2})\rangle=\Sigma_{12}\Sigma_{22}^{-1}m^{\prime}_{2}.$
(5)

$\bar{m}^{\prime}_{1}$ can be seen as the most probable value for each pixel
in zone 1 _individually_ , i.e. what we should get if we know $m^{\prime}_{2}$
and $\xi(\mathbf{r})$. This does not mean that $\bar{m}^{\prime}_{1}$ had the
same statistical properties as the CMB we wanted to reconstruct: it lacked the
deviations between the values in zone 1 inferred from zone 2 only and the
actual values. These deviations that describe the large scales that we cannot
constrain from zone 2 were obtained from the simulation.

Our field $\tilde{m}_{1}$ already contained the appropriate random variations
because we built it with $P(\bm{k})$, but it needed to be offset to have the
same mean as $m_{1}$ should have.

The third step was to compute the residual of the realisation, defined as

$m_{1,r}=\tilde{m}_{1}-\bar{\tilde{m}}_{1}=\tilde{m}_{1}-\Sigma_{12}\Sigma_{22}^{-1}\tilde{m}_{2},$
(6)

which is the field with only the random variations from the most probable zone
1, $\bar{\tilde{m}}_{1}$ (induced from the correlations with zone 2).

The last step was then to add $\bar{m}_{1}$, that is the mean $m_{1}$ should
have:

$m_{1,fill}=m_{1,r}+\bar{m}_{1}=\tilde{m}_{1}+\Sigma_{12}\Sigma_{22}^{-1}(m_{2}-\tilde{m}_{2}).$
(7)

$m_{1,fill}$ is used to fill the zone to mask, it is continuous with the
surrounding pixels and has the expected standard deviation.

###  3.4 Internal linear combination

We had nine point-source-free maps for each cluster, but they were still
contaminated by the tSZ effect (Carlstrom et al. 2002). The maps were written
in column in the from of Eq. 2.

We clearly see where our signal lies and we look for a minimum variance
estimate $\hat{S}$ of the pure lensed CMB map $S$. We used the method provided
by Remazeilles et al. (2011), called constrained Internal Linear Combinations
(ILC) to obtain our minimum variance estimate while nulling the tSZ signal. As
its name suggests, its aim is to recover the signal thanks to a linear
combination of the frequency maps, i.e.

$\hat{S}=\bm{w}^{t}\bm{m}.$ (8)

A derivation of the constraints $\bm{w}^{t}\bm{\alpha}=1$ and
$\bm{w}^{t}\bm{\beta}=0$ gives

$\bm{w}^{t}=\frac{(\bm{\beta}^{t}\hat{\bm{R}}^{-1}\bm{\beta})\bm{\alpha}^{t}\hat{\bm{R}}^{-1}-(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\beta})\bm{\beta}^{t}\hat{\bm{R}}^{-1}}{(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\alpha})(\bm{\beta}^{t}\hat{\bm{R}}^{-1}\bm{\beta})-(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\beta})^{2}},$
(9)

where $\hat{\bm{R}}$ is the empirical covariance matrix of the nine maps. In
their work, Remazeilles et al. (2011) compute $\hat{\bm{R}}$ with an average
on $l$-intervals in spherical harmonics in order to limit the statistical
variations. In our case, because the SPT transfer function is strongly
anisotropic, we could not use the same method. We computed $\hat{\bm{R}}$ on
all the _on_ and _off_ fields and performed the average on individual $\bm{k}$
modes to preserve the maximum of information. We averaged indiscriminately the
_on_ and _off_ maps because the CMB lensing signal is not expected to affect
the ILC reconstruction of the CMB when averaged over all the maps.

This method allowed us to recover the best CMB from the nine frequency maps,
but the noise at small scales was very large in the final map due to the
deconvolution of the beams and transfer functions $\bm{\alpha}$ and
$\bm{\beta}$. We thus applied an effective beam $\tilde{b}$ to the map that
flattened the noise at small scales. This effective beam was taken into
account in the lensing estimator in Sect. 3.5. For the Planck only analysis,
we used a Gaussian beam with a FWHM of 4.9 arcmin, whereas for the SPT-SZ only
analysis we used the value of $\alpha_{\nu_{i}}(k)$ that had the smallest
modulus mode-wise out of the three frequencies. For the joint analysis, the
effective beam was taken as the largest modulus mode-wise between the Planck
effective beam and the SPT-SZ effective beam in Fourier space.

The mean power spectrum of the final map convolved by the effective beam is:

$P_{\hat{S}}=\langle|\tilde{b}\bm{w}^{t}\bm{m}|^{2}\rangle=|\tilde{b}|^{2}\frac{(\bm{\beta}^{t}\hat{\bm{R}}^{-1}\bm{\beta})^{2}(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\alpha})-(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\beta})^{2}(\bm{\beta}^{t}\hat{\bm{R}}^{-1}\bm{\beta})}{\left[(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\alpha})(\bm{\beta}^{t}\hat{\bm{R}}^{-1}\bm{\beta})-(\bm{\alpha}^{t}\hat{\bm{R}}^{-1}\bm{\beta})^{2}\right]^{2}},$
(10)

where $\langle\rangle$ is the ensemble average over various hypothetical sky
realisations. $P_{\hat{S}}$ is used in the definition of the lensing potential
operator in the next section.

This method was used on the maps of the nine frequencies for $|\bm{k}|<175$
(or equivalently $l<6287$), i.e. when still in the original Planck resolution.
For $175\leq|\bm{k}|<278$, where Planck does not provide information, we only
used the three maps from SPT. At $|\bm{k}|\geq 278$, equivalent to $l>10000$,
neither SPT nor Planck provide any reliable information, we thus masked the
corresponding ILC $\bm{k}$ modes.

###  3.5 Lensing potential estimation

We then had to recover the lensing signal from the pure lensed CMB map
obtained in the previous section. This signal is linked to the gravitational
potential, directly related to the mass projected along the line of sight, and
can therefore be used to measure the mass of the cluster. We chose the Hu &
Okamoto (2002) minimum variance lensing estimator to do so. This estimator
uses the fact that the CMB lensing induces correlations between the CMB
$\bm{k}$ modes. They define the estimator in spherical harmonics but we used
the flat-sky equivalent, as done in Melin & Bartlett (2015). The deflection
(or lensing) potential $\hat{\phi}$ was computed in Fourier space as follows:

$\hat{\phi}(\bm{K})=A(\bm{K})\sum_{\bm{k}}\hat{S}^{*}(\bm{k})\hat{S}(\bm{k^{\prime}})F(\bm{k},\bm{k^{\prime}}),$
(11)

where $\bm{K}\equiv\bm{k}-\bm{k^{\prime}}\pmod{n}$, with $n$ the size of the
image in pixels along the x or y axis. The normalisation $A$, that is also the
variance of the reconstructed lensing potential, is defined as

$A(\bm{K})=\left[\sum_{\bm{k}}f(\bm{k},\bm{k^{\prime}})F(\bm{k},\bm{k^{\prime}})\right]^{-1},$
(12)

with the weights $F$ defined to minimise the variance of the estimator
$\hat{\phi}$:

$F(\bm{k},\bm{k^{\prime}})=\frac{f^{*}(\bm{k},\bm{k^{\prime}})}{2P_{\hat{S}}(\bm{k})P_{\hat{S}}(\bm{k^{\prime}})},$
(13)

$P_{\hat{S}}$ being the mean of the ILC power spectrum defined in Eq. 10. The
minimum variance filter $f$ is defined as:

$f(\bm{k},\bm{k^{\prime}})=\tilde{b}^{*}(\bm{k})\tilde{b}(\bm{k^{\prime}})\left[\tilde{C}_{k}\bm{k}.\bm{K}-\tilde{C}_{k^{\prime}}\bm{k^{\prime}}.\bm{K}\right],$
(14)

with $\tilde{b}$ the effective instrumental beam applied to the combination of
the nine maps (Sect. 3.4). $\tilde{C}_{k}$ is the power spectrum of the
unlensed CMB. We obtained it from an all-sky pure CMB Gaussian simulation at
HEALPix resolution $N_{side}=8192$, drawn from the Planck primary CMB power
spectrum, that we projected and apodised as we did for the sky maps (Sect.
3.2).

###  3.6 Matched filter

We then wanted to measure a mass from the reconstructed lensing profile. To
this end, we used the matched filter described in Melin & Bartlett (2015):

$M_{\rm CMBlens}/M_{\rm
fiducial}=\left[\sum_{\bm{K}}\frac{|\Phi(\bm{K})|^{2}}{A(\bm{K})}\right]^{-1}\sum_{\bm{K}}\frac{\Phi^{*}(\bm{K})}{A(\bm{K})}\hat{\phi}(\bm{K}),$
(15)

where $\Phi$ is the lensing potential of a Navarro-Frenk-White (NFW) density
profile (see Navarro et al. 1996) for a given fiducial cluster mass $M_{\rm
fiducial}$. We detailed this profile and its use in Sect. 4.1.

The matched filter returns the amplitude of our mass measurement with respect
to the chosen fiducial mass and can, because of a low S/N, return a negative
amplitude and as a consequence a negative mass. Because the individual error
bars were large, we averaged the individual measurement over the full sample
of 468 clusters. Performing individual measurements also allowed us to check
for redshift and mass dependence of our $M_{\rm CMBlens}/M_{\rm fiducial}$.

##  4 Simulations

We first tested the method on simulated data to study the statistics of the
recovered mass and assess for possible biases. Therefore, we built nine lensed
CMB maps for each cluster, with the instrumental characteristics of the nine
frequencies.

###  4.1 Cluster deflection potential

We chose a NFW density profile (Navarro et al. (1996)) for the cluster,

$\rho(r)=\frac{\delta_{c}\rho_{\rm c}}{(r/r_{s})(1+r/r_{s})^{2}},$ (16)

where $\rho_{\rm c}$ is the critical density at redshift $z$ and
$r_{s}=R_{500}/c_{500}$ the scale radius, with $R_{500}$ the radius inside
which the cluster density is 500 times $\rho_{\rm c}$. $c_{500}$ is the
concentration parameter that we assumed to be constant $c_{500}=3$ for
simplicity, although it is expected to weakly vary with the mass and redshift
of the cluster (e.g. Diemer & Kravtsov 2015; Ludlow et al. 2016). The
characteristic overdensity $\delta_{\rm c}$ is linked to $c_{500}$ as follows:

$\delta_{\rm
c}=\frac{500}{3}\frac{c_{500}^{3}}{\ln(1+c_{500})-c_{500}/(1+c_{500})}.$ (17)

We projected this density profile along the line of sight up to a distance of
$5R_{500}$ from the center of the cluster and obtained a surface density
profile

$\Sigma_{NFW}(r_{\perp})=\int_{los}\rho(r_{\perp},l)\mathop{}\\!\mathrm{d}l.$
(18)

We followed Bartelmann & Schneider (2001) to obtain the lensing profile from
the surface mass density of the cluster. We first got the convergence

$\kappa(\bm{\theta})=\frac{\Sigma_{NFW}(D_{d}\bm{\theta})}{\Sigma_{\rm
crit}},\;\mathrm{where}\;\Sigma_{\rm crit}=\frac{c^{2}}{4\pi
G}\frac{D_{s}}{D_{d}D_{ds}},$ (19)

with $\bm{\theta}$ the angular distance to the center of the cluster.
$\Sigma_{\rm crit}$ is the critical surface mass density above which the
lensing produces multiple images, i.e. it separates weak and strong lensing.
$D_{d}$, $D_{s}$ and $D_{ds}$, are respectively the distances observer-lens
(i.e. observer-cluster), observer-source (i.e. observer-CMB) and lens-source
(i.e. cluster-CMB).

From the convergence $\kappa$, we derived the lensing potential

$\Phi(\bm{\theta})=\frac{1}{\pi}\int_{\mathbb{R}^{2}}\kappa(\bm{\theta^{\prime}})\ln|\bm{\theta}-\bm{\theta^{\prime}}|\mathop{}\\!\mathrm{d^{2}}{\theta^{\prime}}.$
(20)

This potential was used in the simulation of the lensed CMB maps, but also in
the matched filter from Eq. 15. It does not depend on the observed frequency.
We used it for building the lensed maps.

###  4.2 Thermal SZ distortion

We also needed to model the thermal SZ signal from the cluster. This signal is
proportional to the electron pressure of the intracluster gas (Carlstrom et
al. 2002). We used the approach of Nagai et al. (2007) and Arnaud et al.
(2010), with a generalized NFW electron pressure profile (gNFW):

$P(r)=\frac{P_{500}P_{0}}{(r/r_{s})^{\gamma}[1+(r/r_{s})^{\alpha}]^{(\beta-\gamma)/\alpha}},$
(21)

where ($\gamma$, $\alpha$, $\beta$) are the slope of the profile at distances
$r\ll r_{s}$, $r\sim r_{s}$ and $r\gg r_{s}$ respectively. $P_{0}$ is the
normalisation and $P_{500}\propto M_{500}^{5/3}$ is a characteristic pressure
with $M_{500}$ the cluster mass inside $R_{500}$. The SZ signal is used to
detect clusters and to estimate masses in combination with external data
(Bocquet et al. 2019). We used these $M_{\rm SZ}$ masses (rescaled with a
factor 0.8 as discussed in Sect. 2.4) as the fiducial mass $M_{\rm fiducial}$
to simulate our clusters (dark matter NFW, tSZ) and build the lensing profile
for the matched filter.

We used the Arnaud et al. (2010) values for the parameters in Eq. 21 and
derived the Compton parameter $y$:

$y(\theta)=\frac{\sigma_{T}}{m_{e}c^{2}}\int_{los}P(D_{d}\theta,l)\mathop{}\\!\mathrm{d}l,$
(22)

with $\sigma_{T}$ the Thomson cross-section and $m_{e}c^{2}$ the electron rest
mass energy.

The SZ temperature distortion $\Delta T_{SZE}$ at a frequency $\nu$ and a
distance $r$ from the center of the cluster is then given by

$\frac{\Delta T_{SZE}}{T_{CMB}}(\nu,\theta)=j_{\nu}y(\theta),$ (23)

where we see the separated dependence of the tSZ effect in position $\theta$
and frequency $\nu$. We included relativistic corrections to $j_{\nu}$
following Sect. 3.3 of Melin & Pratt (2023). We used the cluster mass-
temperature from Eq. 19 of Bulbul et al. (2019) for each cluster and the
formula of Itoh et al. (1998). The tSZ spectrum including the relativistic
correction was integrated over the spectral bandwidth of each of the nine
frequency bands.

###  4.3 Final maps

For each frequency map, we put the same realisation of CMB $\tilde{S}$,
generated from the CMB power spectrum provided by Planck Collaboration et al.
(2016b). We then used the cluster deflection potential to obtain the lensed
CMB $S$ (Bartelmann & Schneider 2001).

$S=\tilde{S}+\nabla\tilde{S}.\nabla\Phi.$ (24)

This lensed CMB, common between all frequencies was then altered with
frequency or instrumental dependant elements. We added the SZ effect $\Delta
T_{SZE}$ in the maps and convolved them with the corresponding beam, before
adding the appropriate noise, the same way we defined the maps in Eq. 1. For
the simulations, we assumed that the instrumental noise is white with levels
of 40, 18 and 70 $\mu$K.arcmin at 95, 150 and 220 GHz for SPT-SZ, and 38, 31,
46, 1.4$\times 10^{2}$, 1.3$\times 10^{3}$, 5.5$\times 10^{4}$ $\mu$K.arcmin
at 100, 143, 217, 353, 545 and 857 GHz for Planck.

###  4.4 Results on simulations

For each of our 468 SPT-SZ clusters, we simulated nine frequency maps with the
cluster CMB lensing signal ($\Phi$ and $\Delta T_{CMB}$), the _on_ maps. We
also simulated ten sets of nine frequency maps without the CMB lensing signal
as the _off_ maps. By doing _on_ \- $\langle\emph{off}\rangle_{10}$, we
removed spurious signal from instrumental correlated noises, foreground or
background astrophysical sources and any additive systematic bias. In the
simulations, we first created periodic maps, that is where the edges have
continuity with the opposite one. We did not apodise the periodic maps. We
also placed clusters at the center of the maps to avoid any systematic effect
linked to non periodicity. In this scenario, we found no bias in our
measurement, and we did not need to remove the average of the _off_ maps to
estimate the masses. It shows that the method itself does not suffer from any
bias when the true lensing profile is known and when the maps do not contain
any correlated noise, or foreground/background sources.

This case, however, is not fully realistic. Some of our clusters were not
centered in the maps and the CMB is certainly not continuous from one edge to
the opposite one. If a cluster is not centered, it has a non periodic lensing
signal. Since it would otherwise cause problems when performing Fourier
transforms, we still put a periodic signal centered on the cluster position in
the matched filter and we used the _off_ maps to correct the bias it incurs.
In order to create non periodic maps with reduced correlation between their
edges, we created simulated maps five times larger than the final size, and we
kept only the central part.

Results for the joint Planck and SPT-SZ analysis are shown in Fig. 3 and last
line of Table 1. Individual measurements are shown as black dots and have
large errors bars (S/N $\sim 0.1-1$) in blue. Averaging over the full sample
provides $\langle M_{\rm CMBlensing}/M_{\rm SZ}\rangle=1.19\pm 0.18$ (green
stripe) compatible with one as expected.

Figure 3: Recovered CMB lensing mass $M_{\rm CMBlens}$ (divided by the input
fiducial mass $M_{\rm SZ}$) as a function of the input fiducial mass for the
joint Planck and SPT-SZ analysis of non periodic simulations. Black dots are
individual measurements. Error bars are shown in blue but we display only one
out of ten for clarity. The dotted (resp. dashed) red line is the zero (resp.
unity) level. The green stripe is the weighted mean $\langle M_{\rm
CMBlensing}/M_{\rm SZ}\rangle=1.19\pm 0.18$.

In Table 1, we also show the result for Planck only and SPT-SZ only
simulations in line one and two. Both results are compatible with one as
expected. More importantly, the Planck only error is $\sigma_{Planck}=0.27$,
while the SPT-SZ only error is $\sigma_{SPT}=0.25$, close the the Planck
error. The joint analysis provides an error
$\sigma=0.19\sim\sqrt{\frac{1}{1/\sigma_{Planck}^{2}+1/\sigma_{SPT}^{2}}}$
showing that there is indeed some additional information to gain in combining
the two data sets. We investigated farther the complementarity of the two data
sets in Sect. 5.

Table 1: Results on simulations, for a single realisation of a SPT-SZ-like
catalogue. $\begin{array}[]{p{0.5\linewidth}l}\hline\cr\vskip 3.0pt plus 1.0pt
minus 1.0pt\cr Data sets&M_{\rm CMBlensing}/M_{\rm SZ}\\\ \vskip 3.0pt plus
1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr{\em
Planck}&1.42\pm 0.27\\\ SPT-SZ&1.17\pm 0.25\\\ {\em Planck}\ + SPT-SZ&1.19\pm
0.18\\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\hline\cr\end{array}$

These runs on simulation demonstrate that our pipeline is able to extract the
average CMB lensing mass of the cluster sample, and provides a joint error
which is a factor of $\sim$$\sqrt{2}$ better than the error bar of individual
experiments. We also made use of these simulations to assess the impact of
modeling uncertainties on our results in Sect. 6.

##  5 Results

We first present the cluster mass massurements from Planck and SPT-SZ data
separately in Sect. 5.1. We then show the measurement from the joint analysis
in Sect. 5.2 and discuss the improvement with respect to individual analyses.

###  5.1 Masses from individual data sets

The pipeline was first used to measure the masses of the 468 clusters with the
Planck and SPT-SZ data sets separately in order to quantify the gain due to
the combination of both.

####  5.1.1 CMB lensing masses from Planck only

Figure 4: Recovered CMB lensing mass $M_{\rm CMBlens}$ (divided by the SZ mass
$M_{\rm SZ}$) as a function of the SZ mass for the Planck only analysis. Black
dots are individual measurements. Error bars are shown in blue but we display
only one out of ten for clarity. The dotted (resp. dashed) red line is the
zero (resp. unity) level. The green stripe is the weighted mean $\langle
M_{\rm CMBlensing}/M_{\rm SZ}\rangle=1.03\pm 0.27$.

The result for the Planck only data set is shown in Fig. 4. Black dots are
individual measurements, and individual error bars are shown in blue for one
measurement out of ten for clarity. The weighted mean of our measurement is
$M_{\rm CMBlens}/M_{\rm SZ}=1.03\pm 0.27$. This is a $3.7\leavevmode\nobreak\
\sigma$ measurement. It is consistent with unity, showing that the Planck CMB
lensing mass measurement is compatible with the SZ mass measurement $M_{\rm
SZ}$, that is the SPT-SZ mass measurement multiplied by 0.8. The error bars we
obtained are similar to the ones of the simulations, hinting that the Planck
maps are well modelled in the simulations.

####  5.1.2 CMB lensing masses from SPT-SZ only

Figure 5: Recovered CMB lensing mass $M_{\rm CMBlens}$ (divided by the SZ mass
$M_{\rm SZ}$) as a function of the SZ mass for the SPT-SZ only analysis. Black
dots are individual measurements. Error bars are shown in blue but we display
only one out of ten for clarity. The dotted (resp. dashed) red line is the
zero (resp. unity) level. The green stripe is the weighted mean $\langle
M_{\rm CMBlensing}/M_{\rm SZ}\rangle=1.12\pm 0.29$.

The result for the SPT-SZ only data set is shown in Fig. 5. We measured an
average ratio $M_{\rm CMBlens}/M_{SZ}=1.12\pm 0.29$, a
$3.9\leavevmode\nobreak\ \sigma$ measurement. It is in agreement with the
Planck only measurement and the SZ mass $M_{\rm SZ}$. In this case, the error
bars are somewhat larger than on the simulations. This is probably due to the
fact that the noise in the data is more complex than the assumption of simple
white noise in the simulations.

####  5.1.3 Comparison of Planck and SPT-SZ individual error bars

We now compare individual error bars from the Planck only and the SPT-SZ only
measurements. Fig. 6 (top panel) shows the SPT-SZ and the Planck error bars as
a function of redshift. The bottom panel shows the ratio of the two. Planck
provides error bars smaller than SPT-SZ at $z<1$, while SPT-SZ is more
efficient than Planck at $z>1$. However the ratio between the two remains
close to unity with a factor around 1.2 at maximum (reached for the lowest
redshift clusters). We note the strong correlation between the two errors bars
as a function of redshift. Error bars are given by the expression
$\left[\sum_{\bm{K}}\frac{|\Phi(\bm{K})|^{2}}{A(\bm{K})}\right]^{-1/2}$ with
$A(\bm{K})$, the variance of the Hu & Okamoto estimator, being independent of
the cluster properties (mass or redshift). Thus, at fixed mass, the redshift
dependence is given by the redshift dependence of the potential $\Phi$ (or
equivalently $\kappa$). So the redshift dependence of the ratio of the two
errors is driven by the angular diameter distance $D_{d}$ in
$\Sigma_{NFW}(D_{d}\bm{\theta})$.

Figure 6: Comparison of the errors on the SPT and Planck masses as a function
of redshift $z$. The top panel displays the errors while the bottom panel
displays the ratio of the errors. There is a clear correlation between the
ratio of the errors due to the redshift dependance of the lensing potential
$\Phi$.

Fig. 7 shows the SPT-SZ and the Planck error bars as a function of mass, and
the corresponding ratio. Planck outperforms SPT-SZ at high mass ($M_{\rm
SZ}>3\times 10^{14}M_{\odot}$) while SPT-SZ is more efficient at the low mass
end. The dispersion of error bar ratios with respect to mass exceeds that
observed for redshift dependence. This dispersion is driven by the redshift
distribution of clusters at a given mass.

Figure 7: Comparison of the errors on the SPT and Planck masses as a function
of the SZ mass $M_{\rm SZ}$. The top panel displays the errors while the
bottom panel displays the ratio.

###  5.2 Joint extraction of cluster masses

For the Planck and SPT-SZ combined analysis on the 468 clusters, we obtained a
global average ratio $M_{\rm CMBlens}/M_{\rm SZ}=0.92\pm 0.19$, that is a
$4.8\leavevmode\nobreak\ \sigma$ measurement of the lensing signal. As
expected, the joint measurement is also in agreement with the SZ mass $M_{\rm
SZ}$. The individual measurements are shown in Fig. 8. We averaged them in
five equally logarithmically spaced mass bins (red dots and associated error
bars). The global average ratio is shown as a green stripe. The red dots
present a small increasing trend with mass. We show the pull of the
measurements in Fig. 9, as a function of redshift (left) and mass (right). The
pull is close to a normal distribution. However, in the right panel, one can
see the same mass trend as in Fig. 8. This trend may be a hint that the
$M_{\rm SZ}$ masses are overestimated at low mass and underestimated at high
mass, or that there is a residual systematic bias in our measurement as a
function of mass. The pull as a function of redshift (left panel) does not
show any trend.

Figure 8: Recovered CMB lensing mass $M_{\rm CMBlens}$ (divided by the SZ mass
$M_{\rm SZ}$) as a function of SZ mass for the Planck \+ SPT-SZ combined
analysis. Black dots are individual measurements. The red dots are the
weighted means of the five mass bins, spaced logarithmically. The dotted
(resp. dashed) blue line is the zero (resp. unity) level. The green stripe is
the weighted mean within the error $\langle M_{\rm CMBlens}/M_{\rm
SZ}\rangle=0.919\pm 0.190$ (without correcting for the kSZ effect).

Figure 9: Pull of the ratio MCMBlens/MSZ and associated histogram as a
function of redshift (Left) and mass (Right). The blue line in each panel is a
Gaussian centered on 0 with $\sigma=1$ for comparison.

The global improvement from individual (3.7 and 3.9$\sigma$) to the joint
(4.8$\sigma$) measurement is clear. In Fig. 10, we aim at comparing the errors
from Planck only, SPT-SZ only and Planck+SPT-SZ as a function of the spatial
scale $K$ (or equivalently $L$). $A_{i}(K)$ is the variance of the error on
the estimated lensing potential $\hat{\phi}(K)$ for the considered experiment
$i$ ($i$ standing for SPT or Planck). The ratio $A_{\rm
SPT+Planck}(K)/A_{i}(K)$ thus compares the variance of the error of the joint
measurement to the variance of the error from an individual data set
measurement (SPT or Planck). This ratio is lower than unity for SPT and Planck
(blue and red triangles), showing that the joint analysis provides a better
measurement than individual analyses at all scales $K$. The Planck dataset
brings most of the information at large scales ($K<50.9$ or $L<1830$), as
expected, while the SPT-SZ dataset provides information at small scales
($K>145.0$ or $L>5208$). The ranges of scales observed by SPT-SZ and Planck
overlap but are not identical. We also computed the sum $A_{\rm
SPT+Planck}(K)/A_{\rm SPT}(K)+A_{\rm SPT+Planck}(K)/A_{\rm Planck}(K)$ (black
points). It is expected to be equal to unity if the two measurements of
$\hat{\phi}(K)$ from SPT and Planck are independent, that is the inverse
variance of the combination is the sum of the inverse variance of individual
measurements. This is the case for $K=18.9$ ($L=680$). For $K<18.9$ the sum
exceeds unity showing that there is redundant information in SPT-SZ and Planck
data. However, the combined analysis still provides a better measurement than
SPT-SZ or Planck only for these modes as already noticed. More interestingly,
the sum reaches values significantly below unity for $K>18.9$, showing that
the combination provides a better measurement than independent measurements:
the combination takes advantage of correlations between scales $\bm{k}$ from
SPT-SZ and scales $\bm{k^{\prime}}$ from Planck to measure
$\hat{\phi}(K=|\bm{k^{\prime}}-\bm{k}|)$. This is an important result showing
that ongoing and future ground based experiments such as Simons Observatory or
CMB-S4 will still benefit from combining the dataset from their large aperture
telescopes with Planck to improve CMB lensing cluster mass measurements.

Figure 10: Comparison of the variance $A(K)$ of the lensing profile
reconstructed for SPT-SZ and Planck binned in $K=|\bm{k^{\prime}}-\bm{k}|$ (or
equivalently $L=|\bm{l^{\prime}}-\bm{l}|$). The figure displays the ratios
between the combined $A_{\rm SPT+Planck}(K)$ and the single datasets
$A_{i}(K)$, and the sum of both ratios. It shows that the Planck data set
brings most of the information at large scales ($K<50.9$ or $L<1830$), while
the SPT-SZ dataset provides information on the smaller scales ($K>145.0$ or
$L>5208$). The sum is the product between the combined $A_{\rm SPT+Planck}$
term and the inverse variance sum $1/A_{\rm SPT}+1/A_{\rm Planck}$. It is
below unity for $K>18.9$ or $L>680$ demonstrating the significant gain of the
combination with respect to individual measurements.

##  6 Modeling uncertainties

The errors quoted in the previous sections are only statistical. We
investigated the impact of modeling uncertainties on our result in this Sect.
6. We studied the impact of the kSZ in Sect. 6.1. We considered the impact of
the assumed cluster matter profile in the dedicated Sect. 6.2. We studied the
effect of miscentering in Sect. 6.3. In Sect. 6.4 and 6.5, we studied the
impact of the errors on cluster redshift and mass. Finally, we quantified the
impact of the relativistic SZ effect in Sect. 6.6. Results are summarised in
Table 5 and discussed in Sect. 6.7.

###  6.1 Impact of the kinetic SZ

The spectral energy distribution of the kSZ is similar to that of the CMB.
Thus, we cannot separate those two components in the ILC and we want to check
the impact the kSZ has on the final result.

We compare the results on simulated maps without kSZ, and with the kSZ signal
induced by a normal velocity distribution of the galaxy clusters of 300 km/s,
identical for all masses and redshifts. For each cluster, we performed 200
different draws of CMB, instrumental noise and velocity. We then had
$200\times 468$ measurements that we compared with $200\times 468$
measurements with the same CMB and instrumental noise but a null velocity. On
average over these 200 simulations, our result $M_{\rm CMBlens}/M_{\rm SZ}$
was shifted by $-0.060\pm 0.015$. The shift is the same for different CMB
realizations, with a standard deviation (0.015) four times smaller than the
shift (0.060). We decided to correct our final result accordingly. We note
that Melin & Bartlett (2015) did not find a significant shift in their result
due to the kSZ. This is because the shift is much smaller than the statistical
error of their final result and they did not run enough simulations to extract
the kSZ effect from the fluctuations of the statistical noise.

###  6.2 Impact of the assumed lensing profile in the matched filter

Our baseline analysis was performed with a matched filter, assuming a lensing
profile based on a NFW profile truncated at $5\times r_{500}$. The real
lensing profile of individual clusters scatters around the NFW profile, but
one also needs to consider that clusters are not isolated objects. They are
located at the nodes of the cosmic web and are thus connected to filaments and
sheets. Random structures may also be present on the line-of-sight. In order
to assess the impact of our profile assumption on our result, we changed the
truncation radius to a lower ($3\times r_{500}$) and a higher ($7\times
r_{500}$) value. Our result $M_{\rm CMBlens}/M_{\rm SZ}$ changed by $+0.017$
and $-0.019$ respectively (Table 2). The assumed profile thus had an impact on
the result corresponding to about 10% of the statistical error bar. Separating
the effects of filaments, sheets and random structures along the line-of-sight
would require the use of N-body simulations, a work which is beyond the scope
of this article.

Table 2: Results on real maps changing the truncation radius used in the
matched filter. $\begin{array}[]{p{0.5\linewidth}l}\hline\cr\vskip 3.0pt plus
1.0pt minus 1.0pt\cr Truncation radius&\frac{M_{\rm CMBlens}}{M_{\rm SZ}}\\\
\vskip 3.0pt plus 1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr$3\times r_{500}$&0.936\pm 0.202\\\ $5\times r_{500}$
(baseline)&0.919\pm 0.190\\\ $7\times r_{500}$&0.900\pm 0.186\\\ \vskip 3.0pt
plus 1.0pt minus 1.0pt\cr\hline\cr\end{array}$

###  6.3 Impact of the cluster miscentering

In our analysis, we used the cluster positions given by Bocquet et al. (2019).
These positions are affected by an uncertainty. For this test, we used the
simulated maps for which we know the exact positions of the clusters. We drew
at random 468 ”observed” cluster positions from a normal law centered on the
468 ”real” positions, with errors provided by Eq. 16 of Melin & Pratt (2023).
We then performed the analysis on the maps, extracting the signal at the
”real” and ”observed” positions. The results are shown in Table 3. Without
(resp. with) miscentering corresponds to signal extracted at ”real” (resp.
”observed”) positions. We measured a bias $\Delta\frac{M_{\rm CMBlens}}{M_{\rm
SZ}}=+0.005$ much smaller than the error on the result $\sigma=0.181$. We
would need too many simulations to determine if the miscentering is only an
increase in dispersion or if it also adds a systematic shift in the result.

An alternative method allowed us to obtain the same result without
simulations: we used directly the lensing profile modeled for the matched
filter and we computed the matched filter between a modeled profile at the
”real” and at the ”observed” positions. The error used in the matched filter
computation is the one we get from the estimation of the lensing profile on
real maps. We did this calculation 100 times (corresponding to 100 ”observed”
positions) for each of the 468 clusters. With this method, we find a shift of
$\Delta\frac{M_{CMBlens}}{M_{SZ}}=-0.008$, close to the one obtained on
simulated maps but much more precise. There is also a dispersion due to the
velocity dispersion of the clusters of $\Delta\sigma=0.0005$ (calculated as
the standard deviation across the 100 ”observed” positions), negligible with
respect to the shift. With this, we show that miscentering creates a small
negative bias in the result. We also used this alternative method to quantify
the impact of the redshift and mass errors in the next sections.

Table 3: Results on simulated maps with and without miscentering.
$\begin{array}[]{p{0.5\linewidth}l}\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr Miscentering&\frac{M_{\rm CMBlens}}{M_{\rm SZ}}\\\ \vskip 3.0pt plus
1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr
Without&1.192\pm 0.181\\\ With&1.197\pm 0.181\\\ \vskip 3.0pt plus 1.0pt minus
1.0pt\cr\hline\cr\end{array}$

###  6.4 Impact of the redshift errors

Clusters are provided with their redshifts and associated errors. Once again,
we performed the analysis on simulated maps for which we know the ”real”
redshifts. We drew ”observed” redshifts with a normal law centered on the
”real” ones with standard deviation given by the redshift uncertainty in the
catalogue. The beginning of our analysis is not affected by a change in
redshift because no redshift information is used until the modeling of the
lensing profile in the matched filter. Unfortunately, the change in the result
due to the redshift uncertainty is too small to be measured with this method.
We thus used the alternative method: we used the matched filter on a lensing
profile modelled with the 100 ”observed” 468-redshift sets, compared to one
with the ”real” redshifts. The impact of the redshift errors on the matched
filter is a shift of $\Delta\frac{M_{CMBlens}}{M_{SZ}}=+0.0002$ and a
dispersion of about $\Delta\sigma=0.0008$, which are negligible with respect
to our statistical error $\sigma=0.190$.

###  6.5 Impact of the error on the mass $M_{\rm SZ}$

The direct method using simulations did not provide sufficient precision, as
for the redshift. Once again, we used the alternative method, and applied the
matched filter to compare the lensing profiles of the 468 clusters with ”real”
and ”observed” M500. We used 40 draws. We obtained a mean bias of
$\Delta\frac{M_{CMBlens}}{M_{SZ}}=-0.0022$ and a dispersion of
$\Delta\sigma=0.0070$. Both are negligible with respect to the measured
statistical error $\sigma=0.190$.

###  6.6 Impact of the relativistic SZ effect

Our baseline result includes the relativistic SZ effect, that is $j_{\nu_{i}}$
in Eq. 3 is computed including relativistic corrections to the SZ effect.
Details of the implementation are provided in Sect. 4.2. We removed the
relativistic correction from the calculation of $j_{\nu_{i}}$ and re-ran the
full analysis. We obtained a shift $\Delta\frac{M_{CMBlens}}{M_{SZ}}=-0.012$,
of order of 10% of the statistical error bar. In Table 4, we compare our
result on real data with and without including the relativistic SZ effect.

Table 4: Results on real maps with and without the relativistic SZ effect.
$\begin{array}[]{p{0.5\linewidth}l}\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr Relativistic SZ effect&\frac{M_{CMBlens}}{M_{SZ}}\\\ \vskip 3.0pt
plus 1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr
With&0.919\pm 0.190\\\ Without&0.907\pm 0.187\\\ \vskip 3.0pt plus 1.0pt minus
1.0pt\cr\hline\cr\end{array}$

###  6.7 Summary of modeling uncertainties

Our modeling uncertainties are summarised in Table 5. The total error (errors
added in quadrature) is 0.026 corresponding to about 14% of the statistical
error ($0.190$). This demonstrates that the analyses based on the current data
sets are dominated by statistical error driven by the instrumental noise of
the experiments. Note that we did not include the relativistic SZ effect in
our total error for the modeling uncertainties because it is a known effect
which has to be included in the baseline analysis, and we only wanted to
quantify its impact on our final results. This is the reason why it is
included as a separate line at the bottom of the table. Similarly, we only
added the dispersion due to kSZ in the total error as we correct for the kSZ
bias in the final result.

Table 5: Summary of the impact of the considered uncertainties on the
measurement ${M_{\rm CMBlens}}/{M_{\rm SZ}}$. The major uncertainty is the
profile modeling (about 10% of the statistical error). The other uncertainties
are of order of 1%. The considered uncertainties and their modelings are
described in dedicated subsections in Sect. 6.
$\begin{array}[]{p{0.5\linewidth}l}\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr Uncertainty&\Delta\frac{M_{\rm CMBlens}}{M_{\rm SZ}}\\\ \vskip 3.0pt
plus 1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr
Dispersion due to kSZ&\pm 0.015\\\ Profile up to $3\times r_{500}$&+0.017\\\
Profile up to $7\times r_{500}$&-0.019\\\ Miscentering&-0.008\\\ Error on
z&\pm 0.001\\\ Error on $M_{500}$&\pm 0.007\\\ \hline\cr{\bf Total}&{\bf\pm
0.026}\\\ \hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr Bias due to
kSZ&-0.060\\\ No relativistic SZ&-0.012\\\ \hline\cr\end{array}$

##  7 Summary and conclusions

We performed the first CMB-lensing galaxy cluster mass measurement using a
combination of ground and space-based surveys. We used the SPT-SZ cluster
catalogue and extracted its average CMB lensing to SZ mass ratio from the SPT-
SZ and Planck data sets. We found that the joint extraction outperforms the
extraction on single datasets. This was not guaranteed a priori because SPT-SZ
and Planck data observe the same primary CMB anisotropies, so the observations
could have been fully correlated, and the combination unproductive. This is
not the case. We detected the average CMB-lensing signal at $4.8\sigma$ in the
joint analysis compared with $3.7\sigma$ for the Planck only and $3.9\sigma$
for the SPT-SZ only analysis. We measured, before correcting for the kSZ,
$M_{\rm CMBlens}/M_{\rm SZ}=0.92\pm 0.19$ for the combination, $M_{\rm
CMBlens}/M_{\rm SZ}=1.03\pm 0.27$ for the Planck only analysis, and $M_{\rm
CMBlens}/M_{\rm SZ}=1.12\pm 0.29$ for the SPT-SZ only analysis. For the SPT-SZ
only analysis, our measured mass ratio is in agreement with the mass ratio
measured by Baxter et al. (2015) on a similar cluster sample and data set. Our
final combined result, after correcting for the kSZ effect, is

$M_{\rm CMBlens}/M_{\rm SZ}=0.98\pm 0.19\rm{\leavevmode\nobreak\ (stat.)}\pm
0.03\rm{\leavevmode\nobreak\ (syst.)}.$

Our measurement is dominated by the statistical error. We estimated the
modeling uncertainties to be of order of 14% of the statistical error.

We computed $M_{\rm CMBlens}/M_{\rm SZ}$, with $M_{\rm SZ}=0.8\times M_{\rm
SPT}$ where $M_{\rm SPT}$ is the mass given by Bocquet et al. (2019). This
rescales the SPT masses to the level of Planck masses. This shift between
Planck and SPT masses is recognized as a correction of a X-ray measurement
(hydrostatic and instrumental) mass bias $b$, where the SPT analysis finds the
bias to be $b=0.2$. Since we ”re-biased” the measurement, if $M_{\rm SPT}$
were an unbiased estimate of the CMB lensing mass, then our expected
measurement should be $M_{\rm CMBlens}/M_{\rm SZ}=M_{\rm CMBlens}/(0.8\times
M_{\rm SPT})=1.25$ instead of 1. Our measurement is compatible with both
$M_{\rm CMBlens}/M_{\rm SZ}=1.25$ ($b=0.2$) and $M_{\rm CMBlens}/M_{\rm SZ}=1$
($b=0$).

The analysis method is based on an internal linear combination of SPT-SZ and
Planck maps that reconstructs the best primary CMB maps around clusters and
nullifies the SZ effect. An optimal quadratic estimator is then applied to
each CMB map to extract the cluster lensing potential. The last step consists
in applying a matched filter to the lensing potential to measure the lensing
mass and compare its value with the SZ mass. The analysis pipeline has been
fully tested and characterized on simulations before being applied to the SPT-
SZ and Planck data sets. The major technical difficulties of the analysis
consist in taking properly into account the transfer function of the ground-
based experiment and dealing with the masking of point sources in the two data
sets.

We showed that the gain from the ground and space-based combination with
respect to analyses on individual data sets have two origins:

  * •

The spatial scales observed by SPT-SZ and Planck do not fully overlap, so the
two experiments bring different information from different primary CMB scales.
Planck is more efficient at observing large scales while SPT-SZ is better at
small scales.

  * •

The cluster lens correlates scales. Our method allows to extract information
from scales observed by Planck and scales observed by SPT-SZ that are
correlated by the lens. This adds to the efficiency of the combination.

These results demonstrate that both small (of order of arcmin) and large
(greater than 0.5 deg) scales are important for CMB-lensing cluster studies.
CMB-lensing cluster mass measurements with upcoming and future ground-based
large telescopes, such as Simons Observatory and CMB-S4, will therefore
benefit from a combination with Planck to make full use of the lensing
information.

###### Acknowledgements.

We would like to thank Jim Bartlett, Jacques Delabrouille, Boris Bolliet,
Íñigo Zubeldia and Erik Rosenberg for fruitful discussions about the results
presented in this article. We are also grateful to members of the Cosmology
group of the Particule Physics Department of CEA Paris-Saclay for useful
discussions while progressing on this analysis. We would also like to thank
Jim Rich for his helpful comments and suggestions on parts of this article. We
acknowledge the use of the Legacy Archive for Microwave Background Data
Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center
(HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at
the NASA Goddard Space Flight Center. We also acknowledge the use of the
Planck Legacy Archive. We used the HEALPix software (Górski et al. 2005)
available at https://healpix.sourceforge.io and the WebPlotDigitizer (Rohatgi
2022).

## References

  * Abazajian et al. (2019) Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv e-prints, arXiv:1907.04473
  * Abbott et al. (2020) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509
  * Ade et al. (2019) Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Albrecht et al. (2006) Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv e-prints, astro
  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409
  * Allen et al. (2002) Allen, S. W., Schmidt, R. W., & Fabian, A. C. 2002, MNRAS, 334, L11
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, Astronomy and Astrophysics, 517, A92, arXiv:0910.1234 [astro-ph]
  * Bahcall et al. (1997) Bahcall, N. A., Fan, X., & Cen, R. 1997, ApJ, 485, L53
  * Bartelmann & Schneider (2001) Bartelmann, M. & Schneider, P. 2001, Physics Reports, 340, 291, arXiv:astro-ph/9912508
  * Baxter et al. (2015) Baxter, E. J., Keisler, R., Dodelson, S., et al. 2015, ApJ, 806, 247
  * Baxter et al. (2018) Baxter, E. J., Raghunathan, S., Crawford, T. M., et al. 2018, Monthly Notices of the Royal Astronomical Society, 476, 2674
  * Benoit-Lévy et al. (2013) Benoit-Lévy, A., Déchelette, T., Benabed, K., et al. 2013, Astronomy & Astrophysics, 555, A37, arXiv:1301.4145 [astro-ph]
  * Bleem et al. (2015) Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, The Astrophysical Journal Supplement Series, 216, 27, arXiv:1409.0850 [astro-ph]
  * Bocquet et al. (2019) Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, The Astrophysical Journal, 878, 55, publisher: The American Astronomical Society
  * Bulbul et al. (2019) Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50
  * Carlberg et al. (1996) Carlberg, R. G., Yee, H. K. C., Ellingson, E., et al. 1996, ApJ, 462, 32
  * Carlstrom et al. (2002) Carlstrom, J. E., Holder, G. P., & Reese, E. D. 2002, Annual Review of Astronomy and Astrophysics, 40, 643, arXiv:astro-ph/0208192
  * Chan et al. (2023) Chan, V. C., Hložek, R., Meyers, J., & van Engelen, A. 2023, arXiv e-prints, arXiv:2302.13350
  * Chaubal et al. (2022) Chaubal, P. S., Reichardt, C. L., Gupta, N., et al. 2022, ApJ, 931, 139
  * Chown et al. (2018) Chown, R., Omori, Y., Aylor, K., et al. 2018, The Astrophysical Journal Supplement Series, 239, 10, aDS Bibcode: 2018ApJS..239…10C
  * Diemer & Kravtsov (2015) Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108
  * Dodelson (2004) Dodelson, S. 2004, Phys. Rev. D, 70, 023009
  * Geach & Peacock (2017) Geach, J. E. & Peacock, J. A. 2017, Nature Astronomy, 1, 795
  * Geach et al. (2019) Geach, J. E., Peacock, J. A., Myers, A. D., et al. 2019, ApJ, 874, 85
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Gupta & Reichardt (2021) Gupta, N. & Reichardt, C. L. 2021, ApJ, 923, 96
  * Hasselfield et al. (2013) Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, J. Cosmology Astropart. Phys., 2013, 008
  * Hoffman & Ribak (1991) Hoffman, Y. & Ribak, E. 1991, The Astrophysical Journal, 380, L5, aDS Bibcode: 1991ApJ…380L…5H
  * Holder & Kosowsky (2004) Holder, G. & Kosowsky, A. 2004, ApJ, 616, 8
  * Horowitz et al. (2019) Horowitz, B., Ferraro, S., & Sherwin, B. D. 2019, MNRAS, 485, 3919
  * Hu et al. (2007) Hu, W., DeDeo, S., & Vale, C. 2007, New Journal of Physics, 9, 441
  * Hu & Okamoto (2002) Hu, W. & Okamoto, T. 2002, The Astrophysical Journal, 574, 566, arXiv:astro-ph/0111606
  * Itoh et al. (1998) Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Levy et al. (2023) Levy, K., Raghunathan, S., & Basu, K. 2023, arXiv e-prints, arXiv:2305.06326
  * Lewis & Challinor (2006) Lewis, A. & Challinor, A. 2006, Phys. Rep, 429, 1
  * Lewis & King (2006) Lewis, A. & King, L. 2006, Phys. Rev. D, 73, 063006
  * Louis & Alonso (2017) Louis, T. & Alonso, D. 2017, Phys. Rev. D, 95, 043517
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201
  * Ludlow et al. (2016) Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214
  * Madhavacheril et al. (2015) Madhavacheril, M., Sehgal, N., Allison, R., et al. 2015, Phys. Rev. Lett., 114, 151302
  * Madhavacheril et al. (2017) Madhavacheril, M. S., Battaglia, N., & Miyatake, H. 2017, Phys. Rev. D, 96, 103525
  * Madhavacheril & Hill (2018) Madhavacheril, M. S. & Hill, J. C. 2018, Phys. Rev. D, 98, 023534
  * Madhavacheril et al. (2020) Madhavacheril, M. S., Sifón, C., Battaglia, N., et al. 2020, ApJ, 903, L13
  * Maturi et al. (2005) Maturi, M., Bartelmann, M., Meneghetti, M., & Moscardini, L. 2005, A&A, 436, 37
  * Melin & Bartlett (2015) Melin, J.-B. & Bartlett, J. G. 2015, Astronomy & Astrophysics, 578, A21, arXiv:1408.5633 [astro-ph]
  * Melin et al. (2021) Melin, J. B., Bartlett, J. G., Tarrío, P., & Pratt, G. W. 2021, A&A, 647, A106
  * Melin & Pratt (2023) Melin, J. B. & Pratt, G. W. 2023, A&A, 678, A197
  * Nagai et al. (2007) Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, The Astrophysical Journal, 668, 1, arXiv:astro-ph/0703661
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, The Astrophysical Journal, 462, 563, arXiv:astro-ph/9508025
  * Patil et al. (2020) Patil, S., Raghunathan, S., & Reichardt, C. L. 2020, ApJ, 888, 9
  * Planck Collaboration et al. (2016a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A27
  * Planck Collaboration et al. (2016b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A13
  * Planck Collaboration et al. (2016c) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016c, A&A, 594, A24
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A1
  * Planck Collaboration et al. (2020b) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A1
  * Raghunathan et al. (2018) Raghunathan, S., Bianchini, F., & Reichardt, C. L. 2018, Phys. Rev. D, 98, 043506
  * Raghunathan et al. (2019a) Raghunathan, S., Holder, G. P., Bartlett, J. G., et al. 2019a, J. Cosmology Astropart. Phys., 2019, 037
  * Raghunathan et al. (2019b) Raghunathan, S., Patil, S., Baxter, E., et al. 2019b, ApJ, 872, 170
  * Raghunathan et al. (2019) Raghunathan, S., Patil, S., Baxter, E., et al. 2019, Phys. Rev. Lett., 123, 181301
  * Raghunathan et al. (2017) Raghunathan, S., Patil, S., Baxter, E. J., et al. 2017, J. Cosmology Astropart. Phys., 2017, 030
  * Raghunathan et al. (2022) Raghunathan, S., Whitehorn, N., Alvarez, M. A., et al. 2022, ApJ, 926, 172
  * Remazeilles et al. (2011) Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, Monthly Notices of the Royal Astronomical Society, 410, 2481, arXiv:1006.5599 [astro-ph]
  * Rohatgi (2022) Rohatgi, A. 2022, Webplotdigitizer: Version 4.6
  * Ruhl et al. (2004) Ruhl, J., Ade, P. A. R., Carlstrom, J. E., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5498, Z-Spec: a broadband millimeter-wave grating spectrometer: design, construction, and first cryogenic measurements, ed. C. M. Bradford, P. A. R. Ade, J. E. Aguirre, J. J. Bock, M. Dragovan, L. Duband, L. Earle, J. Glenn, H. Matsuhara, B. J. Naylor, H. T. Nguyen, M. Yun, & J. Zmuidzinas, 11–29
  * Seljak & Zaldarriaga (2000) Seljak, U. & Zaldarriaga, M. 2000, ApJ, 538, 57
  * Story et al. (2013) Story, K. T., Reichardt, C. L., Hou, Z., et al. 2013, ApJ, 779, 86
  * Swetz et al. (2011) Swetz, D. S., Ade, P. A. R., Amiri, M., et al. 2011, ApJS, 194, 41
  * Tarrío et al. (2019) Tarrío, P., Melin, J. B., & Arnaud, M. 2019, A&A, 626, A7
  * Vale et al. (2004) Vale, C., Amblard, A., & White, M. 2004, New A, 10, 1
  * Vikhlinin et al. (2009) Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009, ApJ, 692, 1060
  * White et al. (1993) White, S. D. M., Navarro, J. F., Evrard, A. E., & Frenk, C. S. 1993, Nature, 366, 429
  * Yoo & Zaldarriaga (2008) Yoo, J. & Zaldarriaga, M. 2008, Phys. Rev. D, 78, 083002
  * Yoo et al. (2010) Yoo, J., Zaldarriaga, M., & Hernquist, L. 2010, Phys. Rev. D, 81, 123006
  * Zaldarriaga & Seljak (1999) Zaldarriaga, M. & Seljak, U. 1999, Phys. Rev. D, 59, 123507
  * Zubeldia & Challinor (2019) Zubeldia, Í. & Challinor, A. 2019, MNRAS, 489, 401
  * Zubeldia & Challinor (2020) Zubeldia, Í. & Challinor, A. 2020, MNRAS, 497, 5326