© ESO 2020

1. Introduction

Galaxy clusters are cosmological objects that form by hierarchical aggregation of matter under the action of the gravity force. As a consequence of this force the clusters have physical properties that scale to the mass and redshift of the dark matter halo (e.g. Kaiser 1986; Böhringer et al. 2012; Kravtsov & Borgani 2012). This is true not only for their total gravitating mass, but also for observational quantities that depend mostly on the depth of the potential well, like the temperature of the intracluster medium, which indicates how much most of the baryons collapsed in a dark matter halo are heated by the accretion shocks.

Where the physical processes of aggregation and collapse are dominated by the gravity, a set of relations emerges from this self-similar scenario between the global quantities (i.e. integrated over the volume) that describe the observed properties of the galaxy clusters, such as gas temperature, luminosity, mass, and total mass. However, a predictable behaviour is also expected to hold in how these quantities vary with the radial distance from the bottom of the potential well, in particular in regions away from the core. In the cluster’s core, the relative distribution and energetics of the baryons might be affected by feedback from star formation and active galactic nuclei, as well as radiative cooling, making their distribution less predictable, but not regulated by gravity alone. These radial profiles are then defined as universal because they should reproduce the observed profiles of any galaxy cluster, once rescaled by some quantities that are proportional to the mass and the redshift of the halo.

Universal radial profiles of the electron density (Croston et al. 2008), gas temperature (Vikhlinin et al. 2006; Baldi et al. 2012), electron pressure (Nagai et al. 2007; Arnaud et al. 2010), and gas entropy (Pratt et al. 2010) have been obtained recently as average scaled profiles, i.e. by rescaling the measured quantities through the expected global mean value (which affects the normalisation of the profile) and a characteristic physical radius (e.g. R₅₀₀, which defines the radial scale). Ghirardini et al. (2019) present the most recent work in which universal radial profiles of the intracluster medium (ICM) properties are recovered out to R₂₀₀ for the X-COP sample of 12 nearby massive galaxy clusters.

In the present work we investigate how, by assuming a universal radial profile for the gas pressure and a relation for the internal distribution of the cluster dark matter, we can recover the universal integrated properties of a given halo. This semi-analytic approach is then compared with very recent observational results on the radial profiles and the scaling relations between hydrostatic masses, gas masses, gas luminosities, and temperatures obtained for a Planck-selected sample in Lovisari et al. (2020).

The paper is organised as follows. In Sect. 2 we present our assumptions and the predictions for the universal profiles of the thermodynamic quantities. The integrated quantities are described in Sect. 3, where we also discuss an application of the relations presented in Ettori (2015) that account for the deviations from the self-similar scenario in a physical consistent framework. We summarise our main findings and draw our conclusions in Sect. 4.

In the following analysis we refer to radii, R_Δ, and masses, M_Δ, which are the corresponding values estimated at the given overdensity Δ as $M_{\mathrm{\Delta }}=4/3\pi \mathrm{\Delta }{\rho }_{c,z}{R}_{\mathrm{\Delta }}^3$, where ${\rho }_{c,z}=3H_z^2/(8\pi G)$ is the critical density of the universe at the observed redshift z of the cluster, G is the universal gravitational constant, and H_z = H₀ [Ω_Λ+Ω_m(1+z)³]^0.5 = H₀ E_z is the value of the Hubble constant at the same redshift. For the ΛCDM model we adopt the cosmological parameters H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 1 − Ω_Λ = 0.3.

Historically, the characteristic thermodynamic quantities within a given overdensity Δ refer to a singular isothermal sphere with mass M_Δ and radius R_Δ in hydrostatic equilibrium. In this case mass and gas temperature are simply related by the equation (e.g. Voit et al. 2005)

$$\begin{array}{c}\hfill \frac{G{M}_{\mathrm{\Delta }}\mu m_a}{2R_{\mathrm{\Delta }}}=k_{\mathrm{B}}{T}_{\mathrm{\Delta }},\end{array}$$(1)

where μ is the mean molecular weight of the gas, m_a is the atomic mass unit of 1.66 × 10⁻²⁴ g, and k_B is the Boltzmann constant. Associating a mean electron density ${\overline{n}}_{e,\mathrm{\Delta }}=\mathrm{\Delta }{f}_{\mathrm{gas}}{\rho }_{c,z}/({\mu }_e{m}_a)=\mu /{\mu }_e{\overline{n}}_{\mathrm{gas}}$, where f_gas = M_gas(< R_Δ)/M_Δ is the mass gas fraction, ${\overline{n}}_{\mathrm{gas}}$ is the mean gas density, and μ/μ_e = 0.52, with μ_e = 1.17 and μ = 0.61 being respectively the mean molecular weight of electrons and of the gas for a plasma with 0.3 solar abundance tabulated in Anders & Grevesse (1989)¹, we can also write the mean values for pressure, $P_{\mathrm{\Delta }}=k_{\mathrm{B}}{T}_{\mathrm{\Delta }}{\overline{n}}_{e,\mathrm{\Delta }}$, and for entropy, $K_{\mathrm{\Delta }}=k_{\mathrm{B}}{T}_{\mathrm{\Delta }}{\overline{n}}_{e,\mathrm{\Delta }}^{-2/3}$. Some characteristic values of these mean physical scales are quoted in Table 1 with their dependences on mass, redshift, and gas mass fraction.

2. Universal thermodynamic radial profiles

In this section we define the functional forms of the thermodynamic profiles for the ICM that will be integrated to recover the global properties to be compared with the corresponding quantities measured in X-ray observations. To define these functional forms, we assume an ICM in hydrostatic equilibrium within a spherically symmetric dark matter potential. Two simple ingredients are needed: a radial distribution for the mass, and a radial profile for one of the thermodynamic quantities (gas density, temperature, pressure, or entropy).

2.1. The semi-analytic model

To model the mass distribution we adopt a halo concentration–mass–redshift relation, c–M–z, for mass concentration c₂₀₀, global mass value $M_{200}=4/3\pi 200{\rho }_{c,z}{R}_{200}^3$, and redshift of the observed object for a Navarro-Frenk-White (NFW) mass density profile (Navarro et al. 1997) with concentration and radius at Δ = 200 related through the scale radius r_s: c₂₀₀ = R₂₀₀/r_s. We consider the c–M–z relation described in Dutton & Macciò (2014) (hereafter D14), with $log{c}_{200}=A+Blog(M_{200}/{10}^{12}{M}_{\odot }{h}_{100}^{-1})$, B = −0.101 + 0.026z, and A = 0.520 + (0.905 − 0.520)exp(−0.617z^1.21). We also use an alternative relation following the prescriptions described in Bhattacharya et al. (2013) (hereafter B13), which provides lower values of concentration by ∼10% (20%) at z ∼ 0.05 (1) in the mass range considered in the present study (10¹⁴ − 10¹⁵ M_⊙; see also Diemer & Kravtsov 2015, for a detailed comparison between different models as a function of mass and redshift). We note that the quoted relations depend on the assumed cosmological parameters, but less significantly at M >  10¹⁴ M_⊙, where the differences in predicted concentrations are on the order of a few per cent for the values of H₀ and Ω_m adopted here (see e.g. Fig. 9 in D14). We also note that the quoted c–M–z relations refer to the results extracted from dark matter only simulations. It is known that the distribution of baryons in the cluster halo can affect how concentration and total mass relate, causing for instance a steepening of the relation because star formation is fractionally more efficient in low-mass objects, with an overall larger normalisation because this effect is non-vanishing at all masses (see e.g. Fedeli 2012). However, it has been proved that, once a proper selection of N-body simulated systems is done to mimic an observational sample, the predicted c–M–z relation matches the observed one at least for very massive objects (at the 90% confidence level in the CLASH sample, where the total masses were recovered from the gravitational lensing signal; see Merten et al. 2015). Furthermore, considering that the adopted c–M–z relations were also used to infer the total masses in the Early Sunyaev-Zeldovich (ESZ) sample that will be analysed in our work, and that introducing the self-gravity due to the gas (the dominant baryonic component, although accounting for less than 15% of the total mass) would complicate, without much benefit, the calculations presented below, we use the quoted c–M–z relations as input to define the total mass of our systems.

It is worth noting that we are interested here in the global average behaviour of the cluster properties, and do not propagate any error and/or scatter on the mean quantities. The c–M–z relation is known to have an intrinsic scatter of σ _{log c} ∼ 0.16 on the predicted values of the concentration for given mass (e.g. Diemer & Kravtsov 2015) that might be propagated through the relations used in this work to evaluate its impact on the reconstructed distribution of the observed thermodynamic profiles and integrated quantities of the ICM. We postpone to a future work further discussion on the distribution of the input parameters, although we also consider a different set of parameters (both for the c–M–z relation and for the pressure profile) depending on the dynamical state of the systems, such as relaxed cooling core objects that are expected to have higher halo concentration and higher values of pressure in the core (see Sect. 4).

The predictions on the radial behaviour of the interesting physical quantities are then recovered from the inversion of the hydrostatic equilibrium equation (see e.g. Ettori et al. 2013)

$$\begin{array}{c}\hfill n_e=-\frac{d{P}_e}{\mathit{dr}}\frac{r^2}{\mu m_{a}G{M}_{t\mathrm{H}E}},\end{array}$$(2)

where P_e = n_eT_gas in units of keV cm⁻³ is described by a generalised NFW (e.g. Nagai et al. 2007; Arnaud et al. 2010)

$$\begin{array}{cc}& P_e=P_{500}{E}_z^{8/3}{\left(\frac{M_{500}}{3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot }}\right)}^{{\alpha }_M}{h}_{70}^2{P}_r\hfill \\ \hfill & P_r=\frac{P_0}{{(c_{500}x)}^{\mathrm{\Gamma }}{[1+{(c_{500}x)}^A]}^{(B-\mathrm{\Gamma })/A}}\hfill \end{array}$$(3)

with x = r/R₅₀₀, the normalisation $P_{500}\propto \mu /{\mu }_e{f}_{g,500}{M}_{500}^{2/3}=1.65\times {10}^{-3}{(M_{500}/3\times {10}^{14}{h}_{70}^{-1}{M}_{\odot })}^{2/3}$ keV cm⁻³², the parameters (P₀, c₅₀₀, Γ, A, B) equal to (6.41, 1.81, 0.31, 1.33, 4.13), and α_M = 0.12, accounting for the observed deviation from the standard self-similar scaling, set as in Planck Collaboration Int. V (2013). We have also considered an alternative pressure profile from the recent analysis of the joint XMM-Newton and Planck signals of a sample of 12 nearby massive galaxy clusters presented in Ghirardini et al. (2019). Converting the published values to feed Eq. (3), we set (P₀, c₅₀₀, Γ, A, B)) = (5.29, 1.49, 0.43, 1.33, 4.40) and α_M = 0.

The total mass M_tot is defined equal to M_HE/(1 − b), where M_tot is modelled with a NFW profile with parameters set according to the model adopted, and where the factor (1 − b) represents the hydrostatic bias that could affect the estimate of the hydrostatic masses M_HE. We note that the hydrostatic bias is propagated to the shape of the gravitating mass profile through the halo concentration, whereas the input value M_HE is adopted to define M₅₀₀, used for instance to estimate P₅₀₀ and R₅₀₀. Following this procedure, we mimic the observational bias induced from the (biased) measurement of the mass on the normalising factor P₅₀₀ and on the definition of the region over which the physical quantities are integrated. Instead, other quantities that are not directly used in the analysis (e.g. T₅₀₀) are still defined from the bias-corrected M_tot.

We show in Fig. 1 the recovered profiles following these simple prescriptions. We present a few cases of interest, normalised to their average quantity at R₅₀₀: a massive halo of 8 × 10¹⁴ M_⊙ at redshift 0.05, which will be used as reference; the same halo with a different input pressure profile and mass bias of 0.4; a system with 1/4 of the mass at z = 0.05 and 1. While the assumption of the pressure profile has a negligible impact, in particular at r >  0.1 R₅₀₀, large deviations (on the order of ∼20%) are induced by the mass (on all the profiles) and by the bias (on the temperature and density profile and, as a cumulative effect, on the entropy profile). The redshift causes measurable discrepancies, but lower than ∼20%.

Fig. 1. Reconstructed thermodynamic radial profiles for objects with (M500, z) = (8 × 1014 M⊙, 0.05), black solid line; the same, but using a pressure profile from X-COP, red solid line, the B13 c − M − z relation, purple solid line; or assuming a hydrostatic bias b = 0.4, blue solid line; (M500, z) = (2 × 1014 M⊙, 0.05), black dashed line; (2 × 1014 M⊙, 1), red dashed line. The ratio is shown with respect to the reference case, black solid line, M500 = 8 × 1014 M⊙ and z = 0.05. P500, n500, T500, K500 refer to the normalisation values presented in Table 1.

2.2. ESZ sample

To test and validate the predictions of our model, we use the radial profiles and the integrated quantities obtained for the ESZ sample (Lovisari et al. 2017, 2020), which contains the 120 galaxy clusters in the redshift range of 0.059 <  z <  0.546 observed with XMM-Newton and which were originally selected from the Planck Early Sunyaev-Zeldovich (ESZ; Planck Collaboration VIII 2011) sample. As described in Lovisari et al. (2017), these systems, which have mass and redshift distributions representing well the whole ESZ sample of 188 galaxy clusters, are the ones for which R₅₀₀ is completely covered by XMM-Newton observations.

As described in Lovisari et al. (2020), the radial temperature profiles were derived by requiring a S/N >  50 to ensure an uncertainty of ∼10% in the spectrally resolved temperature and a source-to-background count rate ratio higher than 0.6 to reduce the systematic uncertainties in their measurements. On average, the temperature profiles are extracted up to ∼R₅₀₀.

The gas density profiles are presented in Lovisari et al. (2017). They are recovered as the geometrical deprojection of the best-fit results with a double-β model of the surface brightness profile extracted from the background-subtracted vignetting-corrected image in the 0.3–2 keV band defined to maximise the signal-to-noise ratio.

We present in Fig. 2 the comparison between the observed gas density and temperature profiles with those predicted from our model for a mass and redshift equal to the median values in the sample (M₅₀₀ = 5.9 × 10¹⁴ M_⊙, z = 0.193). Both the gas density and temperature profiles predicted from the model, the latter in particular at r >  0.3 R₅₀₀, lie comfortably around the mean, and well within the scatter, of the observed values.

Fig. 2. Comparisons between the stacked profiles from the 120 objects analysed for the ESZ sample in Lovisari et al. (2020) and predictions from different models of (i) the universal pressure profile (Planck, used as reference, and X-COP) and (ii) the c − M − z relation (D14, used as reference, and B13) estimated at the median values in the sample of M500 and redshift. Top panel: electron density profiles recovered from the best-fit parameters of a double-β model. The solid line indicates the median value at each radius; the dotted lines show the 16th and 84th percentile estimated at each radius. The dashed line represents the extrapolation beyond the observational limit for the sample. Bottom panel: stacked temperature profile obtained from the weighted mean of 30 spectral points in each bin. Errors on the mean and dispersion (dotted lines) are overplotted.

3. Integrated quantities

Using the thermodynamic profiles described in the previous section, we can reconstruct other derived profiles (e.g. the gas mass and the hydrostatic mass) and the integrated quantities (e.g. global temperature and luminosity).

The self-similar scenario predicts properties of the ICM that depend on the mass (and redshift) of the halo (see e.g. Voit 2005). In Ettori (2013, 2015) (hereafter E15) we describe the standard self-similar scenario and its possible deviations in terms of physical quantities depending explicitly on the mass. However, the reconstruction of the total mass in clusters is still affected by uncertainties (see e.g. Pratt et al. 2019) that make its use as a variable of reference a possible source of systematic errors. Therefore, in this study we prefer to write the scaling relations with respect to the observed temperature T_spec, which is a direct X-ray observable, also independent from cosmology. For our analysis, T_spec is the gas temperature measured as an integrated spectrum in the radial range 0.15 − 1 R₅₀₀, and estimated from the profiles derived from our model as a “spectroscopic-like” value (see Mazzotta et al. 2004):

$$\begin{array}{c}\hfill T_{\mathrm{spec}}=\frac{\int n_e{(r)}^{2}T{(r)}^{w-0.5}dV(r)}{\int n_e{(r)}^{2}T{(r)}^{w-1.5}dV(r)}.\end{array}$$(4)

Here w = 0.75, T(r) = P_e(r)/n_e(r), and the integrals are performed over the volume V(r) of interest, either up to R₅₀₀ or between 0.15 and 1 R₅₀₀.

The X-ray luminosity is obtained for a given gas density and temperature profile, assuming an apec model in XSPEC with a metallicity fixed to 0.3 times the solar values in Anders & Grevesse (1989) and integrating the emissivity in cylindrical volume extending up to 3 × R₅₀₀ along the line of sight and covering an aperture between 0.15 and 1 R₅₀₀.

In the following we adopt the general description presented in E15. We define the total mass ℳ ≡ E_zM_tot/M₀ (M₀ = 5 × 10¹⁴ M_⊙), the gas mass ℳ_g ≡ E_zM_g/M_g, 0 (M_g, 0 = 5 × 10¹³ M_⊙), and a bolometric luminosity $\mathcal{L}\equiv E_z^{-1}L/L_0$ (L₀ = 5 × 10⁴⁴ erg s⁻¹). We write the normalisations and slopes of the scaling relations that relate these quantities to 𝒯 ≡ k_BT_spec/T₀ (T₀ = 5 keV) as

$$\begin{array}{cc}& \mathcal{M}=k_M{(1-b)}^{-1}{f}_T^{3/2}{\mathcal{T}}^{3/2},\hfill \\ \hfill & {\mathcal{M}}_g=k_M{f}_g{f}_T^{3/2}{\mathcal{T}}^{3/2},\hfill \\ \hfill & \mathcal{L}=k_L{f}_g^2{f}_T^{3/2}{\mathcal{T}}^2,\hfill \end{array}$$(5)

where we have defined the following parameters: f_T = T(R₅₀₀)/T₅₀₀  ×  T₅₀₀/T_spec relates the gas temperature at R₅₀₀ with the observed global value T_spec; the X-ray–measured gas mass fraction f_g = M_g/M_HE = (1 − b)⁻¹M_g/M_tot = ≡C^0.5 f_nc/f_g, 0 is normalised to f_g, 0 = 0.1 and is related to the clumping-free gas mass fraction f_nc through the clumping factor $C=\langle n_{\text{gas}}^2\rangle /{\langle n_{\text{gas}}\rangle }^2$, defined as the ratio of the average squared gas density to the square of the mean gas density and that affects the measurement of the gas density as obtained from the deprojection of the X-ray free-free emission (e.g. Nagai & Lau 2011; Roncarelli et al. 2013; Vazza et al. 2013; Eckert et al. 2015); k_M and k_L (here referring to the case of a bolometric luminosity integrated over a spherical volume up to R₅₀₀) are defined by constants (e.g. the overdensity Δ) and parameters that describe the shape of the thermodynamic profiles and are fully described in Appendix A. We note that the true gas mass fraction, i.e. the unbiased value from the hydrostatic bias and the clumping factor, is written in our notation as f_gas = M_g/M_tot = (1 − b) f_g = (1 − b) C^0.5 f_nc.

As extensively discussed in E15, any deviation from the standard self-similar behaviour can be ascribed to processes that impact the relative distribution of the gas and that, in the present study, we assume to be described as a power-law dependence on the gas temperature measured in the spectroscopic analysis 𝒯 ≡ k_BT_spec/5 keV:

$$\begin{array}{cc}& f_g=f_0{\mathcal{T}}^{f_1},\hfill \\ \hfill & f_T=t_0{\mathcal{T}}^{t_1},\hfill \\ \hfill & (1-b)=h_0.\hfill \end{array}$$(6)

Here f₀ is chosen to be 0.1 and is roughly representative of measured gas mass fractions, while t₀ is set to $\overline{f_T}=0.71$, the median value for our model as derived in Appendix A. We list these parameters, with their descriptions and definitions, in Table 2.

We insert them in Eq. (5) to write these equations in their full dependence on T as

$$\begin{array}{cc}& \mathcal{M}=k_M{h}_0^{-1}{t}_0^{3/2}{\mathcal{T}}^{3/2+3/2t_1}=N_{\mathit{MT}}{\mathcal{T}}^{S_{\mathit{MT}}},\hfill \\ \hfill & {\mathcal{M}}_g=k_M{f}_0{t}_0^{3/2}{\mathcal{T}}^{3/2+3/2t_1+f_1}=N_{\mathit{GT}}{\mathcal{T}}^{S_{\mathit{GT}}},\hfill \\ \hfill & \mathcal{L}=k_L{f}_0^2{t}_0^{3/2}{\mathcal{T}}^{2+3/2t_1+2f_1}=N_{\mathit{LT}}{\mathcal{T}}^{S_{\mathit{LT}}},\hfill \end{array}$$(7)

where we represent with N and S the best-fit estimates of normalisation and slope, respectively, of the corresponding scaling laws that are quoted in Table 3.

Then, by simple algebraic calculations, we can invert Eq. (7) and obtain the normalisations f₀, t₀, h₀ and slopes f₁, t₁:

$$\begin{array}{cc}& f_0=h_0^{-1}{N}_{\mathit{GT}}/N_{\mathit{MT}}=(N_{\mathit{LT}}/k_L){(N_{\mathit{GT}}/k_M)}^{-1},\hfill \\ \hfill & t_0={(h_0{N}_{\mathit{MT}}/k_M)}^{2/3}={(N_{\mathit{GT}}/k_M)}^{4/3}{(N_{\mathit{LT}}/k_L)}^{-2/3},\hfill \\ \hfill & h_0={(N_{\mathit{GT}}/k_M)}^2{(N_{\mathit{MT}}/k_M)}^{-1}{(N_{\mathit{LT}}/k_L)}^{-1},\hfill \\ \hfill & t_1=2/3S_{\mathit{MT}}-1,\hfill \\ \hfill & f_1=S_{\mathit{GT}}-S_{\mathit{MT}}.\hfill \end{array}$$(8)

In Fig. 3, we show the reconstructed relations with the ICM temperature of the gas mass fraction, gas mass, and bolometric luminosity, and the corresponding best-fit relations that allow us to evaluate the normalisations N and slopes S. These values are then used to check the predicted behaviour of some intrinsic properties (e.g. the dependence on T of the gas mass fraction, and of the ratio T(R₅₀₀)/T; see panels at the bottom of Fig. 3).

Fig. 3. Top panels: global properties (empty black dots) and after correction for the Ez factor (solid red dots) recovered from their thermodynamic profiles for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), (2 × 1014 M⊙, 1). The dashed magenta lines identify the best-fit relations of Eq. (7). A core-excised bolometric luminosity is considered. Bottom panels: red dots are the quantities (from left to right: fg, fT, b) estimated in our model for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), and (2 × 1014 M⊙, 1); the blue dotted line indicates Ωb/Ωm = 0.157 (Planck Collaboration XXVII 2016); the dashed magenta lines show the predictions following the relations: fg = 0.124 𝒯0.23, fT = 0.687 𝒯0.07 and 1 − b = h0 = 1, where fg and fT are the absolute values, i.e. multiplied by the assumed normalisations of 0.1 and $\overline{f_T}=0.71$, respectively.

It is worth noting that, in order to compare our predictions with robust estimates of the X-ray luminosity, we integrate the bolometric luminosity within the radial range 0.15 − 1 R₅₀₀, using quantities projected along the line of sight to mimic the observational values. These characteristics of the luminosity L (i.e. if projected or integrated over spherical shells, the radial range of integration, the energy band) define the proper value of the constant k_L. Given the constraints on the L − T relation, we can invert it to recover k_L under the assumption that there is no hydrostatic bias (i.e. (1 − b) = h₀ = 1), so that t₀ = (N_MT/k_M)^2/3, f₀ = N_GT/N_MT, and $k_L=N_{LT}{f}_0^{-2}{t}_0^{-3/2}=0.503$. We adopt this value of k_L in the following analysis.

3.1. Calibration of the model in E15

Our semi-analytic model makes use only of a universal pressure profile and a c–M–z relation, under the assumption that the total mass profile is described by a NFW model. The integrated quantities predicted from this model are used to calibrate the relations presented in E15 and described by Eqs. (5)–(8). By fitting in a robust way the distribution of points plotted in Fig. 3 when b is assumed to be zero, we obtain from Eq. (8) that these integrated quantities relate between them following the scaling laws built from self-similar relations and modified by including the following dependence of the temperature profile f_T and of the gas mass fraction f_g = C^0.5f_nc on T_{spec, 0.15 − 1 R500}:

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_T=0.687(0.711)\times {(T/5\mathrm{keV})}^{0.07(0.10)}\hfill \\ f_g=0.124(0.128)\times {(T/5\mathrm{keV})}^{0.23(0.25)}.\hfill \end{array}\end{array}$$(9)

Here we quote in parentheses the results obtained by assuming a B13 model (where D14 is our model of reference) for the c–M–z relation.

3.2. Role of the hydrostatic bias

As described above, the hydrostatic bias is propagated to the total mass as M_tot = M_HE/(1 − b), where M_tot is modelled with the adopted NFW profile, and the hydrostatic bias modifies the shape of the gravitating mass profile through the halo concentration. This affects all the thermodynamic quantities depending on the dark matter distribution, (e.g. T₅₀₀) that is defined from the bias-corrected M_tot. On the other hand, the input value M_HE is used to rescale the observed properties (e.g. R₅₀₀ and P₅₀₀) to mimic the observational bias once a biased M_HE is measured instead of the true value M_tot.

In Fig. 4, we represent the impact of the hydrostatic bias as the ratio between quantities estimated with a given b value and with b = 0. A clear trend is present, with higher biases producing higher global temperatures, and lower estimates of M_g and L, mostly as consequence of the increased halo mass, associated with the corresponding reduction in the halo concentration and constant value of R₅₀₀. We have also modelled the trend we observe in the physical quantities as a function of the assumed bias b with a functional form

$$\begin{array}{c}\hfill Y/Y_{b=0}={\alpha }_0(1-b)+(1-{\alpha }_0){(1-b)}^{{\alpha }_1},\end{array}$$(10)

Fig. 4. Impact of the hydrostatic bias b on the integrated quantities, represented as the ratio between the quantities estimated for an assigned bias, and with bias equal to zero. From left to right: this ratio for total mass, temperature, gas mass, core-excised bolometric luminosity, and fT = T(R500)/Tspec. Shown is the case of a system with M500 = 5 × 1014 M⊙ at z = 0.05 and 0.5 (the values at different redshifts overlap for most of the quantities).

such that when Y = M_tot, then (α₀, α₁) = (0, −1). For a typical object with M₅₀₀ = 5 × 10¹⁴ M_⊙ at z = 0.05 shown in Fig. 4, this functional form reproduces the observed trends within 2% and provides the following best-fit parameters (α₀, α₁): ( − 0.12, −0.45) for T_spec; (0.34, 0.39) for the gas mass; (0.83, 0.25) for the bolometric luminosity; (0.012, 2.61) for the gas fraction. Very similar results are obtained for systems at lower masses and higher redshifts on T_spec, M_g, and gas mass fraction.

This translates into the following representation in the E15 formalism for the two cases b = 0.2 and b = 0.4 (the latter in brackets), assuming the same set of input values of (M₅₀₀, z)

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_T=0.597(0.500)\times {(T/5\mathrm{keV})}^{0.08(0.10)}\hfill \\ f_g=0.106(0.087)\times {(T/5\mathrm{keV})}^{0.24(0.25)},\hfill \end{array}\end{array}$$(11)

with lower normalisations for f_g and f_T for higher values of b, and no significant change in the dependence on the temperature.

3.3. Results on the ESZ sample

Lovisari et al. (2020) present a study of the X-ray scaling relations for the ESZ sample, based on hydrostatic mass profiles recovered from the B13 c–M–z relation, and using core-excised luminosities and spectroscopic temperatures, as we reproduce in our semi-analytic model. The best-fit relations for the observational data are corrected for the Eddington bias, but not for the Malmquist bias, which is negligible when fitting the X-ray properties of an SZ selected sample (for more details, see Lovisari et al. 2020). The slopes of all the investigated scaling relations are found to deviate significantly from the self-similar predictions, if self-similar redshift evolution is assumed. We re-estimated normalisations and slopes of the scaling relations of interest using the package LIRA (Sereno 2016; see Table 3). In Fig. 5, we overplot these best-fit relations to the data points in the ESZ sample.

Fig. 5. Distribution of the observed values for the ESZ sample with best-fit results from a linear fit in the logarithmic space (using LIRA; Sereno 2016). Predictions from the semi-analytic model (best-fit values in Table 3) are overplotted with dashed lines in magenta (for the D14 c − M − z relation) and blue (B13). Bottom panels: dashed lines indicate the ratios of the model (M) to the observed (O) best-fit relations.

By applying our universal model, we require the self-similar predictions to be corrected by the following dependences on the gas temperature

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_{T,\mathrm{ESZ}}=0.697(\pm 0.103)\times {(T/5\mathrm{keV})}^{0.15(\pm 0.06)}\hfill \\ f_{g,\mathrm{ESZ}}=0.121(\pm 0.045)\times {(T/5\mathrm{keV})}^{0.45(\pm 0.09)}.\hfill \end{array}\end{array}$$(12)

Here we propagate the errors on the best-fit parameters, and (as in the following analysis) we quote the absolute values of f_g and f_T, i.e. multiplied by their normalisation values of f_g, 0 = 0.1 and $\overline{f_T}=0.71$.

Then we compare the observed distribution with the predictions from our semi-analytic model. To mimic the distribution in mass and redshift of the ESZ sample, we simulated ten objects with our semi-analytic code having mass and redshift equal to the median values estimated in bins of 12 clusters each, sorted in redshift. We overplot these predictions to the data points in Fig. 5. We also use a different c–M–z relation to probe the dependence on the assumed model. As previously discussed, the B13 model predicts lower c₂₀₀ at higher redshift, for a fixed halo mass. This behaviour induces differences in the estimates of T, L, and gas mass. For example, a halo with M₅₀₀ = 5 × 10¹⁴ M_⊙ is expected to have c₂₀₀ lower by 5% and 13% at z = 0.05 and 0.5 in B13, inducing a slightly different distribution of the dark matter and consequent reshaping of the thermodynamic profiles that produce T lower by 2% and 4%, M_g higher by 2% and 4% and L higher by 3% and 9%, respectively.

The normalisations are within the 1σ intervals of the observational constraints (the largest tension being the results on N_LT from the D14 model), whereas the slopes of the scaling laws in the ESZ sample tend to be systematically higher than the predicted values by 0.7–1.2σ, 2.6–2.1σ, and 3–2.4σ for the M − T, M_g − T, and L − T relation with D14 and B13 models, respectively. It is important to note, however, that the steepening of the relations might depend on the selection applied and on how the fit is performed (in our case, we consider in LIRA the scatter on both the variables X and Y, i.e. we leave the parameters sigma.XIZ.0 and sigma.YIZ.0 free to vary; for details, see Sereno 2016). For instance, if we fix sigma.XIZ.0 = 0, then the relations become flatter, with slopes of 1.60, 2.00, and 2.82 for the M − T, M_g − T, and L − T relation, respectively. In general, we obtain a better agreement between the observed distributions in the ESZ sample and the simulated data when we apply the B13 model as a consequence of the shift to lower temperatures for a given mass, and conclude that the distribution of the integrated physical quantities follows what is reconstructed from the universal model within 10%, on average.

Using Eq. (8) on these simulated data, and assuming a D14 model, the best-fit values convert into (t₀, t₁) = (0.689, 0.07); (f₀, f₁) = (0.122, 0.22) for h₀ free to vary (best-fit: 1.02) and (t₀, t₁) = (0.680, 0.07); (f₀, f₁) = (0.124, 0.22) for h₀ fixed to one. When B13 is adopted we obtain h₀ = 1.01 once it is left to vary, (t₀, t₁) = (0.698, 0.10); (f₀, f₁) = (0.125, 0.25), and (t₀, t₁) = (0.694, 0.10); (f₀, f₁) = (0.126, 0.25) with h₀ fixed to one. Overall, there is a good consistency with the results shown in Eq. (12), apart from a clear steeper dependence of the gas fraction on T that will be discussed further in the next subsection.

3.4. Including the distribution of Y_SZ

The same plasma responsible for the X-ray emission can also be traced through the Sunyaev-Zeldovich (SZ) effect, generated from the Compton scattering of the photons of the cosmic microwave background on the electrons of the ICM (Sunyaev & Zeldovich 1972). The millimetre wave emission due to the thermal SZ effect is proportional to the integrated pressure of the X-ray emitting plasma along the line of sight and is described, as aperture integrated signal (see e.g. Mroczkowski et al. 2019), by the integrated Compton parameter $Y_{\text{SZ}}{D}_A^2=({\sigma }_T/m_e{c}^2)\int PdV$, where D_A is the angular diameter distance to the cluster, σ_T = 8π/3(e²/m_ec²)² = 6.65 × 10⁻²⁵ cm² is the Thomson cross section, m_e and e are respectively the electron rest mass and charge, c is the speed of light, and P = n_eT is the electron pressure profile.

Using the formalism presented in E15, this signal relates to the gas temperature as

$$\begin{array}{c}\hfill \frac{E_z{Y}_{\mathrm{SZ}}{D}_A^2}{{10}^{-4}{\mathrm{Mpc}}^2}=k_{\mathrm{SZ}}{f}_T^{5/2}{f}_{\mathit{nc}}{\mathcal{T}}^{5/2}=N_{\mathrm{YT}}{\mathcal{T}}^{S_{\mathrm{YT}}},\end{array}$$(13)

where k_SZ is derived in Appendix B.

It is worth noting a few issues concerning the reconstruction of the SZ signal. First, the SZ signal depends linearly on the pressure profile and is sensitive to the upper limit of the integration. We collect the estimates of the spherically integrated SZ flux up to R₅₀₀Y_SZ from the PSZ1 catalogue³ (Planck Collaboration XIV 2015). We refer to Planck Collaboration XIII (2016) for an exhaustive discussion on how Y_SZ is recovered from the integrated Comptonisation Y_5 R500, which represents a nearly unbiased proxy for the total SZ flux within a cylinder of aperture radius 5 R₅₀₀ and needs to be corrected for the degeneracy induced from the signal-size correlation and underlying variations of the pressure profile of reference inducing extra scatter and bias in the extrapolation (see Sects. 5.1–5.3 and Fig. 16 in Planck Collaboration XIII 2016).

A second issue is that the gas fraction f_nc in Eq. (13) refers to the clumping-free value, i.e. a gas mass fraction that does not depend on the gas clumping described by the factor C that appears in Eq. (5) above. By combining now the scaling relations based on SZ data and those using X-ray quantities, and modelling the gas clumping as C = C₀ 𝒯^C₁, we can write $f_g=C^{0.5}{f}_{nc}=C_0^{0.5}{f}_{nc,0}{\mathcal{T}}^{f_{nc,1}+0.5C_1}=f_0{\mathcal{T}}^{f_1}$, where f_nc, 0 and f_nc, 1 indicate the normalisation and slope, respectively, of the gas fraction after the correction for the clumping. After simple calculations we obtain

$$\begin{array}{c}\hfill \{\begin{array}{c}{C}_0^{0.5}=f_0{t}_0^{5/2}{(N_{\mathrm{YT}}/k_{\mathrm{SZ}})}^{-1}\hfill \\ C_1=2f_1-2S_{\mathrm{YT}}+5+5t_1\hfill \\ f_{nc,0}=f_0{C}_0^{-0.5}\hfill \\ f_{nc,1}=f_1-0.5C_1,\hfill \end{array}\end{array}$$(14)

and proceed with the estimates of the gas clumping and the corrected gas fraction using the values of f_T and f_g quoted in Eqs. (9) and (12).

We present the results of our analysis in Fig. 6. By fitting a linear relation in the log space on the sample of ten objects simulated with our semi-analytic model to mimic the observed distribution in mass and redshift of the ESZ sample, we measure (N_YT, S_YT) = (0.33, 2.84) when a D14 c–M–z relation is assumed, and (0.35, 2.92) for the B13 relation. Using the relations imposed by our universal model (see Eq. (14)), we constrain the gas clumping to be (C₀, C₁) = (1.00, 0.12) and (1.04, 0.17) for the D14 and B13 model, respectively. Considering that our model does not include any gas clumping, we can consider these values as indicators of the systematic uncertainties (on the order of a few per cent) that affect our reconstruction of the ICM properties.

Fig. 6. YSZ − T relation for the ESZ sample with the best-fit results from a linear fit in the logarithmic space. The fit on the observed data is performed using LIRA (black dotted line). The predictions from the semi-analytic model overlap (magenta dashed line: D14 model; blue dashed line: B13 model). The dot-dashed orange line is obtained from the propagation of the scaling in the E15 formalism (see Eq. (13)).

When the same analysis is applied to the ESZ sample, we estimate (N_YT, S_YT) = (0.42 ± 0.03, 2.79 ± 0.18). These results are consistent with the best-fit constraints obtained for 62 objects listed in a Planck Early Results work (Planck Collaboration XI 2011, redoing the fit, we measure a normalisation in the adopted units of 0.37 ± 0.01 and slope of 2.97 ± 0.18). Adopting as reference the best-fit values for the ESZ sample, we investigate how we can reconcile the observed differences with our predictions. We note that the normalisation N_YT depends on the gas fraction and f_T. As discussed in Sect. 3.2, any hydrostatic bias induces lower values of f_g and f_T, decreasing the normalisation and enlarging the tension with the observed constraints. By combining the constraints on f_g and f_T with the relations in Eq. (14), we conclude that the ESZ dataset is consistent with a gas clumping C₀ <  1.4 (at the 1σ confidence level), in close agreement with the expected values within R₅₀₀ from hydrodynamical simulations (see e.g. Nagai & Lau 2011; Roncarelli et al. 2013; Vazza et al. 2013; Eckert et al. 2015). From the best-fit results, we can also constrain C₁ ∼ 1.0(±0.5), which represents a further contribution to the dependence on T of the gas fraction, explaining the steeper dependence quoted in Eq. (12) with respect to the predictions from the model.

4. Discussion

Our model combines a universal profile for the gas pressure with the present knowledge on the distribution of the dark matter in galaxy clusters to calibrate properly any deviations from the standard self-similar scenario. In particular, three quantities are introduced to account for these deviations: a temperature-dependent gas mass fraction, f_g; a temperature-dependent ratio between the temperature at given radius and its global value, f_T; a hydrostatic bias, b. To explain self-consistently the scaling relations observed in our model and those recovered in the ESZ sample, we require a significant dependence on T of the gas fraction, and a milder dependence of f_T.

As numerical simulations suggest, and also analyses of the observed gas density distribution as a function of the measured temperature, more massive (hotter) systems tend to have a relatively higher gas density in the cores than groups (see e.g. Croston et al. 2008; Eckert et al. 2016), where feedback from central AGNs can efficiently contrast the attraction from gravity and push gas particles beyond R₅₀₀. The net effect is a reduction of the gas mass fraction in objects at lower masses. This variation of the gas mass fraction with the gas temperature (or total mass) is well documented in past works, mostly based on X-ray selected samples (e.g. Pratt et al. 2009; Eckert et al. 2013, 2016; Lovisari et al. 2015; Ettori 2015), with constraints on normalisations and slopes that are quite similar, but not identical, to the values we obtain in our study and summarise in Eq. (12). For example, by converting to our units, Pratt et al. (2009) estimate f_g = 0.107 × (T/5 keV)^0.36; Lovisari et al. (2015) measure f_g = 0.106 × (T/5 keV)^0.32; Ettori (2015) obtains f_g = 0.107 × (T/5 keV)^0.33; and Eckert et al. (2016) obtain f_g = 0.079 × (T/5 keV)^0.35, which is on the lower side probably due to a bias in the weak-lensing mass measurements (see e.g. Umetsu et al. 2020). Using the estimates within R₅₀₀ of the gas mass, total mass, and temperature obtained for the ESZ sample, we can directly fit the gas fraction–temperature relation and obtain f_g = 0.123(±0.002) × (T/5 keV)^0.39(±0.05), in agreement with our indirect calculations. However, these values of the gas fraction at a given temperature are 15-20% higher than those from the literature cited above. As discussed in Lovisari et al. (2020), part of this offset can be attributed to the nature of the ESZ sample, containing a larger number of disturbed clusters for which higher gas fractions are measured (see also Eckert et al. 2013) as a cumulative effect of gas clumpiness and inhomogeneities (inducing higher gas mass), larger non-thermal contribution to the total pressure (biasing low the total mass), and possible violation of the hydrostatic equilibrium. It is indeed known that SZ selected samples tend to have a larger contribution of dynamically disturbed clusters than X-ray selected ones, as also confirmed for the ESZ sample (see e.g. Lovisari et al. 2017).

We can evaluate the impact of having more relaxed objects on our relations. We repeat the analysis by assuming a pressure profile for cool-core systems, with (P₀, c₅₀₀, Γ, A, B) = (11.82, 0.60, 0.31, 0.76, 6.58) (see Planck Collaboration Int. V 2013), and using the c–M–z relation for relaxed systems in Bhattacharya et al. (2013). Still assuming no hydrostatic bias (b = 0), higher mass concentrations associated with relaxed objects produce higher normalisations of f_g by 5–10%. To compensate for this rise affecting more X-ray selected samples, and to further reduce their estimates of f_g, we are forced to require the presence of a hydrostatic bias of about 0.2. As we discuss in Sect. 3.2, the presence of any hydrostatic bias propagates to all the derived quantities, lowering f_g for a given initial set of (M₅₀₀, z). A bias of 0.2 causes a reduction of 15% or more in the normalisation of f_g, allowing us, by compensating for the increase due to a larger contribution of relaxed clusters, to match the published results based on X-ray selected samples. A similar (or even higher) bias should also be present in the ESZ sample due to the larger contribution of disturbed systems. The presence of more disturbed systems should induce larger inhomogeneities in the gas distribution, which would bias high the reconstructed gas density and the corresponding gas mass fraction (e.g. Roncarelli et al. 2013), so that any correction for it would move f_g to lower values.

On the other hand, the results on f_g in the ESZ sample matches the predictions of our semi-analytic model with b = 0. It is worth noting that the adopted input models describe the ESZ sample best because we consider the average pressure profile recovered for Planck selected systems and use the c–M–z relation estimated for the entire collection of haloes (i.e. also including more disturbed systems). So, it seems that a sort of conspiracy acts to mimic f_g with b = 0, although we are fairly confident that the ESZ sample contains a larger number of disturbed objects than in X-ray selected samples.

As we discuss in Sect. 3.2, the hydrostatic bias impacts all the integrated quantities. We evaluate how it affects the normalisation of f_g (see also Eq. (11)), and model it as

$$\begin{array}{c}\hfill \frac{f_g}{f_{g,b=0}}={\alpha }_0(1-b)+(1-{\alpha }_0){(1-b)}^{{\alpha }_1}\end{array}$$(15)

with (α₀, α₁) = (0.48, 0.44). Assuming that $f_g=C_0^{0.5}{f}_{g,c}\approx 0.123$ is the observed value including the gas clumping, and f_{g, b = 0} = 0.124 (0.126) is set from our models (see end of Sect. 3.3), by adopting the upper limit on the gas clumping obtained from the Y_SZ − T relation (see Sect. 3.4), we conclude that the level of hydrostatic bias allowed in the ESZ sample has to be below 0.24 (0.26) at the 1σ level for the D14 (B13) model (see Fig. 7).

Fig. 7. Dependence of the normalisation of fg upon the bias b (red line; see Eq. (15)). The shaded region (and the region between dashed lines) represents the quantity fnc/fg, b = 0, where $f_{nc}=f_g/C_0^{0.5}$, with fg ≈ 0.123, fg, b = 0 = 0.124 (0.126) and 1 <  C0 <  1.4 being the 1σ upper limit on the level of gas clumping.

The evidence that the ESZ sample includes disturbed systems might then account for the measured f_g by a combination of the presence of hydrostatic bias (b >  0) and gas clumping (C >  1) induced by large gas inhomogeneities.

The other quantity we introduced, f_T = T(R₅₀₀)/T₅₀₀ × T₅₀₀/T_spec = T(R₅₀₀)/T_spec, has never been investigated before and comes from our request to build a proper normalisation for the M − T relation, where the hydrostatic mass depends on the measurement of the temperature at a given radius (T(R₅₀₀)). As we show in Fig. 3, f_T is expected to increase mildly in the hotter systems, with a global effect that might account for some of the deviations observed in the standard self-similar scenario. Moreover, the normalisation decreases with increasing hydrostatic bias b (see Fig. A.1).

5. Conclusions

The gravitational force aggregates matter onto galaxy clusters, driving them to the virialisation and regulating the distribution and the energetic budget of the accreted baryons and the emergent observational properties of the ICM. We show how a universal pressure profile of the ICM, combined with a halo mass – concentration – redshift relation (either the model in D14 or in B13) and the hydrostatic equilibrium equation, allows the radial profiles of the thermodynamic quantities to be reconstructed. Once integrated over some typical scale (e.g. R₅₀₀), these quantities produces observables (such as gas mass, temperature, luminosity, total mass, and SZ Compton parameter), which satisfy universal scaling laws that are the simple combination of self-similar relations, regulated from the mass and redshift of the halo, and further dependences of the gas temperature measured at the radius of reference, T(R₅₀₀), and of the gas mass fraction, f_g = C^0.5f_nc, on the observed spectroscopic measurement, T_{spec, 0.15 − 1 R500}.

We calibrated these dependences both within our framework and using one of the largest samples of X-ray luminous galaxy clusters that has been homogeneously analysed, the ESZ sample (Lovisari et al. 2017, 2020). We demonstrate that self-similar scaling laws hold between the integrated X-ray observables, once f_T = T(R₅₀₀)/T₅₀₀ × T₅₀₀/T_spec and f_g = C^0.5 f_nc are allowed to depend on T ≡ T_spec as f_T = t₀𝒯^t₁ and f_g = f₀𝒯^f₁ (see also Table 2 for a description of these quantities) with the following constraints on the parameters obtained from our semi-analytic model (values for the B13 c–M–z relation in parentheses):

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_T=0.687(0.711)\times {(T/5\mathrm{keV})}^{0.07(0.10)}\hfill \\ f_g=0.124(0.128)\times {(T/5\mathrm{keV})}^{0.23(0.25)}\hfill \end{array}.\end{array}$$(16)

In the ESZ sample, by propagating these dependences to the SZ signal, and interpreting any mismatch between the best-fit results of the observational data and the predictions from our semi-analytic and theoretical models as a measure of the gas clumping, we estimate

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_{T,\mathrm{ESZ}}=0.697(\pm 0.103)\times {(T/5\mathrm{keV})}^{0.15(\pm 0.06)}\hfill \\ f_{g,\mathrm{ESZ}}=0.121(\pm 0.045)\times {(T/5\mathrm{keV})}^{0.45(\pm 0.09)}\hfill \\ C_{\mathrm{ESZ}}=(<1.4)\times {(T/5\mathrm{keV})}^{1.0(\pm 0.5)}\hfill \end{array}.\end{array}$$(17)

This level of gas clumping within R₅₀₀ is consistent with the value measured in hydrodynamical simulations, and supports the evidence that SZ selected samples, like the ESZ one, tend to have a representative contribution of dynamically disturbed clusters (see e.g. Lovisari et al. 2017) with a relative significant presence of gas inhomogeneities, in particular at large scales (see e.g. Roncarelli et al. 2013; Vazza et al. 2013). This upper limit on the gas clumping implies a level of hydrostatic bias b (see Eq. (15) and Fig. 7) below 25%, on average.

We conclude that our semi-analytic model reproduces well the observed properties of galaxy clusters, both resolved spatially and as integrated quantities. By providing the calibration of this physically motivated model in terms of the standard self-similar scaling relations, modified by three components that are able to account completely for the observed deviations, we can also deduce some other interesting properties.

For instance, these explicit forms of the dependences on the observed gas temperature allow us to modify accordingly the characteristic physical quantities that renormalises the observed profiles. In Table 4 we quote how the dependences implied by the universal model propagate and modify some of these quantities. For example, P₅₀₀, which in a self-similar model scales as $M_{500}^{2/3}$, is expected to scale (in parentheses, the values for the B13 model) as T^1.29(1.35) or $M_{500}^{0.81(0.82)}$, which is consistent with the rescaling of $M_{500}^{2/3+0.12}\approx M_{500}^{0.79}$ suggested in Arnaud et al. (2010) (and adopted in Planck Collaboration Int. V 2013).

Moreover, we can write the baryonic depletion parameter Y with its explicit dependence on the gas temperature (or mass) as

$$\begin{array}{c}\hfill Y\frac{{\mathrm{\Omega }}_b}{{\mathrm{\Omega }}_m}=f_g+f_{\mathrm{s}},\end{array}$$(18)

where f_s represents the stellar mass fraction (∼0.015; see e.g. Eckert et al. 2019), and the ratio of the cosmic baryon density to matter density parameters, Ω_b/Ω_m, is equal to 0.157 (Planck Collaboration XXVII 2016). From the results of our semi-analytic model with a D14 (B13) c–M–z relation, we predict an average value of Y = 0.787(0.813) (T/5 keV)^0.23(0.25) + 0.095.

In the same framework, we show that any calibration of the M_tot − T relation with mass measurements that do not rely on the hydrostatic equilibrium equation can be used to constrain the hydrostatic bias b, preserving the use of the self-similar relations in this case as well. On the other hand, by controlling the bias b we can make predictions on the expected variations in the observed properties and, for example, explain the published relations between gas mass fraction and temperature (see Sect. 4) by requiring a hydrostatic bias b ≳ 0.2, as illustrated in Sect. 3.2.

The application of this model to a accurately selected sample of a large number (∼100) of objects analysed homogeneously in their X-ray and lensing signal out to R₅₀₀ and beyond as the one that will soon be available for the XMM-Newton Heritage Galaxy Cluster Project⁴ will allow us to extend the calibration and our understanding of the physical processes that regulate the interplay of baryons and dark matter in galaxy clusters.

Overall, our study demonstrates that the observed properties, both spatially resolved and integrated values, of X-ray luminous galaxy clusters are understood well. This represents a further step in the process to standardise their observables on the basis of a physically motivated model, not only to fully appreciate the phenomena that shape the distribution of baryons and regulate their energetic budget, but also to control biases that could affect their use as cosmological proxies.

Acknowledgments

We thank the anonymous referee for insightful comments that helped in improving the presentation of the work. We acknowledge financial contribution from the contracts ASI-INAF Athena 2015-046-R.0, ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), and from INAF “Call per interventi aggiuntivi a sostegno della ricerca di main stream di INAF”. This research has made use of the SZ-Cluster Database operated by the Integrated Data and Operation Center (IDOC) at the Institut d’Astrophysique Spatiale (IAS) under contract with CNES and CNRS.

References

Appendix A: Derivations of k_M and k_L

In Eq. (5) we introduce some constants that appear in the normalisations of the scaling relations M_tot − T and L − T, namely k_M and k_L. We present here the steps that lead to their definitions.

The hydrostatic mass at the given overdensity Δ is defined as

$$\begin{array}{c}\hfill M_{\mathrm{HE}}=M_{\mathrm{tot}}(1-b)=-\frac{f_T{T}_{\mathrm{spec}}{R}_{\mathrm{\Delta }}}{\mu m_{a}G}\frac{dlogP}{dlogr},\end{array}$$(A.1)

where we have defined the quantity f_T = T(R_Δ)/T_spec that relates the gas temperature at R_Δ with the observed global value T_spec. By definition, the same mass is equal to $M_{\text{HE}}=4/3\pi {\rho }_{c,z}\mathrm{\Delta }{R}_{\mathrm{\Delta }}^3$, where ${\rho }_{c,z}=3H_z^2/(8\pi G)$. By inverting the latter equation, we obtain $R_{\mathrm{\Delta }}^3=2G{M}_{\text{HE}}/(\mathrm{\Delta }{H}_z^2)$ that we use to replace R_Δ in Eq. (A.1). After this substitution, we solve for M_tot and obtain

$$\begin{array}{c}\hfill M_{\mathrm{tot}}={\left(\frac{2G}{H_z^2\mathrm{\Delta }}\right)}^{1/2}{\left(\frac{f_T\beta }{\mu m_{a}G}\right)}^{3/2}{(1-b)}^{-1}{T}_{\mathrm{spec}}^{3/2},\end{array}$$(A.2)

where we define β as the logarithmic slope of the pressure profile at R_Δ (β = −d log P/d log r⁵). Finally, collecting all the constants on the right side of the equation, and recalling that H_z = E_zH₀, we can write

$$\begin{array}{c}\hfill E_z{M}_{\mathrm{tot}}={\left(\frac{f_T\beta }{\mu m_a}\right)}^{3/2}\frac{2^{1/2}{(1-b)}^{-1}}{H_{0}G{\mathrm{\Delta }}^{1/2}}{T}_{\mathrm{spec}}^{3/2},\end{array}$$(A.3)

which allows us to define the constant k_M as the normalisation of the relation

$$\begin{array}{c}\hfill E_z(M_{\mathrm{tot}}/M_0)=k_M{f}_T^{3/2}{(1-b)}^{-1}{(T_{\mathrm{spec}}/T_0)}^{3/2},\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill k_M={\left(\frac{\overline{f_T}\overline{\beta }}{\mu m_a}\right)}^{3/2}\frac{2^{1/2}}{H_{0}G{\mathrm{\Delta }}^{1/2}}\frac{T_0^{3/2}}{M_0}.\end{array}$$(A.5)

Using Δ = 500 and the median values $\overline{\beta }=2.94$ (identical for all the systems being the same as the pressure profile adopted) and $\overline{f_T}=0.71$ (with values in the range of −6%, +2% for b = 0; see Fig. A.1), and the adopted pivot values T₀ and M₀, we measure k_M = 0.883.

Fig. A.1. Distribution of the parameters β, fT, fg, and 1 − b as a function of T and colour-coded according to the assumed hydrostatic bias (black: b = 0; green: b = 0.2; red: b = 0.4). The dashed lines in the panel at the bottom right are recovered through the process described in Eq. (8). Points represent quantities derived for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), with and without a bias b; a further case (2 × 1014 M⊙, 1) with b = 0 is considered.

Consider now the L − T relation. By definition, the X-ray luminosity is equal to the integral over the volume V of interest of the X-ray emissivity ϵ = n_en_pΛ(T, Z):

$$\begin{array}{c}\hfill L=\int ϵdV.\end{array}$$(A.6)

A selection of the X-ray band where the luminosity is measured will affect the estimate of the cooling function Λ(T, Z), which depends on the gas temperature and metallicity. For a bolometric luminosity we can write $ϵ=c_f{n}_e^2{T}_{\mathrm{keV}}^{1/2}$, with c_f = 1.02 × 10⁻²³ erg s⁻¹ cm⁻³ (see discussion in Sect. 2 of E15), and write

$$\begin{array}{c}\hfill L=c_f\frac{f_L}{{\mu }_e^2{m}_a^2}\frac{M_g^2}{V}{T}_{\mathrm{spec}}^{1/2}=c_f\frac{f_L}{{\mu }_e^2{m}_a^2}{f}_g^2\frac{M_{\mathrm{HE}}^2}{V}{T}_{\mathrm{spec}}^{1/2},\end{array}$$(A.7)

where we have defined $f_L=\int {\rho }_{\text{gas}}^{2}dV/{(\int {\rho }_{\text{gas}}dV)}^{2}V$ as the correction needed to consider the gas mass (M_g = ∫ρ_gasdV) instead of the emission integral ($\int n_e^{2}dV$), with ρ_gas = μ_em_an_e, for scaling purposes. Here the integrals are performed over a spherical volume between 0 and R₅₀₀. Using the above relation between total mass and temperature and total mass and volume, we can do the final step and write

$$\begin{array}{c}\hfill L=c_f\frac{f_L}{{\mu }_e^2{m}_a^2}{f}_g^2{f}_T^{3/2}\mathrm{\Delta }{\rho }_{\mathit{cz}}\frac{k_M{M}_0}{E_z{T}_0^{3/2}}{T}_{\mathrm{spec}}^2,\end{array}$$(A.8)

which can be expressed in the form

$$\begin{array}{c}\hfill E_z^{-1}(L/L_0)=k_L{f}_g^2{f}_T^{3/2}{(T_{\mathrm{spec}}/T_0)}^2\end{array}$$(A.9)

with

$$\begin{array}{c}\hfill k_L=k_M\frac{\overline{f_L}{c}_f{f}_{g,0}^2}{{\mu }_e^2{m}_a^2}\frac{3\mathrm{\Delta }{H}_0^2}{8\pi G}\frac{M_0{T}_0^{1/2}}{L_0}.\end{array}$$(A.10)

Using the median value of f_L, $\overline{f_L}=1.90$ (with values in the range −8%, +2%), f_g, 0 = 0.1, and the adopted values of Δ, T₀, M₀, and L₀, we measure k_L = 0.930.

Appendix B: Derivation of k_SZ

In Eq. (13) we introduce the constant k_SZ in the normalisation of the Y − T relation:

$$\begin{array}{c}\hfill \frac{E_z{Y}_{\mathrm{SZ}}{D}_A^2}{{10}^{-4}{\mathrm{Mpc}}^2}=k_{\mathrm{SZ}}{f}_T^{5/2}{f}_{\mathit{nc}}{\mathcal{T}}^{5/2}.\end{array}$$(B.1)

Here f_T and f_g are the values normalised to the median value we observe in our model $\overline{f_T}=0.71$ (see Fig. A.1) and 0.1, respectively, and f_nc is the clumping-free gas mass fraction that is related to the X-ray-measured quantity f_g through the relation f_g = C^0.5 f_nc, where $C=\langle n_{\text{gas}}^2\rangle /{\langle n_{\text{gas}}\rangle }^2$ is the clumping factor.

To define k_SZ, we proceed from the definition of $Y_{\text{SZ}}{D}_A^2=({\sigma }_T/m_e{c}^2)\int P_{e}dV$, where D_A is the angular diameter distance to the cluster, σ_T = 8π/3(e²/m_ec²)² = 6.65 × 10⁻²⁵ cm² is the Thomson cross section, and ∫P_edV = ∫n_eTdV ≈ f_T T_spec M_gC^−0.5/(μ_em_a) = f_g C^−0.5 f_T T_specM_HE/(μ_em_a).

By substituting M_HE from Eq. (A.4), we can write

$$\begin{array}{c}\hfill Y_{\mathrm{SZ}}{D}_A^2=\frac{{\sigma }_T}{m_e{c}^2}\frac{f_{\mathit{nc}}}{{\mu }_e{m}_a}{f}_T^{5/2}\frac{k_M{M}_0{T}_0}{E_z}{\mathcal{T}}^{5/2},\end{array}$$(B.2)

which allows us to define the constant

$$\begin{array}{cc}\hfill k_{\mathrm{SZ}}& =\frac{{\sigma }_T}{m_e{c}^2}{c}_P\frac{2^{1/2}{f}_{g,0}{\overline{\beta }}^{3/2}{\overline{f_T}}^{5/2}}{G{H}_0{\mathrm{\Delta }}^{1/2}{\mu }_e{\mu }^{3/2}{m}_a^{5/2}}\frac{T_0^{5/2}}{Y_0}=0.296,\hfill \end{array}$$(B.3)

where f_g, 0 = 0.1, $\overline{\beta }=2.94$, $\overline{f_T}=0.71$, T₀ = 5 keV, Y₀ = 10⁻⁴ Mpc², and c_P = ∫P_edV × (μ_em_a)/(M_gf_TT_spec) has a typical value of 1.36 (dimensionless) with variations of ±1% in the range of T and redshift studied here.

It should be noted that using Y_X = M_g × T, which is the X-ray analogue of the integrated Compton parameter introduced by Kravtsov et al. (2006), we can write c_P = Y_SZ/(c_XSZY_Xf_T), where c_XSZ = σ_T/(m_ec² μ_em_a)≈1.40 × 10⁻¹⁹ Mpc² $M_{\odot }^{-1}$ keV⁻¹. Therefore, we predict a ratio C = Y_SZ/(c_XSZY_X) = c_P f_T equal to median values of 0.924 (with D14 model) and 0.927 (B13) with relative changes of ±1% in the redshift and mass ranges investigated in this paper, in close agreement with the present observational constraints (see e.g. Arnaud et al. 2010).

All Tables

All Figures

Fig. 1. Reconstructed thermodynamic radial profiles for objects with (M500, z) = (8 × 1014 M⊙, 0.05), black solid line; the same, but using a pressure profile from X-COP, red solid line, the B13 c − M − z relation, purple solid line; or assuming a hydrostatic bias b = 0.4, blue solid line; (M500, z) = (2 × 1014 M⊙, 0.05), black dashed line; (2 × 1014 M⊙, 1), red dashed line. The ratio is shown with respect to the reference case, black solid line, M500 = 8 × 1014 M⊙ and z = 0.05. P500, n500, T500, K500 refer to the normalisation values presented in Table 1.

In the text

Fig. 2. Comparisons between the stacked profiles from the 120 objects analysed for the ESZ sample in Lovisari et al. (2020) and predictions from different models of (i) the universal pressure profile (Planck, used as reference, and X-COP) and (ii) the c − M − z relation (D14, used as reference, and B13) estimated at the median values in the sample of M500 and redshift. Top panel: electron density profiles recovered from the best-fit parameters of a double-β model. The solid line indicates the median value at each radius; the dotted lines show the 16th and 84th percentile estimated at each radius. The dashed line represents the extrapolation beyond the observational limit for the sample. Bottom panel: stacked temperature profile obtained from the weighted mean of 30 spectral points in each bin. Errors on the mean and dispersion (dotted lines) are overplotted.

In the text

Fig. 3. Top panels: global properties (empty black dots) and after correction for the Ez factor (solid red dots) recovered from their thermodynamic profiles for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), (2 × 1014 M⊙, 1). The dashed magenta lines identify the best-fit relations of Eq. (7). A core-excised bolometric luminosity is considered. Bottom panels: red dots are the quantities (from left to right: fg, fT, b) estimated in our model for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), and (2 × 1014 M⊙, 1); the blue dotted line indicates Ωb/Ωm = 0.157 (Planck Collaboration XXVII 2016); the dashed magenta lines show the predictions following the relations: fg = 0.124 𝒯0.23, fT = 0.687 𝒯0.07 and 1 − b = h0 = 1, where fg and fT are the absolute values, i.e. multiplied by the assumed normalisations of 0.1 and $\overline{f_T}=0.71$, respectively.

In the text

Fig. 4. Impact of the hydrostatic bias b on the integrated quantities, represented as the ratio between the quantities estimated for an assigned bias, and with bias equal to zero. From left to right: this ratio for total mass, temperature, gas mass, core-excised bolometric luminosity, and fT = T(R500)/Tspec. Shown is the case of a system with M500 = 5 × 1014 M⊙ at z = 0.05 and 0.5 (the values at different redshifts overlap for most of the quantities).

In the text

Fig. 5. Distribution of the observed values for the ESZ sample with best-fit results from a linear fit in the logarithmic space (using LIRA; Sereno 2016). Predictions from the semi-analytic model (best-fit values in Table 3) are overplotted with dashed lines in magenta (for the D14 c − M − z relation) and blue (B13). Bottom panels: dashed lines indicate the ratios of the model (M) to the observed (O) best-fit relations.

In the text

Fig. 6. YSZ − T relation for the ESZ sample with the best-fit results from a linear fit in the logarithmic space. The fit on the observed data is performed using LIRA (black dotted line). The predictions from the semi-analytic model overlap (magenta dashed line: D14 model; blue dashed line: B13 model). The dot-dashed orange line is obtained from the propagation of the scaling in the E15 formalism (see Eq. (13)).

In the text

	Fig. 7. Dependence of the normalisation of fg upon the bias b (red line; see Eq. (15)). The shaded region (and the region between dashed lines) represents the quantity fnc/fg, b = 0, where $f_{nc}=f_g/C_0^{0.5}$, with fg ≈ 0.123, fg, b = 0 = 0.124 (0.126) and 1 <  C0 <  1.4 being the 1σ upper limit on the level of gas clumping.
In the text

Fig. A.1. Distribution of the parameters β, fT, fg, and 1 − b as a function of T and colour-coded according to the assumed hydrostatic bias (black: b = 0; green: b = 0.2; red: b = 0.4). The dashed lines in the panel at the bottom right are recovered through the process described in Eq. (8). Points represent quantities derived for the following input values for (M500, z): (8 × 1014 M⊙, 0.05), (5 × 1014 M⊙, 0.05), (2 × 1014 M⊙, 0.05), (8 × 1014 M⊙, 0.5), (5 × 1014 M⊙, 0.5), with and without a bias b; a further case (2 × 1014 M⊙, 1) with b = 0 is considered.

In the text