© G. W. Pratt et al. 2022

1. Introduction

In a Λ cold dark matter Universe, halo assembly is driven by the hierarchical gravitational collapse of the dominant dark matter component. To first order, this process is self-similar and scale-free because gravity has no characteristic scale. If the baryon content remains constant, power-law relations link the baryonic observable properties – such as X-ray luminosity L_X, or Sunyaev-Zeldovich effect (SZE) signal Y_SZ, or total optical richness Λ – to the cluster mass and redshift. Each of these observables is sensitive to a different underlying intrinsic physical characteristic (e.g., the distribution of gas or the number of red sequence galaxies above a given threshold). However, the detection of a baryonic observable also depends on the intrinsic properties of the signal itself. For example, while the SZE signal is proportional to the gas density, the X-ray emission is proportional to the square of the gas density. This means that X-ray measurements are very sensitive to the physical conditions in the core regions, while SZE measurements are much less so.

The X-ray luminosity, L_X, is an attractive quantity because it can be measured from very few source counts once the redshift is known. For a virialised galaxy cluster where the intra-cluster medium (ICM) is in hydrostatic equilibrium, L_X depends only on the halo mass, M, redshift, z, and the distribution of the gas in the dark matter potential. Power-law L_X − M relations have indeed been observed (e.g., Maughan 2007; Rykoff et al. 2008; Zhang et al. 2008; Connor et al. 2014; Pratt et al. 2009; Vikhlinin et al. 2009; Schellenberger & Reiprich 2017; Lovisari et al. 2020). However, these relations exhibit a large intrinsic scatter (∼40 per cent), linked to the presence of cool cores and merging activity (e.g., Pratt et al. 2009). Exclusion of the core regions, by measuring L_X in an annulus excluding the cluster centre, significantly reduces the intrinsic scatter (Fabian et al. 1994; Maughan 2007; Pratt et al. 2009; Mantz et al. 2018).

Radial gas density profiles have been a key means of obtaining information about the ICM since the advent of X-ray imaging. The high spatial resolution observations afforded by XMM-Newton and Chandra observations have revealed the complexity of both the core regions, which are strongly affected by non-gravitational processes (Croston et al. 2008; Pratt et al. 2010), and the outskirts, where the gas distribution becomes progressively more inhomogeneous (e.g., Eckert et al. 2015). Key open questions are how the ICM evolves over time in the dark matter potential, and how this connects to the formation and evolution of cool cores.

The advent of SZE surveys has resulted in the detection of large numbers of clusters at z > 0.5 (Hasselfield et al. 2013; Bleem et al. 2015; Planck Collaboration XXVII 2016; Hilton et al. 2021), extending the redshift leverage for studies of how the population changes over time. The suggestion that X-ray-selected and SZE-selected samples may not have the same distribution of dynamical states (e.g., Planck Collaboration IX 2011; Rossetti et al. 2016; Andrade-Santos et al. 2017; Lovisari et al. 2017) has prompted examination of the relationship between the baryon signatures and the true underlying cluster population. Indeed, the dynamical state may well be as fundamental a characteristic as the mass or the redshift (Bartalucci et al. 2019). At the same time, McDonald et al. (2017) found that the new SZE-selected samples suggest that cool cores are in place very early in the history of a cluster, and have not changed in size, density, or total mass up to the present. These authors further suggest that much of what was thought to be cool core growth over time is in fact due to the self-similar evolution of the cluster bulk around this static core.

Here we investigate the properties of the gas density profiles and the core-excised X-ray luminosity, L_Xc, of 31 X-ray-selected clusters at z < 0.2 and 93 SZE-selected clusters at z < 1.13. We describe a universal gas density profile for the SZE-selected objects and quantify its variation with redshift and mass. We quantify the radial variation in scaled profiles with respect to the best-fitting evolution and the best-fitting mass dependence. Outside the cores, the median scaled gas density profiles are remarkably similar, showing no dependence on selection. We obtain the radial variation of the intrinsic scatter in scaled profiles and investigate the evolution of this scatter with redshift. We find that L_Xc, measured in the [0.15−1] R₅₀₀ region¹, is an extremely well-behaved mass proxy that does not depend on cluster selection, and shows little evolution beyond self-similar in the broad redshift range of the sample.

Throughout the paper we assume a flat ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. The Sunyaev-Zeldovich flux in units of square arcminutes is denoted Y_SZ; the quantity $D_{\text{A}}^2{Y}_{\text{SZ}}$, in units of square megaparsecs, is the spherically integrated Compton parameter within R₅₀₀, where D_A is the angular diameter distance of the cluster. Unless stated otherwise, logarithmic quantities, including scatter, are given to base e, and uncertainties are quoted at the 68% confidence level.

2. Data and analysis

2.1. Dataset

The dataset consists of 31 X-ray-selected clusters at 0.05 < z < 0.2 and 93 SZE-selected objects at 0.08 < z < 1.13, with six systems in common. The local X-ray-selected dataset is REXCESS (Böhringer et al. 2007; Croston et al. 2008; Pratt et al. 2009). The SZE-selected systems consist of a local sample comprising a subset of 44 objects at z < 0.5 from the Planck Early SZ sample (ESZ; Planck Collaboration VIII 2011), which Bartalucci et al. (2019) show is representative of the full ESZ; these were complemented by a further 49 distant objects observed by XMM-Newton in a series of three large programmes obtained as part of the M2C project. These cover the redshift ranges 0.5 < z < 0.7 (LP1, ID 069366, 072378), 0.7 < z < 0.9 (LP2, ID 078388), and z > 0.9 (LP3, ID 074440). The LP1 sample was selected from objects detected at a signal-to-noise S/R > 4 in the Planck SZ catalogue, and confirmed by Autumn 2011 to be at z > 0.5. The LP2 sample consists of clusters with estimated masses M₅₀₀ > 5 × 10¹⁴ M_⊙ at 0.7 < z < 0.9 in the PSZ2 catalogue. The LP3 sample is derived from the five highest-mass objects at z > 0.9 from the combined Planck and SPT catalogues (Bartalucci et al. 2017, 2018). Full observation details for the sample can be found in Table A.1.

The left-hand panel of Fig. 1 shows the distribution of the 118 clusters in the redshift-mass plane. Here and in the following, we group the SZE-selected systems into three sub-samples in the redshift ranges z < 0.3 (blue), 0.3 ≤ z < 0.6 (light green), and z ≥ 0.6 (dark green), containing 37, 27, and 29 objects, respectively. The right-hand panel of Fig. 1 shows the stacked mass histogram. This plot makes clear that the REXCESS sample (light blue) has a lower median mass (M₅₀₀ = 2.7 × 10¹⁴ M_⊙) than any of the SZE-selected sub-samples (M₅₀₀ = 6.4, 7.9, and 5.0 × 10¹⁴ M_⊙, respectively). The LP2 sub-sample is subject to significant Eddington bias in the Planck signal, leading in the most extreme case to an estimated mass of only 6.5 × 10¹³ M_⊙ for PSZ2 G208.57−44.31. We show below that this has a negligible effect on our results.

Fig. 1. Sample properties. Left: redshift-mass distribution of the clusters used in this paper. The SZE-selected clusters comprise a subset of 44 systems from the Planck Early SZ sample (Planck Collaboration VIII 2011) at z < 0.5 and a further 49 clusters at z > 0.5. REXCESS (Böhringer et al. 2007) is an X-ray-selected sample of 31 objects at z < 0.25. Right: stacked histogram of the mass distribution. The REXCESS sample has a lower median mass than the SZE-selected samples.

2.2. Analysis

As our aim was to compare the SZE-selected clusters to the low-redshift X-ray-selected REXCESS systems, we followed the X-ray data reduction and analysis procedures described in Croston et al. (2008), Pratt et al. (2009, 2010). Event files were reprocessed with the XMM-Newton Science analysis System v15 and associated calibration files. Standard filtering for clean events (PATTERN < 4 and < 13 for MOS1/2 and pn detectors, respectively, and FLAG = 0) and soft proton flares was applied. The instrumental and particle background was obtained from custom stacked, recast data files derived from observations obtained with the filter wheel in the CLOSED position (FWC), renormalised using the count rate in a high energy band free of cluster emission.

Vignetting-corrected, background-subtracted [0.3–2] keV surface brightness profiles were extracted in annular bins centred on the X-ray peak. Temperature profiles were produced using the procedures described in Pratt et al. (2010). These were extracted in logarithmically spaced annular bins centred on the X-ray peak, with a binning of R_out/R)_in = 1.33−1.5 depending on data quality. After subtraction of the FWC spectra, all spectra were grouped to a minimum of 25 counts per bin. The FWC-subtracted spectrum of the region external to the cluster was fitted with a model consisting of two MEKAL components plus an absorbed power law with a fixed slope of Γ = 1.4. The spectra were fitted in the [0.5−10] keV range using χ² statistics, excluding the [1.4−1.6] keV band (due to the Al line in all three detectors), and, in the pn, the [7.45−9.0] keV band (due to the strong Cu line complex). In these fits the MEKAL models were unabsorbed and have solar abundances, and the temperature and normalisations are free parameters; the powerlaw component is absorbed by the Galactic absorption. Since it has a fixed slope, only its normalisation is an additional free parameter in the fit. This best-fitting model was added as an extra component to the annular spectral fits, with its normalisation rescaled to the ratio of the areas of the extraction regions (corrected for bad pixels, chip gaps, etc). In the annular spectral fits, the temperature and metallicity of the cluster component were left free, and the absorption was fixed to the HI value (Kalberla et al. 2005). The metallicity was fixed to a value of Z = 0.3 Z_⊙ when its relative uncertainty exceeded 30%.

2.2.1. Luminosity

The core-excised X-ray luminosity L_Xc was measured in the [0.15−1] R₅₀₀ region for all objects. Here the only change with respect to the analysis in Pratt et al. (2009) was the use of an updated M₅₀₀ − Y_X relation from Arnaud et al. (2010) to estimate the relevant masses M₅₀₀ and scaled apertures R₅₀₀. We show below that this change has a negligible impact on the results. The core-excised X-ray luminosity L_Xc was calculated both in the bolometric ([0.01−100] keV) and soft ([0.5−2] keV) bands for comparison to previous work. As in Pratt et al. (2009), the luminosities were calculated from the [0.3–2] keV band surface brightness profile count rates, using the best-fitting spectral model estimated in the [0.15−1] R₅₀₀ aperture to convert from count rates to luminosity. In cases where the surface brightness profile did not extend to R₅₀₀ (seven systems), we extrapolated using a power law with a slope measured from the data at large radius. Errors on L_Xc take into account the uncertainties in the spectral model, the count rates, and the value of R₅₀₀, and were estimated from Monte Carlo realisations in which the luminosity calculation was derived for 100 surface brightness profiles, the profiles and R₅₀₀ values each being randomised according to the observed uncertainties. Once obtained, the luminosities were further corrected for point spread function (PSF) effects by calculating the ratio of the observed to PSF-corrected count rates in each aperture (see below). Table B.1 lists the resulting masses and luminosities for the present sample, in addition to the name, redshift, RA, and Dec.

2.2.2. Density profiles

The vignetting-corrected background-subtracted [0.3–2] keV surface brightness profiles were used to obtain the deprojected, PSF-corrected density profiles using the regularised, non-parametric technique described in Croston et al. (2006), and applied to the REXCESS sample in Croston et al. (2008). The surface brightness profiles were converted to gas density by calculating an emissivity profile Λ(θ) in XSPEC, taking into account the absorption and instrumental response, and using a parameterised model of the projected temperature and abundance profiles (see e.g., Pratt & Arnaud 2003). The Croston et al. (2006) method uses the parametric PSF model of Ghizzardi (2001) as a function of the energy and angular offsets, the parameters of which can be found in EPIC-MCT-TN-011² and EPIC-MCT-TN-012³. In Bartalucci et al. (2017), the deprojected density profiles from XMM-Newton observations of a number of clusters obtained using this method were compared to Chandra observations, for which the PSF can be neglected. It was shown that the results obtained with the Croston et al. (2006) method reproduced the deprojected Chandra density profiles accurately down to an effective resolution limit of ∼5 arcseconds (Fig. 6 of Bartalucci et al. 2017). The gas density ρ_gas = n_e × (μ_e m_p), where n_e is the electron density measured in X-rays, m_p is the proton mass, and μ_e = 1.148 is the mean molecular weight per free electron:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{gas}}=1.92\times {10}^{-24}\left(\frac{n_{\mathrm{e}}}{{\mathrm{cm}}^{-3}}\right)\mathrm{g}{\mathrm{cm}}^{-3}.\end{array}$$(1)

3. Gas density profiles

3.1. Model

In the self-similar model, a cluster can be completely defined by only two parameters: its mass, M_V, and its redshift, z. A fundamental property of this model is the cluster overdensity, Δ, with respect to the reference density of the Universe ρ_Uni(z), from which the virial mass $M_{\text{V}}=\mathrm{\Delta }{\rho }_{\text{Uni}}(z)(4\pi /3)R_{\text{V}}{}^3$, and radius R_V, can thereafter be defined. Cluster profiles then exhibit a universal form when the radii are scaled to R_V.

The original, and simplest, self-similar model concerns top-hat spherical collapse in the Standard CDM (Ω = 1) cosmology. Here a cluster at redshift z is represented by a spherical perturbation that has just collapsed, with ρ_Univ(z) being the critical density of the Universe, ρ_c = 3 H²(z)/(8 π G), and Δ = 178. Of course, the hierarchical formation of structure in a ΛCDM Universe is a very complex dynamical process: objects continuously accrete matter along large-scale filaments, and there is no strict boundary that would separate a virialised region from the infall zone. The definitions of the mass and the corresponding overdensity are therefore ambiguous, as is the choice of ρ_Univ(z), as one can use either the critical density or the mean density (see Voit 2005, for a review).

Using numerical simulations, Lau et al. (2015) showed that the structure of the inner part of clusters that is typically covered by X–ray observations is more self-similar when scaling by fixed overdensities with respect to the critical density ρ_c(z). The zone in question corresponds to overdensities of Δ ≳ 200. As a scaling radius, we therefore chose an R₅₀₀ corresponding to Δ = 500, the radius within which the mean matter density is ρ₅₀₀ = 500 ρ_c(z). The corresponding total mass within this radius, M₅₀₀, is

$$\begin{array}{c}\hfill M_{500}=(4\pi /3)R_{500}^3{\rho }_{500},\end{array}$$(2)

with

$$\begin{array}{c}\hfill {\rho }_{500}=9.205\times {10}^{-30}E{(z)}^2\mathrm{g}{\mathrm{cm}}^{-3},\end{array}$$(3)

and where E(z) = [Ω_m(1 + z)³ + Ω_Λ]^1/2 is the evolution of the Hubble parameter with redshift in a flat cosmology. The scaled gas density profile expressed as a function of scaled radius x is then

$$\begin{array}{c}\hfill \rho (x)=\frac{{\rho }_{\mathrm{gas}}(R)}{{\rho }_{500}}\mathrm{where}x=\frac{R}{R_{500}}.\end{array}$$(4)

In the self-similar model, ρ(x) follows a universal shape and its normalisation is independent of mass and redshift. In such a case ⟨ρ_gas⟩∝ρ₅₀₀, where the angle brackets denote the average within R₅₀₀, as expected for gas evolution purely driven by gravitation. Figure 2 shows the scaled density profiles of all 118 systems. If the clusters were perfectly self-similar, they would trace the same locus in this plot. This is clearly not the case. The colour-coding by mass and by redshift highlights that at large radius, the scaled profiles of the higher-mass, higher-redshift systems lie systematically above those of lower-mass, lower redshift-objects. These trends suggest a dependence of the scaling on M₅₀₀ and/or redshift.

Fig. 2. Deprojected, PSF-corrected density profiles for 118 galaxy clusters, normalised by the critical density ρcrit and R500. Fully self-similar clusters would trace the same locus in this plot. The profiles are colour-coded by mass M500 in the left-hand panel, and by redshift z in the right-hand panel. There are clear trends with respect to both quantities.

To better understand this dependence, we fitted the observed scaled profiles with a model consisting of a median analytical profile, the normalisation of which is allowed to vary with z and M₅₀₀, with a radially varying intrinsic scatter. The median profile was expressed as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}(x,z,M_{500})=A(z,M_{500})f(x),\end{array}$$(5)

where f(x) is the function describing the profile shape. Here we adopted a generalised Navarro-Frenk-White (GNFW) model (Nagai et al. 2007):

$$\begin{array}{c}\hfill f(x)=\frac{f_0}{{(x/x_{\mathrm{s}})}^{\alpha }{[1+{(x/x_{\mathrm{s}})}^{\gamma }]}^{(3\beta -\alpha )/\gamma }},\end{array}$$(6)

where x_s is the scaling radius, and the parameters (α, γ, 3β) are the central (x ≪ x_s), intermediate (x ∼ x_s), and outer (x ≫ x_s) slopes, respectively. The case α = 0 and γ = 2 corresponds to the standard β model (Cavaliere & Fusco-Femiano 1976), while the case γ = 2 corresponds to the AB model introduced by Pratt & Arnaud (2002). The latter was used to model the median density profile of the REXCESS sample (Piffaretti et al. 2011).

The normalisation is given by the product f₀ A(z, M₅₀₀), where A(z, M₅₀₀) describes the departure from standard self-similarity in terms of a possible mass and/or redshift dependence of the scaled gas density. For this we assumed a power-law dependence on M₅₀₀ and E(z):

$$\begin{array}{c}\hfill A(z,M_{500})=E{(z)}^{{\alpha }_{\mathrm{z}}}{\left[\frac{M_{500}}{5\times {10}^{14}{M}_{\odot }}\right]}^{{\alpha }_{\mathrm{M}}}.\end{array}$$(7)

The standard self-similar model corresponds to α_M = 0, α_z = 0. We expect α_M > 0, as it is well established that the gas mass fraction of local clusters decreases with decreasing mass due to non-gravitational effects (e.g., Pratt et al. 2010). The model above allows us to disentangle mass dependence and possible evolution.

Equations (6) and (7) translate into a gas mass fraction within R₅₀₀, f_g, 500, which varies with mass and redshift as a function of α_M and α_z:

$$\begin{array}{cc}& f_{\mathrm{g},500}=\frac{M_{\mathrm{gas}}(<R_{500})}{M_{500}}=\frac{1}{3}{\displaystyle {\int }_0^1}\rho (x)x^2\mathrm{d}x\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill & =f_{0}A(z,M_{500})I(x_{\mathrm{s}},\alpha ,\gamma ,3\beta ),\hfill \end{array}$$(9)

where M_gas(< R₅₀₀) is the gas mass within R₅₀₀ and we have used Eqs. (2) and (4). The quantity I(x_s, α, γ, 3β) is the three dimensional integral value for f₀ = 1, which depends solely on the shape parameters.

We introduced a radially varying intrinsic scatter term around the model profile, assuming a log-normal distribution at each radius. Taking into account measurement errors, the probability of measuring a given scaled gas density ρ at given scaled radius x for a cluster of mass M₅₀₀ at redshift z is then

$$\begin{array}{cc}\hfill p(\rho |x,z,M_{500})& =\mathcal{N}[log{\rho }_{\mathrm{m}}(x,z,M_{500}),{\sigma }^2(x)]\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\sigma }^2(x)& ={\sigma }_{\mathrm{int}}^2(x)+{\sigma }_{\mathrm{stat}}^2,\hfill \end{array}$$(11)

where 𝒩 is the log-normal distribution. The variance term, σ²(x), is the quadratic sum of the statistical error, σ_stat on the measured log ρ and of the intrinsic scatter on log ρ_m at radius x, σ_int(x).

We expect the intrinsic scatter to increase towards the centre, as observed in Fig. 2, due to the increasing effect of non-gravitational physics on the density profiles. Ghirardini et al. (2019) studied the intrinsic scatter of massive local clusters in the XCOP sample, modelling the radially varying scatter σ_int(x) with a log-parabola function. However, we found that such an analytical form significantly overestimates the scatter in the inner core, which was not covered by their data. To allow for more freedom we used a non-analytical form for the intrinsic scatter, where σ_int(x) is defined at n equally spaced points in log(x) in the typical observed radial range, [x_min–x_max]. The scatter, σ_int(x) at other radii is computed by spline interpolation. We used n = 7 between x_min = 0.01 and x_max = 1.

The likelihood of a set of scaled density profiles measured for a sample of i = 1, N_c clusters of mass M_500, i and redshift z_i is:

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \prod _{i=0}^{N_{\mathrm{c}}}\prod _{j=0}^{N_{\mathrm{R}}[i]}p({\rho }_{i,j}|x_{i,j},z_i,M_{500,i}),}\end{array}$$(12)

where N_R[i] is the number of points of the profile of cluster i, and the quantity ρ_i, j = ρ_gas[i, j]/ρ₅₀₀(z_i) is the scaled density measured at each scaled radius x_i, j = r[i, j]/R₅₀₀(z_i, M_500, i), with ρ_gas[i, j] and r[i, j] being the physical gas density and radius. The statistical error on log ρ_i, j is σ_{stat, i, j}.

We fitted the data (i.e. the set of observed ρ_i, j) using Bayesian maximum likelihood estimation with Markov chain Monte Carlo (MCMC) sampling. Using the emcee package developed by Foreman-Mackey et al. (2013), we maximised the log of the likelihood, which reads (up to an additive constant)

$$\begin{array}{cc}\hfill ln\mathcal{L}& =-0.5{\displaystyle \sum _{i,j}}[ln{\sigma }_{i,j}^2+\frac{{(log{\rho }_{i,j}-log{\rho }_{\mathrm{m},i,j})}^2}{{\sigma }_{i,j}^2}]\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{m},i,j}& ={\rho }_{\mathrm{m}}(x_{i,j},z_i,M_{500,i})\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\sigma }_{i,j}^2& ={\sigma }_{\mathrm{int}}^2(x_{i,j})+{\sigma }_{\mathrm{stat},i,j}^2.\hfill \end{array}$$(15)

The fit marginalises over a total of fourteen parameters: four describing the shape of the median profile (x_s, α, γ, 3β), a global normalisation, f₀, the slopes α_M and α_z that describe the non-standard mass and evolution dependences, and seven additional parameters describing the intrinsic scatter profile. We used flat priors on all parameters.

3.2. Results

To establish a baseline, we fitted the model described above to the 93 SZE-selected systems. The resulting best-fitting model is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{m}}(x,z,M_{500})=A(z,M_{500})f(x),\end{array}$$(16)

with

$$\begin{array}{c}\hfill A(z,M_{500})=E{(z)}^{2.09\pm 0.02}\times {\left[\frac{M_{500}}{5\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right]}^{0.22\pm 0.01}\end{array}$$(17)

and

$$\begin{array}{c}\hfill f(x)=\frac{f_0}{{(x/x_{\mathrm{s}})}^{\alpha }{[1+{(x/x_{\mathrm{s}})}^{\gamma }]}^{(3\beta -\alpha )/\gamma }},\end{array}$$(18)

where

$$\begin{array}{cc}\hfill f_0& =0.13\pm 0.01,\hfill \\ \hfill x_{\mathrm{s}}& =0.28\pm 0.01,\hfill \\ \hfill \alpha & =0.42\pm 0.06,\hfill \\ \hfill \beta & =0.78\pm 0.03,\mathrm{and}\hfill \\ \hfill \gamma & =1.52\pm 0.16.\hfill \end{array}$$

Figure 3 shows the marginalised posterior likelihood for the parameters of the best-fitting density profile model detailed in Sect. 3.1. All parameters are well constrained: in particular, we note that the mass and evolution parameters α_M and α_z do not show any degeneracies, implying that we clearly separate the mass and redshift effects.

Fig. 3. Marginalised posterior likelihood for the parameters of the best-fitting density profile model detailed in Sect. 3.1.

The left-hand panel of Fig. 4 shows the density profiles of the SZE-selected clusters together with the best-fitting model. The intrinsic scatter term is represented by the orange envelope; the numerical values for this term are given in Table 1. The right-hand panel of Fig. 4 shows the best-fitting model for the SZE-selected systems compared to the profiles from the X-ray-selected REXCESS sample. The agreement is excellent beyond the core; in the inner regions, there is a hint that the X-ray-selected systems may show more dispersion. We will return to this point below in Sect. 5.1.5.

Fig. 4. The universal cluster ICM density profile. Left: scaled density profiles of the SZE-selected clusters (grey points), overplotted with the best-fitting GNFW model with free evolution and mass dependence: ρgas/ρ500(R/R500)∝E(z)αz MαM with αz = 2.09 ± 0.02 and αM = 0.22 ± 0.01 (orange line). The model includes a radially varying intrinsic scatter term (orange envelope). Right: comparison of the best-fitting model, defined on the SZE-selected sample, to the best-fitting model for the X-ray-selected REXCESS sample (light blue). Here, the points with error bars are the REXCESS sample.

4. Luminosity scaling relations

We now turn to the scaling relation between L_Xc and the mass M₅₀₀. The bolometric X-ray luminosity of a cluster can be written (Arnaud & Evrard 1999)

$$\begin{array}{c}\hfill L(T)=f_{\mathrm{gas}}^2(T)[M(T)\mathrm{\Lambda }(T)]\widehat{Q}(T),\end{array}$$(19)

where f_gas = M_gas/M is the gas mass fraction, and Λ(T) is the cooling function. The quantity $\widehat{Q}(T)=⟨{\rho }_{\mathrm{gas}}^2⟩/{⟨{\rho }_{\mathrm{gas}}⟩}^2$ is a dimensionless structure factor that depends on the spatial distribution of the gas density (e.g., clumpiness at small scale, shape at large scale, etc.). Further assuming (i) virial equilibrium of the gas in the dark matter potential [M ∝ T^3/2]; (ii) simple Bremsstrahlung emission [Λ(T)∝T^1/2]; (iii) similar internal structure [$\widehat{Q}(T)=\mathrm{const}.$]; (iv) a constant gas mass fraction [$f_{\text{gas}}^2(T)=\text{const}.$], the standard self-similar relation between bolometric X-ray luminosity and mass, L_X ∝ M^4/3, can be obtained. Similar arguments can be used to obtain the soft-X-ray luminosity-mass relation of L_X ∝ M.

4.1. Fitting method

We fitted the data with a power-law relation of the form

$$\begin{array}{c}\hfill E{(z)}^{n}L/L_0=A_{\mathrm{L}}{(X/X_0)}^{B_{\mathrm{L}}},\end{array}$$(20)

where L₀ = 1 × 10⁴⁴ erg s⁻¹ and 5 × 10⁴⁴ erg s⁻¹ for the soft and bolometric bands, respectively, and X₀ = 5 × 10¹⁴ M_⊙ keV and 4 × 10¹⁴ M_⊙ for Y_X and M, respectively. Fitting was undertaken using linear regression in the log-log plane, taking uncertainties in both variables into account, and including the intrinsic scatter. We fitted the data using a Bayesian maximum likelihood estimation approach with Markov chain Monte Carlo (MCMC) sampling. We write the likelihood as defined by Robotham & Obreschkow (2015)

$$\begin{array}{c}\hfill ln\mathcal{L}={\displaystyle \sum _{i=1}^N}[ln\frac{B_{\mathrm{L}}^2+1}{{\sigma }_i^2}-\frac{{(ln(L_i/L_0)-ln{A}_{\mathrm{L}}-B_{\mathrm{L}}ln(X_i/X_0))}^2}{{\sigma }_i^2}],\end{array}$$(21)

with

$$\begin{array}{c}\hfill {\sigma }_i^2={\sigma }^2+B_{\mathrm{L}}^2{\sigma }_{\mathrm{X},i}^2+{\sigma }_{\mathrm{L},i}^2\end{array}$$(22)

and the intrinsic scatter, σ², as a free parameter. MCMC sampling was undertaken using the emcee package developed by Foreman-Mackey et al. (2013), with flat priors in the ranges [ − 2.0, 3.0] and [0.0, 1.0] for B_L and σ, respectively. The results, reported in Table 2, were compared to those obtained with the LINMIX (Kelly 2007) Bayesian regression package: these were indistinguishable and so are not reported here.

4.2. Results

4.2.1. L_Xc − Y_X

We first fitted the relation between the bolometric L_Xc and the mass proxy Y_X for the SZE-selected systems only. With the evolution term left free, the best-fitting relation, shown in Fig. 5, is

$$\begin{array}{c}\hfill E{(z)}^{-1.79\pm 0.08}(L_{\mathrm{Xc}}/L_0)=(1.49\pm 0.03){(Y_{\mathrm{X}}/Y_0)}^{0.97\pm 0.02},\end{array}$$(23)

Fig. 5. Relation between LXc and YX. The blue envelope is the best-fitting relation given in Eq. (23), and the results from Maughan (2007) are also shown for comparison.

with an intrinsic scatter of σ_LXc|YX = 0.09 ± 0.01. This result yields evolution and mass dependences that are in excellent agreement with the self-similar predictions of −9/5 and 1.0, respectively. It is also in excellent agreement with the REXCESS only results of Pratt et al. (2009) and that of Maughan (2007), the latter of which was estimated from Chandra data and is overplotted on the figure.

4.2.2. L_Xc − M: Low redshift with fixed evolution

We initially fixed the evolution factor, n, to the self-similar values of −2 and −7/3 for the soft and bolometric bands, respectively. We first fitted the bolometric L_Xc–M for the REXCESS data only, using a mass pivot of 2 × 10¹⁴ M_⊙, as used by Pratt et al. (2009). The resulting normalisation, A_L = 1.06 ± 0.04, slope B_L = 1.73 ± 0.05, and intrinsic scatter σ_lnLXc = 0.16 ± 0.03, are in excellent agreement with those found by Pratt et al. (2009) using orthogonal Bivariate Correlated Errors and intrinsic Scatter (BCES; Akritas & Bershady 1996) fitting. Similarly good agreement was found for the soft-band L_Xc–M relation, showing that the scaling relation parameters that are obtained for the X-ray-selected are robust to the change in underlying M − Y_X relation used to estimate the mass, and also to differences in fitting method.

We then fitted the SZE-selected clusters at z < 0.3 (37 systems). The results are given in Table 2 and show that the normalisation and slope for this sub-sample are in agreement within 1σ with those found for REXCESS. This indicates that there is no difference in the scaling relation between the local X-ray and SZE-selected samples, once the core region has been excised. There is a slight hint that the intrinsic scatter of the SZE-selected sample about the best-fitting relation is lower than that for REXCESS, although this is only a ∼2σ effect for the bolometric luminosity, and is less significant for the soft-band. The data and best-fitting relations, including the 1σ scatter envelopes, are shown in the left-hand panel of Fig. 6.

Fig. 6. Relation between the core-excised X-ray luminosity LXc and mass estimated from the YX proxy, for the bolometric and soft-band luminosities of 118 systems. Left: data points with the best fitting relation to the X-ray-selected REXCESS sample with the evolution factor fixed to the self-similar value of n = −2 (grey envelope). The dark grey envelope shows the best fitting relation to the 37 SZE-selected systems at z < 0.3 with n = −2. Right: histogram of the log space residuals from the best fitting relation to the SZE-selected objects at z < 0.3.

4.2.3. L_Xc − M: Free evolution

The right-hand panel of Fig. 6 shows the histogram of the residuals of the full SZE-selected sample (93 systems, 0.08 < z < 1.13) with respect to the best-fitting relation to the systems at z < 0.3. The peak is offset by ΔlnL ≲ 0.1, indicating that some evolution beyond self-similar is in fact needed.

We then fitted the full SZE-selected sample with a power-law relation, including a free evolution factor, n. The results are given in Table 2 and the best-fitting relations are shown in the left-hand panel of Fig. 7; the right-hand panel of Fig. 7 shows the residual histograms. The latter are well-centred on zero. The best-fitting evolution terms, n = −2.23 ± 0.09 for the soft-band and n = −2.50 ± 0.10 for the bolometric luminosity, suggest that stronger than self-similar evolution is significant at the ∼2σ level.

Fig. 7. Relation between the core-excised X-ray luminosity LXc and mass estimated from the YX proxy, for the bolometric and soft-band luminosities of 118 systems. Left: best fitting relation (grey envelope) to the full sample (data points) with the evolution factor n left free to vary. The best-fitting values of n are given in Table 2. Right: histogram of the log space residuals from the best fitting relation. Solid lines show the best-fitting Gaussian distributions with σ corresponding to the best-fitting intrinsic scatter in log space (Table 2).

5. Discussion

5.1. Gas density profiles

5.1.1. Comparison with previous work

Pioneering work on parametric models of scaled density profiles (Neumann & Arnaud 1999) obtained from ROSAT allowed the dispersion in radial slopes to be constrained. Croston et al. (2008) studied the scaled density profiles of the REXCESS sample, obtaining for the first time constraints on the radial dependence of the intrinsic scatter. The scaled density profile of the X-COP sample, obtained assuming a self-similar evolution factor E(z)², was presented in Ghirardini et al. (2019). Figure 8 compares the scaled density profiles and the best-fitting model from our SZE-selected sample to their median scaled density profile and 68% dispersion. The agreement is good out to the maximum X-COP radius of ∼2 R₅₀₀, although with subtle differences in the inner regions (R < 0.1 R₅₀₀). Their profile is less peaked, likely due to their not having corrected the density profiles for PSF effects, and has a smaller dispersion than our sample, which may be linked to their more limited mass coverage. The right-hand panel compares the intrinsic scatter measurements, which also agree quite well, although the scatter of the present sample is better constrained. The best-fitting intrinsic scatter obtained from our sample when the evolution factor is forced to the self-similar value of E(z)² is also shown in grey.

Fig. 8. Scaled density profiles. Left: scaled density profiles (points) and best-fitting model (orange envelope) for SZE-selected systems in our sample compared to the median and 68% dispersion from the X-COP sample (Ghirardini et al. 2019, magenta envelope). Right: comparison of best fitting intrinsic scatter model (blue) with that found by Ghirardini et al. (2019, magenta). The best-fitting intrinsic scatter obtained from our sample when the evolution factor is forced to the self-similar value of E(z)2 is also shown in grey.

We can also compare to the results obtained by Mantz et al. (2016), who modelled the evolution and mass dependence of the scaled density profiles of a morphologically relaxed cluster sample of 40 systems at 0.08 < z < 1.06. While their evolution dependence of 2.0 ± 0.2 is in agreement with our results, they find a mass dependence that is consistent with zero (0.03 ± 0.06). The difference with respect to our results may be due simply to cluster selection. They studied dynamically relaxed, hot systems, for which the mass leverage is more limited. Once scaled, they found a scatter in scaled density of ≲20% at R₂₅₀₀. For the typical mass of the present sample, R₂₅₀₀ ≈ 0.45 R₅₀₀, where our intrinsic scatter measurements are in good agreement with theirs.

5.1.2. Radial dependence of the scaling

Once the best fitting model was obtained, we quantified how well the model represents the data by calculating the variation of scaled density at different scaled radii. To better disentangle redshift and mass evolution, we extracted two sub-samples: one covering a large redshift range at nearly constant mass, and another covering a large mass range at nearly constant redshift. These sub-samples are illustrated in the z − M plane by the green and blue regions in the left-hand panel of Fig. 9. We defined ten radial bins in terms of x = r/R₅₀₀ and measured the gas density in each bin. We then scaled the density by the best-fitting model and fitted a function of the form ρ_gas(x)/ρ_crit, 0 ∝ E(z)^α_z M^α_M. For a self-similarly evolving population, the gas density scales with the critical density ρ_crit and there is no mass dependence, and so α_z = 2 and α_M = 0.

Fig. 9. Deviations from the average scaling with radius. Left: redshift-mass distribution of the SZE-selected sample used in this work. The shaded regions indicate cuts for two sub-samples: a large redshift range at nearly fixed mass, and a large mass range at nearly fixed redshift. Middle: degree to which the radial ICM density profile evolves as a function of redshift at nearly fixed mass. The dotted line shows the self-similar expectation (αz = 2). The dashed line shows the best-fitting evolution, which varies from slower than self-similar in the centre (αz ∼ 0.3) to faster than self-similar around R500 (αz ∼ 2.4). Envelopes show the 1 and 2σ uncertainties. Right: degree to which the radial ICM density profile scales with mass at nearly fixed redshift. The dotted line shows the self-similar expectation (αz = 0). The dashed line shows the best-fitting mass dependence of αz = 0.22. The scaled density at nearly fixed redshift does not depend on radius.

The middle- and right-hand panels of Fig. 9 show the degree to which the two sub-samples vary from self-similar scaling and from the best-fitting model scaling, as a function of scaled radius. Uncertainties are large because of the reduced number of data points and the scatter in the data. At nearly constant mass the overall variation with redshift is slightly greater than self-similar. However, the density evolves differently in the core and in the outer regions. The density evolution is consistent with zero in the core: α_z = 0.28 ± 1.10 at x = 0.01, a value that is in good agreement the result found by McDonald et al. (2017) although with large uncertainties. However, the density in the outer regions appears to evolve more strongly than self-similar: α_z = 2.42 ± 0.22 at x = 1.20, a result that is significant at slightly more than 1σ. At nearly fixed redshift, the density varies with mass in agreement with the α_M = 0.22 ± 0.01 scaling established above. At R₅₀₀, this mass dependence is significant at ∼2σ, but does not depend on radius.

5.1.3. Median and scatter of scaled profiles

We now turn to the ensemble properties of the scaled density profiles. The left-hand panel of Fig. 10 shows the median scaled profile, obtained in the log-log plane, and 68% dispersion for REXCESS X-ray-selected sample and the SZE-selected sample split into three redshift bins. It is clear that once scaled, the four sub-samples are remarkably similar beyond 0.2 R₅₀₀. In the core region, the median central density decreases progressively with redshift. The median scaled central densities of the REXCESS and of the SZE-selected sample at z < 0.3 are virtually indistinguishable. We will return to the central regions in Sect. 5.1.5 below.

Fig. 10. Scaled profiles and scatter. Left: median scaled profiles (solid lines) and 68% scatter (envelopes) for REXCESS and the SZE-selected sample split into three redshift bins. Beyond ∼0.2 R500 the scaled profiles are almost indistinguishable. Right: radial profile of the intrinsic scatter for the various sub-samples. The best-fitting intrinsic scatter model obtained from the SZE-selected sample is also shown. Intrinsic scatter is less than 20% between 0.2 ≲ R500 ≲ 1.0. The gold line shows the model intrinsic scatter profile corrected for the covariance between Mgas and R500 (see Sect. 5.1.4), which results in a suppression of the scatter by a factor of about two at R500.

The radial variation of the intrinsic scatter about these median profiles is quantified in the right-hand panel of Fig. 10, together with that of the best-fitting intrinsic scatter model obtained above (from the SZE-selected clusters only). The intrinsic scatter of this model falls below 20% in the radial range 0.2 ≲ R₅₀₀ ≲ 1.0. There is excellent agreement between the best-fitting intrinsic scatter model and the observed profiles, which all follow broadly the same trend with scaled radius: a steep decrease with a minimum at ∼0.5−0.7 R₅₀₀, followed by an increase towards larger radii.

Within R ≲ 0.2 R₅₀₀ the intrinsic scatter in the scaled gas density profiles climbs steeply towards the centre. This increase is intimately linked to the complex physics of the core regions, dynamical activity, and to the presence or absence of cool core systems in the various samples. In this connection, the sample with the largest intrinsic scatter in the central regions is REXCESS, reflecting the presence of cool core systems in this dataset.

Beyond ∼0.2 R₅₀₀ the relative dispersion of all samples dips below 20%, and at R₅₀₀ the dispersion in profiles is ∼15%. In the SZE-selected sub-samples, there is a clear evolution, in the sense that the low-redshift systems exhibit the lowest intrinsic scatter values while the high-redshift systems show higher values. The scatter will be related to intrinsic cluster-to-cluster variations linked to inhomogeneities that will depend on the mass accretion rate and associated dynamical state, together with a component due to uncertainties in the total mass. It is possible that both of these effects conspire to produce higher intrinsic scatter values for higher-redshift systems: one expects an increase in dynamical activity with redshift, while uncertainties in the cluster mass measurement will also increase in the same sense.

5.1.4. Suppression of scatter due to covariance between M_gas and R₅₀₀

The observed scatter in the density profiles may be suppressed by the use of M_gas in the computation of R₅₀₀ when scaling the radial coordinate. All other things being equal, a cluster with a higher than average ρ_gas (for its mass) at some radius, will have a higher than average M_gas and hence Y_X relative to its mass. Since Y_X is then used to estimate R₅₀₀ assuming a mean scaling relation, one would then overestimate R₅₀₀ for this cluster. The radial scaling for this cluster would then be too large, which would move its density profile back towards the mean profile, reducing the apparent scatter. The reduction in scatter will depend on the slope of the density profile (i.e. the reduction is larger where the profile is steeper), so will be radially dependent.

The possible magnitude of this effect was estimated by generating synthetic cluster density profiles with a known amount of scatter, and then scaling them in radius following the method used for the observed clusters in order to test how much the scatter was changed. In more detail, the best-fitting median profile presented in Sect. 3.2 was normalised to match a 6 keV cluster at z = 0.15 (i.e. when integrated to R₅₀₀, the gas mass and Y_X were consistent with the scaling relations used for the observed clusters). For these reference values, the ‘true’ R₅₀₀ is 1210 kpc.

A large number of realisations of this median profile were then generated by resampling the normalisation from a lognormal distribution with a standard deviation of 0.2 (approximating a constant 20% scatter in ρ_gas at all radii). For each realisation, R₅₀₀ was then computed in the same manner as for the observed clusters; the profile was integrated to compute Y_X (assuming a fixed temperature of 6 keV) and hence R₅₀₀, with the process performed iteratively until R₅₀₀ converged. The profile was then scaled in radius by R₅₀₀ and the process was repeated for each realisation of the density profile.

When the distribution of densities in the realisations was measured at the ‘true’ value of R₅₀₀ = 1210 kpc, the input scatter of 20% was recovered. However, when the profiles were each scaled in radius by the value of R₅₀₀ estimated for each realisation from the M − Y_X relation, the scatter at a scaled radius of unity was found to be 10%; that is, the scatter is suppressed by a factor of two at around R₅₀₀ due to the dependence of R₅₀₀ on M_gas in our analysis. The same factor of two suppression was found for different values of the input scatter. We calculated the reduction in scatter at different scaled radii, obtaining a radial profile of the suppression factor. The intrinsic scatter profile corrected for this suppression factor is plotted in gold in Fig. 10. This method makes a number of simplifying assumptions (e.g., there is no scatter in T, the profile is fixed at the median form), but based on this analysis, we estimate that the scatter measured for the observed clusters is likely to be underestimated by a factor of approximately two at R₅₀₀.

5.1.5. Change of central regions over time

McDonald et al. (2017) showed that while the ICM outside the core regions of their SZE-selected sample, covering the redshift range 0.25 < z < 1.2, evolved self-similarly with redshift, the central absolute median density (i.e. expressed in units of cm⁻³) did not. They interpreted this result as being due to an un-evolving core component embedded in a self-similarly evolving bulk.

Our sample is of comparable size to that of McDonald et al. (2017), but the evolution with mass and redshift has been decoupled and quantified (Sect. 3.1). The effective 5″ resolution of XMM-Newton after PSF correction (Bartalucci et al. 2017) allows us to measure the density of the SZE-selected sample across all redshifts down to a scaled radius of R ∼ 0.05 R₅₀₀. The left-hand panel of Fig. 11 shows a histogram of the resulting values⁴ scaled according to the best-fitting model (Eq. (17)) established in Sect. 3.2. The SZE-selected sample was further divided into three redshift bins to better visualise how the sample changes over time. This histogram is characterised by a strong peak centred on a scaled central density of −0.2 (in log space), which is clearly visible, and coincident, in all three SZE-selected sub-samples. While the histogram of the z > 0.6 sub-sample has no detectable skewness, the histogram of the z < 0.3 sub-sample exhibits a distinct tail to higher scaled central densities, which is characterised by a moderately large positive skewness of G₁ = 0.86 that is significant at > 90% (Doane & Seward 2011). This may indicate the gradual appearance of objects with more peaked scaled central densities towards lower redshifts.

Fig. 11. Central density. Left: histogram of central densities for the SZE-selected systems, scaled according to the best-fitting model (Eq. (17)) derived in Sect. 3.2, measured at 0.05 R500. The solid line is a kernel density plot with a smoothing width of 0.15. Middle: histogram of scaled central densities for the z < 0.3 SZE-selected systems compared to the X-ray-selected sample, measured at 0.015 R500. The solid line is a kernel density plot with a smoothing width of 0.15. Right: histogram of central densities for the SZE-selected sample at 40 kpc.

At z < 0.3, the effective 5″ resolution of XMM-Newton after PSF correction allows us to measure the density down to a scaled radius of R ∼ 0.015 R₅₀₀. The middle panel of Fig. 11 shows a histogram of the scaled central density of the SZE-selected clusters at z < 0.3 compared to that of REXCESS. The positive skewness of the SZE-selected systems is confirmed at greater significance (G₁ = 1.06), while the histogram of the REXCESS sample exhibits two peaks in scaled central density: a main peak that is coincident with the peak of the SZE-selected sub-samples, and a secondary peak at a scaled density of 0.7 (in log space). The latter peak is due to cool core systems and may indicate that centrally peaked systems are over-represented in X-ray-selected samples, as has been argued by Rossetti et al. (2017) from their comparison of the image concentration parameter in Planck clusters to those for X-ray-selected systems, and also by Andrade-Santos et al. (2017).

The right-hand panel of Fig. 11 shows the histogram of the central density of the SZE-selected sample at 40 kpc, measured in physical units (cm⁻³). There is a broad maximum at n_{e, 40 kpc} ∼ 0.01 cm⁻³, and the histograms of the three sub-samples coincide. A Kolmogorov-Smirnov test indicates that all three sub-samples come from the same parent distribution. This result suggests, in agreement with McDonald et al. (2017), that the absolute central density remains constant over the redshift range probed by the current sample.

5.2. The L_Xc − M relations

5.2.1. Comparison with other work

Figure 12 shows the best-fitting bolometric L_Xc − M relation for the SZE-selected clusters in the present sample compared to a number of results from the literature (Maughan 2007; Mantz et al. 2010; Bulbul et al. 2019; Lovisari et al. 2020). With the exception of those obtained by Mantz et al. (2010), these studies generally find slopes that are steeper than the self-similar expectation of 4/3, ranging from 1.63 to 1.92. Studies that put constraints on the evolution with redshift (Mantz et al. 2010; Bulbul et al. 2019; Lovisari et al. 2020) generally find good agreement with self-similar expectations (although with large uncertainties). However, any measurement of the dependence of a quantity on the mass will be affected strongly by the sample selection and on how the mass itself has been measured. Data fidelity and sample sizes are now such that systematic effects are starting to become dominant over measurement uncertainties.

Fig. 12. Comparison of the bolometric LXc − M relation to previous work (Maughan 2007; Mantz et al. 2010; Bulbul et al. 2019; Lovisari et al. 2020).

Concerning the sample selection, the results in Sect. 4.2 show that the L_Xc − M relations of X-ray- and SZE-selected systems are in good agreement, suggesting that once the core regions are excluded, effects due to detection methods relying on the ICM do not have any impact. Similarly, Lovisari et al. (2020) showed that their relaxed and disturbed samples had similar L_Xc − M relation slopes and normalisations. This suggests that the L_Xc − M relation may also be relatively robust to selection effects linked to cluster dynamical state, likely due to the small intrinsic scatter.

A more fundamental issue is the mass measurement itself (Pratt et al. 2019). In the present work we have used Y_X as a mass proxy; the works listed above use variously Y_X (Maughan 2007), the gas mass (Mantz et al. 2010), the SZE signal-to-noise (Bulbul et al. 2019), and the hydrostatic mass (Lovisari et al. 2020), as proxies. In this context, the shallower slope of the Mantz et al. (2010) relation compared to the others can be fully explained by their assumption of a constant gas mass fraction in the mass calculation (e.g., Rozo et al. 2014b).

All of the above mass estimates are derived from ICM observables, and all except Lovisari et al. (2020) use scaling laws that have been calibrated on X-ray hydrostatic mass estimates. Independent mass measurements, such as those available from lensing, galaxy velocity dispersions, or caustic measurements (e.g., Maughan et al. 2016), are critical to making progress on this issue. In this connection, weak-lensing mass measurements for individual clusters have been carried out by several projects, such as the Local Cluster Substructure Survey (LoCuSS; Okabe et al. 2010, 2013; Okabe & Smith 2016), the Canadian Cluster Comparison Project (CCCP; Hoekstra et al. 2012, 2015), the Cluster Lensing And Supernova survey with Hubble (CLASH; Merten et al. 2015; Umetsu et al. 2014, 2016), Weighing the Giants (WtG; von der Linden et al. 2014; Kelly et al. 2014; Applegate et al. 2014), and CHEX-MATE (CHEX-MATE Collaboration 2021). The cluster community is undertaking a major ongoing effort to critically compare various mass estimates, obtained from mass proxies and from direct X-ray, lensing, or velocity dispersion analyses (e.g., Rozo et al. 2014a; Sereno & Ettori 2015, 2017; Sereno 2015; Groener et al. 2016). Ultimately, this effort will help to better constrain the parameters of the scaling relations.

5.2.2. Link between density profile and L_Xc

As noted in for example Maughan et al. (2008), the similarity of the ICM density profiles outside the core implies a low scatter in L_Xc, as is indeed observed here. In order to explore how much of the scatter in L_Xc is due to the variation in density profiles, we computed a ‘pseudo luminosity’ $L_{\text{Xc}}^{\prime }={\displaystyle \int {\rho }_{\text{gas}}{}^2}{T}_{\text{X}}{}^{1/2}\text{d}V$ for each cluster. For this calculation, we used the measured density profile for each cluster and assumed isothermality at the measured core-excised temperature for the cluster. The integral was performed over a cylindrical volume from projected radii of 0.15 R₅₀₀ to R₅₀₀. For 16 of 118 clusters, the density profiles did not reach R₅₀₀, so the integrals were truncated at the maximum observed radius. In all cases the profiles reached to ≈90% of R₅₀₀, and the contribution to the luminosity of the outer parts of the profile is very small, so the effect of this truncation is negligible. The scatter in $L_{\text{Xc}}^{\prime }$ then provides an estimate of the scatter in the bolometric L_Xc due only to the scatter in density profiles.

The intrinsic scatter about the best fitting relation to $L_{\text{Xc}}^{\prime }E{(z)}^{-2.5}$ versus M₅₀₀ was then measured (assuming that the fractional statistical error on $L_{\text{Xc}}^{\prime }$ is the same as that in L_Xc for each cluster), giving a value of 11%. This implies that most or all of the intrinsic scatter in the bolometric L_Xc at fixed mass can be explained by the variation in the ICM density profiles.

In principle, the use of M_Yx to determine R₅₀₀ and hence the aperture within which L_Xc and $L_{\text{Xc}}^{\prime }$ are measured, could introduce additional scatter in the L_Xc − M₅₀₀ relation. If M_Yx were scattered high relative to the true mass of a cluster, then R₅₀₀ would be overestimated and the 0.15−1 R₅₀₀ aperture would be shifted to larger radii. This shift would reduce L_Xc since more of the luminosity comes from the inner edge of the aperture than the outer edge. Hence, if M_Yx were scattered high, then L_Xc and $L_{\text{Xc}}^{\prime }$ would be scattered low, adding to the observed scatter in the $L_{\text{Xc}}^{\prime }-M_{500}$ relation. We examined the impact of this effect by adding 10% scatter to M_Yx when computing $L_{\text{Xc}}^{\prime }$. This increased the measured scatter in $L_{\text{Xc}}^{\prime }$ by less than 1%. The dependence of the luminosity aperture on M_Yx leads to a negligible contribution to the scatter in the L_Xc − M_Yx relation.

We therefore conclude that the measured intrinsic scatter in the bolometric L_Xc at fixed mass is dominated almost entirely by the variation in the ICM density profiles. The results do not depend on the aperture for reasonable assumptions on the scatter between M_Yx and the true mass M₅₀₀. Any residual scatter will come from inhomogeneity and/or substructure in the density distribution, or from the effects of structure in the temperature or metallicity distribution.

5.2.3. Impact of selection bias and covariance

We found above that the bolometric L_Xc − M relation has a slope B_L ∼ 1.7, which is steeper than the self-similar value of 1.3, and has a small scatter of σ_lnLXc ∼ 0.15. A similar difference in slope of Δ B_L ∼ +0.4 is seen relative to the self-similar soft-band L_Xc − M relation. One might wonder whether and to what extent selection effects and covariance may impact these results. Firstly, in survey data the observable that is used for cluster detection, 𝒪_det, can be affected by so-called Malmquist bias. This happens because clusters that are scattered to higher values of 𝒪_det, whether by noise or by the intrinsic scatter between the observable and the mass, will be preferentially detected, leading to a positive bias in the average value of 𝒪_det. This is a particular concern near the detection threshold, and will affect the apparent slope of any scaling relation with 𝒪_det. Secondly, the intrinsic scatter of each quantity around the mean relation may well be correlated, leading to a reduction in the measured scatter with respect to the true underlying value.

Since we showed in Sect. 4.2 that the results for the X-ray-selected sample are compatible within 1σ with those for the local SZE-selected sample, we focus here on the SZE selection. As detailed in Sect. 2.1, the SZE-selected sample is composed of four different sub-samples and with different selections: the ESZ, LP1, LP2 and LP3. Although the ESZ sub-sample was generated from high signal-to-noise ratio (S/N > 6) detections, Andrade-Santos et al. (2021) showed that its selection is dominated by the instrumental and astrophysical scatter, and that the intrinsic scatter in the Y_SZ–M relation has a negligible effect on the recovered scaling parameters. The three other SZE sub-samples (LP1, LP2, LP3) are selected at an S/N lower than that of the ESZ sample, and are thus even more affected by the instrumental and astrophysical scatter. We restricted in consequence our selection study to the ESZ. This will provide upper limits to the effects of selection on the results for the L_Xc − M of the SZE-selected systems in the current study.

Appendix C describes in detail the simulations we used, which were similar to those produced for the Andrade-Santos et al. (2021) study. Simulated clusters, modelled with the Arnaud et al. (2010) pressure profile and drawn from a Tinker mass function (Tinker et al. 2008), were injected into the Planck Early SZ maps in the ‘cosmological’ mask region. The Multi-Matched Filter extraction algorithm (Melin et al. 2006) was then applied, to obtain SZ detections at S/N > 6, corresponding to the threshold for the ESZ sample. We then matched the injected and recovered clusters to produce a mock ESZ catalogue, doing this twenty times to generate 3188 detections in total. To investigate selection and covariance effects, we assumed a Gaussian lognormal correlated distribution for Y_SZ, Y_X, and L_Xc at fixed true mass M, with a covariance matrix to describe correlations between parameters.

A first-order estimate of the expected scatter in the L_Xc − M_Yx relation can be obtained by assuming that selection effects are negligible. In this case,

$$\begin{array}{c}\hfill {\sigma }_{ln{L}_{\mathrm{Xc}}|M_{\mathrm{Yx}}}^2={\sigma }_{ln{L}_{\mathrm{Xc}}|M}^2+\left(B_{\mathrm{L}}^2/B_{\mathrm{Y}}^2\right){\sigma }_{ln{Y}_{\mathrm{X}}|M}^2-2t{\sigma }_{ln{L}_{\mathrm{Xc}}|M}{\sigma }_{ln{Y}_{\mathrm{X}}|M}.\end{array}$$(24)

For a measured intrinsic scatter σ_lnLXc|MYx ∼ 0.13, a covariance between L_Xc and Y_X of t = 0.85 (Farahi et al. 2019), and assuming a scatter between Y_X and true mass M of σ_lnYX|M ∼ 0.16 motivated by recent numerical simulations (Planelles et al. 2014; Le Brun et al. 2017; Truong et al. 2018), the resulting dispersion between core-excised luminosity and true mass is σ_lnLXc|M ∼ 0.22. This suggests that the measured scatter is underestimated by a factor of 0.22/0.13 ∼ 1.7 due to the use of Y_X as a mass proxy.

We generated a series of simulations using the above values for true scatter and covariance between quantities, and used them to investigate the effects of selection, intrinsic scatter, and covariance on the results. We found:

• The Malmquist bias induced on the L_Xc − Y_SZ relation due to intrinsic scatter in the $D_{\text{A}}^2{Y}_{\text{SZ}}-M$ relation is completely negligible, with the slope changing by less than 0.5%. The resulting impact on the slope of the L_Xc − M relation is ΔB_L = −0.01. This is important, because it implies that selection effects due to intrinsic scatter in the Y_SZ observable cannot account for the observed steeper slope of the L_Xc − M relations seen here.

• Use of Y_X as a mass proxy introduces additional intrinsic scatter with respect to the underlying mass. The net impact on the recovered slope of the L_Xc − M relation is minimal, with ΔB_L = −0.04.

• Intrinsic scatter in the L_Xc − Y_X relation has the effect of redressing the slope towards its original value, with ΔB_L = +0.02.

• Covariance between Y_X and L_Xc again changes the slope. For t = 0.85 (Farahi et al. 2019), the impact on the slope of the L_Xc − M relation is ΔB_L = +0.04.

We thus conclude that the slope of the L_Xc − M_Yx relation found here is robust to selection effects due to intrinsic scatter in the Y_SZ and Y_X proxies, and to covariance between quantities. The dispersion, however, is very sensitive to the covariance. We estimate that the measured dispersion in the relation between the core-excised luminosity and true mass, L_Xc − M, is underestimated by a factor of ∼1.7 due to the use of Y_X as a mass proxy.

5.3. Link to f_gas − M

The L_Xc − M relations derived above have a steeper dependence than self-similar, which cannot be explained by selection effects, intrinsic scatter, or covariance. There is evidence for a dependence of the gas content with mass, which has the effect of suppressing the luminosity preferentially in lower-mass systems, leading to the observed steepening of the relation (see discussion in Pratt et al. 2009). Interestingly, with the assumption of a standard dependence of temperature on mass (T ∝ M^2/3, e.g., Arnaud et al. 2005; Mantz et al. 2016), use of Eq. (19) with the observed bolometric L_Xc ∝ M^1.7 relation yields f_gas ∝ M^0.21.

The gas density profile model detailed in Sect. 3.1 yields an alternative method to obtain the dependence of the gas mass fraction with mass, as the model includes a dependence of f_gas = M_gas/M on the total mass M, via Eq. (8). The best-fitting gas density model (Eq. (17)) suggests f_gas ∝ M^{0.22 ± 0.01}. This dependence of the gas content on total mass is therefore in good agreement with the expected relation given the observed L_Xc − M dependence, and with previous findings (e.g., Pratt et al. 2009; Lovisari et al. 2015; Ettori 2015) based primarily on X-ray-selected clusters.

6. Conclusions

We have examined the gas density profiles and the relation between the core excised X-ray luminosity L_Xc and the total mass derived from the Y_X mass proxy for 118 X-ray and SZE-selected objects covering a mass range of M₅₀₀ = [0.5−20]×10¹⁴ M_⊙ and extending in redshift up to z ∼ 1.13. We first examined the scaled density profiles:

• The gas density profiles do not scale perfectly self-similarly, exhibiting subtle trends in mass and redshift.

• Motivated by this finding, we fitted an analytic gas density model to the 93 SZE-selected systems. The analytic model is based on a generalised NFW profile, and correctly reproduces the scaled gas density profile and the radial variation of its intrinsic dispersion. Combined with the empirical mass scaling of the profiles, this analytic model defines the gas density profile of SZE-selected clusters as a function of mass and redshift. This model is given in Eqs. (16)–(18).

• The intrinsic dispersion in scaled profiles is greatest in the central regions, declining to a minimum at ∼0.5−0.7 R₅₀₀, and increasing thereafter. The dispersion is similar for X-ray-selected clusters and for local SZE-selected clusters, except in the centre, where the X-ray-selected systems have a higher dispersion. There is a hint for an evolution of the dispersion with redshift, which may be linked to an increase in perturbed clusters at higher redshifts.

• We investigated the effect of covariance between M_gas and R₅₀₀ due to the use of M_YX as a mass proxy, obtaining a radial profile of the scatter suppression factor. Taking into account this suppression factor, we estimated a scatter in scaled density profiles of approximately 40% at R₅₀₀.

• We quantified deviations from the average scaling with radius. These show no variation with mass, but which show a significant variation with redshift, in the sense that the core regions clearly evolve differently as compared to the bulk.

• We examined the scaled central density measured at R = 0.05 R₅₀₀ for the SZE-selected systems, finding that only the z < 0.3 sample is skewed. This skewness is positive, and may indicate the increased presence of centrally peaked systems at later times.

• We measured the scaled central density at R = 0.015 R₅₀₀ for the X-ray and SZE-selected systems at z < 0.3. The scaled central density of the local X-ray-selected sample exhibits two peaks. The main peak, corresponding to non-cool core systems in the X-ray-selected sample, is slightly offset to higher scaled central density from that of the local SZE-selected sample. The secondary peak in the X-ray-selected sample, corresponding to the cool core systems, is not seen in the SZE-selected sample, although the latter does exhibit a clear tail to higher scaled central density as confirmed by the strongly positively skewed distribution.

• The absolute value of the central density in the SZE-selected sample measured at 40 kpc does not appear to evolve with redshift, consistent with the findings of McDonald et al. (2017).

We then examined the relation between the core excised X-ray luminosity L_Xc and the total mass derived from the Y_{X, 500} mass proxy, M_Yx.

• This relation is extremely tight, with a logarithmic intrinsic scatter of σ_lnLxc|MYx ≲ 0.15 depending on sub-sample and band in which the luminosity is measured. Importantly, at low redshift, the best-fitting parameters of this relation do not depend on whether the sample was selected in X-rays or through the SZE, suggesting that L_Xc is a selection-independent quantity.

• The slope of the bolometric relation fitted to the SZE-selected clusters, B ∼ 1.74 ± 0.02, is significantly steeper than self-similar. When left free to vary, the evolution of n = −2.5 ± 0.09 is in agreement with the self-similar value of −7/3 within < 2σ.

• We thoroughly examined the impact of selection bias and covariance on the relation. We found that the slope of the L_Xc − M_Yx relation is robust to selection effects due to intrinsic scatter in the Y_SZ and Y_X proxies, and to covariance between quantities. The dispersion, however, is very sensitive to the covariance. For reasonable values of covariance, we estimate that the measured dispersion in the L_Xc − M relation is underestimated by a factor of at most ∼1.7 due to the use of Y_X as a mass proxy, implying a true scatter of σ_lnLXc|M ∼ 0.22.

• We show explicitly that the scatter in the L_Xc − M relation can be accounted for almost entirely by object-to-object variations in gas density profiles.

With our study we have examined the mass and redshift dependence of the ICM gas density profile, and made quantitative comparisons between X-ray- and SZE-selected samples. Our overall conclusion is consistent with the view that the ICM bulk evolves approximately self-similarly, with the core regions evolving separately due to cooling and feedback from the central active galactic nucleus. Indeed, it suggests potentially subtle differences in the core regions between X-ray- and SZE-selected systems. It also supports a view where the ICM gas mass fraction depends on mass up to high redshift, with a dependence $f_{\mathrm{gas}}\propto M_{500}^{0.22\pm 0.02}$ for the present sample. Further progress can be undoubtedly be made by bringing to bear fully independent mass estimates, such as those that can be obtained from weak lensing and/or galaxy velocity dispersions. Such studies are one of the goals of the CHEX-MATE project (CHEX-MATE Collaboration 2021).

Acknowledgments

G.W.P., M.A., and J.-B.M. acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP72007-2013) ERC grant agreement no. 340519, and from the French space agency, CNES. B.J.M. acknowledges support from the Science and Technology Facilities Council (grant number ST/V000454/1). The results reported in this article are based on data obtained from the XMM-Newton observatory, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.

References

Appendix A: Sample observation details

Table A.1 contains the sample observation details, including: cluster name(s), redshift, coordinates, column density, exposure time, and XMM-Newton OBSID used for the analysis.

Appendix B: Sample data

Table B.1 contains the name, redshift, RA, Dec., mass from the Y_X proxy, and soft-band and bolometric core-excised X-ray luminosity for the sample.

Appendix C: Sample selection bias and covariance tests

We used similar simulations as those produced for the Andrade-Santos et al. (2021) study. Simulated clusters, modelled with the Arnaud et al. (2010) pressure profile and drawn from a Tinker mass function (Tinker et al. 2008), were injected into the Planck Early SZ maps in the ‘cosmological’ mask region. The Y_SZ value for each object was drawn from the Y_SZ − M relation of Arnaud et al. (2010), with a bias between the X-ray calibrated mass and the true mass of (1 − b) = 0.65, a value used to obtain the observed cluster number counts in the Planck cosmology. As the slope of the Y_SZ − Y_X relation is expected to be close to unity, and we are only interested in slope variations, we assume that the normalisation and slope of the Y_X − M and Y_SZ − M relations are the same.

We draw the Y_SZ, Y_X and L_X quantities associated to each simulated cluster following a correlated Gaussian distribution with covariance matrix. If 𝒬_M is the latent value of 𝒬 at a given mass, obtained if there were no scatter, one can write

$$\begin{array}{cc}\hfill Y_{\mathrm{M}}& =A_{\mathrm{Y}}E{(z)}^{2/3}{M}^{B_{\mathrm{Y}}}\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill L_{\mathrm{M}}& =A_{\mathrm{L}}E{(z)}^2{M}^{B_{\mathrm{L}}}\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill & P(D_{\mathrm{A}}^2{Y}_{\mathrm{SZ}},Y_{\mathrm{X}},L_{\mathrm{Xc}}|Y_{\mathrm{M}},Y_{\mathrm{M}},L_{\mathrm{M}})=\mathcal{N}[(Y_{\mathrm{M}},Y_{\mathrm{M}},L_{\mathrm{M}}),V_{\sigma }],\hfill \end{array}$$(C.3)

where Y_M and L_M are the latent values of Y and L at a given mass (the values obtained if there were no scatter), and $D_{\text{A}}^2{Y}_{\text{SZ}}$, Y_X, and L_Xc are the true values. For these simulations we assume B_Y = 1.79 (Arnaud et al. 2010) and B_L = 1.37 (Table 2). 𝒩 is a Gaussian log-normal correlated distribution at fixed mass, where

$$\begin{array}{c}\hfill V_{\sigma }=\left(\begin{array}{ccc}{\sigma }_{ln{Y}_{\mathrm{SZ}}}^2& r{\sigma }_{ln{Y}_{\mathrm{SZ}}}{\sigma }_{ln{Y}_{\mathrm{X}}}& s{\sigma }_{ln{Y}_{\mathrm{SZ}}}{\sigma }_{ln{L}_{\mathrm{Xc}}}\\ r{\sigma }_{ln{Y}_{\mathrm{SZ}}}{\sigma }_{ln{Y}_{\mathrm{X}}}& {\sigma }_{ln{Y}_{\mathrm{X}}}^2& t{\sigma }_{ln{Y}_{\mathrm{X}}}{\sigma }_{ln{L}_{\mathrm{Xc}}}\\ s{\sigma }_{ln{Y}_{\mathrm{SZ}}}{\sigma }_{ln{L}_{\mathrm{Xc}}}& t{\sigma }_{ln{Y}_{\mathrm{X}}}{\sigma }_{ln{L}_{\mathrm{Xc}}}& {\sigma }_{ln{L}_{\mathrm{Xc}}}^2\end{array}\right).\end{array}$$(C.4)

For the covariance matrix V_σ, we assume σ_lnYSZ = 0.12 (Kay et al. 2012; Le Brun et al. 2017); σ_lnYX = 0.16 (Planelles et al. 2014; Le Brun et al. 2017; Truong et al. 2018); r = 0.4 (Farahi et al. 2019; Nagarajan et al. 2019); s = 0.4 (Farahi et al. 2019); and t = 0.85 (Farahi et al. 2019). As detailed in Sect. 5.2.3, with these assumptions, Eqn. 24 gives a first-order estimate of σ_lnLXc = 0.22 in the absence of selection effects.

The Multi-Matched Filter extraction algorithm (Melin et al. 2006) was then applied, to obtain SZ detections at S/N > 6, corresponding to the threshold for the ESZ sample. We then matched the injected and recovered clusters to produce a mock ESZ catalogue, doing this twenty times, resulting in a total of 3188 detections. Measurement errors were estimated from the maps as described in Melin et al. (2006). Fits to the simulated data were performed as described in Sect. 4.2. Our baseline simulation in the following assumes r = s = t = 0 (i.e. zero covariance between quantities).

We first fitted the Y_M − L_M relation, which assumes that there is zero intrinsic scatter in either of the observables with respect to the mass. The scatter in the Y-axis seen in Fig. C.1 is then entirely due to the observational uncertainties in the SZE measurements. The resulting slope of 0.759 ± 0.001 implies a change in the L_Xc − M relation of Δ B_L = −0.01. Adding intrinsic scatter of σ_lnYSZ = 0.12 (Kay et al. 2012; Le Brun et al. 2017), the slope of the Y_SZ − L_M relation is 0.759 ± 0.001, again implying a negligible change of Δ B_L = −0.01 on the slope of the L_Xc − M relation. These results imply that Malmquist bias in the Y_SZ observable are negligible, and cannot account for the significantly steeper than self-similar slope we find for the L_Xc − M relation.

Fig. C.1. Simulated YSZ reextracted from PSZ2 maps, plotted as a function of LXc. The quantity 𝒬M is the latent variable with respect to the mass M, obtained if there were no scatter in the relation. The Figure shows the effect of progressively adding scatter in YSZ, YX, and LXc with M. The blue points include covariance of t = 0.85 between YX and LXc.

We next studied the robustness of the recovery of the L_Xc − M relation slope to intrinsic scatter in the Y_X proxy, due to Y_X being a scattered estimates of Y_SZ. An intrinsic scatter of σ_lnYX = 0.16 (Planelles et al. 2014; Le Brun et al. 2017; Truong et al. 2018) pushes the slope of the Y_X − L_M relation to 0.741 ± 0.002, in turn changing the L_Xc − M relation slope by Δ B_L = −0.04.

We then added an intrinsic scatter of σ_lnLXc = 0.22, finding that this redresses the slope to 0.775 ± 0.005, implying a change of Δ B_L = +0.02 in the slope of the L_Xc − M relation. Finally, we added a covariance of t = 0.85 between Y_X and L_Xc (Farahi et al. 2019). This pushes the slope to a slightly steeper value of 0.785 ± 0.003, changing the L_Xc − M relation slope by Δ B_L = +0.04. The covariance significantly reduces the dispersion in the Y_X − L_Xc relation, by a factor of two. This can clearly be seen in the difference in dispersion between the red and blue points in Fig. C.1.

In conclusion, we find that selection effects and intrinsic scatter have a negligible effect on the slope of the L_Xc − M relation, and cannot account for the significantly steeper than self-similar value that we find in this work. The dispersion of the L_Xc − M relation that we derive is significantly underestimated, most likely by a factor of ∼1.7, due to the covariance between L_Xc and the mass proxy, Y_X.

All Tables

All Figures

Fig. 1. Sample properties. Left: redshift-mass distribution of the clusters used in this paper. The SZE-selected clusters comprise a subset of 44 systems from the Planck Early SZ sample (Planck Collaboration VIII 2011) at z < 0.5 and a further 49 clusters at z > 0.5. REXCESS (Böhringer et al. 2007) is an X-ray-selected sample of 31 objects at z < 0.25. Right: stacked histogram of the mass distribution. The REXCESS sample has a lower median mass than the SZE-selected samples.

In the text

	Fig. 2. Deprojected, PSF-corrected density profiles for 118 galaxy clusters, normalised by the critical density ρcrit and R500. Fully self-similar clusters would trace the same locus in this plot. The profiles are colour-coded by mass M500 in the left-hand panel, and by redshift z in the right-hand panel. There are clear trends with respect to both quantities.
In the text

	Fig. 3. Marginalised posterior likelihood for the parameters of the best-fitting density profile model detailed in Sect. 3.1.
In the text

Fig. 4. The universal cluster ICM density profile. Left: scaled density profiles of the SZE-selected clusters (grey points), overplotted with the best-fitting GNFW model with free evolution and mass dependence: ρgas/ρ500(R/R500)∝E(z)αz MαM with αz = 2.09 ± 0.02 and αM = 0.22 ± 0.01 (orange line). The model includes a radially varying intrinsic scatter term (orange envelope). Right: comparison of the best-fitting model, defined on the SZE-selected sample, to the best-fitting model for the X-ray-selected REXCESS sample (light blue). Here, the points with error bars are the REXCESS sample.

In the text

	Fig. 5. Relation between LXc and YX. The blue envelope is the best-fitting relation given in Eq. (23), and the results from Maughan (2007) are also shown for comparison.
In the text

Fig. 6. Relation between the core-excised X-ray luminosity LXc and mass estimated from the YX proxy, for the bolometric and soft-band luminosities of 118 systems. Left: data points with the best fitting relation to the X-ray-selected REXCESS sample with the evolution factor fixed to the self-similar value of n = −2 (grey envelope). The dark grey envelope shows the best fitting relation to the 37 SZE-selected systems at z < 0.3 with n = −2. Right: histogram of the log space residuals from the best fitting relation to the SZE-selected objects at z < 0.3.

In the text

Fig. 7. Relation between the core-excised X-ray luminosity LXc and mass estimated from the YX proxy, for the bolometric and soft-band luminosities of 118 systems. Left: best fitting relation (grey envelope) to the full sample (data points) with the evolution factor n left free to vary. The best-fitting values of n are given in Table 2. Right: histogram of the log space residuals from the best fitting relation. Solid lines show the best-fitting Gaussian distributions with σ corresponding to the best-fitting intrinsic scatter in log space (Table 2).

In the text

Fig. 8. Scaled density profiles. Left: scaled density profiles (points) and best-fitting model (orange envelope) for SZE-selected systems in our sample compared to the median and 68% dispersion from the X-COP sample (Ghirardini et al. 2019, magenta envelope). Right: comparison of best fitting intrinsic scatter model (blue) with that found by Ghirardini et al. (2019, magenta). The best-fitting intrinsic scatter obtained from our sample when the evolution factor is forced to the self-similar value of E(z)2 is also shown in grey.

In the text

Fig. 9. Deviations from the average scaling with radius. Left: redshift-mass distribution of the SZE-selected sample used in this work. The shaded regions indicate cuts for two sub-samples: a large redshift range at nearly fixed mass, and a large mass range at nearly fixed redshift. Middle: degree to which the radial ICM density profile evolves as a function of redshift at nearly fixed mass. The dotted line shows the self-similar expectation (αz = 2). The dashed line shows the best-fitting evolution, which varies from slower than self-similar in the centre (αz ∼ 0.3) to faster than self-similar around R500 (αz ∼ 2.4). Envelopes show the 1 and 2σ uncertainties. Right: degree to which the radial ICM density profile scales with mass at nearly fixed redshift. The dotted line shows the self-similar expectation (αz = 0). The dashed line shows the best-fitting mass dependence of αz = 0.22. The scaled density at nearly fixed redshift does not depend on radius.

In the text

Fig. 10. Scaled profiles and scatter. Left: median scaled profiles (solid lines) and 68% scatter (envelopes) for REXCESS and the SZE-selected sample split into three redshift bins. Beyond ∼0.2 R500 the scaled profiles are almost indistinguishable. Right: radial profile of the intrinsic scatter for the various sub-samples. The best-fitting intrinsic scatter model obtained from the SZE-selected sample is also shown. Intrinsic scatter is less than 20% between 0.2 ≲ R500 ≲ 1.0. The gold line shows the model intrinsic scatter profile corrected for the covariance between Mgas and R500 (see Sect. 5.1.4), which results in a suppression of the scatter by a factor of about two at R500.

In the text

Fig. 11. Central density. Left: histogram of central densities for the SZE-selected systems, scaled according to the best-fitting model (Eq. (17)) derived in Sect. 3.2, measured at 0.05 R500. The solid line is a kernel density plot with a smoothing width of 0.15. Middle: histogram of scaled central densities for the z < 0.3 SZE-selected systems compared to the X-ray-selected sample, measured at 0.015 R500. The solid line is a kernel density plot with a smoothing width of 0.15. Right: histogram of central densities for the SZE-selected sample at 40 kpc.

In the text

	Fig. 12. Comparison of the bolometric LXc − M relation to previous work (Maughan 2007; Mantz et al. 2010; Bulbul et al. 2019; Lovisari et al. 2020).
In the text

	Fig. C.1. Simulated YSZ reextracted from PSZ2 maps, plotted as a function of LXc. The quantity 𝒬M is the latent variable with respect to the mass M, obtained if there were no scatter in the relation. The Figure shows the effect of progressively adding scatter in YSZ, YX, and LXc with M. The blue points include covariance of t = 0.85 between YX and LXc.
In the text