© The Authors 2025

1. Introduction

Galaxy clusters constitute the most massive collapsed halos in the Universe. Studying their abundance as a function of redshift and mass provides insights into structure formation history and therefore serves as a powerful tool for constraining cosmological models (e.g., White et al. 1993; Haiman et al. 2001; Vikhlinin et al. 2009; Mantz et al. 2010; Planck Collaboration XIII 2016; Chiu et al. 2023; Bocquet et al. 2024a). The ability to accurately measure cluster masses plays an important role in cluster cosmological studies, as it enables constraints on the rate of cosmic structure growth, the dark energy equation of state, and other cosmological parameters, such as the amplitude of matter fluctuations and the matter density parameter. The development of robust weak lensing (WL) and cosmic microwave background (CMB) lensing has informed mass calibration techniques (Becker & Kravtsov 2011; von der Linden et al. 2014; Dietrich et al. 2019; Zubeldia & Challinor 2019; Grandis et al. 2021; Bocquet et al. 2024b), and the availability of associated high quality WL datasets from, for example, the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP), the Dark Energy Survey (DES), and the Kilo-Degree Survey (KiDS) have set the stage for progress in constraining the standard ΛCDM and wCDM parameters (Costanzi et al. 2019; Abbott et al. 2020; Costanzi et al. 2021; To et al. 2021; Chiu et al. 2023; Bocquet et al. 2024a; Ghirardini et al. 2024) as well as model extensions, including the modification of general relativity and interacting dark matter (e.g., Mantz et al. 2014; Cataneo et al. 2015; Vogt et al. 2024; Mazoun et al. 2024).

These same WL datasets can be employed to study the matter distribution within galaxy clusters. A challenge is that in the existing WL datasets based on large photometric surveys, the matter profiles of individual clusters often have a low signal-to-noise ratio (S/N). Combining WL matter profiles from multiple galaxy clusters provides a way to improve the S/N and reduce the intrinsic variations in the matter distribution from cluster to cluster that arise from their different formation histories. Previous works have employed WL measurements of multiple clusters to constrain cluster masses by combining tangential shear profiles or projected matter profiles of clusters (Oguri & Takada 2011; Umetsu et al. 2016; McClintock et al. 2018; Bellagamba et al. 2019; Giocoli et al. 2021; Lesci et al. 2022). A challenge in this approach is that there are systematic variations in the projected matter profiles of galaxy clusters with cluster mass and redshift. An average WL matter profile therefore reflects the characteristics of the cluster sample, and it depends sensitively on the distribution of the sample in mass and redshift. This approach also requires careful modeling of the spatial distribution and masking of the WL source galaxies on a cluster-by-cluster basis to enable accurate modeling of the multi-cluster matter profile.

If the systematic variations of the matter profiles with mass and redshift can be accurately characterized, then they can also be scaled out, enabling average matter profiles of high S/N that are largely independent of the characteristics of the cluster sample from which they are constructed. In particular, if cluster matter profiles are approximately self-similar in nature – that is, they exhibit similar shapes that vary systematically with mass and redshift – then these systematic trends can be easily removed. In this limit, the need to accurately model the spatial distribution and masking of WL source galaxies is also no longer required.

Approximate self-similarity is a generic prediction of gravitational structure formation (Kaiser 1986). In N-body simulations, cluster halos are well described by so-called Navarro, Frenk, and White (NFW) models (Navarro et al. 1997) that exhibit weak trends in concentration or shape with mass and redshift. In hydrodynamical simulations, self-similar behavior has been seen in cluster gas profiles (Lau et al. 2015) and pressure profiles (Nelson et al. 2014). Observationally, approximate self-similarity has been demonstrated in the intracluster medium (ICM) density, pressure profiles, and temperature profiles (Vikhlinin et al. 2006; Arnaud et al. 2010; Baldi et al. 2012; McDonald et al. 2014), whereas WL studies of cluster matter profiles have tended to focus on whether the NFW model is a good description of the real matter profiles (e.g., Umetsu et al. 2014; Niikura et al. 2015). In particular, Umetsu et al. (2014) studied average matter profiles and found the NFW profile to be an excellent description of real clusters. In contrast, Niikura et al. (2015) analyzed scaled average matter profiles by rescaling with the scale radius of the NFW profile and with the critical density of the Universe and found evidence of a single underlying or universal cluster matter profile.

In this analysis, we use hydrodynamic structure formation simulations and direct WL observations of cluster samples to examine cluster matter profiles, which reveal the remarkable consistency and approximate self-similarity of the simulated and real matter profiles. We exploit this self-similarity to study high S/N cluster matter profiles and then employ them to perform a calibration of observable-mass relations. This mass calibration approach offers a computationally more efficient technique to analyze large cluster WL datasets compared to a cluster-by-cluster approach (Bocquet et al. 2019, 2024b) without loss of information.

The paper is organized as follows. We present the simulated and observed dataset in Sect. 2. The self-similarity of galaxy cluster matter profiles is explored in Sect. 3. The mass calibration method along with the likelihood calculation and the hydrodynamical model is discussed in Sect. 4. In Sect. 5, we validate the analysis method using mock data and present the results using the South Pole Telescope (SPT) clusters and DES-WL data. We conclude with a summary and outlook in Sect. 6. Throughout the paper, we employ a flat ΛCDM cosmology with parameters Ω_m = 0.3 and h = 0.7. All uncertainties are quoted at the 68% credible interval unless otherwise specified.

2. Data

In this section, we first describe the data used in our work: (1) the SPT cluster catalogs and (2) DES (Year 3, hereafter Y3) WL and photometric redshift measurements. Then we summarize the simulation datasets from Magneticum and IllustrisTNG, which are used to explore the impact of baryons on cluster matter profiles.

2.1. South Pole Telescope cluster catalogs

We use a combination of three thermal Sunyaev-Zel’dovich effect (tSZE) selected cluster catalogs that have been extracted from surveys carried out by the South Pole Telescope (Carlstrom et al. 2011) collaboration: SPT-SZ (Bleem et al. 2015; Klein et al. 2024), SPTpol ECS (Bleem et al. 2020) and SPTpol-500d (Bleem et al. 2024). The SPT-SZ survey covers 2500 deg², and the SPTpol ECS survey spans 2770 deg² in the southern sky, while the SPTpol-500d survey pushes to a greater depth within a 500 deg² patch inside the SPT-SZ survey. Galaxy cluster candidates are selected from the mm-wave maps at 90 and 150 GHz using a matched filter technique (Melin et al. 2006), which employs galaxy cluster tSZE models with a range of angular scales (Vanderlinde et al. 2010). Only cluster candidates at redshifts z >  0.25 are considered, because clusters at lower redshifts are larger in angular extent and therefore more strongly impacted by the matched filtering, which is designed to remove atmospheric noise as well as increased noise contributions from the primary CMB. At low redshift the angular scales filtered out overlap with the scales important for the galaxy cluster tSZE, strongly impacting the candidate detection significance and thereby complicating its use as a cluster halo mass proxy. Additionally, we only analyze clusters with a redshift z <  0.95 due to the depth and systematics of the DES WL sample described below.

These cluster candidates are then studied using the Multi-Component Matched Filter cluster confirmation tool (MCMF; Klein et al. 2018). This processing results in a cluster catalog that includes measurements of optical richnesses $\widehat{\lambda }$, sky positions and redshifts. The measured optical richness allows for efficient removal of contaminants from the tSZE candidate list by evaluating the likelihood of each candidate being a random superposition of a physically unassociated optical system with a tSZE noise fluctuation (Klein et al. 2024). The exclusion threshold corresponds to an observed richness threshold that varies with redshift ${\widehat{\lambda }}_{\min }(z)$ (Klein et al. 2019) and has been determined by analyzing the richness distributions along random lines of sight within the survey. The final MCMF-confirmed cluster catalogs have a constant contamination fraction at all redshifts.

The selection threshold in the tSZE detection significance is $\widehat{\zeta }>4.25$ for SPTpol-500d, $\widehat{\zeta }>4.5$ for SPT-SZ and $\widehat{\zeta }>5$ for SPTpol ECS, while the MCMF selection threshold ${\widehat{\lambda }}_{\min }(z)$ is adjusted to maintain a contamination fraction of < 2% in the final MCMF-confirmed cluster lists from both surveys. Figure 1 shows the distribution of observed richness and tSZE detection significance as a function of redshift for the MCMF confirmed SPT sample we study here.

Fig. 1. Distribution of observed tSZE detection significance $\widehat{\zeta }$ and richness $\widehat{\lambda }$ as a function of redshift for the 698 galaxy clusters in the SPT sample that overlap the DES region. The solid line in the top figure shows the detection threshold for the three survey regions. In the bottom figure, the colored lines correspond to the ${\widehat{\lambda }}_{\min }(z)$ detection threshold for each survey.

2.2. DES Y3 lensing

The Dark Energy Survey is a photometric survey in five broadband filters (grizY) which covers an area of ∼5000 deg² in the southern sky. The survey was conducted using the Dark Energy Camera (DECam; Flaugher et al. 2015) at the 4m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO) in Chile. In this work, we use WL data from the first three years of DES observations (DES Y3), which cover the entire 5000 deg² survey footprint.

The DES Y3 shape catalog (Gatti et al. 2021) is constructed from the r, i, z-bands using the METACALIBRATION pipeline (Huff & Mandelbaum 2017; Sheldon & Huff 2017). Other DES Y3 works contain detailed information on the Point-Spread Function modeling (Jarvis et al. 2021), the photometric dataset (Sevilla-Noarbe et al. 2021), and image simulations (MacCrann et al. 2022). After all source selection cuts, the shear catalog consists of roughly 100 million galaxies over an area of 4143 deg². The typical source density is 5 to 6 arcmin⁻², depending on the selection choices of a specific analysis.

Our work follows the selection of lensing source galaxies in four tomographic bins (Fig. 2 shows the three tomographic bins used in our analysis) as employed in the DES 3 × 2pt analysis (Abbott et al. 2022). The selection is defined and calibrated in Myles et al. (2021) and Gatti et al. (2022), where source redshifts are estimated using Self-Organizing Maps Photo-z (SOMPZ). The final calibration accounts for the (potentially correlated) systematic uncertainties in source redshifts and shear measurements. For each tomographic source bin, the mean redshift distribution is provided, and the systematic uncertainties are captured using 1000 realizations of the source redshift distribution. The amplitude of the source redshift distribution is scaled by a factor 1 + m to account for the multiplicative shear bias m. In addition to the tomographic bins and SOMPZ, we used Directional Neighbourhood Fitting (DNF; De Vicente et al. 2016) galaxy photo-z estimates when determining the expected fraction of the lensing source galaxy population in each tomographic bin that is contributed by member galaxies from a particular cluster of interest – the so-called cluster member contamination.

Fig. 2. Dark Energy Survey Y3 lensing source redshift distribution for tomographic bins 2 through 4 that are used in this analysis. The solid line represents the mean and the shaded region depicts the 2σ uncertainties on the redshift distributions.

2.3. Hydrodynamical simulations

In this work, we use the Magneticum Pathfinder suite of cosmological hydrodynamical simulations (Hirschmann et al. 2014; Teklu et al. 2015; Beck et al. 2016; Bocquet et al. 2016; Dolag et al. 2017). We use Box1, which has a box size of 896 h⁻¹ Mpc on a side with 2 × 1526³ particles and the particle mass 1.3 × 10¹⁰ h⁻¹M_⊙ for dark matter particles, and 2.6 × 10⁹ h⁻¹M_⊙ for gas particles. The simulation is run with cosmological parameters (Ω_m = 0.272, Ω_b = 0.0457, H₀ = 70.4, n_s = 0.963, σ₈ = 0.809), which correspond to the WMAP7 constraints for a spatially flat ΛCDM model (Komatsu et al. 2011). From this simulation, we use snapshots at five redshifts z_snap = 0.01, 0.25, 0.47, 0.78, 0.96.

In addition, we also use the data from IllustrisTNG300-1 (Pillepich et al. 2018; Marinacci et al. 2018; Springel et al. 2018; Nelson et al. 2018; Naiman et al. 2018; Nelson et al. 2019). These include 2 × 2500³ resolution elements for a box size of 205 h⁻¹ Mpc on a side. The cosmology corresponds to the Planck2015 constraints for a spatially flat ΛCDM cosmology (Planck Collaboration XIII 2016): Ω_m = 0.3089, Ω_b = 0.0486, σ₈ = 0.8159, n_s = 0.9667, and h = 0.6774. We use snapshots correponding to redshift z_snap ∈ {0.01, 0.24,  0.42,  0.64,  0.95}. From these simulation snapshots, we then extract halos with M_200c >  3 × 10¹³ h⁻¹ M_⊙. Shear maps are generated following Grandis et al. (2021) in a cylinder with a projection depth of 20 h⁻¹ Mpc.

3. Self-similarity in cluster matter profiles

Gravitational lensing is the phenomenon through which photon geodesics are perturbed by gravitational potentials. For a distant galaxy, this causes a distortion in the observed image relative to its true shape (Schneider 2006). In this work, we are interested in WL, where distortions in source galaxy images induced by intervening matter along the line of sight are small. In this regime, the WL signal must be extracted through statistical correlations of source galaxies. The observable of interest in this context is the reduced shear, which is defined as

$$\begin{array}{c}\hfill g=\frac{\gamma }{1-\kappa },\end{array}$$(1)

where γ is the WL shear and κ is the WL convergence (for detailed discussion see Schneider 2006). The ensemble averaged source ellipticity, e, and the shear response, R_γ, are related to the reduced shear as

$$\begin{array}{c}\hfill ⟨g⟩={⟨R_{\gamma }⟩}^{-1}⟨e⟩.\end{array}$$(2)

The term R_γ is the average response of the measured ellipticity to a shear. Due to instrumental and atmospheric effects and noise this shear response typically is less than one. The tangential reduced shear profile induced by an object with a projected mass distribution Σ(R) is related to the critical surface mass density Σ_crit by

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{\mathrm{crit}}{\gamma }_{\mathrm{t}}(R)=⟨\mathrm{\Sigma }(<R)⟩-\mathrm{\Sigma }(R),\end{array}$$(3)

where Σ_crit depends on the geometry of the source-lens system and is defined as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}({\mathrm{z}}_{\mathrm{s}},{\mathrm{z}}_{\mathrm{l}})=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{l}}{D}_{\mathrm{ls}}},\end{array}$$(4)

where z_s and z_l are the source and lens redshifts, respectively, and D_s, D_l, D_ls are the angular diameter distances to the source, lens, and between the source-lens pair. When z_s ≤ z_l, ${\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}$ is defined to be zero.

Analogous to Eq. (3), we introduce an observed quantity

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{reduced}}(R)={\mathrm{\Sigma }}_{\mathrm{crit}}{g}_{\mathrm{t}}(R).\end{array}$$(5)

In this paper, we use the excess surface mass density as defined in the above equation but refer to it as ΔΣ(R). We note that, although Σ(R) is the more fundamental quantity, we use ΔΣ(R) (calculated using Eq. (5)) in our work because it can be directly measured through WL. For convenience we refer to ΔΣ(R) simply as the cluster matter profile.

For a given source redshift distribution P(z_s), we can compute the average lensing efficiency for a given lens as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}({\mathrm{z}}_{\mathrm{l}})={\displaystyle \int {\mathrm{dz}}_{\mathrm{s}}}P({\mathrm{z}}_{\mathrm{s}}){\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}({\mathrm{z}}_{\mathrm{s}},{\mathrm{z}}_{\mathrm{l}}).\end{array}$$(6)

From Eq. (3) we can see that the differential surface mass density at a given projected radius R can be expressed as the difference between the mean enclosed surface mass density and the surface mass density Σ at that projected radius, which is expressed as follows

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)={\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}\chi \rho (\sqrt{R^2+{\chi }^2}),\end{array}$$(7)

$$\begin{array}{c}\hfill ⟨\mathrm{\Sigma }(<R)⟩=\frac{2}{R^2}{\displaystyle {\int }_0^R}\mathrm{d}R\prime R\prime \mathrm{\Sigma }(R\prime ).\end{array}$$(8)

where ρ(r) is the density distribution of the halo and χ is the comoving distance along the line of sight. For the shear signal induced by a halo of mass M, the average excess three-dimensional matter density is given by

$$\begin{array}{c}\hfill \rho (r)={\rho }_{\mathrm{m}}{\xi }_{\mathrm{hm}}(r|M),\end{array}$$(9)

where ρ_m = Ω_m, 0ρ_crit, 0(1 + z)³ is the mean matter density of the universe and ξ_hm(r|M) is the halo-matter correlation function at the halo redshift. At a small radius, ξ_hm(r|M) is dominated by the cluster density profile, and this region is called the “one-halo” region. At a larger radius, most of the contribution comes from correlated structures around the halo, and it is therefore referred to as the “two-halo” region. In this work, we examine both regions but focus on the one-halo region for the cluster mass calibration.

There is a strong theoretical expectation that, barring the impact of baryonic effects, the one-halo region of a halo should be described by the NFW model Navarro et al. (1996, 1997). In this model, the cluster matter profile within the radius r_200c, which encloses a region with a mean density that is 200 times the critical density ρ_crit is well described as

$$\begin{array}{c}\hfill \rho (r)={\delta }_{\mathrm{s}}{\rho }_{\mathrm{crit}}{\left[\frac{r}{c{r}_{200\mathrm{c}}}(1+{\left(\frac{r}{c{r}_{200\mathrm{c}}}\right)}^2)\right]}^{-1},\end{array}$$(10)

where δ_s is a characteristic overdensity depending on c, which is the halo concentration parameter. Such a halo characterized by r_200c has a mass that can be expressed as

$$\begin{array}{c}\hfill M_{200\mathrm{c}}=200{\rho }_{\mathrm{crit}}(z)\frac{4\pi }{3}{r}_{200\mathrm{c}}^3.\end{array}$$(11)

This underlying density profile implies a particular projected ΔΣ matter profile (Bartelmann 1996), whose amplitude scales with the extent of the cluster along the line of sight (i.e., r_200c) and depends on the critical density, which varies with redshift as ρ_crit = 3H²(z)/8πG, where H(z) is the Hubble parameter at redshift z and G is the Gravitational constant. For the projected profiles discussed below, we rename r_200c to be R_200c, corresponding to the projected distance equal to the 3D radius that encloses the mass M_200c.

Within simulations the halo shapes vary from cluster to cluster due to formation history differences, and systematic trends in concentration with mass and redshift have been identified (Bhattacharya et al. 2013; Covone et al. 2014). But in the limit that the systematic variation in c with mass and redshift is small, the average projected matter profiles would have the same shape in the space of R/R_200c. The amplitudes of these projected matter profiles would scale as R_200cρ_crit. This suggests a rescaled projected matter profile, $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$, that would allow for easy exploration of departures from self-similarity:

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(\frac{R}{R_{200\mathrm{c}}},z)=\frac{\mathrm{\Delta }\mathrm{\Sigma }(\frac{R}{R_{200\mathrm{c}}},z)}{R_{200\mathrm{c}}{\rho }_{\mathrm{crit}}(z)}.\end{array}$$(12)

We note here that for any given R_Δ, where Δ represents the overdensity with respect to the critical or mean background densities, one can create a rescaled profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$. If one chooses mean background density ⟨ρ(z)⟩, which corresponds to a scale radius of R_200m rather than R_200c, the rescaling in amplitude would have to be adjusted to follow the correct redshift evolution of the mean background density ⟨ρ(z)⟩ = Ω_mρ_crit(z = 0)(1+z)³.

Finally, because observationally determined cluster centers are not perfect tracers of the true halo center, one must consider also the impact that this mis-centering will have on the observed matter profile. When considering a mis-centering radius R_mis the azimuthal average of the surface mass density can be expressed as

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R,R_{\mathrm{mis}})=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\theta \mathrm{\Sigma }(\sqrt{R^2+R_{\mathrm{mis}}^2-2R{R}_{\mathrm{mis}}\cos \theta }).\end{array}$$(13)

Generally speaking, the mis-centering effects when using optical centers (MCMF adopts the brightest cluster galaxy (BCG) position or the center of mass of the red galaxy distribution if the BCG is significantly offset from that; Klein et al. 2024) or X-ray or even tSZE centers, the impact of mis-centering has only a minor impact outside the inner core region of the cluster. We return to this issue in Sect. 4.4.1.

3.1. Average matter profiles: Hydrodynamical simulations

To enable our study of the average cluster matter profiles, we extracted cluster ΔΣ profiles following the method described in Grandis et al. (2021) for each redshift for both the IllustrisTNG and Magneticum simulations. In total, we extracted 903, 852, 780, 684, and 528 cluster matter profiles at redshifts of 0.01, 0.25, 0.47, 0.78, and 0.96 respectively. In the absence of measurement uncertainties, we constructed average ΔΣ profiles for further study using

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R_j)=\frac{1}{N}{\displaystyle \sum _i}\mathrm{\Delta }{\mathrm{\Sigma }}_i(R_j),\end{array}$$(14)

where the sum is over the N clusters i in the sample and R_j is the radial binning adopted for the cluster matter profiles.

3.1.1. ΔΣ(R) dependence on mass and redshift

To start, we compute the average matter profiles within each of the five redshifts. We compare the average matter profiles ΔΣ(R) at different redshifts using the Magneticum sample in the left panel of Fig. 3. To ensure we are seeing only trends with redshift, we average a narrow range of cluster mass (14.65 <  log(M/M_⊙)< 14.85) for all of the redshifts. The average cluster matter profiles show a significant dependence on redshift.

Fig. 3. Average cluster matter profiles ΔΣ(R) in the Magneticum simulation at five redshifts (color-coded on left) for the same halo-mass bin, 14.65 <  log(M200c/M⊙)< 14.85, and for three mass and redshift bins (color and line-type coded on right). The dependence of the matter profiles on cluster mass and redshift is clearly visible. The thickness of the lines represents the 68% credible region in the average matter profiles.

In the right panel of Fig. 3 we plot average matter profiles ΔΣ(R) at three redshifts of 0.25, 0.47, and 0.78. Each redshift bin is divided into three mass bins (14.30 <  log(M/M_⊙)< 14.45, 14.45 <  log(M/M_⊙)< 14.65 and 14.65 <  log(M/M_⊙)< 14.85). Profiles show a significant dependence on halo mass for all three redshifts. Higher mass clusters have higher amplitude ΔΣ(R) when compared to the low mass clusters for a given redshift, due chiefly to the increased extent of the cluster along the line of sight. Figure 3 makes clear that averaging cluster matter profiles in this way will therefore lead to results that are sensitive to the distribution of the cluster sample in redshift and mass (as well as the spatial distribution and masking of the source galaxies), which complicates the interpretation and characterization of such mean matter profiles.

3.1.2. Evidence for self-similarity in mass and redshift

Motivated by the results from the previous section and the behavior of NFW profiles derived from N-body simulations, we used the same simulated clusters to explore a rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ (see Eq. (12)). This profile would be identical for all samples of clusters if the population were truly self-similar.

To determine the average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ profiles we combine the individual cluster matter profiles ${\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_i$ as

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}{\left(\frac{R}{R_{200\mathrm{c}}}\right)}_j=\frac{1}{N}{\displaystyle \sum _i}\stackrel{\sim }{\mathrm{\Delta }{\mathrm{\Sigma }}_i}{\left(\frac{R}{R_{200\mathrm{c}}}\right)}_j,\end{array}$$(15)

where the summation i is over the N clusters in the sample, and j denotes the radial bin in units of R/R_200c.

In the left panel of Fig. 4 we show average rescaled cluster matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ at five redshifts: 0.01, 0.25, 0.47, 0.78 and 0.96. We average all the clusters for a given redshift and then analyze the redshift trend. The profiles show very small variations and little change in amplitude with redshift as seen in the bottom left panel. This is in contrast to the behavior observed in Fig. 3 where we show ΔΣ(R). The average matter profiles line up for all redshifts from R/R_200c ≈ 0.6 to R/R_200c ≈ 1 with some small, remaining redshift trend at low and high R/R_200c.

Fig. 4. Average cluster matter profiles in the Magneticum simulation rescaled as in Eq. (12) to $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$. In the left plot, we average all the clusters (14.3 <  log(M200c/M⊙)< 15) for a given redshift and analyze the redshift trend. In the right plot, we average clusters over all redshift ranges (0.01 ≤ z ≤ 0.96) within a mass bin and analyze the mass trends. In both panels, the bottom plot shows the deviation of average profiles with respect to the mean of all the profiles. The rescaling dramatically reduces systematic trends in mass and redshift, highlighting the degree of self-similarity in the matter profiles even when baryonic components are included. The thickness of the lines in the upper panels represents the 68% credible region in the average matter profiles.

Similarly, in the right panel of Fig. 4 we combine all the redshift samples and divide them into three mass bins (14.30 <  log(M/M_⊙)< 14.45, 14.45 <  log(M/M_⊙)< 14.65 and 14.65 <  log(M/M_⊙)< 14.85) to study the mass trends in average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$. The profiles show remarkably small variation. This is an indication that even when hydrodynamical effects are included, simulated galaxy clusters over this range of mass and redshift have matter profiles that exhibit strikingly similar shape. As discussed in Sect. 3, the lack of variation in shape is an indication of the self-similarity of cluster matter profiles along dimensions of mass and redshift.

3.1.3. Variations in average matter profiles

In this section we aim to quantify the degree of variation we see in cluster matter profiles. For this, we calculate the fractional scatter in the average matter profiles as a function of redshift and mass. We also compare this to the fractional scatter values obtained when averaging profiles in physical space. The fractional scatter is given by σ_ΔΣ/< ΔΣ> , where σ_ΔΣ is the standard deviation of the sample and ⟨ΔΣ⟩ is the mean of the sample.

We used the Magneticum simulation, which has a larger volume and therefore contains many more halos in comparison to IllustrisTNG. This allows us to measure a higher SNR through averaging profiles, which in turn enables us to better quantify the intrinsic scatter in the average matter profiles. We ignore the uncertainty on the average matter profiles when calculating the fractional scatter with redshift and mass, because the uncertainty is much smaller than the scatter.

Starting with the ΔΣ(R) profiles, we first divide each redshift bin into four mass bins. We calculate the fractional scatter as a function of redshift in a given mass bin and then report the average value of fractional scatter as a function of redshift of four mass bins as a function of radius in the left panel (blue curve) of Fig. 5. The scatter ranges from ≈7% to ≈19% with an average of ≈17%. Similarly, when analyzing the fractional scatter as a function of mass, we compute the scatter at each redshift bin and report the mean value for five redshift bins. The orange curve in the left panel of Fig. 5 shows the mean scatter values as a function of radius R with an average of ≈26%.

Fig. 5. Fractional variation in the average matter profiles versus radius shown for ΔΣ(R) on the left and for average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ on the right. We track the scatter due to redshift variations in blue and mass variations in orange. The average rescaled matter profiles exhibit approximately six and 23 times less variation on average with redshift and mass, respectively, within the one-halo region.

Moving to the average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$, we average all of the clusters for a given redshift and then report the fractional scatter among the five redshift bins. The blue curve in the right panel of Fig. 5 shows the fractional scatter as a function of redshift. The value varies from ≈0.8% to ≈5% with an average value of ≈2.6% and the minimum value is achieved around ≈0.9R/R_200c. The scatter in $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ is reduced by a factor of 6 in comparison to the scatter obtained in ΔΣ(R). Similarly, when analyzing the fractional scatter as a function of mass, we divide each redshift into four mass bins and stack all of the clusters with different redshifts in a given mass bin. The orange curve in Fig. 5 shows the trend of the fractional scatter as a function of mass with scaled radius with an average value of ≈1.1%, which is 23 times smaller relative to the scatter with mass in ΔΣ(R).

3.2. Average matter profiles: Observations

In this section we examine the matter profiles of the 698 SPT tSZE-selected clusters using the WL data from DES in physical and rescaled space. Given the cosmological parameters p (Table 1), the average ΔΣ estimator for the cluster ensemble in a radial bin R_j is a triple sum over clusters k, lensing source galaxy tomographic bins b, and individual lensing source galaxies i as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R_j|\mathit{p})=\frac{{\displaystyle {\sum }_{k,b,i}}\frac{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}{w}_{k,b}{\mathcal{W}}_{k,b,i}^{\mathrm{s}}{e}_{\mathrm{t},k,b,i}}{(1-f_{\mathrm{cl},k,b})}}{{\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\mathcal{W}}_{k,b,i}^{\mathrm{s}}(R_{{\gamma }_{\mathrm{t},i}}+R_{\mathrm{sel}})}.\end{array}$$(16)

Here, ${\mathcal{W}}_{k,b,i}^{\mathrm{s}}$ is the scaled source weight, and w_k, b is the tomographic bin weight

$$\begin{array}{cc}\hfill {\mathcal{W}}_{k,b,i}^{\mathrm{s}}& =w_i^{\mathrm{s}}{\left(\frac{1-f_{\mathrm{cl},k,b}}{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}}\right)}^2\hfill \\ \hfill w_{k,b}& ={\mathrm{\Sigma }}_{\mathrm{crit},k,b}^{-1},\hfill \end{array}$$(17)

where the w_i^s represent individual source weights (defined as the inverse variance in the measured ellipticity). Following Bocquet et al. (2024b), we employed only the tomographic bins 2 to 4 in this analysis. Additionally, we note that we only use the tomographic bins for which the median source redshift is larger than the cluster redshift. Σ_{crit, k, b} is the critical surface density, which depends on cluster redshift and source galaxy redshift distribution, calculated as in Eq. (6). The ellipticity of a source galaxy i from a tomographic bin b and lying in the background of a cluster k is e_{t, k, b, i}, and R_γt, i is the shear response for galaxy i, which is needed to scale the ellipticity to the reduced shear. Additionally, the selection response R_sel accounts for the fact that lensing sources are selected based on their (intrinsically) sheared observations. We use R_sel = −0.0026 for optical centers as measured previously for this sample (Bocquet et al. 2024b). We also scale the ellipticity with a factor of 1/(1 − f_cl, k, b) to correct the profiles for the cluster member contamination, which is measured separately for each tomographic bin (see the discussion in Sect. 4.4).

The corresponding uncertainty in ΔΣ for a radial bin R_j in the average cluster matter profile is calculated as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{\Delta }\mathrm{\Sigma }}^2(R_j|\mathit{p})=\frac{{\displaystyle {\sum }_{k,b,i}}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}{w}_{k,b}{\mathcal{W}}_{k,b,i}^{\mathrm{s}}{\sigma }_{\mathrm{eff},b}}{(1-f_{\mathrm{cl},k,b})}\right)}^2}{{({\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\mathcal{W}}_{k,b,i}^{\mathrm{s}})}^2}.\end{array}$$(18)

Here σ_eff, b² is the effective shape variance for sources in a given tomographic bin, and all other elements are as described in Eq. (16).

Similarly, after accounting for the rescaling described in Eq. (12) and given the cosmological parameters p and the cluster mass M_200c or equivalently the cluster radius R_200c, the inverse variance weighted average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ estimator in a radial bin (R/R_200c)_j is given by

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}{(\frac{R}{R_{200\mathrm{c}}}|{\mathit{M}}_{\mathbf{200}\mathbf{c}},\mathit{p})}_j=\frac{{\displaystyle {\sum }_{k,b,i}}\frac{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}{e}_{\mathrm{t},k,b,i}}{{\rho }_{{\mathrm{crit}}_k}{R}_{{200\mathrm{c}}_k}(1-f_{\mathrm{cl},k,b})}}{{\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}(R_{{\gamma }_{\mathrm{t},i}}+R_{\mathrm{sel}})},\end{array}$$(19)

where ${\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}$ is the re-scaled source weight corresponding to $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$, which is defined as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}=w_i^{\mathrm{s}}{\left(\frac{{\rho }_{{\mathrm{crit}}_k{R}_{{200\mathrm{c}}_k}(1-f_{\mathrm{cl},k,b})}}{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}}\right)}^2,\end{array}$$(20)

and the corresponding uncertainty is given by

$$\begin{array}{c}\hfill {\sigma }_{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}^2{(\frac{R}{R_{200\mathrm{c}}}|{\mathit{M}}_{\mathbf{200}\mathbf{c}},\mathit{p})}_j=\frac{{\displaystyle {\sum }_{k,b,i}}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{crit},k,b}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}{\sigma }_{{\mathrm{eff}}_b}}{{\rho }_{\mathrm{crit},k}{R}_{{200\mathrm{c}}_k(1-f_{\mathrm{cl},k,b})}}\right)}^2}{{({\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}})}^2}.\end{array}$$(21)

The mean estimated scaled radius of a given radial bin j is calculated from the equation

$$\begin{array}{c}\hfill {(\frac{R}{R_{200\mathrm{c}}}|{\mathit{M}}_{\mathbf{200}\mathbf{c}},\mathit{p})}_j=\frac{{\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}{R}_{k,b,i}/R_{{200\mathrm{c}}_k}}{{\displaystyle {\sum }_{k,b,i}}{w}_{k,b}{\stackrel{\sim }{\mathcal{W}}}_{k,b,i}^{\mathrm{s}}},\end{array}$$(22)

where R_k, b, i is the projected separation of ellpiticity i in tomographic bin b from the center of cluster k. A similar expression for the mean radius R_j within a bin pertains, but without the 1/R_200c scaling.

In the left panel of Fig. 6, we show the average cluster matter profile ΔΣ(R), in three redshift bins (for a richness bin, $30<\widehat{\lambda }<90$ containing 426 clusters), and in the right panel, we show the average cluster matter profiles in three richness bins (for a redshift bin, 0.25 <  z <  0.40 containing 123 clusters). While the measurement uncertainties are significant in this sample (color bands represent 68% credible regions), it is still possible to discern variations among the presented average matter profiles.

Fig. 6. SPT cluster average matter profiles ΔΣ(R) for three redshift bins in a given richness bin. In the right panel we show the average matter profiles ΔΣ(R) for different richness bins in a given redshift bin. The color bands encode the 68% credible region for each profile. The profiles show variation with redshift and richness that is consistent with that shown in Fig. 3 for the simulated clusters in redshift and mass.

To study the average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles, we need a robust value of R_200c for each cluster, which then implies we need good mass constraints for each system (Eq. (11)). For this purpose, we adopt the observable-mass relation (ζ-mass and λ-mass – see discussion in Sect. 4.1) posteriors from Bocquet et al. (2024a), where the mass calibration constraints were obtained using the same DES WL data as part of the cosmological cluster abundance analysis of the sample (Bocquet et al. 2024b). We used the full sample while analyzing clusters in redshift or richness bins. The average rescaled cluster matter profiles exhibit less evidence for variation than the ΔΣ(R) profiles in Fig. 7. In other words, the tSZE selected clusters show indications of self-similarity with redshift and richness similar to those presented above for the clusters from hydrodynamical simulations with redshift and mass.

Fig. 7. Average rescaled SPT cluster matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ constructed using the mass calibration posteriors from Bocquet et al. (2024a). In the left panel, we show $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles for three redshift bins, and in the right panel we show $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles for different richness bins. The average rescaled matter profiles show regularity similar to that seen in Fig. 4 for the cluster simulations.

This simplicity in the average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ of the simulated and observed galaxy cluster population offers some advantages. It allows us to combine large ensembles of clusters with a wide halo mass and redshift range, creating higher SNR cluster matter profiles, which can be used to test different models of structure formation. Moreover, the modeling of average cluster matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ for observed cluster samples becomes more straightforward, because the final average profile is insensitive to the redshift and mass distribution of the cluster sample (and to the spatial distribution and masking of the source galaxies). We employ this simplicity in matter profiles in Sect. 4, where we present a new method of galaxy cluster mass calibration that exploits the approximate self-similarity of galaxy clusters.

4. Mass calibration method

Measurements of the WL signal induced by foreground galaxy clusters can be used to robustly estimate the cluster mass. However, because the S/N of the WL signal is low for individual clusters, it is practical to perform mass calibration using the lensing signal averaged over many clusters (e.g., Umetsu et al. 2014, 2016; Okabe et al. 2010, 2013). Because the average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ are particularly simple to model, they offer the possibility to improve upon previous average matter profile based WL mass calibration methods.

The method presented below involves (1) building an ensemble of average rescaled matter profiles – one for each bin in cluster observable, (2) extracting a likelihood of these matter profiles given a model profile and then (3) iterating with a Markov Chain Monte Carlo method to characterize the posteriors of the model parameters that describe the galaxy cluster observable-mass relations discussed in Sect. 4.1. The mean posteriors describe the parameters for which the observed and model average matter profiles are consistent over the full range of cluster observables. The average matter profile model is discussed in Sect. 4.2, and the mass calibration likelihood is presented in Sect. 4.3. A discussion of the systematic effects and their correction then appears in Sect. 4.4.

4.1. Observable-mass relations

In our analysis, each confirmed cluster has four associated observable quantities. These include the tSZE detection significance $\widehat{\zeta }$, the MCMF obtained richness $\widehat{\lambda }$ and redshift z and the WL mass M_WL that is derived using the WL shear and photometric redshift measurements of the background, lensed source galaxies. The WL masses are measured using average matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$, and these masses are used to constrain the so-called cluster observable–mass relations (e.g., Mohr & Evrard 1997; Mohr et al. 1999; Finoguenov et al. 2001; Chiu et al. 2016a, 2018; Bulbul et al. 2019) that describe the redshift dependent statistical binding between the observables (i.e., detection significance, richness and WL mass) and the underlying halo mass, which for this analysis we take to be M_200c.

4.1.1. tSZE detection significance $\widehat{\zeta }$

As described in an early SPT analysis (Vanderlinde et al. 2010), the tSZE detection significance or signal-to-noise ratio $\widehat{\zeta }$ is related to the unbiased significance ζ as

$$\begin{array}{c}\hfill P(\widehat{\zeta }|\zeta )=\mathcal{N}(\sqrt{{\zeta }^2+3},1),\end{array}$$(23)

where 𝒩 denotes a Gaussian distribution. This relationship accounts for the maximization bias in $\widehat{\zeta }$ caused during the cluster matched filter candidate selection (Melin et al. 2006), which has three free parameters (cluster sky location and cluster model filter scale). The normal distribution models the impact of the unit noise in the appropriately rescaled mm-wave maps. The mean unbiased detection significance is modeled as a power-law relation in mass and redshift:

$$\begin{array}{c}\hfill ⟨ln\zeta |M_{200\mathrm{c}},z⟩=ln{\zeta }_0+{\zeta }_{\mathrm{M}}ln\left(\frac{M_{200\mathrm{c}}}{M_{\mathrm{piv}}}\right)+{\zeta }_{\mathrm{z}}ln\left(\frac{E(z)}{E(z_{\mathrm{piv}})}\right),\end{array}$$(24)

where lnζ₀ is the normalization, ζ_M is the mass trend, ζ_z is the redshift trend, E(z) is the dimensionless Hubble parameter, M_piv = 3 × 10¹⁴h⁻¹ M_⊙ is the pivot mass and z_piv = 0.6 is the pivot redshift, which are chosen to reflect the median mass and redshift of our confirmed cluster sample. To account for the variable depth of the SPT survey fields, we rescaled the amplitude ζ₀ on a field-by-field basis

$$\begin{array}{c}\hfill {\zeta }_{0,i}={\gamma }_i{\zeta }_0,\end{array}$$(25)

where γ_i is obtained from simulated maps (Bleem et al. 2015, 2020, 2024). This approach allows us to combine the full SPT cluster sample when empirically modeling the ζ-mass relation. We model the intrinsic scatter in ζ at fixed mass and redshift as log-normal σ_lnζ. This single scatter parameter has been shown to be sufficient to model the SPT tSZE-selected cluster sample ζ-mass relation (Bocquet et al. 2019, 2024a). We return to this question with new validation tools in Sect. 5.2.1.

4.1.2. Cluster richness $\widehat{\lambda }$

The observed cluster richness $\widehat{\lambda }$ is related to the intrinsic richness λ as

$$\begin{array}{c}\hfill P(\widehat{\lambda }|\lambda )=\mathcal{N}(\lambda ,\sqrt{\lambda }),\end{array}$$(26)

which models the Poisson sampling noise in the limit of a normal distribution where the dispersion is $\sigma =\sqrt{\lambda }$. The mean intrinsic richness is modeled as a power law in mass and redshift

$$\begin{array}{c}\hfill ⟨ln\lambda |M_{200\mathrm{c}},z⟩=ln{\lambda }_0+{\lambda }_{\mathrm{M}}ln\left(\frac{M_{200\mathrm{c}}}{M_{\mathrm{piv}}}\right)+{\lambda }_{\mathrm{z}}ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\end{array}$$(27)

where λ₀ is the normalization, λ_M is the mass trend, λ_z is the redshift trend and, as above, M_piv = 3 × 10¹⁴h⁻¹ M_⊙ and z_piv = 0.6. The intrinsic scatter of the intrinsic richness λ at fixed mass and redshift is modeled as a log-normal distribution with the parameter σ_lnλ. This scatter is the same for all redshifts and masses, which has been shown to be adequate for modeling the λ-mass relation of the SPT selected cluster sample (Bocquet et al. 2024a). We return to this question also with new validation tools in Sect. 5.2.1.

4.1.3. Weak lensing mass M_WL

In addition to the ζ-mass and λ-mass relations, we also include a mapping between the so-called WL mass M_WL, which is the mass one would infer by fitting a model profile to an individual cluster matter profile, and the halo mass M_200c. This follows the approach adopted in previous work (Becker & Kravtsov 2011; Dietrich et al. 2019; Grandis et al. 2021; Bocquet et al. 2024b) and is a mechanism for incorporating corrections for systematic biases that may arise from the interpretation of the average matter profiles and for marginalizing over the remaining systematic uncertainties in those bias corrections. For example, the uncertainties associated with hydrodynamical simulations and the subgrid physics they incorporate can be modeled with this M_WL-mass relation, incorporating an effective systematic floor in the accuracy of the final, calibrated masses. We characterize this relation as

$$\begin{array}{cc}\hfill \langle \ln \left(\frac{M_{\mathrm{WL}}}{M_{\mathrm{piv}}}\right)\rangle & =ln{M}_{{\mathrm{WL}}_0}+M_{{\mathrm{WL}}_{\mathrm{M}}}\ln \left(\frac{M_{200\mathrm{c}}}{M_{\mathrm{piv}}}\right)\hfill \\ \hfill & +M_{{\mathrm{WL}}_{\mathrm{z}}}\ln \left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$(28)

where lnM_WL0 is the logarithmic bias at M_200c = M_piv and M_WLM and M_WLz are the mass and redshift trends, respectively, of this bias. For symmetry with the other observable-mass relations, we explicitly include the redshift trend parametrization, whereas in previous analyses (see Bocquet et al. 2024b) the relation has been defined at specific redshifts where the required simulation outputs are available.

The WL mass M_WL estimated from individual clusters exhibits a mass dependent log-normal scatter σ_lnWL about the mean relation, which we model as

$$\begin{array}{cc}\hfill ln{\sigma }_{ln\mathrm{WL}}^2& =ln{\sigma }_{ln{\mathrm{WL}}_0}^2+{\sigma }_{ln{\mathrm{WL}}_{\mathrm{M}}}^{2}ln\left(\frac{M_{200\mathrm{c}}}{M_{\mathrm{piv}}}\right)\hfill \\ \hfill & +{\sigma }_{ln{\mathrm{WL}}_{\mathrm{z}}}^2\ln \left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$(29)

where lnσ_lnWL0² is the logarithm of the variance of M_WL around M_200c at M_piv and z_piv and σ_lnWLM² and σ_lnWLz² are the mass and redshift trends, respectively, of this variance. For an average rescaled matter profile that is produced using N clusters, the effective scatter of the extracted M_WL about the true mass would scale down as $1/\sqrt{N}$, reducing the stochasticity associated with the estimate of the underlying halo mass and reducing the importance of possible correlations between the scatter in M_WL and other observables.

The parameter posteriors on these relations are extracted through a M_WL calibration exercise carried out on hydrodynamical simulations of clusters output over a range of redshifts (Grandis et al. 2021). This calibration exercise employs a model profile or set of model profiles as discussed in the next section and characterizes the biases and scatter associated with that model. In addition, systematic uncertainties on photometric redshifts, the multiplicative shear bias, the cluster member contamination model and the cluster mis-centering model are also incorporated into the posteriors on these parameters, making it straightforward to marginalize over all critical systematic uncertainties in the mass calibration analysis. We return to this in Sect. 4.4.

4.2. Average rescaled matter profile model ${\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}}$

We used the IllustrisTNG and the Magneticum simulations at five redshifts between 0 to 1 (as described in the Sect. 2.3) to create an average WL model for use in mass calibration. Because there is little variation in the average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ with mass and redshift, we could adopt a single average matter profile at all redshifts and masses for the model used in mass calibration. We could then correct for the small biases introduced by this assumption of perfect self-similarity using the M_WL-mass relations (Eq. (28) and (29)). However, an examination of the average matter profiles presented in Fig. 4 provides clear evidence for small departures from self-similarity with redshift, while showing no convincing evidence of departures from self-similarity in mass. Therefore, we adopted a model $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profile that varies with redshift, while assuming perfect self-similarity in mass. This approach sets the mean logarithmic bias lnM_WL0 in the M_WL-mass relation to zero over all masses and redshifts.

To construct an average WL model that represents the typical behavior across both sets of simulations, we select the same number of simulated clusters from both the Magneticum and IllustrisTNG simulations. Since the subgrid physics differ between the two, the average matter profiles may differ. Given that IllustrisTNG has fewer halos, we randomly choose an equal number of halos from Magneticum. Specifically, we select 301, 284, 260, 228, and 176 halos corresponding to the redshifts 0.01, 0.25, 0.47, 0.78, and 0.96, respectively.

For each halo, we extracted three mis-centered cluster matter profiles (using the method described in Grandis et al. 2021) and average all of the clusters at a given redshift from both simulations. We follow the mis-centering distribution model as described in Sect. 4.4.1. Because mis-centering depends on the richness of the cluster, we assign each halo a richness value using the richness-mass relation (Eq. (27)) using the parameters obtained in a previous cluster cosmology analysis (Chiu et al. 2023). These individual ΔΣ(R) profiles are then rescaled into $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles and averaged following Eq. (15).

Given that the two sets of simulations do not have outputs at exactly the same redshifts (e.g., 0.42 versus 0.47 and 0.64 versus 0.78), we quantified the differences in the average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles at these redshifts, verifying that this induces a negligible uncertainty. This is achieved by interpolating the profiles from both simulations as a function of redshift separately and comparing for each simulation the differences in the profiles at both the redshifts, and finding them to be very small (percentage error of ≈ 0.4 %, see Fig. A.1). We therefore adopted the mean value of the redshift (in case the redshifts are different) while combining the profiles from both simulations. The combined profiles are then interpolated as a function of redshift to capture the slight differences we see with redshift in the average $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ profiles.

Once an average rescaled matter profile model has been chosen, it is used to characterize the bias and scatter in the M_WL estimates with respect to the true underlying halo masses, determining posteriors of the parameters in Eqs. (28) and (29). As part of this calibration process, uncertainties on the other crucial systematics (uncorrelated large-scale structure covariance, cluster mis-centering, cluster member contamination of the source galaxy sample, and hydrodynamical uncertainties on the model) are also included. We note that our model uses Eq. (3) to calculate ΔΣ, while we use Eq. (5) to calculate the observed profiles. This introduces a 1% bias in our modeling; however, given that our analysis is dominated by statistical uncertainties, this has no significant impact on our conclusions.

4.3. Mass calibration likelihood

In this section we present the mass calibration likelihood that we employ with the average rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$. The lowest level observational constraint from weak gravitational lensing is a tangential reduced shear profile (Eq. (1)) constructed for each of a series of tomographic bins within which the shear galaxy sample is organized. A complication with using the tangential shear profiles (see, e.g., Bocquet et al. 2024b), is that the profiles from the different bins have amplitudes that depend on Σ_crit, which in turn depends on the redshift distributions of the background galaxies (Eq. (4)). The matter profile ΔΣ(R) (Eq. (3)) is simpler in that the profiles for each tomographic bin are all estimators of the same underlying projected matter density of the cluster (see, e.g., McClintock et al. 2018). The observable we adopt here $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ (Eq. (12)) offers additional simplicity, because this profile is approximately the same for all clusters, independent of their mass and redshift. However, the matter profile ΔΣ(R) and rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ are no longer pure observables. They both have dependences on cosmological parameters that impact the distance-redshift relation, and the rescaled matter profile also has dependences on the masses and redshifts of the constituent clusters. This dependence has to be considered within the likelihood, as outlined in the next subsection.

4.3.1. Likelihood of the rescaled matter profile

The lensing likelihood for an average rescaled matter profile is given by a product of the independent Gaussian probabilities of obtaining the observed matter profile given the model within each radial bin. Because the rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(R/R_{200\mathrm{c}})$ depends on cosmological parameters and the masses and radii of the constituent clusters, the likelihood has to be altered to account for these dependencies. The likelihood transformation for a data vector t_θ which is some function of a data vector y (independent of model parameters) and model parameters θ is given by (Severini 2004)

$$\begin{array}{c}\hfill P(t_{\theta }|\theta )\mathrm{d}\theta ⟶P(t_{\theta }|\theta )|\frac{\partial t_{\theta }}{\partial \mu }|\mathrm{d}\mu .\end{array}$$(30)

Here, μ denotes a function of data y that has the same dimension as t_θ and is independent of the model parameters θ. In the case where μ cannot be expressed with the same dimension as t_θ, the likelihood transformation is then given by

$$\begin{array}{c}\hfill P(t_{\theta }|\theta )\mathrm{d}\theta ⟶P(t_{\theta }|\theta )|\frac{\partial t_{\theta }}{\partial \mu }(\frac{\partial t_{\theta }}{\partial \mu }{)}^{\mathrm{T}}{|}^{1/2}\mathrm{d}\mu .\end{array}$$(31)

Using Eq. (31), we can write the transformed Gaussian likelihood (lensing likelihood) for a rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}({\mathit{M}}_{\mathbf{200}\mathbf{c}},\mathit{p})$ (given by Eq. (19)) with j radial bins as

$$\begin{array}{c}\hfill P(\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}({\mathit{M}}_{\mathbf{200}\mathbf{c}},\mathit{z},\mathit{p})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},\mathit{z})={\displaystyle \prod _j}{P}_{\mathrm{G},j}|\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial {\mathit{e}}_{t,j}}(\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial {\mathit{e}}_{t,j}}{)}^{\mathrm{T}}{|}^{1/2},\end{array}$$(32)

where ${\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}}$ is the model profile and the second factor in the above equation is the transformation calculation for all the radial bins, e_t, j = [e_t¹, e_t², e_t³…e_tⁿ] is the vector containing the n source ellipticities in a given radial bin j and P_G, j is the Gaussian likelihood for that bin

$$\begin{array}{c}\hfill P_{\mathrm{G},j}=(\sqrt{2\pi }{\sigma }_{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }},j}{)}^{-1}exp[-\frac{1}{2}(\frac{{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j-{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod},\mathrm{j}}}{{\sigma }_{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }},j}}{)}^2],\end{array}$$(33)

where ${\sigma }_{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }},j}$ is the rescaled shape noise as described in Eq. (21). The likelihood transformation is a one-dimensional partial derivative matrix with length n and is given by

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial {\mathit{e}}_{t,j}}=[\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^1},\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^2},\dots ,\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^n}].\end{array}$$(34)

The transformation factor can then be expressed as

$$\begin{array}{cc}\hfill |\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial {\mathit{e}}_{t,j}}(\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial {\mathit{e}}_{t,j}}{)}^T{|}^{1/2}=& ({\left(\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^1}\right)}^2+{\left(\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^2}\right)}^2+...\hfill \\ \hfill & +{\left(\frac{\partial {\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_j}{\partial e_{t,j}^n}\right)}^2{)}^{1/2},\hfill \end{array}$$(35)

where we just compute the derivative of Eq. (19) with respect to the measured ellipticities.

4.3.2. Likelihood of single cluster rescaled matter profile

The intrinsic scatter in the observable mass relations, the measurement noise on the observables, and the posteriors of the observable-mass relation parameters and the cosmological parameters all contribute to create the posterior mass distribution $P(M|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})$ for a given cluster. Even in the limit of perfect knowledge of the observable-mass relation and cosmological parameters, this posterior distribution has some characteristic width determined by the mass trends in each observable together with the sources of scatter mentioned above. Moreover, even in the case of a perfect match between the observed and model rescaled matter profiles, the resulting WL mass estimate is a biased and scattered estimator of the true halo mass M_200c as described by the M_WL-halo mass relations (Eq. (28) and (29)).

The cluster mass uncertainty represented by this mass posterior and any biases and scatter in the WL mass estimate have to be accounted for. Therefore, the single cluster lensing likelihood for cluster k with the WL mass posterior $P_k(M_{\mathrm{WL}}|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})$ and rescaled matter profile ${\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_k$ is written as

$$\begin{array}{cc}\hfill P({\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_k|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},z)=\int \mathrm{d}{M}_{\mathrm{WL}}& P_k(M_{\mathrm{WL}}|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})\hfill \\ \hfill & P({\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_k(M_{\mathrm{WL}})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},z),\hfill \end{array}$$(36)

where we are explicit with subscript k to emphasize that this expression represents a weighted likelihood for a single cluster. We note that the M_WL dependence of ${\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_k$ is due to the cluster radius (Eq. (12)), which is mass dependent as in Eq. (11).

We calculate the mass posterior using three observables $\widehat{\zeta }$, $\widehat{\lambda }$ and z (we neglect the cluster photometric redshift uncertainty because it is too small relative to other sources of scatter to be important). According to Bayes’ theorem, the expression for the mass posterior of a cluster that accounts for intrinsic and measurement scatter in the observables is

$$\begin{array}{cc}& P(M_{\mathrm{WL}}|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})=\hfill \\ \hfill & \frac{\int \int \int \mathrm{d}M\mathrm{d}\lambda \mathrm{d}\zeta P(\widehat{\lambda }|\lambda )P(\widehat{\zeta }|\zeta )P(\zeta ,\lambda ,M_{\mathrm{WL}}|M,z,\mathit{p})P(M|z,\mathit{p})}{P(\widehat{\lambda },\widehat{\zeta }|z,\mathit{p})},\hfill \end{array}$$(37)

where the measurement noise is represented by $P(\widehat{\lambda }|\lambda )$ and $P(\widehat{\zeta }|\zeta )$, the intrinsic scatter and any bias in the observable about mass by P(ζ, λ, M_WL|M, z, p), P(M|z, p) is the halo mass function factor which allows us to account for Eddington bias due to the selection, and $P(\widehat{\lambda },\widehat{\zeta }|z,\mathit{p})$ is just the numerator integrated over M_WL.

In the context of multiple observables, the single cluster mass calibration likelihood ℒ_single can be written (assuming $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ is uncorrelated with other observables) as a product of the single cluster lensing likelihood and the likelihood of the observables

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{single}}=P({\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_k|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},z)P(\widehat{\lambda }|\widehat{\zeta },z,p),\end{array}$$(38)

where the second component is the likelihood of the observed richness $\widehat{\lambda }$ given the observed tSZE detection significance $\widehat{\zeta }$ and redshift z. It can be calculated using Bayes’ theorem accounting for intrinsic scatter in the observables

$$\begin{array}{cc}& P(\widehat{\lambda }|\widehat{\zeta },z,\mathit{p})=\hfill \\ \hfill & \frac{\int \int \int \mathrm{d}M\mathrm{d}\lambda \mathrm{d}\zeta P(\widehat{\lambda }|\lambda )P(\widehat{\zeta }|\zeta )P(\zeta ,\lambda |M,z,\mathit{p})P(M|z,\mathit{p})}{P(\widehat{\zeta }|z,\mathit{p})},\hfill \end{array}$$(39)

where $P(\widehat{\zeta }|z,\mathit{p})$ is just the normalization that comes from integrating the numerator over all $\widehat{\lambda }$, including importantly the ${\widehat{\lambda }}_{\min }(z)$ selection threshold, which is crucial for accounting for Malmquist bias. Additionally, in this study, we assume that there is no correlated scatter between ζ and λ, so P(ζ, λ|M, z, p) can be further simplified as P(ζ|M, z, p)P(λ|M, z, p).

4.3.3. Likelihood of multi-cluster average rescaled matter profile

For a given $\widehat{\zeta }-\widehat{\lambda }-z$ bin containing n clusters, the rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ is notionally calculated as in Eq. (19). However, as discussed above for the single cluster rescaled matter profile, the mass posteriors of the clusters must be included. Rather than extracting an average likelihood by marginalizing over the WL mass posterior $P(M_{\mathrm{WL}}|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})$ as in the single cluster case (Eq. (36)), in the multi-cluster case we adopt a Monte Carlo integration approach that allows us to efficiently marginalize over the mass posteriors of all n clusters simultaneously. In effect, we rebuild the average matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ for the cluster ensemble many times and use those profiles to extract likelihoods and then estimate the average likelihood of the rescaled matter profile.

Following the likelihood for a single cluster matter profile in Eq. (38), we write the WL mass calibration likelihood for an ensemble of n clusters with associated observables ${\widehat{\lambda }}_i$, ${\widehat{\zeta }}_i$ and z_i as

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{bin}}=⟨P(\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(\widehat{\mathit{\zeta }}\mathbf{,}\widehat{\mathit{\lambda }}\mathbf{,}\mathit{z},\mathit{p})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},\mathit{z})⟩\prod _{i=1}^{n}P({\widehat{\lambda }}_i|{\widehat{\zeta }}_i,z_i,\mathit{p}),\end{array}$$(40)

where $⟨P(\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(\widehat{\mathit{\zeta }}\mathbf{,}\widehat{\mathit{\lambda }}\mathbf{,}\mathit{z},\mathit{p})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},\mathit{z})⟩$ is the average lensing likelihood of the average rescaled matter profile built from the ensemble. The observable vectors $\widehat{\mathit{\zeta }}$, $\widehat{\mathit{\lambda }}$, and z each contain the measurements for the n clusters in the ensemble. For an n cluster ensemble, it takes the form

$$\begin{array}{cc}\hfill ⟨P& (\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(\widehat{\mathit{\zeta }}\mathbf{,}\widehat{\mathit{\lambda }}\mathbf{,}\mathit{z},\mathit{p})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},\mathit{z})⟩=\hfill \\ \hfill & \mathrm{d}{M}_{{\mathrm{WL}}_1}\dots \mathrm{d}{M}_{{\mathrm{WL}}_n}P(M_{{\mathrm{WL}}_1}|{\widehat{\zeta }}_1,{\widehat{\lambda }}_1,z_1,\mathit{p})\times \dots \hfill \\ \hfill & P(M_{{\mathrm{WL}}_n}|{\widehat{\zeta }}_n,{\widehat{\lambda }}_n,z_n,\mathit{p})\hfill \\ \hfill & P(\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}(M_{{\mathrm{WL}}_1},\dots ,M_{{\mathrm{WL}}_n},\mathit{p})|{\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}}_{\mathrm{mod}},\mathit{z}),\hfill \end{array}$$(41)

where we note that the M_WL is needed to build the rescaled matter profile (using Eq. (19) with M_WL instead of M_200c) and is calculated from observables ($\widehat{\lambda }$, $\widehat{\zeta }$ and z) and the M_WL-mass relation using Eq. (28).

The final likelihood for m$\widehat{\zeta }-\widehat{\lambda }-z$ bins can be written as the product of the likelihood of individual bins

$$\begin{array}{c}\hfill \mathcal{L}=\prod _{\mathrm{bin}=1}^m{\mathcal{L}}_{\mathrm{bin}}.\end{array}$$(42)

4.4. Modeling and correcting for systematic effects

4.4.1. Cluster mis-centering distribution

For each of the clusters in our sample, we have two measurements of the cluster center. The first is the tSZE center as measured by the SPT and the second is the optical center extracted using the MCMF algorithm. MCMF adopts the BCG as the center if it is within 250 kpc of the cluster position determined by SPT; otherwise, the position of the peak of the galaxy density map is used. We only make use of MCMF centers for our analysis. As the observationally determined center is not a perfect tracer of the true halo center, the effect of this mis-centering must be taken into account when modeling the cluster matter profile. We adopt the mis-centering model and the parameters from the recent work by Bocquet et al. (2024b). The mis-centering distribution for the tSZE and optical centers is modeled using the double Rayleigh distribution.

$$\begin{array}{cc}\hfill P_{\mathrm{offset}}(r)& =\rho \text{Rayl}(r,{\sigma }_0)+(1-\rho )\text{Rayl}(r,{\sigma }_1),\hfill \\ \hfill {\sigma }_i& ={\sigma }_{i,0}{\left(\frac{\lambda }{60}\right)}^{1/3}\text{for}i\in \{0,1\}.\hfill \end{array}$$(43)

The double Rayleigh distribution is a good description of the mis-centering of the optical center with respect to the true halo center. The constraints on the mis-centering parameters (ρ, σ₀, σ₁) are obtained by simultaneously fitting for SPT and optical centers. A large fraction of clusters (ρ ≈ 0.89) are well centered, and the two scatter parameters are σ₀ ≈ 0.007 h⁻¹ Mpc and σ₁ ≈ 0.18 h⁻¹Mpc (for additional details, see Bocquet et al. 2024b).

Crucial for our mass calibration analysis is to include the effects of uncertainties on the mis-centering distribution. We do this as part of the M_WL-mass relation calibration. However, given the radial range we adopt for mass calibration (R >  500 h⁻¹ kpc), the miscentering itself and the uncertainties on the mis-centering have little impact on our results.

4.4.2. Cluster member contamination

Cluster galaxies are generally included in the WL source galaxy sample in the case where the cluster redshift lies within the source galaxy redshift distribution associated with a tomographic bin (see Fig. 2). These cluster galaxies are not sheared by their host cluster halo, and therefore their inclusion biases the lensing signal toward being low. To correct for this shear bias, we follow the methodology described in Paulus (2021), developed for the DES Y1 WL dataset, and in Bocquet et al. (2024b), where the method was extended for application to the DES Y3 WL tomographic bin based dataset. The method follows notionally the work by Varga et al. (2019) but more explicitly accounts for varying cluster redshift z_cl. Here we extend this method again by modeling the contamination within each tomographic bin separately.

In our analysis, the fractional contamination by cluster members f_cl, b in tomographic source bin b is extracted by modeling the probability density function of galaxies in redshift along the line of sight toward the cluster as the weighted sum of the cluster member distribution P_cl, b(z) and the average source galaxy field distribution in tomographic bin b as

$$\begin{array}{c}\hfill P_b(z)=f_{\mathrm{cl},b}{P}_{\mathrm{cl},b}(z)+(1-f_{\mathrm{cl},b})P_{\mathrm{field},b}(z),\end{array}$$(44)

where P_field, b(z) corresponds to the field component, which has been divided into three tomographic bins b ∈ 2, 3, 4. For this analysis the individual source galaxy DNZ photo-z’s (De Vicente et al. 2016) are employed.

To determine both the cluster and field components as described, the source galaxies from the DES Y3 shear catalog associated with each cluster are divided into nine logarithmically spaced bins in projected radius, ranging from 0.7h⁻¹ Mpc to 10h⁻¹ Mpc from each cluster center. Because the projected cluster galaxy population falls off rapidly with radial distance from the cluster, the outermost two bins are dominated by the field source distribution. As previously reported in Paulus (2021), the depth inhomogeneities and masking variations in the DES WL source galaxy catalog lead to a field component surface density associated with each tomographic bin being relatively homogeneous on the scale of a cluster but varying significantly over the survey. Thus, we used a local measure of the field surface density around each cluster when modeling the contamination.

The redshift distribution of the cluster component P_cl, b(z) is modeled as a Gaussian distribution in the space of z − z_cl with a redshift dependent characteristic width of σ_z, b(z) and redshift offset parameter z_off, b(z) as

$$\begin{array}{c}\hfill P_{\mathrm{cl},b}(z)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{z},b}(z)}{e}^{-\frac{{(z-z_{\mathrm{cl}}-z_{\mathrm{off},b}(z))}^2}{2{\sigma }_{\mathrm{z},b}^2(z)}}.\end{array}$$(45)

The characteristic width of the cluster members in redshift reflects the typical photo-z uncertainties, and the offset in redshift away from the cluster redshift can result from, for example, the selection applied in dividing source galaxies into tomographic bins.

The spatial distribution of the cluster component is modeled as a projected NFW profile f_NFW(R, c_λ) whose amplitude varies with observed cluster richness as ${\widehat{\lambda }}^{{\mathrm{B}}_{\lambda }}$ and whose redshift variation is extracted directly through measurements within independent redshift bins (Paulus 2021). For convenience, this projected NFW spatial model is normalized to a value of 1.0 at a projected radius of 1 h⁻¹ Mpc, and the NFW concentration 10^c_λ is modeled as a function of richness as

$$\begin{array}{c}\hfill r_s=\frac{{(\widehat{\lambda }/60)}^{1/3}}{{10}^{c_{\lambda }}}{h}^{-1}\mathrm{Mpc}.\end{array}$$(46)

This approach allows adequate freedom to describe the galaxy populations around tSZE selected clusters over our redshift range of interest (Hennig et al. 2017). We note that we measure the contamination around the same cluster centers used for the matter profile analysis, and so mis-centering effects are automatically included. In our analysis, we exclude the core region of the cluster where blending is more common. Moreover, we apply no correction for WL magnification bias (e.g., Chiu et al. 2016b), which is strongest in the excluded cluster core region.

The fractional cluster member contamination extracted is used to apply a radially dependent correction 1/(1 − f_cl, b(R)) to the amplitude of each cluster matter profile ΔΣ(R) derived from each tomographic bin b before these matter profiles are averaged. It can be expressed as

$$\begin{array}{cc}\hfill \frac{1}{1-f_{\mathrm{cl},b}(R)}& =1+e^{A_{\mathrm{eff},b}(z_{\mathrm{cl}})}{\left(\frac{\widehat{\lambda }}{60}\right)}^{B_{\lambda ,b}}{f}_{\mathrm{NFW}}(R,c_{\lambda ,b}),\hfill \end{array}$$(47)

where

$$\begin{array}{c}\hfill A_{\mathrm{eff},b}(z_{\mathrm{cl}})=A_{\infty }+{\displaystyle \sum _i}{A}_i{e}^{-\frac{1}{2}\frac{{(z_{\mathrm{cl}}-z_i)}^2}{{\rho }_{\mathrm{corr}}^2}},\end{array}$$(48)

is a normalization factor that is dependent on the cluster redshift and the extracted amplitudes A_i that are extracted within redshift bins z_i, where the redshift centers of the bins are z_i ∈ {0.2, 0.28, 0.36, 0.44, 0.52, 0.6, 0.68, 0.76, 0.84, 0.92, 1.0}. In this expression we have integrated over the cluster member redshift distribution. For each tomographic bin b, the parameters A_∞, A_i, ρ_corr, B_λ, b, c_λ, b, z_off, b(z) and σ_z, b(z) are obtained by fitting to the entire cluster population. The redshift dependence of the redshift width and offset in Eq. (45) are assumed to be simple linear functions in redshift around a pivot redshift z = 0.5 as σ_z(z)=σ_z0 + σ_zz(z − 0.5) and z_off(z)=z_off0 + z_offz(z − 0.5), respectively (again, within each tomographic bin b).

To solve for the parameters of the cluster member contamination model, we iterate over all clusters in the sample comparing our model (Eq. (44)) to the observed surface density of source galaxies as a function of redshift and projected separation from the cluster. Following Paulus (2021), we apply a regularization term to the likelihood with a correlation length in redshift ρ_{corr − z} to the pairs of neighboring amplitudes e^A_i, b, which then prefers a solution with smooth variation in redshift (as noted explicitly in Eqs. (24) through (27) in Bocquet et al. 2024b). With larger cluster samples these regularization terms would no longer be important. The parameter constraints for each tomographic bin are given in Table C.1. Because the method applied here is similar to that in Bocquet et al. (2024b), we direct the reader to Fig. 9 in Bocquet et al. (2024b) for validation tests of the model.

Figure 8 shows the fractional cluster member contamination at a projected radius of 1.4h⁻¹ Mpc as a function of cluster redshift for the three source galaxy tomographic bins. The mean model is shown in color coded solid lines corresponding to the model parameters as listed in Table C.1 over the redshift range for which each tomographic bin is used. The dashed lines show the 68% credible intervals. For a cluster at a given redshift, the total fractional contamination would be a weighted sum of the fractional contamination within each tomographic bin, where we apply the weight ${\mathrm{\Sigma }}_{\mathrm{crit}}^{-1}$ appropriate for each bin. This weighted contamination is shown as the black black line with associated dashed lines corresponding to the 68% credible interval. The radius 1.4h⁻¹ Mpc is the characteristic radius for the cluster fitting when considering our radial fitting range and the increase of the number of source galaxies with radius. Thus, the typical contamination for a $\widehat{\lambda }=60$ cluster varies from ∼1% at the lowest cluster redshifts to ∼6% at the highest, and the contamination varies as ${\widehat{\lambda }}^{B_{\lambda }}$ where B_λ = 0.78, 0.60 and 0.53 for tomographic bins 2, 3 and 4, respectively.

Fig. 8. Fractional cluster member contamination for a cluster with $\widehat{\lambda }=60$ as a function of redshift for the three source galaxy tomographic bins (see Fig. 2) and their weighted combination. The solid (dashed) lines show the mean (68% credible interval) of the contamination. The colored lines extend over the redshift ranges for which each tomographic bin is employed in constructing the cluster matter profiles. The 1.4h−1 Mpc radius is chosen because that is a characteristic radius from which the WL constraints are coming, given the adopted radial fitting range and the increasing number of source galaxies with radius.

5. Results and discussion

In this section we first present validation tests of our new calibration method applied to mock data in Sect. 5.1. In Sect. 5.2 we show the resulting constraints from the real data, discuss the choices made in this analysis, such as which systematics are included, and present a validation of our adopted form for the observable-mass relations using the real data. We then present the average cluster matter profile out to larger radii (including the two-halo dominated region) and compare it with the simulation in Sect. 5.3. Section 5.4 contains a comparison of our results with those reported in the recent literature.

5.1. Validation of method using mock dataset

We validate our analysis method using a realistic SPT+MCMF mock catalog with DES-like mock WL data. For this purpose, we create mock clusters and lensing data following the approach taken in the recent SPT×DES analysis (Bocquet et al. 2024b).

The first step in creating a mock SPT+MCMF cluster catalog within a fiducial cosmological model is to calculate the expected number of halos as a function of redshift and mass; for this we adopt the Tinker et al. (2008) halo mass function scaled by the surveyed volume as a function of redshift, imposing a mass range 10¹³M_⊙ <  M_200c <  10¹⁶ M_⊙ and a redshift range 0.25 < z < 0.94. We draw a Poisson realization of this sample, and then for each halo we assign cluster observables using the observable-mass relations presented in Sect. 4.1. We then apply survey selection cuts in tSZE detection significance $\widehat{\zeta }$ and optical richness $\widehat{\lambda }$ consistent with those used to produce the real sample. These realistic mocks follow the SPT and MCMF related survey depth geometry and produce a mock sample similar to the real dataset (see Fig. 1) that is fully consistent with the form of the observable-mass relations that we use to analyze the real dataset.

To generate mock WL data for each selected cluster in our mock catalog, we first estimate the total source galaxies associated with the cluster by calculating its area on the sky (corresponding to the radial distance to cluster center R = 10h⁻¹Mpc) and assuming a source galaxy density of 6 arcmin⁻². We assume the same source redshift distribution as the DES Y3 data for each tomographic bin for the mock source galaxies. For each source, we assign a weight by randomly drawing source weights from real DES data. We divide the total source galaxies equally among the three tomographic bins. We then assign a radial distance to each source galaxy by randomly drawing distance ∝ R. We then sample the amplitude of the rescaled matter profile $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ at the distance R_i/R_200c for each galaxy i, given the cluster radius R_200c and using our rescaled matter profile model (see Sect. 4.2). We add tomographic bin-dependent cluster member contamination, consistent with our measurement using the DES data (see Sect. 4.4.2) and then convert the $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ amplitude for each galaxy to g_t using ⟨Σ_crit, b⟩ for the appropriate tomographic bin b. We apply scatter to each g_t measurement by drawing g_t from a normal distribution with σ_eff = 0.3 (effective shape noise for DES data). This process produces realistically noisy and biased tangential shear data for each cluster. Specifically, these shear data include all the known systematic and stochastic effects needed to model cluster shear profiles in DES data.

We create several statistically independent mock catalogs to assess the performance of our likelihood model and the software. For the analysis, we divide our data into 3 × 3 × 3 $\widehat{\zeta }-\widehat{\lambda }-z$ observable bins, and the likelihood calculation is performed following the formalism outlined in Sect. 4.3. Figure 9 shows the posteriors for three mock catalogs with the same set of input observable-mass relation and cosmological parameters but with different uniform random deviate seeds. As expected, all the mock catalogs contain a number of clusters that is similar to the real SPT sample. The corresponding sets of lensing data are also generated with different random number seeds. Mock-1 and mock-2 are generated with a shape noise (σ_eff = 0.3), similar to DES Y3, while mock-3 is created with a shape noise value of 0.05, which is equivalent to scaling up the lensing source galaxy density by a factor of 36, therefore providing a more stringent test of the software.

Fig. 9. Posterior constraints for average cluster mass calibration of three different mock SPT samples in a fiducial flat ΛCDM cosmology. The gray dashed line shows the input parameters used to generate the mock samples, which are recovered within the uncertainties.

To effectively and efficiently sample the high dimensional parameter space, we used the Markov Chain Monte Carlo algorithm MULTINEST (Feroz et al. 2009, 2019) for our likelihood analysis. As is clear in Fig. 9, the posteriors are in good agreement with the input parameters (plotted as dashed lines). These validation tests show no signs of biases.

Our $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ profile analyses of both mock and real datasets typically converge in a factor of five less time on similar computing resources than in the case of the cluster-by-cluster analysis. Interestingly, the time required for a single iteration of the $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$ likelihood is similar to that for the cluster-by-cluster analysis, but the number of iterations required for convergence is typically five times less. This faster convergence seems to be due to the difference in SNR of the average profiles as compared to the individual cluster matter profiles, which influences the stability of the likelihood far away from the best fit parameter values. In future analyses using the average profile method, we plan to present further efficiency improvements. Ongoing testing indicates that these approximate methods reduce the time required for a single likelihood evaluation by more than an order of magnitude.

5.2. SPT×DES analysis

With the validation of the code, we move on to apply the analysis method outlined above to the SPT×DES sample. Following Grandis et al. (2021), we restrict our analysis to the radial range 0.5 <  R/(h⁻¹Mpc)< 3.2/(1 + z_cl). This radial cut allows us to restrict the analysis to the one-halo region while simultaneously avoiding the central region of the cluster, which is most affected by cluster member contamination, mis-centering, and baryonic processes. Throughout the analysis, we use DNF redshifts (De Vicente et al. 2016), which are used to calibrate cluster member contamination. The shape noise per tomographic bin for the DES Y3 data is taken from Amon et al. (2022), which is in good agreement with our bootstrap error estimates

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{eff},b}=\{\begin{array}{cc}0.262\hfill & b=2\hfill \\ 0.259\hfill & b=3\hfill \\ 0.301\hfill & b=4\hfill \end{array}.\end{array}$$

We analyzed the SPT clusters in the redshift range 0.25 ≤ z <  0.94, which contains 698 clusters with DES WL data. As with the mock validation, we divide the data into 3 × 3 × 3 $\widehat{\zeta }-\widehat{\lambda }-z$ bins, leading to 27 independent average rescaled matter profiles $\stackrel{\sim }{\mathrm{\Delta }\mathrm{\Sigma }}$. Because the WL signal in our sample has a higher SNR for lower redshift clusters, the highest redshift bin is chosen to be wider so that it has a sufficient SNR to approximately equal the SNR in the two lower redshift bins. Each redshift bin is further divided into 9 bins (3×3 $\widehat{\zeta }-\widehat{\lambda }$). Here, also, we choose bin boundary values such that each bin has a similar SNR. The observable bins for each redshift range are as follows

$$\begin{array}{cc}\hfill 0.25\le z& <0.33\hfill \\ \hfill & 0\le \widehat{\lambda }<60,60\le \widehat{\lambda }<120,120\le \widehat{\lambda }<250\hfill \\ \hfill & 4.25\le \widehat{\zeta }<5,5\le \widehat{\zeta }<7,7\le \widehat{\zeta }<50\hfill \\ \hfill 0.33\le z& <0.43\hfill \\ \hfill & 0\le \widehat{\lambda }<75,75\le \widehat{\lambda }<120,120\le \widehat{\lambda }<250\hfill \\ \hfill & 4.25\le \widehat{\zeta }<6.5,6.5\le \widehat{\zeta }<8.5,8.5\le \widehat{\zeta }<50\hfill \\ \hfill 0.43\le z& <0.94\hfill \\ \hfill & 0\le \widehat{\lambda }<65,65\le \widehat{\lambda }<110,110\le \widehat{\lambda }<250\hfill \\ \hfill & 4.25\le \widehat{\zeta }<5,5\le \widehat{\zeta }<7,7\le \widehat{\zeta }<50.\hfill \end{array}$$

Parameter priors for our run are listed in Table 1. In summary, we fixed the sum of neutrino masses to the minimum allowed value of 0.06 eV. Nearly all other cosmological parameters were fixed to their mean Planck values (Planck Collaboration VI 2020), except for Ω_m, which has a Gaussian prior (our results are unaffected when using a wide flat prior on Ω_m). The observable-mass relation parameters are assigned a wide flat prior. Moreover, we assume no correlated scatter between ζ and λ.

Rather than precisely following the form of the M_WL − M_200c relation presented in Eqs. (28) and (29), our analysis follows the approach adopted in Bocquet et al. (2024b) where they analyze the same SPT×DES sample. We follow the redshift variation by interpolating between the relations determined through mock observations of specific simulation outputs at redshifts z ∈ {0.252, 0.470, 0.783, 0.963}. We set the amplitude and mass trend of the bias to 0 and 1, respectively, because we have adopted the simulation matter profiles as our model. We adopt the scatter parameter priors from the mock analysis of the simulations at four different redshifts, and we interpolate linearly to obtain the expectation at a given redshift. Additionally, we include the redshift dependent uncertainty on the amplitude of the bias with Gaussian random deviates whose values are scaled using the parameters σ_lnMWL0, 1 and σ_lnMWL0, 2.

The resulting posteriors inferred from applying the new mass calibration software to the SPT×DES data are shown in Fig. 10. The mean posterior, along with the corresponding 68% credible intervals, are listed in Table 2. We find that the mass trend ζ_M of the ζ-mass relation has a value of 1.797 ± 0.195, which is close to the 5/3 scaling one would expect for the tSZE measured within the cluster virial region (e.g., r <  r_200c). However, given that the angular filtering in the tSZE cluster detection removes more flux from larger, more massive clusters, this measured mass trend is likely evidence for a steeper than self-similar relation. The redshift trend ζ_z is consistent with 0. We find that the redshift trend λ_z of the λ-mass relation and the mass trend λ_M, are statistically consistent with 0 and 1, respectively.

Fig. 10. Contour plot showing the posterior constraints for cluster mass calibration using average matter profiles (in red) of the SPT cluster sample in a ΛCDM model. The blue contour shows cluster-by-cluster WL-only mass calibration results from Bocquet et al. (2024a) for the same SPT sample. The dashed blue line shows the prior boundary used in Bocquet et al. (2024a) for two parameters which is smaller than ours.

We note that we marginalize over all the crucial systematic errors in our analysis by adopting the M_WL − M_200c relation. We have repeated the analysis without marginalizing over the systematics, and we do not notice any significant difference. This is expected, because our analysis is shape-noise dominated. Bocquet et al. (2024a) find the same as shown in their Fig. 3.

The blue contours in Fig. 10 show the constraints from Bocquet et al. (2024a) WL analysis of the same SPT×DES sample with MCMF center and DNF redshifts. They perform the analysis on a cluster-by-cluster basis compared to our average profile approach and they use the same radial fitting range as ours. Our analysis is done with a prior on Ω_m unlike in Bocquet et al. (2024a), where they use a flat prior. We use a wider prior on σ_lnζ and ζ_z compared to their work as can be seen in Fig. 10. The slight differences in intrinsic scatter posteriors can be attributed to the different modeling of the observational Poisson noise on the richness in this and the Bocquet et al. (2024a) analyses. In general, the results from our analysis show very good agreement with their single cluster analysis method and every parameter agrees on average at the ≈0.4σ level.

5.2.1. Goodness of fit

To assess the goodness of fit of the data to our model, we compare the average matter profile model to the average cluster matter profile in 3 × 2 × 2 $\widehat{\zeta }-\widehat{\lambda }-z$ bins and perform a χ² fit to all bins. Figure 11 shows 12 average matter profiles (observations represented with black data points), each for a given $\widehat{\zeta }-\widehat{\lambda }-z$ bin along with the model and its corresponding 2σ region (shown in shaded blue region). The profiles are extracted corresponding to the mean parameter obtained from our posterior. We obtain a chi-squared value of χ² = 97.83 from 84 data points, effectively constraining 8 parameters. This corresponds to a reduced chi-squared value of χ_red² = 1.27 and a probability of exceeding the observed χ² of p=0.047.

Fig. 11. Average SPT cluster matter profiles corresponding to the mean recovered parameters in twelve $\widehat{\zeta }-\widehat{\lambda }-z$ bins are shown with black data points with 1σ error bars. The WL model is shown in the red line and the shaded blue region represents 2σ error region on the model.

5.2.2. Observable-mass relation validation

To further validate our mass calibration results, we perform a series of tests on the richness and tSZE observable-mass relations. In these tests we are determining whether the data are consistent with our model description of the observable-mass relations. In Fig. 12 we show the $\widehat{\lambda }$-mass (top) and $\widehat{\lambda }$-redshift (bottom) relations. To analyze the $\widehat{\lambda }$-mass relation, we calculate the richnesses at the pivot redshift by simply dividing $\widehat{\lambda }$ by ((1 + z)/(1 + z_piv))^λ_z. This factor removes the measured redshift trend and allows us to study only the mass trend. On the x-axis, we show mean mass posteriors obtained using dual-observables, given by the following equation:

Fig. 12. Observed optical richness $\widehat{\lambda }$ of SPT clusters as a function of the cluster halo mass (top) and redshift (bottom). The clusters are shown with filled circles, where the error bars also capture the error in the observable-mass relations and the estimated cluster halo mass. The top and bottom plots shows the richness $\widehat{\lambda }$ normalized at the pivot redshift zpiv = 0.6 and the pivot mass Mpiv = 3 × 1014h−1 Mpc, respectively. The intrinsic model (Eq. (27)) is shown in blue. The light and dark-shaded regions in both the panels represent 68% and 95% credible intervals of the mean model respectively.

$$\begin{array}{cc}& P(M_{200\mathrm{c}}|\widehat{\zeta },\widehat{\lambda },z,\mathit{p})=\hfill \\ \hfill & \frac{{\displaystyle \int \int \mathrm{d}}\lambda \mathrm{d}\zeta P(\widehat{\lambda }|\lambda )P(\widehat{\zeta }|\zeta )P(\zeta ,\lambda |M_{200\mathrm{c}},z,\mathit{p})P(M_{200\mathrm{c}}|z,\mathit{p})}{{\displaystyle \int \int \int \mathrm{d}}M\mathrm{d}\lambda \mathrm{d}\zeta P(\widehat{\lambda }|\lambda )P(\widehat{\zeta }|\zeta )P(\zeta ,\lambda |M_{200\mathrm{c}},z,\mathit{p})P(M_{200\mathrm{c}}|z,\mathit{p})}.\hfill \end{array}$$(49)

where of course the parameters of the observable-mass relations have been constrained using the WL mass calibration described above.

Similarly, we study the richness-redshift relation by normalizing the richnesses with (M_200c/M_piv)^λ_M. In both panels the solid black line shows the intrinsic observable-mass relation (Eq. (27)) corresponding to the mean posterior values. The dark and light-shaded blue bands represent 68% and 95% credible intervals on the intrinsic observable-mass relation model. The gray error bar on the data points represents the statistical errors and uncertainties in the observable-mass relations and the estimated cluster halo masses.

In Fig. 13 we show $\widehat{\zeta }$-mass and $\widehat{\zeta }$-redshift relation. On the y-axis, we plot $(\sqrt{{\widehat{\zeta }}^2-3})/\gamma$, where γ is a scale factor that is used to correct for the different depths of fields in the SPT survey. Again, we normalize the y-axis with a factor of (E(z)/E(z_piv))^ζ_z while analyzing the relation with mass and with a factor of (M_200c/M_piv)^ζ_M while analyzing its relation with redshift. The intrinsic observable-mass relation model (Eq. (24)) with mean posterior parameters is shown with a black line and the error region is shown in blue bands as in the previous figure.

Fig. 13. Observed debiased detection significance $\widehat{\zeta }$ of SPT clusters as a function of the cluster halo mass (top) and redshift (bottom). The intrinsic model (Eq. (24)) is shown in blue. The plotting scheme is the same as Fig. 12. The effects of Eddington bias and selection on $\widehat{\zeta }$ can be seen in the above plots.

Compared to the intrinsic observable-mass relations (defined by blue bands) the observables in our analysis suffer from selection effects, which can be seen in Figs. 12 and 13 (we further dicuss this in Appendix B). The Eddington bias is clearly visible in the left side (low mass portion) of the top panels in Figs. 12 and 13 from the points lying above the best-fitting line in $\widehat{\zeta }$-mass and below the relation in $\widehat{\lambda }$-mass. These points are preferentially scattered away from the mean expected observables $\widehat{\zeta }$ and $\widehat{\lambda }$. A particularly strong feature is the two extensions of the observable $\widehat{\zeta }$ to lower mass away from the mean relation. These two features are created by the $\widehat{\zeta }$ selection thresholds in the two SPT survey fields, and these are often referred to as Malmquist bias.

To understand whether the data points are behaving consistently with expectation, we create a mock sample 100 times larger than SPT sample using the mean parameter values from the measured posteriors and apply the same selection in $\widehat{\zeta }$ and $\widehat{\lambda }$ that was applied in creating the real cluster sample. The comparison of the mock and the real data then allows for robust validation.

In Fig. 14 we plot the difference between the observed $\widehat{\lambda }$ and the mean intrinsic observable-mass relation model in log space, within two different mass bins (top panels) and two different different redshift bins (bottom panels). The orange line represents the 100 times larger mock sample and the green histogram is the real SPT sample. Poisson uncertainties are shown as error bars for each real data bin. We see good agreement between the mocks and the real sample in bins of high (χ² = 20.68, p = 0.04) and low mass (χ² = 15.15, p = 0.17) and high (χ² = 26.92, p = 0.002) and low redshift (χ² = 16.86, p = 0.05), indicating that our observable-mass relation model (including intrinsic and observed scatter components) provides a reasonably good description of the real dataset. The p value in the high redshift deviation plot indicates a ∼3σ tension, which is notable. This tension could likely be reduced by introducing a redshift dependent intrinsic scatter component in the λ-mass relation, but we leave that for a future discussion. We note that the deviations in the top panels have the redshift trend scaled out, and the deviations in the bottom panels have the mass trend scaled out, similar to what we show in Fig. 12.

Fig. 14. Distributions of deviations in $\widehat{\lambda }$ about the intrinsic observable-mass relation shown for the SPT sample (green histogram) and a mock dataset (100 times larger) drawn from that same intrinsic relation (orange line). The top (bottom) panels show the deviations in low and high mass (redshift) with the measured redshift (mass) trend removed as in Fig. 12.