<figure><img src="image_1.png" /></figure>

# Joint measurement of the galaxy cluster pressure profile with Planck and SPT-SZ

J.-B Melin, G.W Pratt

► To cite this version:

J.-B Melin, G.W Pratt. Joint measurement of the galaxy cluster pressure profile with Planck and SPT-SZ. Astronomy and Astrophysics - A&A, 2023, 678, pp.A197. 10.1051/0004-6361/202346690. hal-04092607

HAL Id: hal-04092607

https://hal.science/hal-04092607

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Joint measurement of the galaxy cluster pressure profile with Planck and SPT-SZ

J.-B. Melin¹ and G. W. Pratt²

¹ Université Paris-Saclay, CEA, Département de Physique des Particules, 91191 Gif-sur-Yvette, France
e-mail: jean-baptiste.melin@cea.fr

² Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

Received 18 April 2023 / Accepted 8 August 2023

## ABSTRACT

We measured the average Compton profile of 461 clusters detected jointly by the South Pole Telescope (SPT) and Planck. The number of clusters included in this analysis is about one order of magnitude larger than in previous analyses. We propose an innovative method developed in Fourier space to combine optimally the Planck and SPT-SZ data, allowing us to perform a clean deconvolution of the point spread and transfer functions while simultaneously rescaling by the characteristic radial scale $R_{500}$ with respect to the critical density. The method additionally corrects for the selection bias of SPT clusters in the SPT-SZ data. We undertake a generalised Navarro-Frenk-White (gNFW) fit to the profile with only one parameter fixed, allowing us to constrain the other four parameters with excellent precision. The best-fitting profile is in good agreement with the universal pressure profile based on REXCESS in the inner region and with the Planck intermediate Paper V profile based on Planck and the XMM-Newton archive in the outer region. We investigate trends with redshift and mass, finding no indication of redshift evolution but detecting a significant difference in the pressure profile of the low- versus high-mass subsamples, in the sense that the low mass subsample has a profile that is more centrally peaked than that of the high mass subsample. We also scaled the average Compton profile by the mean Universe density ($R_{200m}$) and provide the best-fitting gNFW profile. Using the profiles scaled by both the critical ($R_{500}$) and the mean Universe density ($R_{200m}$), we studied the outskirt regions by reconstructing the average Compton parameter profile in real space. These profiles show multiple pressure drops at $\theta > 2\theta_{500}$, but these cannot clearly be identified with the accretion shocks predicted by hydrodynamical simulations. This is most probably due to our having reached the noise floor in the outer parts of the average profile with the current data sets.

**Key words.** galaxies: clusters: general – galaxies: clusters: intracluster medium – methods: statistical

## 1. Introduction

The Sunyaev-Zeldovich (SZ, Sunyaev & Zeldovich 1970, 1972) effect is the inverse Compton scattering of cosmic microwave background (CMB) photons by hot electrons in the Universe. It allows the electronic pressure of the hot gas in the cosmic web to be probed, as seen for example in the all-sky map obtained by Planck (Planck Collaboration XXII 2016). The SZ effect depends on the Compton parameter

$$y = \frac{\sigma_T}{m_e c^2} \int_{\text{los}} P_e dl \qquad (1)$$

where $\sigma_T$ is the Thomson cross-section, $m_e$ is the electron rest mass, c is the speed of light and $P_e = n_e \times T_e$ is the electronic pressure (product of electronic density and temperature). The integral is performed along the line-of-sight (los).

The cosmic web is composed of filaments, sheets, and nodes. Clusters of galaxies are located at the nodes. They constitute the densest regions of the cosmic web, dominated by gravitational processes. As a consequence, they do not show a preferred spatial scale and are expected to be self similar: at fixed cosmic time, any cluster can be matched to another according to a single parameter (the mass). This property allows for statistical analyses of clusters after rescaling, which we take advantage of in this work. In this article we use as our main parameter the cluster mass $M_{500}$, corresponding to the mass enclosed in a

sphere of radius $R_{500}$, whose mean density is 500 times the critical density of the Universe at the cluster redshift z. We note $\theta_{500} = R_{500}/D_{\text{ang}}(z)$ the angular radius of the cluster ($D_{\text{ang}}(z)$ being the angular distance of the cluster). We also use $M_{200m}$, $R_{200m}$ and $\theta_{200m}$, the equivalent quantities defined with respect to 200 times the mean density of the Universe.

Galaxy clusters are filled with a hot, gaseous, intracluster medium (ICM), the radial pressure profile of which decreases from the centre (radial coordinate $r = 0$) to the outskirts ($r \ge R_{500}$, $R_{500}$). The SZ effect is ideal to probe the pressure profile $P_e(r)$ of the ICM in galaxy clusters, and hence constrain models of structure formation for different cosmological scenarios. First measurements of pressure profiles based on X-ray data were performed by Nagai et al. (2007) and Arnaud et al. (2010) using Chandra and XMM-Newton respectively on relatively small cluster samples (a few to a few tens of objects). Nagai et al. (2007) fitted their data with a generalised Navarro-Frenk-White (gNFW) profile (Eq. (7) hereafter, see also Zhao 1996) which has been widely adopted by the cluster community since then. The goal of these studies was to support the imminent first SZ observations by providing a generic profile to improve the efficiency of cluster detection algorithms, to help in the analysis and interpretation of SZ data, and to improve in the theoretical modelling of the ICM (Appendix A of Nagai et al. 2007). The Nagai et al. (2007) and Arnaud et al. (2010) studies were extended to the low mass scale by Sun et al. (2011) using galaxy
groups and *Chandra* data. These three data sets are compatible. In the following, we use the pressure profile from Arnaud et al. (2010) based on the REXCESS sample (Böhringer et al. 2007) to compare with our results. We refer to this as the universal pressure profile (UPP). We also use the profiles extracted from the REXCESS Cool-core and REXCESS Morphologically disturbed subsamples.

First SZ profiles were observed with the Sunyaev–Zel’dovich Array (Mroczkowski et al. 2009, based on three clusters; Bonamente et al. 2012, 25 clusters), the South Pole Telescope (Plagge et al. 2010, 15 clusters), and the Atacama Cosmology Telescope (Sehgal et al. 2011, nine clusters). Planck Collaboration Int. V (2013) studied the average SZ profile of 62 clusters detected by Planck which were part of the XMM-Newton archive, and Sayers et al. (2013) the average profile of 45 X-ray selected clusters observed with Bolocam. These six studies all found profiles in good agreement with the UPP, but the latter two suggested some possible deviations at small ($r < 0.1R_{500}$) and large ($r > R_{500}$) radii. In particular, the average profiles from Planck Collaboration Int. V (2013; hereafter PIPV) and Sayers et al. (2013) showed higher pressure than the UPP at $r > R_{500}$.

The advent of high resolution SZ observations combined with Planck and/or XMM-Newton has allowed for further progress. Ruppin et al. (2017) studied the SZ profile of PSZ1 G045.85+57.71 with NIKA and Planck, and Ruppin et al. (2018) studied PSZ2 G144.83+25.11 with NIKA2 and Planck, demonstrating the high efficiency of the NIKA and NIKA2 Kinetic Inductance Detector experiments to map profiles at high (~12 arcsec) resolution. Recently, attention has focused on the synergy in spatial scales between different experiments for the measurement of SZ profiles. Romero et al. (2017) combined data from MUSTANG and Bolocam for 14 cluster profiles, and Romero et al. (2018) combined the data from MUSTANG, NIKA, Bolocam, and Planck to study the profile of CLJ 1226.9+3332 at redshift $z = 0.89$. Oppizzi et al. (2023) performed a measurement of the pressure profiles of six clusters observed with Planck and SPT-SZ. Similarly Pointecouteau et al. (2021) used the Planck and ACT data jointly to measure the average pressure profile of 31 clusters, demonstrating that space-based and ground-based data can be efficiently combined to improve constraints on the cluster pressure profile. Finally, Sayers et al. (2023) combined Bolocam, Planck and Chandra data to study the redshift evolution and mass dependence of the average pressure profile of 40 clusters.

The quest for high resolution individual SZ profile observations at high redshift is now ongoing. Using NIKA2, Kéruzoré et al. (2020) presented the profile of ACT-CL J0215.4+0030 at $z \sim 0.9$, and Ruppin et al. (2020) combined NIKA2 and Chandra to study MOO J1142+1527 at $z = 1.2$. These high redshift SZ profiles are in agreement with the UPP, indicating probable self similar evolution of the average cluster pressure profile with time. However, the studies with MUSTANG-2 and Chandra of IDCS J1426.5+3508 at $z = 1.75$ (Andreon et al. 2021), and JKCS 041 at $z = 1.8$ with MUSTANG-2 (Andreon et al. 2023) have indicated disagreement between the pressure profiles of these high $z$ clusters and the UPP.

In parallel to this observational effort, impressive progress has been made to understand the impact of astrophysical processes on the pressure profile of galaxy clusters using numerical simulations (Molnar et al. 2009; Shaw et al. 2010; Battaglia et al. 2012; McCarthy et al. 2014). He et al. (2021) recently combined simulations with XMM-Newton and Chandra data to better understand the hydrostatic X-ray mass

bias and its impact on the average pressure profile. Other recent simulation studies focused on cluster outskirts (Aung et al. 2021; Baxter et al. 2021; Zhang et al. 2021). These have showed that a pressure shock is expected at $2 < r < 2.5R_{200m}$ in the average pressure profile of relaxed clusters.

The observation of cluster outskirts at large radius in SZ is of growing interest (Ghirardini et al. 2018; Hurier et al. 2019). Pratt et al. (2021) stacked the Planck y-maps of ten nearby galaxy groups and studied the resulting average profile, finding a $3\sigma$ feature at $2R_{500}$ that they attributed to a signature of internal accretion shocks. Recently, Anbajagane et al. (2022) studied the stacked profile of the SPT-SZ clusters in the combined SPT-SZ and Planck y-maps (Bleem et al. 2022) and found a pressure deficit at $1.1R_{200m}$ ($3.1\sigma$) and a decrease in pressure at $4.6R_{200m}$ ($2.0\sigma$). They attributed these features to a shock-induced thermal non-equilibrium between electrons and ions for the former, and accretion shocks for the latter.

In this work, we average the SZ profiles of 461 South Pole Telescope (SPT) clusters in the Planck and SPT-SZ data sets. We compare the average profile extracted from Planck only, SPT-SZ only, and from the two data sets jointly, and provide the best-fitting gNFW model fit. We then compare our average profile to previous work, study cluster cores, and investigate possible trends with mass and redshift. We then perform the averaging of the SZ profiles in rescaling with $R_{200m}$, in a first attempt to obtain constraints on the profile in the cluster outskirts. In Sect. 2, we present the data sets we use for our study. We detail our analysis method in Sect. 3. We present our results in Sect. 4 and discuss them in Sect. 5. In Sect. 6, we summarise our results and conclude. Throughout the article we assume a flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km}^{-1} \text{ s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_m = 1 - \Omega_\Lambda = 0.3$.

## 2. Data sets

We present the cluster sample used in the analysis in Sect. 2.1 and the data sets in Sects. 2.2 and 2.3. The Planck and SPT-SZ data sets along with their formatting are the same as in Melin et al. (2021). We provide the main characteristics below, but encourage the reader to consult the above article for more detailed information.

### 2.1. Cluster sample

We used the fiducial SPT-SZ 2500d cluster catalogue published by the SPT collaboration (Bleem et al. 2015). It contains 677 candidates, out of which 620 are outside the SPT-SZ boundary mask provided with the public data<fnref n="1" />. We selected among them the 477 detections with redshift. We discarded five objects because they are on the boundary of the SPT-SZ footprint so we cannot extract the background statistics: less than 50% of the pixels are observed in the tangential maps centered on these clusters. Finally, we removed eleven clusters with $Y^{0.75} \text{ SPT} \ge 2 \times 10^{-4} \text{ arcmin}^2$, for which we showed that the filter transfer function (FTF) is not accurate (Fig. B.3. of Melin et al. 2021). Our final sample contains 461 clusters. $M_{500}$ is provided with the catalogue. Throughout this work, we renormalise the values of $M_{500}$ by a factor 0.8 to make them compatible with the Arnaud et al. (2010) mass definition, as was done in Tarrío et al. (2019) and Melin et al. (2021). All the quantities derived from the mass (e.g., $R_{500}$) are based on these rescaled masses. We only used the original SPT masses at four specific places in the

<footnote n="1">¹ https://lambda.gsfc.nasa.gov/product/spt/spt ProD_table.cfm</footnote>
article: in Sect. 3.3 for the temperature-mass relation; in Sect. 3.4 for the ‘unbiased SPT-significance’-mass relation; in Table 1 for naming the low-mass/high-mass subsamples; and in Sect. 5.5.2 for the scaling with $\theta_{200m}$. These choices are explained in each respective section.

## 2.2. Planck

We use the Planck full mission maps from the High Frequency Instrument. The instrument observed the all sky in six frequency bands: 100, 143, 217, 353, 545, and 857 GHz. We assume that the beams are Gaussian with full width at half maximum (FWHM) of 9.659, 7.220, 4.900, 4.916, 4.675, 4.216 arcmin respectively as done in Planck Collaboration XXVII (2016). The maps are provided in galactic coordinates at HEALPix resolution $N_{side} = 2048$ (1.72 pixel$^{-1}$) in the Planck Legacy Archive<fnref n="2" /> with their associated products, including the frequency response of each band. We upgraded the maps to $N_{side} = 8192$ (0.43 pixel$^{-1}$) by zero padding in harmonic space, that is, by adding new modes at high multipole with zero power. We then rotated the coordinate system from Galactic to equatorial, to match the SPT-SZ data format.

## 2.3. SPT-SZ

We work with the “SPT Only Data Maps”. SPT-SZ observed in three frequency bands: 90, 150 and 220 GHz. The maps are all provided at Gaussian resolution with FWHM of 1.75 arcmin in HEALPix format at resolution $N_{side} = 8192$ (0.43 pixel$^{-1}$) and in equatorial coordinates. The FTF are provided in the public archive with the maps. We retrieved the frequency responses of the instrument from Fig. 10 of Chown et al. (2018) using the WebPlotDigitizer<fnref n="3" />.

## 2.4. Consistency between Planck and SPT-SZ

Melin et al. (2021) performed extensive consistency tests between Planck and SPT-SZ; we summarise here the main results. Fluxes of SPT-SZ clusters extracted from Planck maps are fully compatible with published SPT fluxes. This result is discussed in Sect. 5.1 and illustrated in Fig. 1 of Melin et al. (2021). We note that Fig. 1 of Melin et al. (2021) shows that the SPT-SZ fluxes are subject to Malmquist bias, as expected. It is therefore important to correct for this bias when performing a joint Planck and SPT-SZ analysis of SPT-SZ selected clusters as we do in Sect. 3.4.

The signal-to-noise ratio of SPT-SZ clusters extracted from public SPT-SZ maps is 10% lower than the signal-to-noise ratio published by the SPT collaboration. This is expected because the public SPT-SZ maps are shallower than the private SPT-SZ maps. The cluster sizes extracted from the public SPT-SZ maps are in good agreement with the cluster sizes published by the SPT collaboration. Finally, cluster fluxes extracted from SPT-SZ public maps are also fully compatible with the flux published by the SPT-SZ collaboration except for clusters with $Y^{0.75} SPT \ge 2 \times 10^{-4}$ arcmin<fnref n="2" />, due to inaccuracy of the FTF. This is the reason for having excluded eleven clusters from our sample in Sect. 2.1. All of these consistency tests are presented in Appendix B of Melin et al. (2021).

# 3. Pressure profile in harmonic space

We give important conventions for the cluster profile rescaling and fitting in Sect. 3.1. We discuss the advantages of working in harmonic space in Sect. 3.2. We then introduce the internal linear combination to extract SZ maps in Sect. 3.3. In Sect. 3.4, we detail the Malmquist bias correction of the SPT-SZ signal. We describe in Sect. 3.5 the averaging procedure of the SZ maps. Then we present the likelihood used to perform the gNFW fits in Sect. 3.6.

## 3.1. Conventions

Following conventions in Arnaud et al. (2010), we defined the cluster scaled pressure profile as

$$p(x) = \frac{P_e(r)}{P_{e,500}} \qquad (2)$$

with the scaled radius

$$x = \frac{r}{R_{500}} \qquad (3)$$

and the characteristic electronic pressure

$$P_{e,500} = \frac{\mu}{\mu_e} P_{g,500}. \qquad (4)$$

Here $\mu = 0.6$ is the mean molecular weight and $\mu_e = \frac{5}{1+2/\mu} \approx 1.15$ the mean molecular weight per free electron, and

$$P_{g,500} = 1.65 \times 10^{-3} E(z)^{8/3} \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{2/3} h_{70}^2 \text{keV cm}^{-3} \qquad (5)$$

is the characteristic gas pressure (Eq. (3) of Nagai et al. 2007; Eq. (5) of Arnaud et al. 2010), with $E(z)$ the Hubble parameter at redshift $z$ normalised to its present value and $h_{70}$ the Hubble constant normalised to $70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Arnaud et al. (2010) noticed that $p(x)$ depends on $M_{500}$ and proposed to define (their Eq. (8))

$$p(x, M_{500}) = \mathbb{P}(x) \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{0.12}, \qquad (6)$$

where $\mathbb{P}(x)$ is the average scaled profile (gNFW) that we constrain in this work:

$$\mathbb{P}(x) = \frac{P_0}{(c_{500}x)^{\gamma} [1 + (c_{500}x)^{\alpha}]^{(\beta-\gamma)/\alpha}} \qquad (7)$$

Above, $P_0$ is the normalisation (independent of the mass), $c_{500} = R_{500}/r_s$ is the concentration parameter, $r_s$ is the scale radius, and $\alpha, \beta$ and $\gamma$ are the intermediate ($r \sim r_s$), outer ($r \gg r_s$) and inner ($r \ll r_s$) slopes of the profile, respectively. The average scaled profile $\mathbb{P}(x)$ thus has five parameters: $\Theta = [P_0, c_{500}, \gamma, \alpha, \beta]$.

We can now define the characteristic Compton parameter

$$y_{500} = \frac{\sigma_T}{m_e c^2} \times (2R_{500}) \times P_{e,500} \times \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{0.12} \qquad (8)$$

In Sect. 3.5, we renormalise the Compton parameter map of each cluster by the above quantity.

<footnote n="2">² https://pla.esac.esa.int</footnote>

<footnote n="3">³ https://automeris.io/WebPlotDigitizer/</footnote>
## 3.2. Advantages of working in harmonic space

The observed clusters were convolved by the instrumental point spread functions (PSF) of the experiments, and also by the FTF for SPT-SZ. The averaging process requires a radial rescaling of the profiles by $\theta_{500}$, with $\theta_{500}$ ranging from 0.88 to 16.5 arcmin. The effect of the PSF and FTF will therefore be different on the rescaled pressure profile from one cluster to another. To undertake the cleanest possible PSF and FTF deconvolution across the large range of cluster sizes, we performed our analysis in Fourier space. The PSF and FTF convolutions in real space translate into multiplications in Fourier space, making the deconvolutions easier. We adopted the flat sky approximation which allows us to associate the amplitude of a Fourier mode $||k||$ to the multipole $l$ in spherical harmonics. The conversion from $||k||$ to $l$ has the advantage to present results independent of the specific characteristics of the tangential maps. As a consequence, the measurements are presented as $y_l/y_{500}/\theta_{500}^2$ versus $l\theta_{500}$ (harmonic space, $l$ being the multipole) instead of $y(\theta)/y_{500}$ versus $\theta/\theta_{500}$ (real space, $\theta$ being the angular radius). Although the measurement is actually performed in Fourier space and then converted to harmonic space, we will use the term 'harmonic space' for the remainder of the article.

Working in harmonic space has another major advantage. To a first approximation, the instrumental and astrophysical noise in harmonic space can be considered as independent from one multipole to another. This is rigorously the case for instrumental white noise and un lensed primary CMB anisotropies. However, real data can contain some correlated instrumental noise in Fourier space, and also some lensing of the primary anisotropies. In the following, we therefore assume that the measurements in harmonic space are correlated from one multipole to another, and we use a covariance matrix $S$, measured on the data with bootstrap resampling, to take this correlation into account. In particular, we used this matrix to write the likelihood, detailed in Sect. 3.6. Full details on the calculation of $S$ are given in Appendix A.

The method to extract the average cluster pressure profile is now presented in detail in the following Sections. We first extracted an SZ map for each cluster using an internal linear combination (Sect. 3.3). An important aspect, detailed in Sect. 3.4, is to correct for the Malmquist bias of the SPT clusters in the SPT-SZ maps. We then renormalised, rescaled, and azimuthally averaged each map, and finally bin averaged them (Sect. 3.5).

## 3.3. Internal linear combination

We followed closely the method developed in Melin et al. (2021) for extracting clusters from Planck and SPT-SZ data jointly, with two main differences. First, we did not use a matched multi-filter (which assumes a cluster profile) but instead applied an internal linear combination (ILC) to reconstruct a Compton parameter $y$-map for each cluster. Secondly, we did not apply the method blindly, but centred the extraction at the location of the SPT-SZ clusters. We recall briefly below how we prepare the data and then provide details on the ILC.

We first cut 10 deg $\times 10$ deg tangential maps centred on the 461 clusters of our sample for each Planck and SPT-SZ map, yielding nine frequency maps in total for each cluster with pixel size of about 0.43 arcmin. As in Melin et al. (2021), we started by detecting and masking bright point sources ($S/N > 8$) with a single frequency matched filter, applied to individual channel maps. Map areas corresponding to masked point sources were

set to zero. This step is important to avoid contamination of the cluster signal by bright sources, and to lower the astrophysical noise in the data. We note that fourteen out of the 461 clusters are overlapping by less than $3 \times \theta_{500}$. We left these in the analysis, but checked their impact on the average profile (see Sect. 4.1). We write the tangential maps as

$$ \mathbf{m}(\mathbf{x}) = y_o \mathbf{t}_{\theta_s}(\mathbf{x}) + \mathbf{n}(\mathbf{x}) \quad (9) $$

where $\mathbf{m}(\mathbf{x})$ is the column vector corresponding to the tangential maps, $y_o \mathbf{t}_{\theta_s}(\mathbf{x})$ the cluster component, and $\mathbf{n}(\mathbf{x})$ the instrumental and astrophysical noise. The quantity $y_o \mathbf{t}_{\theta_s}(\mathbf{x})$ can be decomposed into the central Compton parameter $y_o$, and the cluster profile $\mathbf{t}_{\theta_s}(\mathbf{x})$, whose $i$th component is $f_v(v_i)[b_i * T_{\theta_s}](\mathbf{x})$. The quantity $T_{\theta_s}$ is the two dimensional SZ profile$^4$, convolved by $b_i$ (defined hereafter) and multiplied by $f_v$, the SZ frequency spectrum integrated over the frequency response. We included the relativistic corrections to the calculation of $f_v$ using the formula from Itoh et al. (1998) and the cluster temperature-mass relation$^5$ from Eq. (19) of Bulbul et al. (2019). The quantity $b_i$ is the PSF for Planck frequencies, and the PSF convolved by the FTF for SPT-SZ frequencies, respectively. As stated above, we used the flat sky approximation. We adopted the fast Fourier transform to perform all the calculations, and only translated to harmonic space as a final step. Equation (9) can be written in Fourier space as

$$ \mathbf{m}(\mathbf{k}) = \mathbf{a}(\mathbf{k}) y_o T_{\theta_s}(\mathbf{k}) + \mathbf{n}(\mathbf{k}) \quad (10) $$

with $\mathbf{a}(\mathbf{k})$ the mixing vector whose $i$th component is $f_v(v_i)b_i(\mathbf{k})$. The ILC provides an estimator of $y_o T_{\theta_s}$ (e.g., Remazeilles et al. 2011)

$$ \widehat{y_o T_{\theta_s}}(\mathbf{k}) = [\mathbf{a}^t(\mathbf{k}) \cdot \mathbf{P}^{-1} \cdot \mathbf{a}(\mathbf{k})]^{-1} \mathbf{a}^t(\mathbf{k}) \cdot \mathbf{m}(\mathbf{k}), \quad (11) $$

with $\mathbf{P}(\mathbf{k}) = \langle \mathbf{n}^*(\mathbf{k})\mathbf{n}(\mathbf{k}) \rangle$ being the power spectrum matrix of the noise across frequency channels. In practice $\mathbf{P}(\mathbf{k})$ was estimated directly from the maps by masking a circular region of $4\theta_{500}$ radius to exclude the studied cluster. The error on the estimator is given by

$$ \sigma(\mathbf{k}) = [\mathbf{a}^t(\mathbf{k}) \cdot \mathbf{P}^{-1} \cdot \mathbf{a}(\mathbf{k})]^{-1/2}. \quad (12) $$

We then divided the profiles $\widehat{y_o T_{\theta_s}}(\mathbf{k})$ and associated error $\sigma(\mathbf{k})$ by the pixel window function $w(\mathbf{k})$ to correct our profile for the 0.43 arcmin pixelisation. We rejected modes $\mathbf{k}$ corresponding to multipoles $l > 10 000$ in SPT-SZ data, for which the FTF is not provided (Chown et al. 2018).

## 3.4. Malmquist bias correction of the SPT-SZ signal

An important point of the analysis is to correct for the Malmquist bias of the SZ signal in the SPT-SZ data. Our cluster sample being selected in the SPT-SZ data, we expect the signal to be overestimated in the SPT-SZ maps, especially at low detection significance. We calculated a correction factor for each cluster to be applied to the SPT-SZ maps before combining them in the ILC. We did not apply a correction to the Planck data, since the sample has not been selected from it. We computed the

<footnote>$^4$ That is, the three dimensional SZ profile integrated along the line-of-sight.</footnote>

<footnote>$^5$ We note that we have not recalibrated the masses $M_{500}$ by 0.8 to compute the temperature because the relation has been derived with the SPT mass definition.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Effect of the Malmquist bias (MB) correction on the measured average SPT-SZ rescaled harmonic transform Compton parameter profile. The correction lowers the SPT profile by about 10%, which is larger than the statistical error in the range $3 \times 10^3$ arcmin $< l\theta_{500} < 2 \times 10^4$ arcmin.</figcaption></figure>

correction factor using the formalism in Bleem et al. (2015). We adopted the ‘unbiased SPT-significance’-mass relation

$$ \zeta = A_{SZ} \left( \frac{M_{500}}{3 \times 10^{14} M_{\odot} h^{-1}} \right)^{B_{SZ}} \left( \frac{E(z)}{E(0.6)} \right)^{C_{SZ}} \quad (13) $$

with $A_{SZ} = 4.14$, $B_{SZ} = 1.44$, $C_{SZ} = 0.59$, and a log-normal scatter, $D_{SZ} = 0.22$. h is the Hubble constant normalised to $100 \text{ km s}^{-1} \text{ Mpc}^{-1}$. The significance $\xi$ is related to the unbiased significance $\zeta$ via

$$ \zeta = \sqrt{\langle\xi\rangle^2 - 3} \quad (14) $$

for $\zeta > 2$. For each published SPT-SZ cluster Malmquist bias corrected mass $M_{500}$, and significance $\xi$, we drew $N = 1000000$ values of $\zeta_{\text{true}}$ following Eq. (13) including the log-normal scatter. We then computed $\xi_{\text{true}} = \sqrt{\zeta_{\text{true}}^2 + 3}$ for the $N$ values. We added a noise with normal distribution $N(0, 1)$ to $\xi_{\text{true}}$ to obtain observed significances $\xi_{\text{obs}}$. We selected the values that verify $\xi_{\text{obs}} > 4.5$ corresponding to the selection threshold of the SPT-SZ cluster catalogue, and $\zeta_{\text{true}} > 2$ corresponding to the validity range of Eq. (14). We computed the average of $\xi_{\text{obs}}$ in bins of $\xi_{\text{true}}$ for the selected values. We thus obtained a $\xi_{\text{obs}}-\xi_{\text{true}}$ relation for the considered cluster. We interpolated this relation at the published observed significance $\xi$ to obtain the average true significance $\xi_t$. We adopted the ratio $\xi_t/\xi$ as the correction factor to be applied to each SPT-SZ map. The correction is negligible for $\xi > 7$ clusters but is important for lower significance clusters, reaching 0.4 at the detection threshold $\xi = 4.5$. We note that we have not recalibrated the masses $M_{500}$ by 0.8 to compute the Malmquist bias correction factor, since Eq. (13) has been derived by the SPT collaboration with their mass definition.

The impact of the Malmquist bias correction on the average SPT-SZ rescaled harmonic transform Compton parameter profile is shown in Fig. 1. The uncorrected SPT-SZ profile is about 10% higher than the corrected profile. The correction is larger than the statistical error in the range $3 \times 10^3$ arcmin $< l\theta_{500} < 2 \times 10^4$ arcmin (see Sect. 3.5 below for definition of $l\theta_{500}$).

### 3.5. Averaging procedure

We renormalised each Fourier map obtained after the ILC $\hat{y}(k) = \hat{y}_0 \hat{T}_{\theta_s}(k)/w(k)$ by dividing it by $\theta_{500}^2 y_{500}$. We then bin averaged it in $\|k\|\theta_{500}$ bins. This procedure performs the spatial rescaling of the profile in $\theta_{500}$ and the azimuthal averaging at the same time. The averaging in bins was performed with inverse variance weights ($\sigma^2(k)/w^2(k)/\theta_{500}^4/\hat{y}_{500}^2$). The profiles binned in $\|k\|\theta_{500}$ were averaged together also using inverse variance weights. The final measurement is an estimate of the average of the two dimensional SZ profile of the cluster sample and its associated errors. This measurement is obtained in Fourier space $y_k/y_{500}/\theta_{500}^2$ vs. $k\theta_{500}$ and we converted it to harmonic space $y_l/y_{500}/\theta_{500}^2$ vs. $l\theta_{500}$. Errors on $\tilde{y}_l = y_l/y_{500}/\theta_{500}^2$ are the inverse weighted average errors. We denote these $\tilde{\sigma}_l$. We chose to adopt a logarithmic binning with 25 points per decade.

The link between the profile shape in harmonic space and its shape in real space is not intuitive for all. In Appendix B, we provide a visual representation of the correspondence between three dimensional pressure profiles in real space, the two dimensional projected SZ profiles in real space, and the two dimensional projected SZ profiles in harmonic space, to help the reader pass from one space to the other.

### 3.6. Likelihood for gNFW profile fits

We write the log likelihood as $\ln \mathcal{L}(\Theta) = -\frac{\chi^2(\Theta)}{2}$ with

$$ \chi^2(\Theta) = \sum_{l\theta_{500\text{bin}}, l'\theta_{500\text{bin}}}^{<l\theta_{500\text{max}}} [\tilde{y}_l - \tilde{y}_{\text{mod},l}(\Theta)] S_{ll'}^{-1} [\tilde{y}_{l'} - \tilde{y}_{\text{mod},l'}(\Theta)] \quad (15) $$

where $\Theta$ is the set of the five gNFW profile parameters. The quantity $\tilde{y}_{\text{mod},l}$ is the model, that is, the harmonic transform of the SZ profile constructed from the integrated gNFW profile along the line-of-sight. $S$ is the covariance matrix, estimated as detailed in Appendix A. The quantity $\tilde{y}_l$ is the average profile measurement obtained as described in the previous section. The sum is performed over the binned multipoles $l\theta_{500\text{bin}}, l'\theta_{500\text{bin}}$ lower than $l\theta_{500\text{max}}$.

For the fit to the full sample, we fixed $l\theta_{500\text{max}} = 8 \times 10^4$ arcmin, a value chosen to avoid the inclusion of systematics at high $l\theta_{500}$. Figure 2 (top) shows the number of pixels (resp. clusters) contributing to the measurement in a given $l\theta_{500\text{bin}}$ as the red (resp. blue) histogram. At low $l\theta_{500}$ values (around $10^2$ arcmin), the number of contributing clusters is around 100. As $l\theta_{500}$ increases, the number of contributing clusters increases together with the number of pixels. The number of pixels increases both because of the increasing number of clusters in each bin, and because of the logarithmic binning in $l\theta_{500}$ (the annulus contains more pixels when $l\theta_{500}$ increases). Around $l\theta_{500} = 10^3$ arcmin, all the clusters of the sample contribute, which explains the flattening of the histogram of clusters. However, the number of pixels continues to increase due to the logarithmic binning. Around $l\theta_{500} = 2 \times 10^4$ arcmin, the number of clusters decreases along with the number of pixels. At $l\theta_{500} = 8 \times 10^4$ arcmin (marked as the dashed black vertical line in Fig. 2), the number of clusters continues to decrease but the number of pixels increases again. This effectively means that the reconstructed profile is based on more pixels but fewer clusters. The information contained in the average is thus subject to possible systematic effects in the remaining about ten individual cluster profiles included in the averaging procedure. To minimise the risk of including systematic effects, we therefore exclude $l\theta_{500}$
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Top: number of pixels and clusters contributing to measurement of the harmonic transform Compton parameter profile in a given $l\theta_{500}$ bin. See text in Sect. 3.6 for a description of the two histograms and their interplay. Bottom: relative error on the profile. The solid red line is the inverse variance weighted errors used throughout the analysis (shown as red error bars in Figs. 3 and 5). The dotted-dashed green line corresponds to bootstrap errors. The dashed blue line is the raw scatter in a given $l\theta_{500}$ bin divided by the square root of the number of clusters contributing to this bin.</figcaption></figure>

bins above $l\theta_{500\text{max}} = 8 \times 10^4$ arcmin in our likelihood for the full sample. Using the same criteria to fix $l\theta_{500\text{max}}$ for the subsamples, we find the same value of $l\theta_{500\text{max}} = 8 \times 10^4$ arcmin for the high mass and low z subsamples.

We maximised the likelihood using Bayesian Monte Carlo Markov chain (MCMC) sampling. We used the emcee package (Foreman-Mackey et al. 2013) to marginalise over the gNFW parameters, adopting the same flat priors for all considered samples and subsamples except for the low mass subsample: $\mathcal{U}_{[-3,2]}$ for $\log(P_0)$, $\mathcal{U}_{[0.01,3]}$ for $c_{500}$ (or $c_{200}$), $\mathcal{U}_{[0.01,2.5]}$ for $\gamma$ and $\mathcal{U}_{[1,20]}$ for $\beta$. For the low mass sample, we also adopted flat priors but modified the upper bounds adopting $\mathcal{U}_{[-3,1.5]}$ for $\log(P_0)$, $\mathcal{U}_{[0.01,2]}$ for $c_{500}$, $\mathcal{U}_{[0.01,1.5]}$ for $\gamma$ and $\mathcal{U}_{[1,75]}$ for $\beta$. The Arnaud et al. (2010) UPP parameters were used as a starting point of the fit, shortening the burn-in phase, and the final converged chains contained ~150 000 iterations, which were subsequently thinned by a factor of five.

# 4. Results

We present results for the full cluster sample in Sect. 4.1. In Sect. 4.2, we provide the best gNFW fit. In Sect. 4.3, we study the intrinsic scatter around the profile.

## 4.1. Full sample profile

The blue diamonds in Fig. 3 show the average rescaled Compton parameter profile joint measurement in harmonic space for the full sample (461 clusters). We overplot the measurement based on Planck data only with black diamonds, and that based on SPT-SZ data only with red diamonds. For clarity, in this figure we multiplied (divided) the Planck (SPT-SZ) points by a factor of two. The Compton parameter profile corresponding to PIPV is plotted in green. We also plot the PIPV profile multi-

plied (divided) by two as long dashed (short dashed) green lines. The Planck data only and SPT-SZ data only profiles are compatible. The bottom panel shows the ratio between the Planck, SPT-SZ, Planck+SPT-SZ measurements and the PIPV profile on a logarithmic scale. It is clear from the figure that Planck provides better constraints than SPT-SZ on large scales ($l\theta_{500} < 2500$ arcmin). Conversely, on small scales ($l\theta_{500} > 6500$ arcmin), SPT-SZ outperforms Planck. The joint analysis therefore allows us to combine the two data sets optimally at all scales. Removing the fourteen clusters overlapping within $3 \times \theta_{500}$ from the average modifies only weakly the intermediate and external part of the profile ($l\theta_{500} < 10^4$ arcmin), and always within the statistical errors.

## 4.2. Best-fitting gNFW profile

We now discuss the best-fitting parameters and associated errors obtained from the likelihood described in Sect. 3.6. As found in other studies, leaving the five gNFW parameters free led to strong degeneracies between the three slope parameters $\gamma$, $\alpha$, $\beta$, and did not allow meaningful constraints to be put on all five parameters. We therefore fixed the intermediate slope to the UPP value ($\alpha = 1.05$) and left the other four parameters (amplitude $\log(P_0)$, concentration parameter $c_{500}$, internal slope $\gamma$ and external slope $\beta$) free. Results are shown in Fig. 4. Fixing $\alpha$ allowed the chains to converge to a small region of the parameter space and led to closed contours for all four free parameters. The best-fitting parameters are displayed with the dotted line; their numerical values (which maximise the likelihood given by the MCMC chains) and associated errors (at 68% confidence) are provided in the first row of Table 1. The $\chi^2$ of the fit is good, with $\chi^2 = 92.96$ and 68 degrees of freedom (72 data points minus four parameters), leading to a reduced $\chi^2$/d.o.f. = 1.37.

In the left hand panel of Fig. 5, we show the joint measurement (blue diamonds, identical to Fig. 3) with the best-fitting gNFW model overplotted (solid blue line). As in Fig. 3, the bottom panel shows the ratio with respect to the PIPV profile, but this time in a linear instead of a logarithmic y-scale. The best gNFW fit is similar to the UPP in the inner part and to the PIPV profile in the outer part. We will further discuss the comparison to previous work in Sects. 5.1 and 5.2. The right hand panel of Fig. 5 shows the 3D pressure profiles in real space corresponding to the 2D Compton parameter profiles shown in the left hand panel of the figure. The best-fitting gNFW parameters provided in Table 1 are sufficient to reproduce the average gNFW profile but are strongly correlated, as can be seen in Fig. 4. They are therefore not sufficient to fully describe the uncertainties on the best-fitting profile (that is, the error envelope). To compute the error envelope, we used the MCMC chains. At each r (in real space) and each $l\theta_{500}$ (in harmonic space), we used the 16 and 84 percentiles of the distribution of the profiles provided by the MCMC chains. The resulting 68% C.L. envelope based on the MCMC chains is delimited with dashed blue lines in both panels of Fig. 5. We provide the tabulated best-fitting profile and associated error envelope in harmonic and real space for the full sample in Table E.1. We graphically compare the error envelopes of the various subsamples in Fig. E.2.

## 4.3. Intrinsic scatter around the average profile

The average profile is computed as the inverse variance weighted mean of the individual profiles. The error bars are statistical. We computed 1000 bootstraps of the full analysis, where a single bootstrap corresponded to a draw with replacement of a sample
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Average rescaled harmonic transform Compton parameter profile measurements from *Planck* only, SPT-SZ only and joint (*Planck*+SPT-SZ) in black, red and blue respectively. We multiplied (resp. divided) the *Planck* (resp. SPT-SZ) measurements by a factor two for clarity. The joint analysis allows to take advantage of the two data sets at all scales. The solid green line is the PIPV profile. The long dashed (resp. short dashed) green line is the PIPV profile multiplied (resp. divided) by a factor of two.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Marginalised posterior likelihood for the parameters of the best-fitting gNFW model (Eq. (7)) with fixed $\alpha = 1.05$. Contours represent the 68% and 95% confidence regions. Blue contours show the results for the fit to the full data set; the dotted blue lines show the best-fitting value of each parameter as reported in Table 1.</figcaption></figure>
<table><caption>**Table 1.** Best-fitting parameters (given by the maximum of the likelihood) for the full sample when the profiles are rescaled with $\theta_{500}$.</caption><thead><tr><th></th><th>log(P<sub>0</sub>)</th><th>c<sub>500</sub></th><th>γ</th><th>β</th><th>α</th><th>χ<sup>2</sup>/d.o.f.</th></tr></thead><tbody><tr><td><strong>All</strong></td><td><strong>0.23</strong><sup>+0.16</sup><sub>-0.21</sub></td><td><strong>0.61</strong><sup>+0.14</sup><sub>-0.17</sub></td><td><strong>0.71</strong><sup>+0.09</sup><sub>-0.06</sub></td><td><strong>6.32</strong><sup>+0.88</sup><sub>-0.91</sub></td><td><strong>1.05</strong></td><td><strong>1.37</strong></td></tr><tr><td>z &lt; 0.6</td><td>0.59<sup>+0.17</sup><sub>-0.20</sub></td><td>1.07<sup>+0.22</sup><sub>-0.34</sub></td><td>0.55<sup>+0.12</sup><sub>-0.08</sub></td><td>4.90<sup>+0.49</sup><sub>-0.53</sub></td><td>1.05</td><td>1.19</td></tr><tr><td>z &gt; 0.6</td><td>1.11<sup>+0.13</sup><sub>-0.63</sub></td><td>1.46<sup>+0.02</sup><sub>-0.82</sub></td><td>0.19<sup>+0.53</sup><sub>-0.16</sub></td><td>5.06<sup>+0.89</sup><sub>-0.61</sub></td><td>1.05</td><td>1.37</td></tr><tr><td>M<sub>500</sub> &lt; 3.6 &times; 10<sup>14</sup> M<sub>&odot;</sub></td><td>-0.74<sup>+0.49</sup><sub>-0.13</sub></td><td>0.03<sup>+0.03</sup><sub>-0.01</sub></td><td>0.62<sup>+0.09</sup><sub>-0.15</sub></td><td>82.83<sup>+5.8</sup><sub>-63.93</sub></td><td>1.05</td><td>1.56</td></tr><tr><td>M<sub>500</sub> &gt; 3.6 &times; 10<sup>14</sup> M<sub>&odot;</sub></td><td>0.45<sup>+0.19</sup><sub>-0.19</sub></td><td>0.82<sup>+0.19</sup><sub>-0.22</sub></td><td>0.62<sup>+0.09</sup><sub>-0.10</sub></td><td>5.55<sup>+0.57</sup><sub>-0.61</sub></td><td>1.05</td><td>1.59</td></tr></tbody></table>

**Notes.** Errors are the 16 and 84 percentiles. In the fitting procedure, the gNFW parameter $\alpha$ is fixed to the UPP value of $\alpha = 1.05$. In this table, $M_{500}$ is the original SPT mass, uncorrected for the 0.8 factor (see Sect. 2.1 for details). log is the decimal logarithm. The final column gives the reduced $\chi^2$ for the corresponding fit.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Average pressure profile. Left: average rescaled harmonic transform Compton parameter joint measurement and best fitting gNFW model, with the REXCESS and PIPV profiles overplotted for comparison. The best fit (solid blue line) is close to the UPP at large $l\theta_{500}$ (the central regions), and to the PIPV profile at small $l\theta_{500}$ (the outer regions). Right: corresponding 3D pressure profiles. In both panels the dashed blue lines delineate the 68% C.L. uncertainties of our best fitting gNFW profile.</figcaption></figure>

of 461 clusters from the original sample of 461 clusters. We performed the full analysis again for each bootstrap sample. The bootstrap errors were then computed as the standard deviation of the 1000 recovered profiles. The result is shown in the bottom panel of Fig. 2, where the bootstrap errors are displayed in green. They are of the same order of magnitude as the statistical inverse variance weighted errors (shown in red) at low and high $l\theta_{500}$, but they are larger than the statistical errors in the intermediate $l\theta_{500}$ range ($2 \times 10^3$ arcmin $< l\theta_{500} < 2 \times 10^4$ arcmin). This effect is an indication of some intrinsic scatter ($\sigma_{int}$) in the profiles in the intermediate $l\theta_{500}$ range. We tried to measure it using the approach described in Sect. 5.3 of Planck Collaboration X (2011). We computed the raw scatter $\sigma_{raw}$ as the standard deviation of the 461 individual profiles. We show this quantity divided by $\sqrt{N}$ (the number of clusters contributing to the $l\theta_{500}$ bin under consideration) as the dashed blue line in the bottom panel of Fig. 2. We also computed $\sigma_{stat}$ from the square root of the mean of the variances of the individual profiles. Unfortunately, as this figure shows, $\sigma_{raw}$ is of the same order of magnitude as $\sigma_{stat}$ across the full $l\theta_{500}$ range, and the estimated difference $\sigma_{int}^2 = \sigma_{raw}^2 - \sigma_{stat}^2$ is compatible with zero. This estimator of the intrinsic scatter $\sigma_{int}$ is simple, robust, and has well-characterised uncertainties, but it is suboptimal when the variances of the

individual profiles are very different with respect to each other, which is the case here.

### 4.4. Impact of positional uncertainty

Before quantitative comparison to X-ray selected samples, it is prudent to consider what effect the significantly reduced positional uncertainty of the SZ experiments has on accuracy of the recovery of the inner slope of the gNFW profile. We estimated the impact of cluster positional uncertainty on the inner slope of the profile by injecting simulated SPT clusters at random locations in the SPT footprint. We adopted the UPP parameters to describe the pressure profile of the simulated objects. We first extracted clusters at the exact location of injection and ran the gNFW fitting analysis including the MCMC chains. We recovered the parameters of the UPP. We then extracted the clusters but with randomly shifted positions, with the positional uncertainty described by a Gaussian distribution with FWHM given by (arXiv v1 of Story et al. 2011; Song et al. 2012; Zenteno et al. 2020)

$$ \sqrt{\theta_{\text{beam}}^2 + (k\theta_c)^2 / \xi}. \qquad (16) $$
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Average rescaled harmonic transform Compton parameter joint measurement and its best fitting gNFW model compared to other gNFW profile fits published in other works. See text for discussion of the differences.</figcaption></figure>

Here, $\theta_{\text{beam}} = (1.1 + 1.6)/2$ arcmin is the mean of the FWHM of the 95 and 150 GHz original beams, $k = 2$, $\theta_c$ is the cluster core radius, and $\xi$ is the signal-to-noise ratio. The positional uncertainty varies between 0 and about 1 arcmin, with a median around 0.1 arcmin. Posterior contours for the best-fitting gNFW models are presented in Appendix D, and show no significant difference between the average profile extracted from the true and the randomly shifted positions. We therefore conclude that positional uncertainty has a negligible effect on our result.

# 5. Discussion

We now compare our results to previous work (Sect. 5.1). We discuss the consequence of our result on cluster cores in Sect. 5.2. In Sects. 5.3 and 5.4, we study the evolution of the profile with redshift and mass. We discuss the outskirts in Sect. 5.5.

## 5.1. Comparison to previous work

In Fig. 6, we show the joint measurement of the average rescaled harmonic transform Compton parameter (blue diamonds, similar to Figs. 3 and 5), with gNFW fits from various previous works overplotted. As already noted, the UPP (Arnaud et al. 2010), corresponding to the average profile of all 31 X-ray selected REXCESS clusters, is in excellent agreement with our measurement in the inner part of the profile (at high $l\theta_{500}$) but is slightly below in the outer part (at low $l\theta_{500}$). The second and third columns of Fig. B.1 in Appendix B help the reader to pass from harmonic to real space for the inner and outer part of the profile. Conversely, the PIPV profile (Planck Collaboration Int. V 2013), based on 62 nearby massive clusters observed by Planck, is in excellent agreement with our measurement in the outer part, but is below in the inner part. Our measurement is thus very well reproduced by the UPP in the inner part and by the PIPV profile in the outer part.

Sayers et al. (2013) averaged the signal of 45 massive clusters observed with Bolocam. The corresponding gNFW fit (solid blue line) is above our data points in the inner and in the outer

parts but passes through our data points in the central range $7 \times 10^3$ arcmin $< l\theta_{500} < 4 \times 10^4$ arcmin. Sayers et al. (2016) used Planck and Bolocam data for 47 massive clusters to constrain the amplitude $P_0$ and outer slope $\beta$ of the profile, fixing the other shape parameters to the UPP. This profile (triple dot dashed blue line) is in excellent agreement with our data in the inner region, but has an outer slope that is considerably steeper than our measurement. Pointecouteau et al. (2021) measured the profile of 31 clusters in Planck and ACT data. The gNFW fit is shown as the triple dot dashed red line in the figure. In contrast to Sayers et al. (2016), their profile is in excellent agreement with our measurement in the outer part, but lies significantly below in the inner part.

Our study comprises about a factor of ten more clusters than these previous works, and used an optimal treatment of the PSF and cluster rescaling thanks to the novel analysis in harmonic space. This allows us to significantly improve the constraints on the average pressure profile based on SZ measurements only. In particular, we only fixed the intermediate shape ($\alpha$) and determined the other four parameters ($\beta, \gamma, c_{500}, \log(P_0)$). Our contours (Fig. 4) are consequently reduced with respect to Fig. 5 of Sayers et al. (2013), for which one parameter was also fixed. Contours from Fig. 6 of Pointecouteau et al. (2021) are smaller than ours, but were obtained by fixing two parameters instead of just one. Anbajagane et al. (2022) also used the SPT-SZ cluster sample for their study. They used the y-map from Bleem et al. (2022), constructed from Planck and SPT-SZ data. They stacked the SPT-SZ clusters but did not provide a gNFW fit to the average profile. They focused instead on the outer regions to search for pressure shocks. One important missing piece in their analysis is the selection bias of the SPT-SZ sample in SPT-SZ data (Sect. 3.4), which is expected to impact their result.

## 5.2. Cluster cores

We now focus on the inner part of the average profile in Fig. 5. We recall that the positional uncertainty has a negligible impact on the recovered gNFW parameters, as discussed in Sect. 4.4. We also note that the constraint on the inner part of the profile originates primarily from the low redshift objects, as can be seen in Fig. 7. The joint measurement is in good agreement with the UPP (solid black line, based on all REXCESS clusters) in the inner part of the profile ($l\theta_{500} > 10^4$ arcmin). It is lower than the REXCESS cool-core profile and higher than the REXCESS morphologically disturbed and the PIPV profiles in this high $l\theta_{500}$ region. This indicates that the central slope in real space is closer to that of the UPP, shallower than that of the REXCESS cool-cores, and steeper than that of the REXCESS morphologically disturbed and the PIPV samples. This is an indication that the REXCESS and the SPT-SZ samples may contain similar fractions of cool-core and morphologically disturbed clusters. Rossetti et al. (2016) found that the Planck sample contains more disturbed clusters than three X-ray selected samples, including REXCESS. If the SPT-SZ selection were similar to that of Planck, we would have expected the SPT-SZ sample to contain more disturbed clusters than REXCESS. This would lead to us finding a shallower inner slope than the UPP, closer to that seen for the REXCESS morphologically disturbed profile, such as the PIPV profile. This is not what we find. Zenteno et al. (2020) also found no evidence that SPT-SZ selected clusters include a higher fraction of mergers than X-ray selected clusters. The only way to reconcile the results obtained by Rossetti et al. (2016) and Zenteno et al. (2020) and our work is to consider that
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Redshift subsamples. Left: average pressure profile for $z < 0.6$ (cyan circles), and $z \ge 0.6$ (orange stars). The solid blue line is the best-fitting gNFW model for the full sample (Table 1); the ratio plot on the bottom panel is calculated with respect to the best fit for the full sample. Right: marginalised posterior likelihood for the parameters of the best-fitting gNFW model of the low $z$ (cyan) and high $z$ (orange) subsamples, compared to the results for the full sample (blue). Contours represent the 68% and 95% confidence regions.</figcaption></figure>

the Planck and SPT-SZ selections differ, in the sense that they do not select the same fraction of disturbed clusters.

### 5.3. Trends with redshift

We built low and high redshift subsamples by splitting our sample in two at $z = 0.6$, a value close to the median redshift of the sample ($z = 0.56$). The low $z$ subsample contains 266 clusters, and the high $z$ subsample 195 objects. The resulting Compton parameter profile measurements are shown in Fig. 7. The low $z$ subsample is well measured across all $l\theta_{500}$, while the high $z$ subsample lacks data points at $l\theta_{500} > 2 \times 10^4$ arcmin (small scales) due to the small angular size of the high $z$ clusters with respect to the resolution of the two instruments.

We fitted the subsample profiles as described in Sect. 3.6, limiting the fit to $l\theta_{500\text{max}} < 2 \times 10^4$ arcmin for the high $z$ subsample owing to the more limited range of this profile. The corner plot for the parameters of the two samples is shown in Fig. 7 and the best fit values are given in Table 1. The fits for the low and high $z$ subsamples are both in agreement with the best fit of the full sample. We thus find no strong indication for evolution of the mean cluster pressure profile between $z = 0.40$ (median redshift of the low $z$ sample) and $z = 0.74$ (median redshift of the high $z$ sample).

We provide the tabulated best-fitting profile and associated error envelope in harmonic and real space for the low and high redshift subsamples in Tables E.2 and E.3. Furthermore, the low and high redshift subsample error envelopes are graphically compared to that of the full sample in Fig. E.2.

### 5.4. Trends with mass

We next divided our sample into two mass subsamples at $M_{500} = 3.6 \times 10^{14} M_{\odot}$. This mass is close to the median mass of the sample ($M_{500} = 3.58 \times 10^{14} M_{\odot}$). The low mass subsample contains 233 clusters with a median mass of $M_{500} = 3.15 \times 10^{14} M_{\odot}$.

The high mass subsample contains 228 clusters with a median mass of $M_{500} = 4.46 \times 10^{14} M_{\odot}$. The resulting subsample profiles are shown in Fig. 8. The profile of the low mass sample is significantly below the profile of the high mass sample over the whole $l\theta_{500}$ range. The profile of the high mass sample is slightly above the full sample gNFW fit.

We fitted the subsample profiles as described in Sect. 3.6, limiting the fit to $l\theta_{500\text{max}} < 6 \times 10^4$ arcmin for the high mass subsample owing to the more limited range of this profile. The best fitting gNFW models for the two subsamples are shown in Fig. 8 and the best fitting parameter values are given in Table 1. These fits confirm the visual impression that the low and high mass samples are not compatible. In Fig. 8 the 95% C.L. contours of the low mass subsample do not overlap for any of the parameters. In particular, for a given $\log(P_0)$, the profile of the low mass subsample has a lower $\gamma$ than the profile of the high mass subsample. This shows that the profile is less peaked for the low mass than for the high mass subsample at fixed $\log(P_0)$. This can be physically interpreted as an increase in the impact of non-gravitational processes on the pressure profile at low masses, or alternatively, by the faction of cool-core clusters being smaller for the low mass systems selected by SPT. On the other hand, at fixed $\gamma$, $\log(P_0)$ is lower for the low mass than for the high mass subsample. If one assumes that there is no change of the pressure profile with mass, this result could be due to an overestimated $P_{\text{e},500}$ value which would drive the fit to a lower $P_0$. This might point towards a possible under-correction of the Malmquist bias of the SPT data for the low mass sample. At this stage, it is not possible to decide between the two interpretations (less peaked profiles for low mass clusters, or residual systematics in the Malmquist bias correction of the SPT subsample).

We provide the tabulated best-fitting profile and associated error envelope in harmonic and real space for the low and high mass subsamples in Tables E.4 and E.5. We also compare the low and high mass subsample error envelopes to that of the full sample in Fig. E.2.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Mass subsamples. Left: average pressure profile for $M_{500} < 3.6 \times 10^{14} M_{\odot}$ (cyan circles), and $M_{500} \ge 3.6 \times 10^{14} M_{\odot}$ (orange stars). The solid blue line is the best-fitting gNFW model for the full sample (Table 1); the ratio plot on the bottom panel is calculated with respect to the best fit for the full sample. Right: marginalised posterior likelihood for the parameters of the best-fitting gNFW model of the low mass (cyan) and high mass (orange) subsamples, compared to the results for the full sample (blue). Contours represent the 68% and 95% confidence regions. Note that some parameters are plotted with logarithmic axes.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

## 5.5. Outskirts

### 5.5.1. External slope

We now focus on the outer part of the average profile of the full sample shown in Fig. 5. For $l\theta_{500} < 10^3$ arcmin, we showed that the joint measurement is in good agreement with the PIPV profile, and therefore has a shallower slope than all the profiles based on REXCESS. This result confirms the excess of pressure in the outer part of the cluster ($r > R_{500}$) with respect to Arnaud et al. (2010) found by Planck Collaboration Int. V (2013). For reference, the third row of Fig. B.1 shows harmonic and real space representations for the outer part.

The best-fitting gNFW model to the average profile of the full sample yields $\beta = 6.32_{-0.91}^{+0.88}$, a value which is higher than $\beta = 5.49$, the value for the UPP, although the slope for our profile is shallower in the outer part than the UPP. This result is due to degeneracies between parameters. In fact, the external slope is not determined only by $\beta$, but by a combination of $c_{500}$ and $\beta$, as can be seen in the first column of Fig. B.1. The degeneracy between $c_{500}$ and $\beta$ for our profile can clearly be seen in Fig. 4. Imposing the UPP value $c_{500} = 1.18$ would lead to a best-fitting value for $\beta$ lower than 5.49, as expected for a shallower slope.

We now look for a possible dependence of the outer slope with cluster mass, focussing on the $c_{500}$ versus $\beta$ plot in the right hand panel of Fig. 8. The contours for the low mass and high mass subsamples are not compatible, indicating a dependence of the external slope on mass. This can also be seen to some extent in the left hand panel of Fig. 8, but the data points of the low mass subsample have large uncertainties at $l\theta_{500} < 10^3$ arcmin.

### 5.5.2. Scaling in $\theta_{200m}$

Lau et al. (2015) found that the self-similarity in the outskirts is better preserved if the profiles are averaged with respect to the mean instead of the critical density of the Universe. We thus performed a new analysis but we rescaled the cluster

profiles with $\theta_{200m}$ instead of $\theta_{500}$. We converted SPT masses $M_{500}$ to $M_{200m}$, assuming a NFW profile and using $c_{200c}$ from Diemer & Kravtsov (2015). We did not include the renormalisation factor 0.8 on $M_{500}$ since the profile rescaled with $R_{200m}$ will not be compared to other works with different mass definition. We renormalised the profiles by $y_{200m}\theta_{200m}^2$ instead of $y_{500}\theta_{500}^2$, and we rescaled the radial quantities by $\theta_{200m}$ instead of $\theta_{500}$. Details of the derivation of $y_{200m}$ and associated characteristic pressure $P_{c,200m}$ are provided in Appendix C. The average rescaled harmonic transform Compton parameter profile is shown in Fig. 9 and the values for the best gNFW fit are given in Table 2. The fit is acceptable, with a reduced $\chi^2$/d.o.f. = 1.42.

We provide the tabulated best-fitting profile and associated error envelope in harmonic and real space for the profile scaled by $\theta_{200m}$ in Table E.6.

### 5.5.3. Can we detect an accretion shock?

To investigate possible accretion features in the outskirts, we reconstructed the 2D real space Compton profile directly from the data points of the left hand panel of Fig. 5. This was undertaken by drawing a Gaussian realisation with a mean corresponding to the data points and correlations from the covariance matrix S. We then projected this realisation into a two dimensional map and performed the inverse Fourier transform. This provides one realisation of a possible 2D real space profile. We performed 1000 such random realisations, projections and inverse Fourier transforms. The mean of the 2D real space profiles is shown as a black line in the left hand panel of Fig. 10, and the region encompassing the raw standard deviation of the realisations is shown by the red envelope. The solid green line is the PIPV profile and the blue solid and dashed lines show the best fitting gNFW model and its 68% confidence envelope.

The reconstructed profile follows closely the gNFW fit up to $\theta \sim 2\theta_{500}$, and does not show a single feature. Instead, multiple pressure drops are visible at $\theta > 2\theta_{500}$. To test if these multiple
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Pressure profile scaled by $\theta_{200m}$ and best fitting model. Left: same as Fig. 5 Left but averaged according to $\theta_{200m}$ instead of $\theta_{500}$. Right: same as Fig. 4 but for the profile averaged according to $\theta_{200m}$. Values for the gNFW best-fitting parameters are given in Table 2.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Table 2.** Best-fitting parameters (given by the maximum of the likelihood) for the full sample when the profile is rescaled by $\theta_{200m}$.</figcaption></figure>

<table><thead><tr><th>log(P<sub>0</sub>)</th><th>c<sub>200m</sub></th><th>γ</th><th>β</th><th>α</th></tr></thead><tbody><tr><td><b>0.35</b><sup>+0.21</sup><sub>-0.32</sub></td><td><b>0.91</b><sup>+0.19</sup><sub>-0.39</sub></td><td><b>1.07</b><sup>+0.07</sup><sub>-0.05</sub></td><td><b>6.86</b><sup>+1.28</sup><sub>-1.40</sub></td><td><b>1.05</b></td></tr></tbody></table>

**Notes.** Errors are the 16 and 84 percentiles. $\alpha = 1.05$ is fixed. log is the decimal logarithm.

pressure drops are a significant feature of the reconstructed 2D Compton profile, we performed exactly the same reconstruction but for simulated SPT clusters injected at random locations. For this test, we adopted the UPP, and performed three injections of the full simulated SPT sample. Results for the simulated clusters are shown in the right hand panel of Fig. 10 right, where the three injections are shown in red, orange and cyan. The test shows that the UPP is well reconstructed up to $\theta = 2\theta_{500}$. However, at higher $\theta$, the reconstructed profiles also display pressure drops. These features are not simulated in the injected profiles, which are smooth. They are thus most probably due to the noise in the outer part of the profile in the current data sets.

According to recent simulation studies (Aung et al. 2021; Baxter et al. 2021; Zhang et al. 2021), a single accretion shock is expected in the average profile of relaxed clusters in the radial range $2 < r < 2.5R_{200m}$. Zenteno et al. (2020) studied the dynamical state of 288 SPT-SZ clusters in the Dark Energy Survey footprint and provided a list of 43 disturbed and 41 relaxed clusters. Among the 43 disturbed (41 relaxed) clusters, all (39) are included in our sample. We measured the average profile for each of the subsamples. The two profiles are unfortunately compatible within their error bars due to the small size of the subsamples.

In conclusion, we do not see any obvious single accretion shock in the outskirts of the average pressure profile. This is most probably due to the fact that we are close to the noise floor in the outer parts of the profile with the current data sets. The fact that

we cannot separate the clusters in our sample according to their dynamical state makes any such study even more difficult. In future studies, it will thus be crucial to work with larger samples, including clusters with well characterised dynamical states, and with deeper data.

## 6. Summary and conclusions

We measured the average pressure profile of 461 SPT clusters in the Planck and SPT-SZ data jointly. The size of the sample is about ten times larger than in previous studies (Planck Collaboration Int. V 2013; Sayers et al. 2013, 2016; Pointecouteau et al. 2021), which used tens of clusters. The optimal combination of space- and ground-based data allowed to constrain simultaneously the pressure profile at small and large scales. For the first time, we proposed to work in Fourier space, and provided results on the average profile in harmonic instead of real space. This approach greatly simplifies PSF deconvolutions and the rescaling of individual profiles before averaging. We obtained very competitive constraints on the cluster pressure profile based on SZ data only, without the need for including X-ray data to constrain the inner part or simulations to constrain the outer part. Our main results are summarised below.

- The Malmquist bias of the SPT clusters in the SPT-SZ maps must be corrected for before combining Planck and SPT-SZ maps.

- The average pressure profile measured from the joint Planck+SPT-SZ dataset is in excellent agreement with the UPP (Arnaud et al. 2010) in the inner part and with the PIPV (Planck Collaboration Int. V 2013) profile in the outer part. A gNFW profile provides a good fit to the data. Fixing $\alpha$ to the UPP value, we were able to constrain four of the gNFW parameters ($\gamma$, $c_{500}$, $\beta$, $P_0$), for the first time in harmonic space from SZ data only.

- We split the sample in two low and high redshift subsamples. The resulting profiles are compatible, thus showing no indication of evolution with redshift.
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Reconstructed Compton profiles (2D real space). Left: reconstructed Compton profile (2D real space) from the data points of Fig. 5 left. The 68% C.L. envelope is shown in red and the mean in black. The solid green line is the PIPV profile. The solid blue line is the best gNFW fit. The dashed blue lines delineate the 68% C.L. of the gNFW fit. In the bottom panel, the ratio is computed with respect to the best fit. Right: same as the left panel, but with simulated UPP profiles injected into actual data, extracted and reconstructed using the same method as on the real clusters. We performed three injections of simulated SPT clusters at different locations. The reconstructed profiles are in red, orange and cyan. The solid bright blue line is the UPP.</figcaption></figure>

- Splitting the sample in two low and high mass subsamples, we find inconsistent profiles. This inconsistency may be of physical origin, with lower mass clusters being on average less peaked in the centre, or may be a residual systematic effect of the analysis where the Malmquist bias being underestimated for the low mass subsample.

- We reconstructed the two dimensional Compton parameter profile in real space to study the outskirts. We found multiple pressure drops at $\theta > \theta_{500}$. Testing this result by injection of smooth simulated cluster samples into real maps, we conclude that these pressure drops are in fact most probably due our being close to the noise floor in the outer parts of the profile with the current data sets. We lack both sensitivity and information on the dynamical state of the clusters in the SPT sample to probe this further.

For practical applications, the average scaled pressure profile from the combination of *Planck* and SPT-SZ observations in units of $R_{500}$ is given by Eq. (7), with best-fitting parameters given in bold in Table 1. For clusters of given mass $M_{500}$ and redshift $z$, the physical pressure profile can then be derived from Eq. (5). The equivalent gNFW parameters for the scaled pressure profile in units of $R_{200m}$ are given in Table 2, which can be used in conjunction with Eq. (C.2) to obtain the physical pressure profile as a function of $M_{200m}$ and $z$. The MCMC chains corresponding to the best fitting model to the average profiles, scaled as a function of $R_{500}$ and $R_{200m}$, are available on request to the authors. Finally, the tabulated best-fitting profiles and associated error envelopes are provided in Appendix E.

The harmonic method presented in the work is promising. Upcoming and future ground-based experiments, such as the Simons Observatory (SO, Ade et al. 2019) and CMB-S4 (Abazajian et al. 2019), will detect around 15 000 and 100 000 objects, respectively. Averaging the cluster profiles in bins of redshift and/or mass in SO/CMB-S4+Planck will allow the study of the cluster pressure profile and its evolution to unprecedented precision.

**Acknowledgements.** We thank the A&A referee, Shutaro Ueda, for useful suggestions and corrections, which helped to clarify and improve this article. We thank Monique Arnaud for valuable discussions at the beginning of this work,

and for comments on the final draft of this article. We thank Rémi Adam, Stefano Ettori, Fabio Gastaldello and Mariachiara Rossetti for useful suggestions while presenting a first version of this work at the NIKA2 conference in June 2021, and Colin Hill, Nick Battaglia and Daisuke Nagai for useful suggestions while presenting a second version of this work at the SZ workshop (Center for Computational Astrophysics, Flatiron Institute) in June 2022. We also thank Etienne Pointecouteau for clarifications and discussions about the PACT 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 by Ankit Rohatgi. GWP acknowledges financial support from CNES, the French space agency.

## References

- Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv e-prints [arXiv:1907.04473]

- Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J Cosmol. Astropart. Phys., 2019, 056

- Anbajagane, D., Chang, C., Jain, B., et al. 2022, MNRAS, 514, 1645

- Andreon, S., Romero, C., Castagna, F., et al. 2021, MNRAS, 505, 5896

- Andreon, S., Romero, C., Aussel, H., et al. 2023, MNRAS, 522, 4301

- Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

- Aung, H., Nagai, D., & Lau, E. T. 2021, MNRAS, 508, 2071

- Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 75

- Baxter, E. J., Adhikari, S., Vega-Ferrero, J., et al. 2021, MNRAS, 508, 1777

- Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

- Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2022, ApJS, 258, 36

- Böhringer, H., Schuecker, P., Pratt, G. W., et al. 2007, A&A, 469, 363

- Bonamente, M., Hasler, N., Bulbul, E., et al. 2012, New J. Phys., 14, 025010

- Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50

- Chown, R., Omori, Y., Aylor, K., et al. 2018, ApJS, 239, 10

- Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108

- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

- Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7

- Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

- He, Y., Mansfield, P., Rau, M. M., Trac, H., & Battaglia, N. 2021, ApJ, 908, 91

- Hurier, G., Adam, R., & Keshet, U. 2019, A&A, 622, A136

- Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7

- Kéruzoré, F., Mayet, F., Pratt, G. W., et al. 2020, A&A, 644, A93

- Lau, E. T., Nagai, D., Avestruz, C., Nelson, K., & Vikhlinin, A. 2015, ApJ, 806, 68
McCarthy, I. G., Le Brun, A. M. C., Schaye, J., & Holder, G. P. 2014, MNRAS, 440, 3645

Melin, J. B., Bartlett, J. G., Tarrío, P., & Pratt, G. W. 2021, A&A, 647, A106

Molnar, S. M., Hearn, N., Haiman, Z., et al. 2009, ApJ, 696, 1640

Mroczkowski, T., Bonamente, M., Carlstrom, J. E., et al. 2009, ApJ, 694, 1034

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Oppizzi, F., De Luca, F., Bourdin, H., et al. 2023, A&A, 672, A156

Plagge, T., Benson, B. A., Ade, P. A. R., et al. 2010, ApJ, 716, 1118

Planck Collaboration X. 2011, A&A, 536, A10

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration XXVII. 2016, A&A, 594, A27

Planck Collaboration Int. V. 2013, A&A, 550, A131

Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73

Pratt, C. T., Qu, Z., & Bregman, J. N. 2021, ApJ, 920, 104

Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 410, 2481

Romero, C. E., Mason, B. S., Sayers, J., et al. 2017, ApJ, 838, 86

Romero, C., McWilliam, M., Macías-Pérez, J. F., et al. 2018, A&A, 612, A39

Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515

Ruppin, F., Adam, R., Comis, B., et al. 2017, A&A, 597, A110

Ruppin, F., Mayet, F., Pratt, G. W., et al. 2018, A&A, 615, A112

Ruppin, F., McDonald, M., Brodwin, M., et al. 2020, ApJ, 893, 74

Sayers, J., Czakon, N. G., Mantz, A., et al. 2013, ApJ, 768, 177

Sayers, J., Golwala, S. R., Mantz, A. B., et al. 2016, ApJ, 832, 26

Sayers, J., Mantz, A. B., Rasia, E., et al. 2023, ApJ, 944, 221

Sehgal, N., Trac, H., Acquaviva, V., et al. 2011, ApJ, 732, 44

Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452

Song, J., Zenteno, A., Stalder, B., et al. 2012, ApJ, 761, 22

Story, K., Aird, K. A., Andersson, K., et al. 2011, ApJ, 735, L36

Sun, M., Sehgal, N., Voit, G. M., et al. 2011, ApJ, 727, L49

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Comment. Astrophys. Space Phys., 2, 66

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comment. Astrophys. Space Phys., 4, 173

Tarrío, P., Melin, J. B., & Arnaud, M. 2019, A&A, 626, A7

Zenteno, A., Hernández-Lang, D., Klein, M., et al. 2020, MNRAS, 495, 705

Zhang, C., Zhuravleva, I., Kravtsov, A., & Churazov, E. 2021, MNRAS, 506, 839

Zhao, H. 1996, MNRAS, 278, 488
# Appendix A: Correlated noise in harmonic space

<figure><img src="image_15.png" /><figcaption>Fig. A.1. Correlation matrix $\mathcal{R}$ between the $l\theta_{500}$ bins of the average cluster pressure profile estimated by bootstrap resampling. The top left half-matrix is estimated at the location of the clusters, the bottom right half is estimated on blank fields.</figcaption></figure>

As discussed in Sect. 3.2, the $l\theta_{500}$ bins are expected to be uncorrelated to first approximation. We tested this hypothesis by computing the covariance matrix $S$ between the bins using bootstrap resampling. In practice, we performed 1,000 bootstraps, drawing with replacement 1,000 SPT-like samples (461 clusters) from the original 461 SPT clusters. For each of the samples, we performed the same analysis as for the original data set. We then estimated the covariance between the $l\theta_{500}$ bins of the 1,000 samples. The result for the correlation matrix $\mathcal{R}$ is shown as the top left half-matrix in Fig. A.1. We recall the definition

$$R_{ll'} = \frac{S_{ll'}}{\sqrt{S_{ll}} \sqrt{S_{ll'}}} \qquad (A.1)$$

In fact, there is a significant contribution of the off-diagonal terms for $10^3 < l\theta_{500} < 10^4$ arcmin, with correlation factors larger than 0.5. To further investigate these correlations, we performed the same analysis on blank fields. The results are shown as the bottom right half-matrix in Fig. A.1. The off-diagonal correlations in the range $10^3 < l\theta_{500} < 10^4$ arcmin are not visible for the blank fields. We notice however that some correlations exist both on cluster and blank fields in the range $5 \times 10^4 < l\theta_{500} < 10^5$ arcmin but are not very significant. As significant correlations in the range $10^3 < l\theta_{500} < 10^4$ arcmin

are only present at the locations of clusters and not on blank fields, this points toward a cluster-related origin of these correlations, such as lensing of the primary CMB anisotropies, or sources behind the clusters.

<figure><img src="image_16.png" /><figcaption>Fig. A.2. Ratio between the diagonal terms of the covariance matrix $\sqrt{S_{ll}}$ and the statistical error $\hat{\sigma}_l$ as a function of $l\theta_{500}$.</figcaption></figure>

We also compared the diagonal terms of the covariance matrix $\sqrt{S_{ll}}$ to the statistical error $\hat{\sigma}_l$ estimated by our averaging procedure detailed in Sect. 3.5. Results are shown in Fig. A.2. As expected, the $\sqrt{S_{ll}}$ are in average higher than the statistical $\hat{\sigma}_l$.

We thus adopted the full covariance matrix $S$ in the likelihood (Sect. 3.6) to be conservative in our gNFW fit.

# Appendix B: The gNFW profile in real and harmonic space

We produced illustrative plots based on variations of the five parameters of the gNFW around the fiducial values of the UPP (Arnaud et al. 2010). These are shown in Fig. B.1. The first column corresponds to the three dimensional pressure profile in real space, the second column to the two dimensional Compton profile in real space that is the three dimensional pressure profile integrated along the line-of-sight, the third column is the two dimensional Fourier transform of the Compton profile expressed in harmonic space. From top to bottom, rows correspond to variations around $P_0, \gamma, \beta, \alpha$ and $c_{500}$.
<figure><img src="image_17.png" /><figcaption>**Fig. B.1.** Pressure profiles. *From left to right:* 3D real space, 2D real space, 2D harmonic space. *From top to bottom:* the gNFW parameters $P_0$, $\gamma$, $\beta$, $\alpha$, $c_{500}$ are varied.</figcaption></figure>
## Appendix C: Self-similar quantities

Following Appendix A of Arnaud et al. (2010), we define the characteristic temperature $kT_{200m} = \mu m_p GM_{200m}/(2R_{200m})$ where $\mu$ is the mean molecular weight, $m_p$ the proton mass and G the gravitational constant. We note $n_{e,200m} = \rho_{g,200m}/(\mu_e m_p)$ the characteristic electronic density where $\rho_{g,200m} = 200 f_B \rho_m(z)$ is the characteristic gas density with $\mu_e$ the mean molecular weight per free electron, $f_B$ the baryon fraction in the Universe and $\rho_m(z)$ the mean density of the Universe at redshift $z$. We define the characteristic pressure in 200m as

$$P_{e,200m} = n_{e,200m} kT_{200m} \qquad (C.1)$$

which can be written as a function of z and $M_{200m}$

$$P_{e,200m} = \frac{3}{8\pi} \left[ \frac{200 G^{-1/4} \Omega_m(z) H(z)^2}{2} \right]^{4/3} \frac{\mu}{\mu_e} f_B M_{200m}^{2/3} \qquad (C.2)$$

Then we define

$$y_{200m} = \frac{\sigma_T}{m_e c^2} \times (2R_{200m}) \times P_{e,200m} \qquad (C.3)$$

as in Eq. 8 but we intentionally do not include any mass dependent factor $M_{500}^{0.12}$ because it has not been constrained for overdensity 200m.

## Appendix D: Impact of positional uncertainty on recovered gNFW parameters

Fig D.1 shows the impact of the positional uncertainty on recovered gNFW parameters for clusters injected in the data assuming the universal pressure profile. Details of the procedure are provided in Sect. 4.4. The shift of the posteriors is negligible with respect to their size.

<figure><img src="image_18.png" /><figcaption>Fig. D.1. Marginalised posterior likelihood for the parameters of the best-fitting gNFW model for clusters injected in the data assuming the universal pressure profile. Blue contours assume the positions of the injected clusters are perfectly known. Orange contours include random positional uncertainties. The shift of the posteriors is negligible with respect to their size.</figcaption></figure>
# Appendix E: Tabulated best-fitting profiles and error envelopes

For each sample and subsample, we provide tabulated best-fitting profiles and associated 68% and 95% C.L. envelopes (upper and lower bounds). The profiles are provided in 3D real space for the pressure (first six columns of each table), in 2D real space for the Compton parameter (following six columns) and in the 2D harmonic space for the harmonic transform of the

Compton parameter (final six columns). The best-fitting pressure profile for each sample and subsample, extracted from the first six columns of each table, is shown in Fig. E.1. The different error envelopes are compared in Fig. E.2.

The complete set of profiles for the full sample are provided in Table E.1, for the low redshift sample in Table E.2, for the high redshift sample in Table E.3, for the low mass sample in Table E.4, for the high mass sample in Table E.5, and for the full sample scaled in $\theta_{200m}$ in Table E.6.

<figure><img src="image_19.png" /><figcaption>**Fig. E.1.** Best fitting pressure profiles in real space (solid lines) and associated 68% (dashed lines) and 95% (dashed dotted lines) C.L. envelopes. From left to right and top to bottom: full sample scaled in $\theta_{500}$, low redshift subsample, low mass subsample, full sample scaled in $\theta_{200m}$, high redshift subsample, high mass subsample.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. E.2.** Best fitting pressure profiles in real space (solid lines) and associated 68% (dark) and 95% (light) C.L. envelopes from the MCMC chains. Left: Full sample (dark blue) compared to the low (cyan) and high (orange) mass subsamples. The bottom panel shows the ratio with respect to the full sample. Right: Same but for the low (cyan) and high (orange) redshift subsamples.</figcaption></figure>