© The Authors 2023

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

$$\begin{array}{c}\hfill y=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{\displaystyle {\int }_{\mathrm{los}}}{P}_{\mathrm{e}}\mathrm{d}l\end{array}$$(1)

where σ_T is the Thomson cross-section, m_e is the electron rest mass, c is the speed of light and P_e = n_e × 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₅₀₀, corresponding to the mass enclosed in a sphere of radius R₅₀₀, whose mean density is 500 times the critical density of the Universe at the cluster redshift z. We note θ₅₀₀ = R₅₀₀/D_ang(z) the angular radius of the cluster (D_ang(z) being the angular distance of the cluster). We also use M_200m, R_200m and θ_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 ≥ R₅₀₀, R₅₀₀). 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₅₀₀) and large (r > R₅₀₀) 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₅₀₀.

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 ∼ 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σ feature at 2R₅₀₀ 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σ) and a decrease in pressure at 4.6R_200m (2.0σ). 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 ΛCDM cosmology with H₀ = 70 km⁻¹ s⁻¹ Mpc⁻¹ and Ω_m = 1 − Ω_Λ = 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¹. 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 tangental maps centered on these clusters. Finally, we removed eleven clusters with Y^0.75 SPT ≥ 2 × 10⁻⁴ arcmin², 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₅₀₀ is provided with the catalogue. Throughout this work, we renormalise the values of M₅₀₀ 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₅₀₀) are based on these rescaled masses. We only used the original SPT masses at four specific places in the 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 θ_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$\stackrel{\prime }{.}$72 pixel⁻¹) in the Planck Legacy Archive² with their associated products, including the frequency response of each band. We upgraded the maps to N_side = 8192 (0$\stackrel{\prime }{.}$43 pixel⁻¹) 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$\stackrel{\prime }{.}$43 pixel⁻¹) 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³.

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 ≥ 2 × 10⁻⁴ arcmin², 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

$$\begin{array}{c}\hfill p(x)=\frac{P_{\mathrm{e}}(r)}{P_{\mathrm{e},500}}\end{array}$$(2)

with the scaled radius

$$\begin{array}{c}\hfill x=\frac{r}{R_{500}}\end{array}$$(3)

and the characteristic electronic pressure

$$\begin{array}{c}\hfill P_{\mathrm{e},500}=\frac{\mu }{{\mu }_{\mathrm{e}}}{P}_{g,500}.\end{array}$$(4)

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

$$\begin{array}{c}\hfill 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{\mathrm{keV}\mathrm{cm}}^{-3}\end{array}$$(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₇₀ the Hubble constant normalised to 70 km s⁻¹ Mpc⁻¹. Arnaud et al. (2010) noticed that p(x) depends on M₅₀₀ and proposed to define (their Eq. (8))

$$\begin{array}{c}\hfill p(x,M_{500})=\mathbb{p}(x){\left[\frac{M_{500}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{0.12},\end{array}$$(6)

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

$$\begin{array}{c}\hfill \mathbb{p}(x)=\frac{P_0}{{(c_{500}x)}^{\gamma }{[1+{(c_{500}x)}^{\alpha }]}^{(\beta -\gamma )/\alpha }}.\end{array}$$(7)

Above, P₀ is the normalisation (independent of the mass), c₅₀₀ = R₅₀₀/r_s is the concentration parameter, r_s is the scale radius, and α, β and γ are the intermediate (r ∼ r_s), outer (r ≫ r_s) and inner (r ≪ r_s) slopes of the profile, respectively. The average scaled profile $\mathbb{p}(x)$ thus has five parameters: Θ = [P₀, c₅₀₀, γ, α, β].

We can now define the characteristic Compton parameter

$$\begin{array}{c}\hfill y_{500}=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}\times (2R_{500})\times P_{\mathrm{e},500}\times {\left[\frac{M_{500}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{0.12}.\end{array}$$(8)

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

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 θ₅₀₀, with θ₅₀₀ 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θ₅₀₀ (harmonic space, l being the multipole) instead of y(θ)/y₅₀₀ versus θ/θ₅₀₀ (real space, θ 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 unlensed 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 𝒮, 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 𝒮 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 × 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 × θ₅₀₀. We left these in the analysis, but checked their impact on the average profile (see Sect. 4.1). We write the tangental maps as

$$\begin{array}{c}\hfill \mathit{m}(\mathit{x})=y_{\mathrm{o}}{\mathit{t}}_{{\mathit{\theta }}_{\mathbf{s}}}(\mathit{x})+\mathit{n}(\mathit{x})\end{array}$$(9)

where m(x) is the column vector corresponding to the tangental maps, y_ot_θs(x) the cluster component, and n(x) the instrumental and astrophysical noise. The quantity y_ot_θs(x) can be decomposed into the central Compton parameter y_o, and the cluster profile t_θs(x), whose ith component is f_ν(ν_i)[b_i * T_θs](x). The quantity T_θs is the two dimensional SZ profile⁴, convolved by b_i (defined hereafter) and multiplied by f_ν, the SZ frequency spectrum integrated over the frequency response. We included the relativistic corrections to the calculation of f_ν using the formula from Itoh et al. (1998) and the cluster temperature–mass relation⁵ 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

$$\begin{array}{c}\hfill \mathit{m}(\mathit{k})=\mathit{a}(\mathit{k})y_{\mathrm{o}}{T}_{{\theta }_{\mathrm{s}}}(\mathit{k})+\mathit{n}(\mathit{k})\end{array}$$(10)

with a(k) the mixing vector whose ith component is f_ν(ν_i)b_i(k). The ILC provides an estimator of y_oT_θs (e.g., Remazeilles et al. 2011)

$$\begin{array}{c}\hfill \widehat{y_{\mathrm{o}}{T}_{{\theta }_{\mathrm{s}}}}(\mathit{k})={[{\mathit{a}}^t(\mathit{k})·{\mathit{P}}^{-1}·\mathit{a}(\mathit{k})]}^{-1}{\mathit{a}}^t(\mathit{k})·\mathit{m}(\mathit{k}),\end{array}$$(11)

with P(k) = ⟨n^*(k)n(k)⟩ being the power spectrum matrix of the noise across frequency channels. In practice P(k) was estimated directly from the maps by masking a circular region of 4θ₅₀₀ radius to exclude the studied cluster. The error on the estimator is given by

$$\begin{array}{c}\hfill \sigma (\mathit{k})={[{\mathit{a}}^t(\mathit{k})·{\mathit{P}}^{-1}·\mathit{a}(\mathit{k})]}^{-1/2}.\end{array}$$(12)

We then divided the profiles $\widehat{{\mathit{y}}_{\mathrm{o}}{T}_{{\theta }_{\mathrm{s}}}}(\mathit{k})$ and associated error σ(k) by the pixel window function w(k) to correct our profile for the 0.43 arcmin pixelisation. We rejected modes 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 correction factor using the formalism in Bleem et al. (2015). We adopted the ‘unbiased SPT-significance’-mass relation

$$\begin{array}{c}\hfill \zeta =A_{\mathrm{SZ}}{\left(\frac{M_{500}}{3\times {10}^{14}{M}_{\odot }{h}^{-1}}\right)}^{B_{\mathrm{SZ}}}{\left(\frac{E(z)}{E(0.6)}\right)}^{C_{\mathrm{SZ}}}\end{array}$$(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 km s⁻¹ Mpc⁻¹. The significance ξ is related to the unbiased significance ζ via

$$\begin{array}{c}\hfill \zeta =\sqrt{{⟨\xi ⟩}^2-3}\end{array}$$(14)

for ζ > 2. For each published SPT-SZ cluster Malmquist bias corrected mass M₅₀₀, and significance ξ, we drew N = 1 000 000 values of ζ_true following Eq. (13) including the log-normal scatter. We then computed ${\xi }_{\mathrm{true}}=\sqrt{{\zeta }_{\mathrm{true}}^2+3}$ for the N values. We added a noise with normal distribution 𝒩(0, 1) to ξ_true to obtain observed significances ξ_obs. We selected the values that verify ξ_obs > 4.5 corresponding to the selection threshold of the SPT-SZ cluster catalogue, and ζ_true > 2 corresponding to the validity range of Eq. (14). We computed the average of ξ_obs in bins of ξ_true for the selected values. We thus obtained a ξ_obs − ξ_true relation for the considered cluster. We interpolated this relation at the published observed significance ξ to obtain the average true significance ξ_t. We adopted the ratio ξ_t/ξ as the correction factor to be applied to each SPT-SZ map. The correction is negligible for ξ > 7 clusters but is important for lower significance clusters, reaching 0.4 at the detection threshold ξ = 4.5. We note that we have not recalibrated the masses M₅₀₀ 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 × 10³ arcmin < lθ₅₀₀ < 2 × 10⁴ arcmin (see Sect. 3.5 below for definition of lθ₅₀₀).

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 × 103 arcmin < lθ500 < 2 × 104 arcmin.

3.5. Averaging procedure

We renormalised each Fourier map obtained after the ILC $\widehat{\mathit{y}}(\mathit{k})=\widehat{{\mathit{y}}_{\mathrm{o}}{T}_{{\theta }_{\mathrm{s}}}}(\mathit{k})/\mathit{w}(k)$ by dividing it by ${\theta }_{500}^2{y}_{500}$. We then bin averaged it in ∥k∥θ₅₀₀ bins. This procedure performs the spatial rescaling of the profile in θ₅₀₀ and the azimutal averaging at the same time. The averaging in bins was performed with inverse variance weights (${\sigma }^2(\mathit{k})/w^2(k)/{\theta }_{500}^4/y_{500}^2$). The profiles binned in ∥k∥θ₅₀₀ 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θ₅₀₀ and we converted it to harmonic space $y_l/y_{500}/{\theta }_{500}^2$ vs. lθ₅₀₀. Errors on ${\stackrel{\sim }{\mathit{y}}}_l={\mathit{y}}_l/{\mathit{y}}_{500}/{\theta }_{500}^2$ are the inverse weighted average errors. We denote these ${\stackrel{\sim }{\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 correspondance 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}(\mathrm{\Theta })=-\frac{{\chi }^2(\mathrm{\Theta })}{2}$ with

$$\begin{array}{c}\hfill {\chi }^2(\mathrm{\Theta })={\displaystyle \sum _{{l{\theta }_{500}}_{\mathrm{bin}},{l^{\prime }{\theta }_{500}}_{\mathrm{bin}}}^{<{l{\theta }_{500}}_{\max }}}[{\stackrel{\sim }{y}}_l-{\stackrel{\sim }{y}}_{\mathrm{mod},l}(\mathrm{\Theta })]{\mathcal{S}}_{ll″}^{-1}[{\stackrel{\sim }{y}}_{l^{\prime }}-{\stackrel{\sim }{y}}_{\mathrm{mod},l^{\prime }}(\mathrm{\Theta })]\end{array}$$(15)

where Θ is the set of the five gNFW profile parameters. The quantity ${\stackrel{\sim }{\mathit{y}}}_{\mathrm{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. 𝒮 is the covariance matrix, estimated as detailed in Appendix A. The quantity ${\stackrel{\sim }{\mathit{y}}}_l$ is the average profile measurement obtained as described in the previous section. The sum is performed over the binned multipoles lθ_500bin, l′θ_500bin lower than lθ_500max.

For the fit to the full sample, we fixed lθ_500max = 8 × 10⁴ arcmin, a value chosen to avoid the inclusion of systematics at high lθ₅₀₀. Figure 2 (top) shows the number of pixels (resp. clusters) contributing to the measurement in a given lθ_500bin as the red (resp. blue) histogram. At low lθ₅₀₀ values (around 10² arcmin), the number of contributing clusters is around 100. As lθ₅₀₀ 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θ₅₀₀ (the annulus contains more pixels when lθ₅₀₀ increases). Around lθ₅₀₀ = 10³ 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θ₅₀₀ = 2 × 10⁴ arcmin, the number of clusters decreases along with the number of pixels. At lθ₅₀₀ = 8 × 10⁴ 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θ₅₀₀ bins above lθ_500max = 8 × 10⁴ arcmin in our likelihood for the full sample. Using the same criteria to fix lθ_500max for the subsamples, we find the same value of lθ_500max = 8 × 10⁴ arcmin for the high mass and low z subsamples.

Fig. 2. Top: number of pixels and clusters contributing to measurement of the harmonic transform Compton parameter profile in a given lθ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θ500 bin divided by the square root of the number of clusters contributing to this bin.

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: 𝒰_{[ − 3, 2]} for log(P₀), 𝒰_[0.01, 3] for c₅₀₀ (or c₂₀₀), 𝒰_{[0.01, 2.5]} for γ and 𝒰_[1, 20] for β. For the low mass sample, we also adopted flat priors but modified the upper bounds adopting 𝒰_{[ − 3, 1.5]} for log(P₀), 𝒰_[0.01, 2] for c₅₀₀, 𝒰_{[0.01, 1.5]} for γ and 𝒰_[1,75] for β. 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 multiplied (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θ₅₀₀ < 2500 arcmin). Conversely, on small scales (lθ₅₀₀ > 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 × θ₅₀₀ from the average modifies only weakly the intermediate and external part of the profile (lθ₅₀₀ < 10⁴ arcmin), and always within the statistical errors.

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.

4.2. Best-fitting gNFW profile

We now discuss the best-fitting p arameters 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 γ, α, β, and did not allow meaningful constraints to be put on all five parameters. We therefore fixed the intermediate slope to the UPP value (α = 1.05) and left the other four parameters (amplitude log(P₀), concentration parameter c₅₀₀, internal slope γ and external slope β) free. Results are shown in Fig. 4. Fixing α 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 χ² of the fit is good, with χ² = 92.96 and 68 degrees of freedom (72 data points minus four parameters), leading to a reduced χ²/d.o.f. = 1.37.

Fig. 4. Marginalised posterior likelihood for the parameters of the best-fitting gNFW model (Eq. (7)) with fixed α = 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.

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θ₅₀₀ (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.

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θ500 (the central regions), and to the PIPV profile at small lθ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.

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 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θ₅₀₀, but they are larger than the statistical errors in the intermediate lθ₅₀₀ range (2 × 10³ arcmin < lθ₅₀₀ < 2 × 10⁴ arcmin). This effect is an indication of some intrinsic scatter (σ_int) in the profiles in the intermediate lθ₅₀₀ range. We tried to measure it using the approach described in Sect. 5.3 of Planck Collaboration X (2011). We computed the raw scatter σ_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θ₅₀₀ bin under consideration) as the dashed blue line in the bottom panel of Fig. 2. We also computed σ_stat from the square root of the mean of the variances of the individual profiles. Unfortunately, as this figure shows, σ_raw is of the same order of magnitude as σ_stat across the full lθ₅₀₀ range, and the estimated difference ${\sigma }_{\text{int}}^2={\sigma }_{\text{raw}}^2-{\sigma }_{\text{stat}}^2$ is compatible with zero. This estimator of the intrinsic scatter σ_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)

$$\begin{array}{c}\hfill \sqrt{{\theta }_{\mathrm{beam}}^2+{(k{\theta }_{\mathrm{c}})}^2}/\xi .\end{array}$$(16)

Here, θ_beam = (1.1 + 1.6)/2 arcmin is the mean of the FWHM of the 95 and 150 GHz original beams, k = 2, θ_c is the cluster core radius, and ξ 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θ₅₀₀) but is slightly below in the outer part (at low lθ₅₀₀). 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.

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.

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 × 10³ arcmin < lθ₅₀₀ < 4 × 10⁴ arcmin. Sayers et al. (2016) used Planck and Bolocam data for 47 massive clusters to constrain the amplitude P₀ and outer slope β 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 (α) and determined the other four parameters (β, γ, c₅₀₀, log(P₀)). 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θ₅₀₀ > 10⁴ 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θ₅₀₀ 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 the Planck and SPT-SZ selections differ, in the sense that they do not select the same fraction of disturbed clusters.

Fig. 7. Redshift subsamples. Left: average pressure profile for z < 0.6 (cyan circles), and z ≥ 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.

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θ₅₀₀, while the high z subsample lacks data points at lθ₅₀₀ > 2 × 10⁴ 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θ_500max < 2 × 10⁴ 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₅₀₀ = 3.6 × 10¹⁴ M_⊙. This mass is close to the median mass of the sample (M₅₀₀ = 3.58 × 10¹⁴ M_⊙). The low mass subsample contains 233 clusters with a median mass of M₅₀₀ = 3.15 × 10¹⁴ M_⊙. The high mass subsample contains 228 clusters with a median mass of M₅₀₀ = 4.46 × 10¹⁴ M_⊙. 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θ₅₀₀ range. The profile of the high mass sample is slightly above the full sample gNFW fit.

Fig. 8. Mass subsamples. Left: average pressure profile for M500 < 3.6 × 1014 M⊙ (cyan circles), and M500 ≥ 3.6 × 1014 M⊙ (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.

We fitted the subsample profiles as described in Sect. 3.6, limiting the fit to lθ_500max < 6 × 10⁴ 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₀), the profile of the low mass subsample has a lower γ 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₀). 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 γ, log(P₀) 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_e, 500 value which would drive the fit to a lower P₀. 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.

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θ₅₀₀ < 10³ 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₅₀₀) 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 β = 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 β, but by a combination of c₅₀₀ and β, as can be seen in the first column of Fig. B.1. The degeneracy between c₅₀₀ and β for our profile can clearly be seen in Fig. 4. Imposing the UPP value c₅₀₀ = 1.18 would lead to a best-fitting value for β 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₅₀₀ versus β 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θ₅₀₀ < 10³ arcmin.

5.5.2. Scaling in θ_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 θ_200m instead of θ₅₀₀. We converted SPT masses M₅₀₀ 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₅₀₀ since the profile rescaled with R_200m will not be compared to other works with different mass definition. We renormalised the profiles by $y_{\text{200m}}{\theta }_{\text{200m}}^2$ instead of $y_{500}{\theta }_{500}^2$, and we rescaled the radial quantities by θ_200m instead of θ₅₀₀. Details of the derivation of y_200m and associated characteristic pressure P_e, 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 χ²/d.o.f. = 1.42.

Fig. 9. Pressure profile scaled by θ200m and best fitting model. Left: same as Fig. 5 Left but averaged according to θ200m instead of θ500. Right: same as Fig. 4 but for the profile averaged according to θ200m. Values for the gNFW best-fitting parameters are given in Table 2.

We provide the tabulated best-fitting profile and associated error envelope in harmonic and real space for the profile scaled by θ_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 𝒮. 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.

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.

The reconstructed profile follows closely the gNFW fit up to θ ∼ 2θ₅₀₀, and does not show a single feature. Instead, multiple pressure drops are visible at θ > 2θ₅₀₀. To test if these multiple 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 θ = 2θ₅₀₀. However, at higher θ, 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 sates, 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 α to the UPP value, we were able to constrain four of the gNFW parameters (γ, c₅₀₀, β, P₀), 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.

• 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 θ > θ₅₀₀. 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₅₀₀ is given by Eq. (7), with best-fitting parameters given in bold in Table 1. For clusters of given mass M₅₀₀ 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₅₀₀ 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.

Acknowledgments

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

Appendix A: Correlated noise in harmonic space

As discussed in Sect. 3.2, the lθ₅₀₀ bins are expected to be uncorrelated to first approximation. We tested this hypothesis by computing the covariance matrix 𝒮 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θ₅₀₀ bins of the 1,000 samples. The result for the correlation matrix ℛ is shown as the top left half-matrix in Fig. A.1. We recall the definition

$$\begin{array}{c}\hfill {\mathcal{R}}_{l{l}^{\prime }}=\frac{{\mathcal{S}}_{l{l}^{\prime }}}{\sqrt{{\mathcal{S}}_{\mathit{ll}}}\sqrt{{\mathcal{S}}_{l^{\prime }{l}^{\prime }}}}.\end{array}$$(A.1)

Fig. A.1. Correlation matrix ℛ between the lθ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.

In fact, there is a significant contribution of the off-diagonal terms for 10³ < lθ₅₀₀ < 10⁴ 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³ < lθ₅₀₀ < 10⁴ 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 × 10⁴ < lθ₅₀₀ < 10⁵ arcmin but are not very significant. As significant correlations in the range 10³ < lθ₅₀₀ < 10⁴ 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.

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

Fig. A.2. Ratio between the diagonal terms of the covariance matrix $\sqrt{{\mathcal{S}}_{\mathit{ll}}}$ and the statistical error ${\stackrel{\sim }{\sigma }}_l$ as a function of lθ500.

We thus adopted the full covariance matrix 𝒮 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₀, γ, β, α and c₅₀₀.

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 P0, γ, β, α, c500 are varied.

Appendix C: Self-similar quantities

Following Appendix A of Arnaud et al. (2010), we define the characteristic temperature kT_200m = μm_pGM_200m/(2R_200m) where μ is the mean molecular weight, m_p the proton mass and G the gravitational constant. We note n_e, 200m = ρ_g, 200m/(μ_em_p) the characteristic electronic density where ρ_g, 200m = 200f_Bρ_m(z) is the characteristic gas density with μ_e the mean molecular weight per free electron, f_B the baryon fraction in the Universe and ρ_m(z) the mean density of the Universe at redshift z. We define the characteristic pressure in 200m as

$$\begin{array}{c}\hfill P_{e,200\mathrm{m}}=n_{e,200\mathrm{m}}k{T}_{200\mathrm{m}}\end{array}$$(C.1)

which can be written as a function of z and M_200m

$$\begin{array}{c}\hfill P_{e,200\mathrm{m}}=\frac{3}{8\pi }{\left[\frac{200G^{-1/4}{\mathrm{\Omega }}_{\mathrm{m}}(z)H{(z)}^2}{2}\right]}^{4/3}\frac{\mu }{{\mu }_e}{f}_{\mathrm{B}}{M}_{200\mathrm{m}}^{2/3}.\end{array}$$(C.2)

Then we define

$$\begin{array}{c}\hfill y_{200\mathrm{m}}=\frac{{\sigma }_{\mathrm{T}}}{m_e{c}^2}\times (2R_{200\mathrm{m}})\times P_{e,200\mathrm{m}}\end{array}$$(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.

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.

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.

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 θ500, low redshift subsample, low mass subsample, full sample scaled in θ200m, high redshift subsample, high mass subsample.

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.

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 θ_200m in Table E.6.

All Tables

All Figures

	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 × 103 arcmin < lθ500 < 2 × 104 arcmin.
In the text

Fig. 2. Top: number of pixels and clusters contributing to measurement of the harmonic transform Compton parameter profile in a given lθ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θ500 bin divided by the square root of the number of clusters contributing to this bin.

In the text

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.

In the text

	Fig. 4. Marginalised posterior likelihood for the parameters of the best-fitting gNFW model (Eq. (7)) with fixed α = 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.
In the text

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θ500 (the central regions), and to the PIPV profile at small lθ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.

In the text

	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.
In the text

Fig. 7. Redshift subsamples. Left: average pressure profile for z < 0.6 (cyan circles), and z ≥ 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.

In the text

Fig. 8. Mass subsamples. Left: average pressure profile for M500 < 3.6 × 1014 M⊙ (cyan circles), and M500 ≥ 3.6 × 1014 M⊙ (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.

In the text

	Fig. 9. Pressure profile scaled by θ200m and best fitting model. Left: same as Fig. 5 Left but averaged according to θ200m instead of θ500. Right: same as Fig. 4 but for the profile averaged according to θ200m. Values for the gNFW best-fitting parameters are given in Table 2.
In the text

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.

In the text

	Fig. A.1. Correlation matrix ℛ between the lθ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.
In the text

	Fig. A.2. Ratio between the diagonal terms of the covariance matrix $\sqrt{{\mathcal{S}}_{\mathit{ll}}}$ and the statistical error ${\stackrel{\sim }{\sigma }}_l$ as a function of lθ500.
In the text

	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 P0, γ, β, α, c500 are varied.
In the text

	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.
In the text

	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 θ500, low redshift subsample, low mass subsample, full sample scaled in θ200m, high redshift subsample, high mass subsample.
In the text

	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.
In the text