© The Authors 2024

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 Λ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 Ω_m and σ₈ (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 XXIV 2016; Boc-quet 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¹ and the Vera Rubin Observatory². 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 2009). The major uncertainty for cosmology with clusters at ɀ < 1 will, therefore, be overcome in the next decade. Extending cosmological studies with clusters atredshifts of ɀ > 1 will not be possible with the weak gravitational lensing on galaxies due to the lack of background sources at these red-shifts. 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 ɀ ~ 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 of 0 < ɀ < 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 datasets for which the lensing signal became detectable only appeared in the early 2010s (Ruhl et al. 2004; Swetz et al. 2011; Planck Collaboration VI 2020). The signal-to-noise ratio (S/N) for CMB lensing mass measurement is very low in these datasets. 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 XXIV 2016). 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. 2019; Gupta & Reichardt 2021; Levy et al. 2023; Chan et al. 2023; Saha et al. 2024). The motivation is to be ready to analyse data from the new generation of CMB instruments such as the Simons Observatory (SO, Ade et al. 2019) and CMB-S4 (Abazajian et al. 2019), which will provide datasets that will enable individual mass measurements for the first time. 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). In particular, Madhavacheril et al. (2020) performed the CMB cluster lensing extraction on the ACT and Planck data jointly. These measurements were performed using the intensity (i.e. temperature) maps of current datasets, but the polarised signal is expected to significantly improve the precision of the measurement in future datasets. In particular, Raghunathan et al. (2019c) recently succeeded in making the first detection of the CMB cluster lensing in the SPTpol data using only the polarised signal. In parallel to these efforts, the first cluster cosmological analyses using CMB cluster lensing to float the mass scale were performed (Planck Collaboration XXIV 2016; 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 dataset, 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 datasets provide. The combined S/N 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 datasets 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 summarise and conclude in Sect. 7.

Throughout the article, we work in the flat Λ CDM model with parameters H₀ = 70 km/s/Mpc, Ω_m = 0.3.

2 Datasets

We used two sets of sky maps in our analysis, SPT-SZ and Planck. Our analysis was 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 an ~2500 deg² area of the southern sky between 2008 and 2011, in three frequency bands centred on 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, which was 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 Sect. 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, which 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 co-ordinates, and in the 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 SPT-SZ band from Fig. 10 of Chown et al. (2018), which are useful for computing the frequency dependence 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 I 2020) 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 centred 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 XXVII 2016). 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 co-ordinates, 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 co-ordinate 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 datasets

The two datasets 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.5deg) 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 datasets.

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₅₀₀ (mass enclosed in a sphere of radius R₅₀₀ within which the average density is 500 times the critical density at the cluster’s redshift) are estimated using the SZ S/N versus mass relation for a flat ΛCDM model (h = H₀/(100 km/s/Mpc) = 0.7, Ω_m = 0.3, σ₈ = 0.8), which is calibrated from external X-ray and weak lensing datasets. 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_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 red-shift plane is shown in Fig. 1. Clusters span a large redshift range ɀ ∈ [0, 1.7] above a mass threshold that varies slowly with redshift M_sz≳2× 10¹⁴M_⊙.

Fig. 1 Mass-redshiſt distribution of the 468 SPT-SZ clusters used in our analysis. We display in blue the error bars for the clusters with Msz > 7 × 1014M⊙ 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_CMBlens/M_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 centred on 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 centred on 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: on − 〈off〉. 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 }_{on}^2+{\sigma }_{off}^2}=\sqrt{{\sigma }_{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 $$\{\begin{array}{l}{m}_{{\nu }_0}(k)={\alpha }_{{\nu }_0}(k)S(k)+{\beta }_{{\nu }_0}(k)y(k)+n_{{\nu }_0}(k)\hfill \\ m_{{\nu }_1}(k)={\alpha }_{{\nu }_1}(k)S(k)+{\beta }_{{\nu }_1}(k)y(k)+n_{{\nu }_1}(k),\hfill \\ \dots \hfill \end{array}$$(1)

or in a more compact way $$m(k)=\alpha (k)S(k)+\beta (k)y(k)+n(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 ν_i. We also have $$\{\begin{array}{l}{\alpha }_{{\nu }_i}(k)=b_{{\nu }_i}(k)\times t_{{\nu }_i}(k)\hfill \\ {\beta }_{{\nu }_i}(k)={\alpha }_{{\nu }_i}(k)\times j_{{\nu }_i}\hfill \end{array},$$(3)

where $b_{{\nu }_i}\text{and\hspace{0.17em}}{t}_{{\nu }_i}$ are respectively the isotropic beam and the transfer function of the instrument in the frequency band ν_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} \text{and} {\beta }_{{\nu }_i}$ depend on the vector 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 ν_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 10x 10 deg² map centred on it, when possible. This size matches the size used for the cluster extraction in the Planck official analysis (Planck Collaboration XXVII 2016). 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 (~40 pixels), we shifted the cluster as close as possible to the centre 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 deg² 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 artefacts when using the (two-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 centred on 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 nine maps – six Planck and three SPT-SZ - that we masked separately. We first computed the S/N 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 disc 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 m from the two-dimensional power spectrum P(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 annu-lus 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 disc and annulus in Fig. 2.

We now had a zone to mask with a field m₁ and a constraining zone with a field m₂ in the original map, and the same for the Gaussian realisation (m̃₁ and m̃ ₂).

The second step involved computing the 2-dimensional correlation function of m, that we noted ξ(r) with r the distance between two pixels. ξ(r) is the inverse Fourier transform of P(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 Σ₂₂ and the other one the covariance matrix between the masked zone and the constraining one Σ₁₂. 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 Σ₂₂ accurately and, since the values of m₁ are the ones we have to get rid of, we cannot use only them to compute Σ₁₂.

The conditional probability distribution function of a field $m_1^{\prime }$ constrained by $m_2^{\prime }$ is $$𝒫\left(m_1^{\prime }\mid m_2^{\prime }\right)=\frac{𝒫\left(m_2^{\prime }\mid m_1^{\prime }\right)𝒫\left(m_1^{\prime }\right)}{𝒫\left(m_2^{\prime }\right)},$$(4)

and it is a Gaussian centred on $${\overline{m}}_1^{\prime }=\langle 𝒫\left(m_1^{\prime }\mid m_2^{\prime }\right)\rangle ={\Sigma }_{12}{\Sigma }_{22}^{-1}{m}_2^{\prime }.$$(5)

${\overline{m}}_1^{\prime }$ can be seen as the most probable value for each pixel in zone 1 individually; in other words, what we should get if we know $m_2^{\prime }$and ξ(r). This does not mean that ${\overline{m}}_1^{\prime }$ 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 m̃₁ already contained the appropriate random variations because we built it with P(k), but it needed to be offset to have the same mean as m1 should have.

The third step was to compute the residual of the realisation, defined as $$m_{1,r}={\tilde{m}}_1-{\overline{\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, ${\overline{\tilde{m}}}_1$ (induced from the correlations with zone 2).

The last step was then to add ${\overline{m}}_1$, that is the mean m₁ should have: $$m_{1,fill}=m_{1,r}+{\overline{m}}_1={\tilde{m}}_1+{\Sigma }_{12}{\Sigma }_{22}^{-1}\left(m_2-{\tilde{m}}_2\right).$$(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.

Fig. 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 ri is 3 (resp. 10) arcmin for SPT (resp. Planck) and the outer radius ro is twice the inner one.

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 Ŝ 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; in other words, $$\widehat{S}=w^{t}m.$$(8)

A derivation of the constraints w^tα = 1 and w^tβ = 0 gives $$w^t=\frac{\left({\beta }^t{\widehat{R}}^{-1}\beta \right){\alpha }^t{\widehat{R}}^{-1}-\left({\alpha }^t{\widehat{R}}^{-1}\beta \right){\beta }^t{\widehat{R}}^{-1}}{\left({\alpha }^t{\widehat{R}}^{-1}\alpha \right)\left({\beta }^t{\widehat{R}}^{-1}\beta \right)-{\left({\alpha }^t{\widehat{R}}^{-1}\beta \right)}^2},$$(9)

where R̂ is the empirical covariance matrix of the nine maps. In their work, Remazeilles et al. (2011) compute 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 R̂ on all the on and off fields and performed the average on individual 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 a and ß. We thus applied an effective beam 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_{\widehat{S}}=\langle {\left|\tilde{b}{w}^{t}m\right|}^2\rangle =|\tilde{b}{|}^2\frac{{\left({\beta }^t{\widehat{R}}^{-1}\beta \right)}^2\left({\alpha }^t{\widehat{R}}^{-1}\alpha \right)-{\left({\alpha }^t{\widehat{R}}^{-1}\beta \right)}^2\left({\beta }^t{\widehat{R}}^{-1}\beta \right)}{{\left[\left({\alpha }^t{\widehat{R}}^{-1}\alpha \right)\left({\beta }^t{\widehat{R}}^{-1}\beta \right)-{\left({\alpha }^t{\widehat{R}}^{-1}\beta \right)}^2\right]}^2},$$(10)

where 〈〉 is the ensemble average over various hypothetical sky realisations. P_Ŝ 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 |k| < 175 (or equivalently l < 6287); that is, when still in the original Planck resolution. For 175 ≤ |k| < 278, where Planck does not provide information, we only used the three maps from SPT. At |k| ≥ 278, equivalent to l > 10000, neither SPT nor Planck provide any reliable information, we thus masked the corresponding ILC 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 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 $\widehat{\varphi }$ was computed in Fourier space as follows: $$\widehat{\varphi }(K)=A(K){\displaystyle \sum _k{\widehat{S}}^{*}}(k)\widehat{S}\left(k^{\prime }\right)F\left(k,k^{\prime }\right),$$(11)

where K ≡ k − k′ (mod 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(K)={\left[{\displaystyle \sum _{k}f}\left(k,k^{\prime }\right)F\left(k,k^{\prime }\right)\right]}^{-1},$$(12)

with the weights F defined to minimise the variance of the estimator $\widehat{\varphi }$: $$F\left(k,k^{\prime }\right)=\frac{f^{*}\left(k,k^{\prime }\right)}{2P_{\widehat{S}}(k)P_{\widehat{S}}\left(k^{\prime }\right)},$$(13)

P_Ŝ being the mean of the ILC power spectrum defined in Eq. (10). The minimum variance filter f is defined as: $$f\left(k,k^{\prime }\right)={\tilde{b}}^{*}(k)\tilde{b}\left(k^{\prime }\right)\left[{\tilde{C}}_{k}k\cdot K-{\tilde{C}}_{k^{\prime }}{k}^{\prime }\cdot 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 lens-ing profile. To this end, we used the matched filter described in Melin & Bartlett (2015): $$M_{\text{CMBlens}}/M_{\text{fiducial }}={\left[{\displaystyle \sum _K\frac{|\text{Φ}(K){|}^2}{A(K)}}\right]}^{-1}{\displaystyle \sum _K\frac{{\text{Φ}}^{*}(K)}{A(K)}}\widehat{\varphi }(K),$$(15)

where Φ is the lensing potential of a Navarro-Frenk-White (NFW) density profile (see Navarro et al. 1996) for a given fiducial cluster mass M_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_CMBlens / M_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 }_{\text{c}}}{\left(r/r_s\right){\left(1+r/r_s\right)}^2},$$(16)

where ρ_c is the critical density at redshift ɀ and r_s = R₅₀₀/c₅₀₀ the scale radius, with R₅₀₀ the radius inside which the cluster density is 500 times ρ_c. c₅₀₀ is the concentration parameter that we assumed to be constant c₅₀₀ = 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 δ_c is linked to c₅₀₀ as follows: $${\delta }_{\text{c}}=\frac{500}{3}\frac{c_{500}^3}{\ln \left(1+c_{500}\right)-c_{500}/\left(1+c_{500}\right)}.$$(17)

We projected this density profile along the line of sight up to a distance of 5R₅₀₀ from the centre of the cluster and obtained a surface density profile $${\Sigma }_{NFW}\left(r_{\perp }\right)={\displaystyle {\int }_{los}\rho }\left(r_{\perp },l\right)\text{d}l.$$(18)

We followed Bartelmann & Schneider (2001) to obtain the lens-ing profile from the surface mass density of the cluster. We first got the convergence $$\kappa (\theta )=\frac{{\Sigma }_{NFW}\left(D_d\theta \right)}{{\Sigma }_{\text{crit }}}\text{, where }{\Sigma }_{\text{crit }}=\frac{c^2}{4\pi G}\frac{D_s}{D_d{D}_{ds}}\text{, }$$(19)

with θ the angular distance to the centre of the cluster. Σ_crit is the critical surface mass density above which the lensing produces multiple images; that is, 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 κ, we derived the lensing potential $$\Phi (\theta )=\frac{1}{\pi }{\displaystyle {\int }_{{ℝ}^2}\kappa }\left({\theta }^{\prime }\right)\ln \left|\theta -{\theta }^{\prime }\right|{\text{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 intr-acluster gas (Carlstrom et al. 2002). We used the approach of Nagai et al. (2007) and Arnaud et al. (2010), with a generalised NFW electron pressure profile (gNFW): $$P(r)=\frac{P_{500}{P}_0}{{\left(r/r_s\right)}^{\gamma }{\left[1+{\left(r/r_s\right)}^{\alpha }\right]}^{(\beta -\gamma )/\alpha }},$$(21)

where (γ, α, β) are the slope of the profile at distances r << r_s, r ~ r_s and r >> r_s respectively. P₀ is the normalisation and $P_{500}\propto M_{500}^{5/3}$ is a characteristic pressure with M₅₀₀ the cluster mass inside R₅₀₀. 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_SZ masses (rescaled with a factor 0.8 as discussed in Sect. 2.4) as the fiducial mass M_flducial 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}{\displaystyle {\int }_{los}P}\left(D_d\theta ,l\right)\text{d}l,$$(22)

with σ_T the Thomson cross-section and m_ec² the electron rest mass energy.

The SZ temperature distortion ∆T_SZE at a frequency ν and a distance r from the centre 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 θ and frequency ν. We included relativistic corrections to j_ν 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 XIII (2016). We then used the cluster deflection potential to obtain the lensed CMB S (Bartelmann & Schneider 2001). $$S=\tilde{S}+\nabla \tilde{S}.\nabla \text{Φ}.$$(24)

This lensed CMB, common between all frequencies was then altered with frequency or instrumental dependant elements. We added the SZ effect Δ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 µK.arcmin at 95, 150 and 220 GHz for SPT-SZ, and 38, 31, 46, 1.4×10², 1.3×10³, 5.5×10⁴ µK.arcmin at 100, 143, 217, 353, 545 and 857 GHz for Planck.

4.4 Results of simulations

For each of our 468 SPT-SZ clusters, we simulated nine frequency maps with the cluster CMB lensing signal (Φ and Δ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 - 〈off〉₁₀, we removed spurious signal from instrumental correlated noises, foreground or background astrophysi-cal 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 centre 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, foreground or background sources.

This case, however, is not fully realistic. Some of our clusters were not centred in the maps and the CMB is certainly not continuous from one edge to the opposite one. If a cluster is not centred, it has a non-periodic lensing signal. Since it would otherwise cause problems when performing Fourier transforms, we still put a periodic signal centred 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 ~ 0.1-1) in blue. Averaging over the full sample provides 〈M_CMBlensing/M_sz〉 = 1.19 ± 0.18 (green stripe) compatible with one as expected.

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 σ_Planck = 0.27, while the SPT-SZ only error is σ_SPT = 0.25, close to the Planck error. The joint analysis provides an error $\sigma =0.19~\sqrt{\frac{1}{1/{\sigma }_{\text{Planck }}^2+1/{\sigma }_{SPT}^2}}$ showing that there is indeed some additional information to gain in combining the two datasets. We investigated farther the complementarity of the two datasets in Sect. 5.

These tests on simulations 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 $~\sqrt{2}$ better than the error bar of individual experiments. We also made use of these simulations to assess the impact of modelling uncertainties on our results in Sect. 6.

Fig. 3 Recovered CMB lensing mass MCMBlens (divided by the input fiducial mass Msz) 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 〈MCMBlensing/Msz〉= 1.19 ±0.18.

Fig. 4 Recovered CMB lensing mass MCMBlens (divided by the SZ mass Msz) 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 〈MCMBlensing/Msz〉 = 1.03 ± 0.27.

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 datasets

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

5.1.1 CMB lensing masses from Planck only

The result for the Planck only dataset 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_CMBlens/M_sz = 1.03 ± 0.27. This is a 3.7σ measurement. It is consistent with unity, showing that the Planck CMB lensing mass measurement is compatible with the SZ mass measurement Λf_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

The result for the SPT-SZ only dataset is shown in Fig. 5. We measured an average ratio M_CMBlens/M_sz = 1.12 ± 0.29, a 3.9σ measurement. It is in agreement with the Planck only measurement and the SZ mass M_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.

Fig. 5 Recovered CMB lensing mass MCMBlens (divided by the SZ mass Msz) 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 〈MCMBlensing/Msz〉 = 1.12 ± 0.29.

Fig. 6 Comparison of the errors on the SPT and Planck masses as a function of redshift ɀ. 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 dependence of the lensing potential Φ.

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. Figure 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 ɀ < 1, while SPT-SZ is more efficient than Planck at ɀ > 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[{\displaystyle \sum _K\frac{|\text{Φ}(\mathit{K}){|}^2}{A(\mathit{K})}}\right]}^{-1/2}$ with A(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 Φ (or equivalently κ). So the redshift dependence of the ratio of the two errors is driven by the angular diameter distance, D_d, in Σ_NFW(D_dθ).

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_sz > 3 × 10¹⁴M_⊙) 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.

Fig. 7 Comparison of the errors on the SPT and Planck masses as a function of the SZ mass MSZ. 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_CMBlens/M_sz = 0.92 ± 0.19, that is a 4.8σ measurement of the lensing signal. As expected, the joint measurement is also in agreement with the SZ mass M_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_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.

The global improvement from individual (3.7 and 3.9σ) to the joint (4.8σ) 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 $\widehat{\varphi }(K)$ for the considered experiment i (i standing for SPT or Planck). The ratio A_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 dataset measurement (Planck or SPT). This ratio is lower than unity for Planck and SPT (blue and red triangles), showing that the joint analysis provides a better measurement than individual analyses at all scales K. The Planck dataset provides the lowest error on the lensing potential at large scales (K < 50.9 or L < 1830), as expected, while the SPT-SZ dataset performs better at small scales (K > 145.0 or L > 5208). However, the gain on the error at L > 5000 does not mean that there is much information at these scales (see Appendix A). The ranges of scales observed by SPT-SZ and Planck overlap but are not identical. We also computed the sum A_SPT+Planck(K)/A_SPT(K) + A_SPT+Planck(K)/A_Planck(K) (black dots). It is expected to be equal to unity if the two measurements of $\widehat{\varphi }(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 k from SPT-SZ and scales k' from Planck to measure $\widehat{\varphi }\left(K=\left|k^{\prime }-k\right|\right)$. 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 (LATs) with Planck or their small aperture telescopes to improve CMB lensing cluster mass measurements from temperature data. Indeed, the transfer functions of the LATs might be similar to the ones of SPT-SZ, causing the loss of large scale modes that we proved to be relevant.

Fig. 8 Recovered CMB lensing mass MCMBlens (divided by the SZ mass MSZ) 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 〈MCMBlens/MSZ〉 = 0.919 ± 0.190 (without correcting for the kSZ effect).

Fig. 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 centred on 0 with σ = 1 for comparison.

6 Modelling uncertainties

The errors quoted in the previous sections are only statistical. We investigated the impact of modelling 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 miscentring in Sect. 6.3. In Sects. 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.

Fig. 10 Comparison of the variance A(K) of the lensing profile reconstructed for SPT-SZ and Planck binned in K = |k′ − k| (or equivalently L = |l′ − l|). The figure displays the ratios between the combined ASPT+Planck(K) and the single datasets Ai(K), and the sum of both ratios. It shows that the Planck dataset provides the lowest error on the lensing potential at large scales (K < 50.9 or L < 1830), while the SPT-SZ dataset performs better at small scales (K > 145.0 or L > 5208). The sum is the product between the combined ASPT+Planck term and the inverse variance sum l/ASPT + 1/APlanck- It is below unity for K > 18.9 or L > 680 demonstrating the significant gain of the combination with respect to individual measurements.

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 × 468 measurements that we compared with 200 × 468 measurements with the same CMB and instrumental noise but a null velocity. On average over these 200 simulations, our result M_CMBlens/M_SZ was shifted by −0.060 ± 0.015. The shift is the same for different CMB realisations, 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 × r₅₀₀. 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 × r₅₀₀) and a higher (7 × r₅₀₀) value. Our result M_CMBlens/M_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.

6.3 Impact of the cluster miscentring

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 centred 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) miscentring corresponds to signal extracted at real (resp. observed) positions. We measured a bias $\text{Δ}\frac{M_{\text{CMBlens}}}{M_{\text{SZ}}}=+0.005$ much smaller than the error on the result σ = 0.181. We would need too many simulations to determine if the miscentring 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 modelled for the matched filter and we computed the matched filter between a modelled 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 $\text{Δ}\frac{M_{\text{CMBlens }}}{M_{\text{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 ∆σ = 0.0005 (calculated as the standard deviation across the 100 observed positions), negligible with respect to the shift. With this, we show that miscentring 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.

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 centred 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 modelling 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 $\text{Δ}\frac{M_{\text{CMBlens}}}{M_{\text{SZ}}}=+0.0002$and a dispersion of about Δσ = 0.0008, which are negligible with respect to our statistical error σ = 0.190.

6.5 Impact of the error on the mass M_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 M₅₀₀. We used 40 draws. We obtained a mean bias of $\text{Δ}\frac{M_{\text{CMBlens }}}{M_{\text{SZ}}}=-0.0022$ and a dispersion of Δσ = 0.0070. Both are negligible with respect to the measured statistical error σ = 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 $\text{Δ}\frac{M_{\text{CMBlens }}}{M_{\text{SZ}}}=-0.012$, of the 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.

6.7 Summary of modelling uncertainties

Our modelling 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 datasets are dominated by statistical error driven by the instrumental noise of the experiments. We note that we did not include the relativistic SZ effect in our total error for the modelling 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.

7 Summary and conclusions

We used the SPT-SZ cluster catalogue and extracted its average CMB lensing to SZ mass ratio from the SPT-SZ and Planck datasets. 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σ in the joint analysis compared with 3.7σ for the Planck only and 3.9σ for the SPT-SZ only analysis. We measured, before correcting for the kSZ, M_CMBlens/M_SZ = 0.92 ± 0.19 for the combination, M_CMBlens/M_SZ = 1.03 ± 0.27 for the Planck only analysis, and M_CMBlens/M_SZ = 1.12 ± 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 dataset. Our final combined result, after correcting for the kSZ effect, is $$M_{\text{CMBlens }}/M_{\text{SZ}}=0.98\pm 0.19\text{ (stat.) }\pm 0.03\text{ (syst.). }$$

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

We computed M_CMBlens/M_SZ, with M_SZ = 0.8 × M_SPT where M_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 recognised 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_SPT were an unbiased estimate of the CMB lensing mass, then our expected measurement should be M_CMBlens /M_SZ = M_CMBlens/(0.8 × M_SPT) = 1.25 instead of 1. Our measurement is compatible with both M_CMBlens/M_SZ = 1.25 (b = 0.2) and M_CMBlens/M_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 characterised on simulations before being applied to the SPT-SZ and Planck datasets. 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 datasets.

We showed that the gain from the ground and space-based combination with respect to analyses on individual datasets 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 the 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 or their small telescopes to make full use of the lensing information in the temperature maps.

Acknowledgements

We thank the anonymous referee for his careful reading of the manuscript and for raising important points regarding the interpretation of our results. We would also 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 Web-PlotDigitizer software (Rohatgi et al. 2018) available at https://automeris.io.

Appendix A Inverse variance of the lensing estimator as a function of scale

Following Fig. 1 of Saha et al. (2024), we looked at the weight of each Fourier mode in the lensing signal. Fig. A.1 shows K|ϕ(K)|² /A(K) with respect to K, that is the contribution of each Fourier mode K to the lensing estimator S/N. The black curves show the expected result of the combination if the two datasets were independent. The three panels cover three different angular diameter sizes of lensing potential. They were chosen to be the smallest, median and largest profiles, computed within 5 × R₅₀₀ of the clusters in our sample. The method to model them is described in Sect. 4.1. The variance A(K) of the quadratic estimator is the one given in Eq. 12. This figure shows how Planck brings not only information at large scales as expected, but also more lensing signal at small scales, for all cluster sizes: the combination (dashed green line) is higher than the sum (black dots) of both datasets for L > 650.

We produced the same figure with Planck and CMB-S4 LATs simulations instead of Planck and SPT-SZ data. In this case, we assumed that CMB-S4 LATs have transfer functions similar to SPT-SZ, and we used the following instrumental characteristics:

• Frequencies: 27, 39, 93, 145, 225, 278 GHz.

• Respective FWHMs: 7.4, 5.1, 2.2, 1.4, 1.0, 0.9 arcmin.

• Respective noise levels: 21.34, 11.67, 1.89, 2.09, 6.9, 16.88 µK.arcmin.

The results are shown in Fig A.2. Even when the temperature data have much lower noise at small scales than the Planck data, Planckis still useful to measure the large scales (not observed by LATs) to better recover the lensing signal at all scales.

The gain in signal is not only due to a better cleaning of the tSZ. We computed the lensing estimator error in the case of a simple ILC, which does not explicitly separate the tSZ and simply treats it like other foregrounds and noise. We plot our results in Fig. A.3, for Planck and SPT-SZ in the same way as in Fig A.1. We see that there is still a significant gain (dashed green line) with respect to the case where Planck and SPT-SZ datasets are considered as independent (black dots).

Fig. A.1 Contribution of each Fourier mode K to the lensing estimator S/N. Four curves are drawn: SPT-SZ (blue triangles), Planck (red downward triangles), the sum of the two first curves (black dots) and the combination of both datasets (dashed green line). The three panels cover the range of ϕ(K) that are used, corresponding to the smallest (left), median (middle) and largest (right) cluster of our sample in term of angular diameter size.

Fig. A.2 Same as Fig A.1, but for simulations of CMB-S4 LATs (instead of SPT-SZ) and of Planck.

Fig. A.3 Same as Fig A.1, but without including the tSZ.

References

All Tables

All Figures

	Fig. 1 Mass-redshiſt distribution of the 468 SPT-SZ clusters used in our analysis. We display in blue the error bars for the clusters with Msz > 7 × 1014M⊙ and of one out of 20 under this threshold for clarity.
In the text

	Fig. 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 ri is 3 (resp. 10) arcmin for SPT (resp. Planck) and the outer radius ro is twice the inner one.
In the text

Fig. 3 Recovered CMB lensing mass MCMBlens (divided by the input fiducial mass Msz) 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 〈MCMBlensing/Msz〉= 1.19 ±0.18.

In the text

	Fig. 4 Recovered CMB lensing mass MCMBlens (divided by the SZ mass Msz) 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 〈MCMBlensing/Msz〉 = 1.03 ± 0.27.
In the text

	Fig. 5 Recovered CMB lensing mass MCMBlens (divided by the SZ mass Msz) 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 〈MCMBlensing/Msz〉 = 1.12 ± 0.29.
In the text

	Fig. 6 Comparison of the errors on the SPT and Planck masses as a function of redshift ɀ. 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 dependence of the lensing potential Φ.
In the text

	Fig. 7 Comparison of the errors on the SPT and Planck masses as a function of the SZ mass MSZ. The top panel displays the errors while the bottom panel displays the ratio.
In the text

Fig. 8 Recovered CMB lensing mass MCMBlens (divided by the SZ mass MSZ) 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 〈MCMBlens/MSZ〉 = 0.919 ± 0.190 (without correcting for the kSZ effect).

In the text

	Fig. 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 centred on 0 with σ = 1 for comparison.
In the text

Fig. 10 Comparison of the variance A(K) of the lensing profile reconstructed for SPT-SZ and Planck binned in K = |k′ − k| (or equivalently L = |l′ − l|). The figure displays the ratios between the combined ASPT+Planck(K) and the single datasets Ai(K), and the sum of both ratios. It shows that the Planck dataset provides the lowest error on the lensing potential at large scales (K < 50.9 or L < 1830), while the SPT-SZ dataset performs better at small scales (K > 145.0 or L > 5208). The sum is the product between the combined ASPT+Planck term and the inverse variance sum l/ASPT + 1/APlanck- It is below unity for K > 18.9 or L > 680 demonstrating the significant gain of the combination with respect to individual measurements.

In the text

Fig. A.1 Contribution of each Fourier mode K to the lensing estimator S/N. Four curves are drawn: SPT-SZ (blue triangles), Planck (red downward triangles), the sum of the two first curves (black dots) and the combination of both datasets (dashed green line). The three panels cover the range of ϕ(K) that are used, corresponding to the smallest (left), median (middle) and largest (right) cluster of our sample in term of angular diameter size.

In the text

	Fig. A.2 Same as Fig A.1, but for simulations of CMB-S4 LATs (instead of SPT-SZ) and of Planck.
In the text

	Fig. A.3 Same as Fig A.1, but without including the tSZ.
In the text