11institutetext: INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via
Piero Gobetti 93/3, 40129 Bologna, Italy 22institutetext: INFN, Sezione di
Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy 33institutetext: Center
for Astrophysics $|$ Harvard $\&$ Smithsonian, 60 Garden Street, Cambridge, MA
02138, USA 44institutetext: Department of Astronomy, University of Geneva, ch.
d’Ecogia 16, 1290 Versoix, Switzerland

#  i(cm)z, a semi-analytic model for the thermodynamic properties in galaxy
clusters: Calibrations with mass and redshift,
and implication for the hydrostatic bias

S. Ettori 1122 L. Lovisari 1133 D. Eckert 44

###### Abstract

Context. In the self-similar scenario for galaxy cluster formation and
evolution, the thermodynamic properties of the X-ray emitting plasma can be
predicted in their dependencies on the halo mass and redshift only. However,
several departures from this simple self-similar scenario have been observed.

Aims. We show how our semi-analytic model i(cm)z, which modifies the self-
similar predictions through two temperature-dependent quantities, the gas mass
fraction $f_{g}=f_{0}T^{f_{1}}E_{z}^{f_{z}}$ and the temperature variation
$f_{T}=t_{0}T^{t_{1}}E_{z}^{t_{z}}$, can be calibrated to incorporate the mass
and redshift dependencies.

Methods. We used a published set of 17 scaling relations to constrain the
parameters of the model. We were subsequently able to make predictions as to
the slope of any observed scaling relation within a few percent of the central
value and about one $\sigma$ of the nominal error. Contextually, the evolution
of these scaling laws was also determined, with predictions within $1.5\sigma$
and within 10 percent of the observational constraints. Relying on this
calibration, we have also evaluated the consistency of the predictions on the
radial profiles with some observational datasets. For a sample of high-quality
data (X-COP), we were able to constrain a further parameter of the model, the
hydrostatic bias $b$.

Results. By calibrating the model versus a large set of X-ray scaling laws, we
have determined that (i) the slopes of the temperature dependence are
$f_{1}=0.403\;(\pm 0.009)$ and $t_{1}=0.144\;(\pm 0.017)$; and that (ii) the
dependence upon $E_{z}$ are constrained to be $f_{z}=-0.004\;(\pm 0.023)$ and
$t_{z}=0.349\;(\pm 0.059)$. These values, which are inserted in the scaling
laws that propagate the mass and redshift dependence to the integrated
quantities, permit one to estimate directly how the normalizations of a given
quantity $Q_{\Delta}$ changes as a function of the mass (or temperature) and
redshift halo in the form $Q_{\Delta}\sim M^{a_{M}}\,E_{z}^{a_{z}}\sim
T^{a_{T}}\,E_{z}^{a_{Tz}}$, which is in very good agreement with the current
observational constraints. When applied to the best spatially resolved data,
we obtained estimates of the hydrostatic bias $b$ that are lower than, but
still comparable with, the results obtained by other, more standard, means. We
conclude that the calibrated semi-analytic model i(cm)z is able to make
valuable predictions on the slope and redshift evolution of the X-ray scaling
laws, and on the expected radial behavior of the thermodynamic quantities,
including any possible hydrostatic bias.

###### Key Words.:

galaxy: clusters: general – galaxy: clusters: intracluster medium – X-rays:
galaxy: clusters

##  1 Introduction

Galaxy groups and clusters are overdensities in the cosmic field where
gravity, mostly produced from a still-unknown dark component, rules and
regulates most of their appearance. The baryons flow across these halos, and
the largest part of those is heated up during this circulation process
reaching temperatures that make them shine in X-ray and affect the
distribution of the cosmic microwave background (CMB) photons at millimeter
wavelengths. Being gravitationally bound and dark matter dominated, these
massive halos have observed properties that depend, in general and at first
approximation, on the mass and redshift of the halo only (see e.g., Voit,
2005; Allen et al., 2011).

However, several departures from this simple, self-similar scenario have been
observed and are thought to be consequences of the time and mass-dependent
feedback from active galactic nuclei (AGNs) and star formation, with a greater
impact on systems where this feedback contributes in a non-negligible way to
the overall energy budget. This seems to occur more regularly and
systematically in halos at lower masses ($\sim 10^{13}-10^{14}M_{\odot}$; see
e.g., Eckert et al., 2021; Lovisari et al., 2021; Oppenheimer et al., 2021),
where this extra energy becomes comparable to the total binding one. The
action of this feedback over the cosmic life of these structures impacts the
overall baryon distribution in such a way that the ideal condition of an
isolated “closed box,” where the total baryon mass fraction matches the
cosmological value, is reached well beyond any nominal virial radius (at
$>6R_{500}$; see e.g., Angelinelli et al., 2022).

In the following analysis, we refer to masses $M_{\Delta}$ that are associated
with a given overdensity $\Delta$ as $M_{\Delta}=4/3\,\pi\,\Delta\,\rho_{\rm
c,z}R_{\Delta}^{3}$, where $\rho_{\rm 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_{0}\,\left[\Omega_{\Lambda}+\Omega_{\rm
m}(1+z)^{3}\right]^{0.5}=H_{0}\,E_{z}$ is the value of the Hubble constant at
the same redshift. For the $\Lambda$ Cold Dark Matter (CDM) model, we adopt
the cosmological parameters $H_{0}=70$ km s-1 Mpc-1 and $\Omega_{\rm
m}=1-\Omega_{\Lambda}=0.3$.

##  2 The i(cm)z model

The semi-analytic model i(cm)z has been presented in Ettori et al. (2020)
(hereafter E20). It generates radial profiles of the thermodynamic quantities
of galaxy groups and clusters by putting a universal pressure profile in
hydrostatic equilibrium with a gravitational potential described via a
concentration-mass relation. In detail, the procedure combines the following
steps:

  * $\bullet$

as input, a total mass ($M_{500}$ in our analysis) and a redshift $z$ are
provided;

  * $\bullet$

a concentration-mass-redshift relation is assumed,

$\log c_{200}=A+B\log(M_{200}/10^{12}M_{\odot}h_{100}^{-1}),$ (1)

with $B=-0.101+0.026\,z$ and $A=0.520+(0.905-0.520)\exp(-0.617z^{1.21})$
(Dutton & Macciò, 2014); alternative forms can be considered (e.g.,
Bhattacharya et al., 2013; Ludlow et al., 2016); the conversion from $M_{200}$
to $M_{500}$ is obtained through iteration assuming a Navarro, Frenk & White
(Navarro et al., 1997) profile;

  * $\bullet$

a functional form is adopted for the universal electronic pressure profile,

$P_{\rm e}=P_{500}P_{r},$ (2)

with $P_{500}=1.65\times 10^{-3}\left(M_{500}/3\times
10^{14}h_{70}^{-1}M_{\odot}\right)^{2/3+\alpha_{M}}$ keV
cm${}^{-3}\,E_{z}^{8/3}\,h_{70}^{2}$ and
$P_{r}=P_{0}\,(c_{500}x)^{-\gamma}\,\left[1+(c_{500}x)^{\alpha}\right]^{(\beta-\gamma)/\alpha}$
with $x=r/R_{500}$, the parameters $(P_{0},c_{500},\gamma,\alpha,\beta)$ equal
to $(8.403,1.177,0.308,1.051,5.491)$, and $\alpha_{M}=0.12$, accounting for
the observed deviation from the standard self-similar scaling, set as in
Arnaud et al. (2010) (alternative values for the same model or different
models are available in, e.g., Planck Collaboration et al., 2013; Ghirardini
et al., 2019a, b; Sayers et al., 2022);

  * $\bullet$

the hydrostatic equilibrium equation is considered and inverted to recover the
electron density

$n_{\rm e}=-\frac{dP_{\rm e}}{dr}\frac{r^{2}}{\mu\,m_{a}\,G\,M_{\rm HE}},$ (3)

where $\mu\approx 0.6$ is the mean molecular weight of the gas, $m_{a}$ is the
atomic mass unit of $1.66\times 10^{-24}$ g, and $G$ is the universal
gravitational constant;

  * $\bullet$

the gas temperature profile is obtained using the perfect gas law and dividing
the electron pressure in eq. 2 and the electron density in eq. 3

$T_{\rm gas,3D}=P_{\rm e}/n_{\rm e};$ (4)

  * $\bullet$

in the case a hydrostatic bias $b=1-M_{\rm HE}/M$ is considered (see Sect. 6),
the input mass is corrected by a factor $1/(1-b)$ only for eq. 1, leaving the
nominal, uncorrected value $M_{\rm HE}$ in all the remaining cases (e.g., for
the determination of $R_{500}$ or for the normalization of $P_{500}$) to mimic
what observers would do.

From these radial profiles, the integrated quantities are recovered
numerically over the volume of interest (see E20 for further details).

Figure 1: Differences between the predicted and observed values of the slope
(filled circle) and the redshift dependence (cross) of the listed scaling
relations. (Top panel) Deviations in terms of $\sigma=(O-M)/\epsilon_{O}$,
where $O$ is the estimate constrained from the observations with error
$\epsilon_{O}$ and $M$ is the predicted value; (bottom panel) systematic
deviation $\delta=(O-M)/M$. Shaded regions represent $[-3,+3]\,\sigma$ and
$[-0.3,+0.3]\,\delta$ in the top and bottom panel, respectively. Labels
indicate the following: luminosity, either bolometric ($L_{\rm b}$) or in
bands (0.5-2 keV, $L_{05}$; 0.1-2.4 keV, $L_{01}$); total ($M$) and gas
($M_{\rm g}$) mass; and temperature $T$ and X-ray $Y_{x}=M_{\rm g}\,T$.

###  2.1 Departures from self-similarity: The temperature (mass) and redshift
dependence

In gravity-dominated systems, the so-called self-similar relations between the
integrated quantities, such as luminosity, gas mass, and temperature, hold.
However, departures from this self-similar scenario are expected, and they are
observed now systematically, going down in mass, where radiative and feedback
processes impact the intra-cluster medium (ICM) distribution and the action of
the gravity alone, inducing deviations in the gas distribution both as a
function of radius and in the single power-law behavior between integrated
quantities.

As described in our previous work (see e.g., Ettori, 2015; Ettori et al.,
2020), we account for these departures by introducing temperature-
(mass-)dependent quantities such as the gas mass fraction $f_{g}$, and the
variation between the global spectroscopic temperature $T_{\rm spec}$ and its
value at $R_{500}$, $f_{T}=T(R_{500})/T_{500}\times T_{500}/T_{\rm
spec}=T(R_{500})/T_{\rm spec}$. These quantities allow us to express the
corrections needed to accommodate for the observed deviations from the self-
similar scenario, both radially and in the scaling relations in the form of a
power-law dependence on either the gas temperature or the halo mass.

Following E20, we use the following notation hereafter:

$\begin{array}[]{lp{0.8\linewidth}}M&total mass used in
eq.\leavevmode\nobreak\ \ref{eq:cmz}\\\ M_{\rm HE}&hydrostatic mass used in
eq.\leavevmode\nobreak\ \ref{eq:puniv}, \ref{eq:mhe}\\\ b=1-M_{\rm
HE}/M&hydrostatic mass bias\\\ f_{T}=T(R_{500})/T_{\rm
spec}=t_{0}T^{t_{1}}&variation in temperature\\\ f_{g}=M_{g}/M_{\rm
HE}=f_{0}T^{f_{1}}&gas mass fraction dependence on $T$.\\\ \end{array}$

This will allow us to write a modified, more flexible, version of the standard
self-similar scaling laws. We refer readers to E20 for a further detailed
treatment of the normalizations $t_{0}$ and $f_{0}$. For the sake of clarity,
as we also discuss in E20, we note that the X-ray-measured gas mass fraction
$f_{g}=M_{g}/M_{\rm HE}$ is related to the “true” gas mass fraction (i.e., the
one unbiased from the hydrostatic bias and the clumping factor) by the
relation $f_{\rm
gas,true}=M_{g}/M_{tot}=(1-b)^{2/3}\,f_{g}=(1-b)^{2/3}\,C^{0.5}\,f_{nc}$,
where we assume that $M_{g}$ is proportional to $R_{\Delta}$ and
$R_{\Delta,\rm HE}=(1-b)^{1/3}R_{\Delta,tot}$, $f_{nc}$ is the clumping-free
gas mass fraction, and $C=<n_{\rm gas}^{2}>/<n_{\rm gas}>^{2}$ is the clumping
factor.

In E20 (see also Ettori, 2015), we developed all the equations in function of
the gas temperature and quote the rescaling also as a function of the total
mass only for the normalizations of the thermodynamic quantities (see Table 4
in Ettori et al., 2020). For the sake of completeness and because it is useful
for the new calibrations that we discuss in the present work, we have
rewritten those relations, focusing only on their dependency on $T$ (or $M$):

$\displaystyle\mathcal{M}$
$\displaystyle\sim(1-b)^{-1}\,f_{T}^{3/2}\,\mathcal{T}^{3/2}\sim\mathcal{T}^{3/2+3/2t_{1}}$
$\displaystyle\mathcal{M}_{g}$ $\displaystyle\sim
f_{g}\,f_{T}^{3/2}\,\mathcal{T}^{3/2}\sim\mathcal{T}^{3/2+3/2t_{1}+f_{1}}$
$\displaystyle\mathcal{L}$ $\displaystyle\sim
f_{g}^{2}\,f_{T}^{3/2}\,\mathcal{T}^{2}\sim\mathcal{T}^{3/2+c+3/2t_{1}+2f_{1}}$
$\displaystyle\mathcal{T}$ $\displaystyle\sim\mathcal{M}^{2/3\times
1/(1+t_{1})}.$ (5)

Here, we have defined the following quantities: the total mass
$\mathcal{M}\equiv E_{z}M_{tot}/M_{0}$ ($M_{0}=5\times 10^{14}M_{\odot}$), the
gas mass $\mathcal{M}_{g}\equiv E_{z}M_{g}/M_{g,0}$ ($M_{g,0}=5\times
10^{13}M_{\odot}$), a bolometric luminosity $\mathcal{L}\equiv
E_{z}^{-1}L/L_{0}$ ($L_{0}=5\times 10^{44}$ erg s-1), and a gas temperature
$\mathcal{T}\equiv k_{B}T_{\rm spec}/T_{0}$ ($T_{0}=$ 5 keV). The exponent $c$
indicates the temperature dependence of the cooling function and is equal to
$0.44$ (for the pseudo-bolometric 0.01-100 keV band), $-0.11$ (0.1–2.4 keV
band), and $-0.09$ (0.5–2 keV band), with a slight dependence on the
metallicity and on the abundance table for the reference assumed (0.3 times
the solar value in Asplund et al. 2009 in our case; see also Ettori, 2015;
Lovisari et al., 2021, and some further considerations in Sect. 7).

###  2.2 Evolution with redshift

So far, the redshift evolution has been considered only in the part related to
the self-similar model through the quantity $E_{z}$. We have expand this
approach by also allowing a redshift evolution to the $f_{g}$ and $f_{T}$
parameters in the following form: $f_{g}=f_{0}T^{f_{1}}E_{z}^{f_{z}}$ and
$f_{T}=t_{0}T^{t_{1}}E_{z}^{t_{z}}$, with two extra parameters $f_{z}$ and
$t_{z}$.

By propagating the dependence through the scaling laws, we obtain

$\displaystyle M$
$\displaystyle\sim\mathcal{T}^{3/2+3/2t_{1}}\;E_{z}^{3/2t_{z}-1}$
$\displaystyle M_{g}$
$\displaystyle\sim\mathcal{T}^{3/2+3/2t_{1}+f_{1}}\;E_{z}^{f_{z}-3/2t_{z}-1}$
$\displaystyle L$
$\displaystyle\sim\mathcal{T}^{3/2+c+3/2t_{1}+2f_{1}}\;E_{z}^{2f_{z}+3/2t_{z}+1}$
$\displaystyle T$ $\displaystyle\sim M^{2/3\times
1/(1+t_{1})}\;E_{z}^{(3/2-t_{z})/(1+t_{1})},$ (6)

where we indicate, on the right side, the entire contribution assigned to the
evolution in the form of $E_{z}$; that is to say, the quantities on the left
side are the ones actually measured and are the same as the ones in eq. 5 once
the dependence on $E_{z}$ is moved to the right side of the equation (e.g.,
$M=\mathcal{M}E_{z}^{-1}$, $L=\mathcal{L}E_{z}$).

##  3 Calibration of the model parameters with observed scaling laws

The i(cm)z model assumes some power-law relations among the integrated
quantities. These relations are presented in eq. 5, and their redshift
evolution in eq. 6. Any new relation can be obtained by propagation of those,
allowing us to predict the expected values of the slope (and redshift
evolution) for any scaling relation as a function of only two parameters,
$f_{1}$ and $t_{1}$, for the slope, and two others for the $z$ evolution,
$f_{z}$ and $t_{z}$.

To calibrate these values, we have collected from recent literature the best-
fit results for the slopes and redshift evolution of several scaling
relations. To make full use of the entire set of considered relations and of
the statistical uncertainties on their best-fit constraints, we proceeded as
follows:

  1. 1.

we randomly drew, from a flat distribution between 0 and 1, two numbers
representing the exploring set $i$ of values $\\{f_{1i},t_{1i}\\}$;

  2. 2.

for any scaling relation $j$ considered, we assigned a probability
$P_{ij}=\exp{(-(O_{j}-M_{j})^{2}/\epsilon_{O_{j}}^{2}}$ to each set $i$ of
values, where $O_{j}$ is the observational constraint on the slope with an
associated error $\epsilon_{O_{j}}$ and $M_{j}$ is the predicted value
provided from the set $\\{f_{1i},t_{1i}\\}$;

  3. 3.

for each set $i$ of values, we built a cumulative probability
$P_{i}=\prod_{j}P_{ij}$, where the product was done over the entire collection
of scaling relations $j$ under consideration that are considered independent
of each other here;

  4. 4.

we repeated the procedure $N$ (with $N\sim 10^{6}$) times and estimate the
mean and the dispersion of the values as $\bar{f_{1}}=\sum_{i}w_{i}f_{1i}$ and
$\sigma_{f_{1}}=\left(\sum w_{i}(f_{1i}-\bar{f_{i}})^{2}\right)^{1/2}$, with
the weights $w_{i}=P_{i}/\sum P_{i}$. The same calculation was performed to
estimate $\bar{t_{1}}$ and $\sigma_{t_{1}}$.

We selected any scaling law where the normalization, slope, and redshift
evolution were considered as free parameters. If relations for more subsamples
are available, we would use the one with the largest coverage in mass and
redshift. We used a total number of 17 scaling laws extracted from three
different samples:

Lovisari et al. (2020)

discuss the integrated properties of 120 galaxy clusters in the Planck Early
Sunyaev-Zel’dovich (ESZ) sample spanning the mass range $2.4-17.6\times
10^{14}M_{\odot}$ in the redshift interval $0.059<z<0.546$;

Bahar et al. (2022)

present the analysis of scaling relations for a sample of 265 clusters
extracted from the eROSITA Final Equatorial-Depth Survey (eFEDS) field with a
contamination level of $<10$% (including AGNs and spurious fluctuations) and
covering ranges in total mass of $6.9\times 10^{12}M_{\odot}<M_{500}<7.8\times
10^{14}M_{\odot}$ (with 68 systems at masses below $10^{14}M_{\odot}$
representing one of the largest group samples detected uniformly to date) and
in redshift of $0.02<z<0.94$;

Pratt et al. (2022)

studied the scaling laws of 93 SZ-selected objects with XMM-Newton exposures
spanning a mass range of $M_{500}=0.5-20\times 10^{14}M_{\odot}$ and lying at
redshifts $0.05<z<1.13$.

We noticed that no specific information on the dynamical state of the objects
of these samples is available. However, given the SZ-selected samples of
Lovisari et al. (2020) and Pratt et al. (2022), we expect a contribution from
cool-core (more relaxed) objects to be lower than in samples selected through
their X-ray emission (in particular, when the core contribution is not
compensated for; see e.g., Eckert et al., 2011), and to be around $\sim$
30-50% of the total (see e.g., Rossetti et al., 2016; Nurgaliev et al., 2017;
Andrade-Santos et al., 2017). For what concerns the sample in Bahar et al.
(2022), the investigation of the morphological properties in a more extended
eFEDS sample of 325 objects presented in Ghirardini et al. (2022) indicates
that about 30-40% of the clusters are dynamically relaxed. For the sake of
simplicity, but not affecting the overall efficiency in reproducing the
observed datasets, we have not considered any dependency on the dynamical
state in the current parametrization of the i(cm)z model (see more comments in
Sect. 7).

The final outcomes of this process are the following constraints on $f_{1}$
and $t_{1}$:

$\displaystyle\bar{f_{1}}$ $\displaystyle=$ $\displaystyle
0.403\;(\sigma_{f_{1}}=0.009)$ $\displaystyle\bar{t_{1}}$ $\displaystyle=$
$\displaystyle 0.144\;(\sigma_{t_{1}}=0.017),$ (7)

with an associated mild correlation coefficient $\rho={\rm
cov}(f_{1},t_{1})/(\sigma_{f_{1}}\sigma_{t_{1}})$ of 0.44.

When the analysis was repeated for each of the three samples independently, we
measured the following: $f_{1}=0.45\;(0.04)$ and $t_{1}=0.05\;(0.03)$, only
using the seven relations from Lovisari et al. (2020); $f_{1}=0.47\;(0.03)$
and $t_{1}=0.25\;(0.07)$, only using the three relations from Pratt et al.
(2022); and $f_{1}=0.22\;(0.02)$ and $t_{1}=0.40\;(0.04)$, with the seven
scaling laws from Bahar et al. (2022). The best-fit results presented in
equation 7 represent a sort of weighted mean of these values. We noticed that
SZ-selected samples (such as the ones in Lovisari et al. 2020 and Pratt et al.
2022) seem to prefer $f_{1}>t_{1}$, whereas the reverse can be observed in the
X-ray-selected sample of Bahar et al. (2022). With our goal being to provide a
first, exhaustive calibration of the physically meaningful parameters of the
i(cm)z model $f_{g}$ and $f_{T}$, we have made use of the entire collection of
scaling relations, assigning them equal weight (and with a statistical weight
defined by their corresponding error) once combined to obtain the best-fit
results presented in equation 7.

Figure 2: Predicted dependence on mass and redshift of the quantities $f_{g}$
(blue contours) and $f_{T}$ (red contours) renormalized to 1 at
$(M_{500},z)=(5\times 10^{14}M_{\odot},0.05)$ (purple diamond).

This procedure has also been validated by comparing the derived and the
observed constraints of the slopes of the scaling laws considered. In Fig. 1,
we plot the deviations between the observed constraints and the best-fit
results of eq. 7. We show both the statistical deviations, as the difference
between the central values in terms of the quoted error,
$\sigma=(O-M)/\epsilon_{O}$, and the systematic one as the difference between
the central values divided by the model prediction, $\delta=(O-M)/M$. We
measured a median (and first and third quartiles) of
$\sigma=-0.01\,(-1.59,1.02)$ and of $\delta=0.0\,(-0.05,0.03)$ on the
measurements of the slope.

Using the observed limits on the evolution of the scaling relations, together
with the constraints obtained above on $f_{1}$ and $t_{1}$, we put constraints
on the redshift dependency of $f_{g}$ and $f_{T}$ (see Sect. 2.2):

$\displaystyle\bar{f_{z}}$ $\displaystyle=$
$\displaystyle-0.004\;(\sigma_{f_{z}}=0.023)$ $\displaystyle\bar{t_{z}}$
$\displaystyle=$ $\displaystyle 0.349\;(\sigma_{t_{z}}=0.059),$ (8)

with a correlation coefficient of 0.63. Adding this contribution to the
evolution improved the modeling of the observed constraints. We measured
$\sigma=0.16\,(-1.79,0.67)$ and $\delta=-0.08\,(-0.43,0.17)$ when these extra
parameters were not included and the entire evolution is described by the
modified self-similar model. When we also propagated $f_{z}$ and $t_{z}$, we
obtained $\sigma=-0.15\,(-0.92,1.39)$ and $\delta=-0.03\,(-0.30,0.07)$.

The best-fit results presented in equations 7 and 8 provide a full description
of the behavior of the physically meaningful parameters $f_{g}$ and $f_{T}$ of
the i(cm)z model. As we show in Fig. 2, we can recover their expected
dependencies on the halo mass and redshift, suggesting that the gas fraction
$f_{g}$ does not change with cosmic time, whereas a mild positive evolution is
predicted for $f_{T}$, and that $f_{g}$ decreases significantly with $M_{500}$
(by 50% going down in mass by an order of magnitude), and more severely than
$f_{T}$.

Table 1: Dependences of the characteristic physical scales on the temperature,
mass, and redshift (as $E_{z}^{\alpha_{z}}$) in the self-similar and i(cm)z
models.

Quantity | $f(M,z)$ | $f(f_{\rm gas})$ | $f(T,z)$ | $f(M,z)$
---|---|---|---|---
| self-similar | i(cm)z
$T_{\Delta}$ | $M^{2/3}\,E_{z}^{2/3}$ | $f_{\rm gas}^{0}$ | $T^{1+t_{1}}\,E_{z}^{t_{z}}$ | $M^{2/3}\,E_{z}^{2/3}$
$n_{\Delta}$ | $E_{z}^{2}$ | $f_{\rm gas}^{1}$ | $T^{f_{1}}\,E_{z}^{2+f_{z}}$ | $M^{2/3\,f_{1}/(1+t_{1})}\,E_{z}^{2+f_{z}+f_{1}(2/3-t_{z})/(1+t_{1})}$
$P_{\Delta}$ | $M^{2/3}\,E_{z}^{8/3}$ | $f_{\rm gas}^{1}$ | $T^{1+t_{1}+f_{1}}\,E_{z}^{2+f_{z}+t_{z}}$ | $M^{2/3\,+2/3f_{1}/(1+t_{1})}\,E_{z}^{8/3+f_{z}+f_{1}(2/3-t_{z})/(1+t_{1})}$
$K_{\Delta}$ | $M^{2/3}\,E_{z}^{-2/3}$ | $f_{\rm gas}^{-2/3}$ | $T^{1+t_{1}-2/3f_{1}}\,E_{z}^{-4/3-2/3f_{z}+t_{z}}$ | $M^{2/3\,-4/9f_{1}/(1+t_{1})}\,E_{z}^{-2/3-2/3f_{z}-2/3f_{1}(2/3-t_{z})/(1+t_{1})}$
111The basic equations are $T_{\Delta}=f_{T}T\sim
T^{1+t_{1}}\,E_{z}^{t_{z}}=(E_{z}M)^{2/3}=(E_{z}^{3}R^{3})^{2/3}$ and
$n_{\Delta}\sim\Delta\rho_{cz}\sim f_{g}E_{z}^{2}\sim
T^{f_{1}}E_{z}^{2+f_{z}}$. All the other relations were obtained by
combinations of those.

Table 2: Best-fit results for
$Q_{\Delta}\sim\,M^{a_{M}}\,E_{z}^{a_{z}}\sim\,T^{a_{T}}\,E_{z}^{a_{Tz}}$.

$Q_{\Delta}$ | $a_{M}$ | $a_{z}$ | $a_{T}$ | $a_{T,z}$
---|---|---|---|---
$T_{\Delta}$ |  $2/3$ $[2/3]$ |  $2/3$ $[2/3]$ |  $1.14\;(0.02)$ $[1]$ |  $0.35\;(0.06)$ $[0]$
$n_{\Delta}$ |  $0.23\;(0.01)$ $[0]$ |  $2.11\;(0.03)$ $[2]$ |  $0.40\;(0.01)$ $[0]$ |  $2.00\;(0.02)$ $[2]$
$P_{\Delta}$ |  $0.90\;(0.01)$ $[2/3]$ |  $2.78\;(0.03)$ $[8/3]$ |  $1.55\;(0.02)$ $[1]$ |  $2.35\;(0.06)$ $[2]$
$K_{\Delta}$ |  $0.51\;(0.01)$ $[2/3]$ |  $-0.74\;(0.02)$ $[-2/3]$ |  $0.88\;(0.02)$ $[1]$ |  $-0.98\;(0.06)$ $[-4/3]$
222Mean and dispersion for the parameters $a_{M}$, $a_{z}$, $a_{T}$, and
$a_{Tz}$ obtained from the constraints on $f_{g}\sim T^{f_{1}}E_{z}^{f_{z}}$
and $f_{T}\sim T^{t_{1}}E_{z}^{t_{z}}$ in equations 7 and 8. In the square
brackets, we quote the self-similar predictions. We note that, in combination
with the relations in Tab. 1, other rescaling can be obtained. For example,
the emission measure $EM\approx\int n^{2}dl$, which can be directly
reconstructed from the observed X-ray surface brightness, is expected to scale
as
$M^{1/3+4/3f_{1}/(1+t_{1})}E_{z}^{10/3+2f_{z}+2f_{1}(2/3-t_{z})/(1+t_{1})}\sim
M^{0.80(0.01)\,[1/3]}E_{z}^{3.55(0.06)\,[10/3]}$ and as
$T^{1/2+2f_{1}+t_{1}/2}E_{z}^{3+2f_{z}+t_{z}/2}\sim
T^{1.38(0.02)\,[1/2]}E_{z}^{3.17(0.05)\,[3]}$.

Figure 3: Radial profiles of the relative uncertainties ($\epsilon/Y-1$) on
(left) gas density, (center) 3D gas temperature, and (right) projected
spectroscopic-like temperature due to (i) the errors on the assumed parameters
that describe the mass dependence, (ii) the scatter in the $c-M-z$ relation,
and (iii) the intrinsic scatter on the “universal” pressure profile. Red
(green) dotted lines indicate the median values of the 16th and 84th
percentiles (and of the medians) estimated at each radius. The case refers to
the input values $(M_{500},z)=(4.6\times 10^{14}M_{\odot},0.6)$, which
represents the median values of the S18 sample.

Figure 4: Differences between the predicted and observed gas density and
temperature profiles using the datasets from S18 (four panels on the left) and
G19 (four panels on the right). (Top panels) Deviations in terms of
$\sigma=(O-M)/\epsilon$, where $O$ is the observed profile with error
$\epsilon_{O}$ and $M$ is the predicted value with an estimated uncertainty of
$\epsilon_{M}$ (see Fig. 3) that combines in quadrature with $\epsilon_{O}$ to
produce the total error $\epsilon$; and (bottom panels) systematic deviation
$\delta=(O-M)/M$. Grey lines represent the single object; solid, dashed, and
dotted lines indicate mean, median, and scatter, respectively. In the four
panels on the left, blue and red lines correspond to the properties of the
subsamples at $z<0.6$ and $z>0.6$, respectively.

##  4 Radial profiles: Calibration in mass and redshift

Characteristic thermodynamic quantities defined within an overdensity
$\Delta$, $Q_{\Delta}$, were used to renormalize the radial profiles of the
given quantity only as a function of the halo properties (i.e., mass and
redshift; see Tab. 1). Effects of the deviation from self-similarity propagate
to these values in a way that can be represented in the form

$Q_{\Delta}\sim\,M^{a_{M}}\,E_{z}^{a_{z}}\sim\,T^{a_{T}}\,E_{z}^{a_{Tz}}.$ (9)

Self-similar and modified values, after the calibration in the context of the
i(cm)z model presented in Eq. 7 and 8, are quoted in Table 2.

Recently, Pratt et al. (2022) have estimated the evolution and mass dependence
of the “universal” gas density profile for a sample of 93 SZ-selected galaxy
clusters in the mass range $M_{500}=[0.5-20]\times 10^{14}M_{\odot}$ at
$0.05<z<1.13$. They measured $\bar{a}_{M,n}=0.22\pm 0.01$ and
$\bar{a}_{z,n}=2.09\pm 0.02$. Using the best-fit values in equations 7 and 8,
we obtained with the i(cm)z model $a_{M,n}=0.23(\pm 0.01)$ and
$a_{z,n}=2.11(\pm 0.03)$, which are in remarkable agreement with both the
dependencies.

A constraint of the universal pressure profile has been presented in Arnaud et
al. (2010): $\bar{a}_{M,P}=2/3+a_{P}+a_{P1}\approx 0.88$ with $a_{P}=0.12$ and
$a_{P1}=0.10-(a_{P}+0.10)\times(x/0.5)^{3}/(1+(x/0.5)^{3})$ and $x\approx
0.1$, where a self-similar evolution of $\bar{a}_{z,P}=8/3=2.67$ has been
assumed. We obtained $a_{M,P}=0.90(\pm 0.01)$ and $a_{z,P}=2.77(\pm 0.03)$,
which does indeed suggest a deviation in the redshift dependence from the
self-similar model as well. These predicted values of $a_{M,P}$ and $a_{z,P}$
are consistent with the looser constraints we obtained by fitting the gas
density and temperature profiles directly from Sanders et al. (2018) and
Ghirardini et al. (2019a) (see Fig. 9).

Pratt et al. (2010) constrained the mass dependence of the entropy profiles at
a given overdensity. For $\Delta=500$, they measured $\bar{a}_{M,K}=0.62\pm
0.17$, still assuming a self-similar behavior in redshift
($\bar{a}_{z,K}=-2/3=-0.67$). With the calibrated i(cm)z model, we predict
$a_{M,K}=0.51(\pm 0.01)$ and $a_{z,K}=-0.74(\pm 0.02)$.

The predicted and observed profiles were compared and the free parameters
($a_{m},a_{z}$) were constrained through a minimization of the following
$\chi^{2}_{m}$:

$\chi^{2}_{m}=\sum_{i}\left[\frac{\left(d_{i}-m_{i}\right)^{2}}{\epsilon_{d_{i}}^{2}+\sigma_{m_{i}}^{2}}+\ln{\left(\epsilon_{d_{i}}^{2}+\sigma_{m_{i}}^{2}\right)}\right],$
(10)

where the index $m$ refers to the observable under investigation (e.g.,
surface brightness, $k=x$; gas density, $k=n$; and gas temperature, $k=T$);
$m_{i}$ are the predicted profiles for the given observable from i(cm)z
estimated for each object $i$ for the assumed mass and redshift; and $d_{i}$
and $\epsilon_{d_{i}}$ are the values, and the relative error, of the
investigated quantities. The quantity $\sigma_{m_{i}}$ accounts for the
systematic uncertainties affecting the adopted model $m_{i}$. We have
estimated those through a Monte-Carlo process where we propagated the errors
on the assumed mass dependence and the scatter on both the $c-M-z$ relation
and on the universal pressure profiles. For the latter ones, we adopted values
of 0.16 in $\log_{10}c$ for a given mass (e.g., Diemer & Kravtsov, 2015), and
0.10 in $\log_{10}P$ (e.g., Arnaud et al., 2010). We show in Fig. 3 an example
of the estimates of $\sigma_{m_{i}}$ we recovered. For the sake of simplicity,
we propagated the median relative errors to eq. 10.

Depending on the data available, more methods can be combined by adding the
corresponding $\chi_{m}^{2}$. For instance, when the profiles of the gas
density and spectral temperature are in use, the total $\chi^{2}$ is estimated
as $\chi_{n}^{2}+\chi_{T}^{2}$.

We compared the predictions from the model i(cm)z with the radial profiles
available in literature for a few recent datasets:

Sanders et al. (2018) (S18)

The sample of Sanders et al. (2018) includes 83 galaxy clusters detected at
the South Pole Telescope for their SZ signal and studied in X-ray with Chandra
exposures. They span a range in $M_{500}$ between $1.2$ and $17.6\times
10^{14}M_{\odot}$ at redshift $0.28-1.22$, with a median value of $4.6\times
10^{14}M_{\odot}$ and $z=0.6$. The authors recovered the thermodynamic
profiles through a publicly available forward-modeling projection code,
MBPROJ2, that combines the information from exposure-corrected and background-
subtracted X-ray surface brightness profiles extracted in ten independent
energy bands between 0.5 and 7 keV, assuming hydrostatic equilibrium. We used
the publicly available density and temperature profiles.

Ghirardini et al. (2019a) (G19)

This work describes the universal profiles recovered for the 12 massive
($M_{500}>3\times 10^{14}M_{\odot}$; median value, $5.7\times
10^{14}M_{\odot}$), local ($0.05<<0.1$; median value, $0.065$) clusters
constituting the XMM Cluster Outskirts Project (X-COP) sample (Eckert et al.,
2017), that is a very large program (VLP) on XMM-Newton that targeted 12 the
most significant Sunyaev-Zel’dovich (SZ) sources in the Planck survey in order
to combine X-ray and SZ information out to the virial radius.

In Fig. 4, we present the reconstructed profiles for the S18 and G19 samples.
The deviations in $\sigma$ are, on average, within 1 over the radial range
$0-1\;R_{500}$ both in gas density and in temperature, with no evident
dependency on the redshift. For the gas density, a systematic mean deviation
$\delta$ from –30% to +30% moving outward in S18, and of about –10% in G19 is
observed; for the temperature profile, $\sim$ +20/30% ($<$10%) was measured in
S18 (G19), with the scatter around these central values, however, encompassing
the zero corresponding to the perfect match.

##  5 Integrated quantities: Predictions and comparison with the modified
scaling laws

As detailed in E20, the i(cm)z model provides an estimate of the global
quantities by integrating the thermodynamic profiles over the volume of
interest. In Fig. 5, we compare these estimates for a set of representative
quantities (gas mass fraction, gas temperature, and luminosity) with the
relations that should stand among them (see Sect. 2 and equations 6) once the
revised laws have been calibrated as described in Sect. 3.

Figure 5: Comparison between some interesting global values $X$ (hydrostatic
mass, spectroscopic-like temperature, bolometric luminosity, and gas mass
fraction) estimated from the i(cm)z model and the predicted revised scaling
laws (see equations 6). For each quantity $X$, we plot the value normalized at
$(M_{500},z)=(5\times 10^{14}M_{\odot},0.05)$ and estimated at $\Delta=500$
estimated for the grids of $M_{500}/10^{14}M_{\odot}=[0.15,0.5,1,2,5,8]$ and
redshift $[0.05,0.2,0.6,1,1.5]$. Filled dots indicate the estimates at the
extreme of the redshift distribution, with the smallest dots corresponding to
the highest $z$. The crosses show the predicted values from the revised
scaling laws (orange, $f_{g}$; purple, $T_{\rm sl}$; and green, $L_{\rm
bol}$).

In general, the agreement is remarkable; for instance, the gas mass fraction
was recovered within 1% over the entire range of mass
$(10^{13}-10^{15}M_{\odot})$ and redshift $(0-1.5)$ investigated here. In
addition, the spectroscopic temperature matches the predictions from the
modified scaling laws, for a given mass and redshift, within 1% (median value;
with the range of the first and third quartile being –7% and 9%). Larger
tension between the estimated and predicted values is present on the
(bolometric) luminosity: the median (mean) differences are of 8% (26%), with
the largest deviations (¿50%) associated with the haloes at $M_{500}<5\times
10^{13}M_{\odot}$ (at higher masses, the median deviation is –4% with first
and third quartile of –13% and 9%). We obtained a better agreement by
considering the X-ray luminosity in a soft band (e.g., 0.5–2 keV), with a the
median deviation of about 5% (first and third quartile of 0 and 15%).

##  6 Constraints on the hydrostatic bias

The procedure described in Sect. 2 allows the hydrostatic and the total mass
to be treated separately, with the consequent possibility of being to
constrain the hydrostatic bias $b=1-M_{\rm HE}/M$. Overall, the method is
sensitive to two effects: (i) the rescaling of the radius with $R_{500}$, with
the latter being an estimate derived from the hydrostatic mass $M_{\rm
HE}\equiv M_{500}$ defined as input; and (ii) the actual depth of the
potential well, and thus the measurement of the gas temperature. We show in
Fig. 6 the predicted impact of different levels of $b$: with no strong
dependence on the halo mass, the gas density should vary by about 10-15% in
the radial range $0.1-1R_{500}$ for $b=0.2$, and between 30% and 50% for
$b=0.6$; the gas temperature should be higher by $\sim$ 10% when $b=0.2$ and
up to 90% with respect to the expected value evaluated with no bias when
$b=0.6$. To constrain the hydrostatic bias, the technique looks for a minimum
$\chi^{2}$ for a grid of values of $b$ adopted as input, together with the
other parameters of the i(cm)z model, the hydrostatic mass $M_{500}$, and the
redshift.

Figure 6: Predicted electron density (left) and temperature (right) profiles
of two halos at $z=0.05$ and a hydrostatic mass $M_{500}$ of $5\times
10^{13}M_{\odot}$ (blue lines) and $5\times 10^{14}M_{\odot}$ (red lines),
respectively. Different line styles indicate the level of hydrostatic bias:
$b=0$ (solid lines); $b=0.2$ (dotted line); and $b=0.6$ (dashed line). Bottom
panels: differences $\delta=(M_{b}-M_{0})/M_{0}$ between the profiles with
respect to the case with $b=0$.

We have applied this technique to the X-COP sample, for which independent
measurements of $b$ have been obtained either by comparing different mass
estimators (Ettori et al., 2019), or by imposing a universal gas mass fraction
and various models of nonthermal pressure support (Eckert et al., 2019; Ettori
& Eckert, 2022). We show the best-fit constraints for all 12 objects in Fig.
13. Table 3 and Fig. 7 summarize the limits from the published and current
work. The technique proposed in this study does indeed seem to identify the
objects more affected from a hydrostatic bias; although, they do suggest
values on the lower end of the distribution overall.

Figure 7: Predicted $1\sigma$ upper limits on the value of the hydrostatic bias $b$ here represented as the ratio between the nonthermal ($P_{\rm NT}$) and total ($P_{\rm tot}$) pressure in the ICM ($b=P_{\rm NT}/P_{\rm tot}$; see e.g., Eckert et al., 2019; Ettori & Eckert, 2022). The dashed line represents the upper limit when only the temperature profile was used in the likelihood; when the new universal pressure profile also formed X-COP data (see Sect. B), the upper limits indicated with a dotted line were obtained. We note that the high level of nonthermal pressure support required for A2319 is explained by the ongoing merger described in Ghirardini et al. (2018).  Table 3: Hydrostatic bias $b(R_{500})$ in X-COP clusters. Cluster |  $\alpha_{500}$ (Eck19) |  $\alpha_{500}$ (Ett22) |  $b$ (this work) / only $\chi^{2}_{T}$
---|---|---|---
A1644 | $<0.10$ | $<0.00$ |  $0.05_{-0.04}^{+0.04}$ / $0.09_{-0.05}^{+0.04}$
A1795 | $0.02_{-0.02}^{+0.06}$ | $0.07_{-0.02}^{+0.04}$ |  $<0.01$ / $<0.01$
A2029 | $0.06_{-0.06}^{+0.06}$ | $0.11_{-0.04}^{+0.05}$ |  $0.06_{-0.05}^{+0.05}$ / $0.03_{-0.03}^{+0.06}$
A2142 | $0.16_{-0.05}^{+0.05}$ | $0.20_{-0.02}^{+0.03}$ |  $0.09_{-0.03}^{+0.04}$ / $0.11_{-0.05}^{+0.04}$
A2255 | $0.06_{-0.06}^{+0.07}$ | $0.11_{-0.02}^{+0.08}$ |  $<0.05$ / $0.18_{-0.07}^{+0.06}$
A2319 | $0.44_{-0.04}^{+0.04}$ | $0.54_{-0.07}^{+0.12}$ |  $0.12_{-0.02}^{+0.02}$ / $0.19_{-0.02}^{+0.02}$
A3158 | $0.09_{-0.06}^{+0.06}$ | $0.13_{-0.06}^{+0.06}$ |  $0.04_{-0.04}^{+0.04}$ / $0.15_{-0.06}^{+0.05}$
A3266 | $<0.11$ | $<0.00$ |  $<0.02$ / $<0.01$
A644 | $0.03_{-0.03}^{+0.06}$ | $<0.15$ |  $<0.01$ / $<0.01$
A85 | $0.10_{-0.06}^{+0.05}$ | $0.15_{-0.03}^{+0.03}$ |  $0.01_{-0.01}^{+0.04}$ / $0.04_{-0.04}^{+0.04}$
RXC1825 | $0.05_{-0.05}^{+0.05}$ | $<0.29$ |  $0.17_{-0.03}^{+0.03}$ / $0.22_{-0.03}^{+0.04}$
ZW1215 | $<0.12$ | $<0.00$ |  $0.11_{-0.06}^{+0.06}$ / $<0.01$

We further considered the role of the universal pressure profile, in
particular, by estimating a profile through a joint fit of the X-COP
nonparametric deprojected points from the recent analysis in Eckert et al.
(2022) and including the corrections due to the nonthermal pressure to
$M_{500}$ (and, consequently, to $R_{500}$) from the analysis in Eckert et al.
(2019) (see Sect. B for details). The purpose of this exercise was to evaluate
how the universal pressure profile changes once the nonthermal pressure
contribution is included and propagated to the rescaling in the radius and
normalization. The new radial profile is latter than the one adopted here from
Arnaud et al. (2010) (see Fig. 12), implying, for a given mass, lower (higher)
values of gas density (temperature). The overall net effect is the reduction
of the estimated mass bias.

A larger effort as to the calibration of the universal pressure profile, also
relying on careful and possibly bias-corrected mass estimates, is needed to
strengthen the constraints through this method. Reversely, the method outlined
here can potentially be used to assess the robustness of the assumption on the
universal gas fraction that is the imposed to constrain the level of
nonthermal pressure in Eckert et al. (2019) and Ettori & Eckert (2022).

##  7 Limitations of the i(cm)z model

The i(cm)z model, introduced in E20 and extended in the present work, relies
on the following: (i) some relations between the total mass and halo
concentration and redshift as well as a universal pressure profile put in
hydrostatic equilibrium, and (ii) some simple modifications of the scaling
relations that account from the departure from the self-similar predictions
(see Equations 5 and 6) also for what concerns the normalizations of the
universal profiles (see Tab. 1). In E20, we discuss what the impact is of
assuming alternative $c-M-z$ relations and some different sets of parameters
describing the universal pressure profile. The i(cm)z model has been
calibrated here with some relations among X-ray-integrated quantities
published in recent work. As we discuss in Sect. 3, each analyzed sample
provides constraints on the parameters of the model (primarily $f_{1}$ and
$t_{1}$, and then $f_{z}$ and $t_{z}$) that differ between them and that might
depend on the relations used and on the characteristic of the sample selected.
With the purpose of providing a first, general calibration of our model, we
relied on all the scaling laws published recently for samples of hundreds of
objects over a large range of masses and redshifts.

Some quantities are more prone than others to the few assumptions made in
reconstructing them. For instance, the X-ray luminosity depends on the assumed
metallicity (0.3 times the solar value as tabulated in Asplund et al., 2009).
We have also investigated what the impact is in assuming a different value for
the mean metallicity or a different table of reference. Using a metallicity of
1 solar, the tension between the predicted and recovered luminosity presented
in Fig. 5 increases by more than a factor of 4. On the contrary, still
assuming a value of 0.3, but with respect to the abundance table in Anders &
Grevesse (1989), the median differences increase to about 15% for the X-ray
luminosity in 0.5–2 keV (first and third quartile of 8 and 27%) and to 17% for
the pseudo-bolometric one (0 and 66%). This assumption has a more subtle
impact on the exponent of the cooling function (see also e.g., Lovisari et
al., 2021). While the assumption on the solar table does not significantly
affect (at fixed metallicity) the slope ($c\approx 0.43-0.44$ in the pseudo-
bolometric band 0.01-100 keV; $c\approx-0.1$ in the soft X-ray bands); changes
in metallicity have a larger impact, with $c$ decreasing to about $0.37$ and
$-0.18$ for the pseudo-bolometric and soft bands, respectively, for a metal
abundance of 1 solar, which induce a small change (of about 10%) on the
estimate of $f_{1}$.

Furthermore, we decided to calibrate the only few parameters that keep the
model simple and, at the same time and more importantly, effective. A possible
extension of the model would be the inclusion of a parameter representing the
dynamical state (i.e., relaxed and disturbed, see e.g., Campitiello et al.,
2022), which, for instance, might affect the choice of the $c-M-z$ relation
and of the universal pressure profile, both showing slight differences between
the most relaxed systems and the majority of the population.

##  8 Summary and conclusions

We present a detailed calibration in temperature (mass) and redshift of the
i(cm)z model we have introduced in Ettori et al. (2020). The calibration was
done using both a large sample of recently evaluated scaling relations (see
Lovisari et al., 2020; Bahar et al., 2022; Pratt et al., 2022) and spatially
resolved thermodynamic quantities (from Sanders et al., 2018; Ghirardini et
al., 2019a) at halo masses down to $M_{500}\approx 10^{13}M_{\odot}$ and up to
redshift 1.2.

Our main findings are as follows.

  1. (i)

We effectively modified the self-similar scenario by introducing temperature-
(mass-)dependent quantities such as the gas mass fraction $f_{g}$, the
variation between its value at $R_{500}$, and the global spectroscopic
temperature $T_{\rm spec}$, $f_{T}=T(R_{500})/T_{500}\times T_{500}/T_{\rm
spec}=T(R_{500})/T_{\rm spec}$; by expressing the dependencies on the gas
temperature and redshift as $f_{g}\sim T^{f_{1}}E_{z}^{f_{z}}$ and $f_{T}\sim
T^{t_{1}}E_{z}^{t_{z}}$, we recovered how the slope and redshift evolution of
the (mostly X-ray) scaling relations have to be modified (see Table 1, and
Sections 2.1 and 2.2).

  2. (ii)

Using a large dataset of published slopes and redshift evolution of X-ray
scaling laws, in Sect. 3, we have constrained $\bar{f_{1}}=0.403\;(\pm
0.009)$, $\bar{t_{1}}=0.144\;(\pm 0.017)$, $\bar{f_{z}}=-0.004\;(\pm 0.023)$,
and $\bar{t_{z}}=0.349\;(\pm 0.059)$. These constraints are represented as
predictions of the values of $f_{g}$ and $f_{T}$ as a function of the halo
mass and redshift in Fig. 2. We expect that the gas fraction $f_{g}$ does not
change with the cosmic time, whereas a mild positive evolution is predicted
for $f_{T}$; on the contrary, $f_{g}$ decreases significantly with $M_{500}$
(by 50% going down in mass by an order of magnitude), and definitely more
drastically than $f_{T}$.

  3. (iii)

We verified that these values are well in agreement with the (looser)
constraints provided by the observed thermodynamic profiles of large samples
of galaxy clusters (Sanders et al., 2018; Ghirardini et al., 2019a).

  4. (iv)

By propagating the modified laws that account for a dependency on both the
mass and the redshift, the i(cm)z model also allows one to modify any “virial”
quantity accordingly and to evaluate how it affects the rescaling of any
universal profile (see Sect. 4). We used our best-fit values in equations 7
and 8 in the scaling relations detailed in Table 1, and quote in Table 2 the
expected variations on the slope and redshift dependence of a given quantity
$Q$ estimated within an overdensity $\Delta$ (see eq. 9), finding remarkable
agreement with the current observational constraints.

  5. (v)

In Sect. 5, we have verified the self-consistency of the i(cm)z model by
evaluating the integrated quantities from the thermodynamic profiles (as
detailed in E20). The integrated values of some representative and fundamental
quantities, such as the gas mass fraction and global temperature, match the
predictions obtained from the modified and calibrated scaling laws within a
few percent. Only the bolometric luminosity shows differences up to 100%, but
only in low mass ($M_{500}<5\times 10^{13}M_{\odot}$) systems. When a
luminosity in a soft X-ray band is considered, the tension in these objects
decreases by more than a factor of two; part of this tension on the recovered
luminosity can be ascribed to the assumptions made as to the plasma
metallicity (see discussion in Sect. 7).

  6. (vi)

In Sect. 6, we have validated the i(cm)z model as a reliable indicator of the
level of hydrostatic bias present in the total mass reconstruction. By
applying the technique to the X-COP objects, a sample that offers high-quality
data from which a detailed analysis of their hydrostatic mass profile has been
produced (Ettori et al., 2019; Eckert et al., 2022), we assessed a level of
hydrostatic bias quite in agreement with the current available constraints
based on completely different assumptions (e.g., an universal gas mass
fraction of reference; see Eckert et al., 2019; Ettori & Eckert, 2022).

Figure 8: Distribution in the mass-redshift plane of the objects (X-COP,
Ghirardini et al. 2019; SPT, Sanders et al. 2018) used in the present work for
the calibration of the i(cm)z model, compared with the samples that will be
available in the near future (CHEX-MATE, CHEX-MATE Collaboration et al. 2021;
X-GAP).

Furthermore, we speculate that the i(cm)z model can be used to solve several
variegated issues as follows.

  * $-$

A natural application is the construction of the expected observational
quantities (from gas density to gas temperature and pressure profiles, from
X-ray surface brightness profiles in an assigned observational band to the
integrated X-ray luminosity, temperature, and gas mass) for each object in a
given catalog of haloes with an assigned mass and redshift.

  * $-$

Its predictions can be used to validate and calibrate the outcomes of
hydrodynamical simulations.

  * $-$

By fitting the observed thermodynamic radial profiles, it can provide an
estimate of the relevance of the hydrostatic bias in a given object.

  * $-$

The calibrated model predicts that the gas mass fraction depends on the gas
temperature, but it does not evolve with redshift, allowing for the systematic
uncertainties that affect the use of $f_{\rm gas}=M_{\rm gas}/M_{\rm tot}$ as
a cosmological ruler to be defined formally (see e.g., Ettori et al., 2009;
Mantz et al., 2022).

  * $-$

Once the thermodynamic profiles are calibrated, we predict that, for example,
the effective polytropic index $\gamma=\partial\log P/\partial\log n$ in the
range $0.8-1.2R_{500}$ should decrease with increasing halo masses and at
higher redshifts, spanning values within 5% (f.i. from 1.27 in halos with
$M_{500}=10^{13}M_{\odot}$ at $z=0.05$ to 1.14 for unlikely objects with
$M_{500}=10^{15}M_{\odot}$ at $z=2$) of the canonical value of 1.2 (see, e.g.,
Ghirardini et al., 2019b).

  * $-$

The most extreme objects in their observed properties can be identified and
their (eventual) tension quantified; for example, from equations 5 and 6 and
the best-fit results in equations 7 and 8, we predict that, in a $\Lambda$CDM
universe, a galaxy cluster with a global $T$ of 5 keV for example is expected
to be 4.5 times more luminous (bolometric luminosity; 3.5, if we refer to
$L_{500}$ in the 0.5-2 keV band) than a 3 keV system, and 1.6, 2.3, and 5.0
times more luminous at $z=[0.6,1,2]$ than at $z=0.05$. Significantly higher
(or lower) values, and outside the intrinsic scatter of the population
investigated, would represent a conundrum for the cosmological formation of
these systems.

Samples of spatially resolved thermodynamic properties out to $R_{500}$ at a
high signal-to-noise ratio for a large number of objects will be soon
available through numerous dedicated large XMM-Newton programs, such as CHEX-
MATE333http://xmm-heritage.oas.inaf.it/, a Planck SZ selected sample of 118
objects covering masses higher than $M_{500}\sim 2\times 10^{14}M_{\odot}$ up
to $z\approx 0.6$ (CHEX-MATE Collaboration et al., 2021) and
X-GAP444https://www.astro.unige.ch/xgap/, where 49 local galaxy groups with
$M_{500}$ between $10^{13}$ and $10^{14}M_{\odot}$ will be observed
homogeneously to measure $f_{\rm gas}(<R_{500})$ with statistical
uncertainties of about 10% (see Fig. 8). For these data, we will expand our
modelization by constraining the adopted parameters at higher confidence, and
introducing even new parameters (for instance related to the dynamical state),
if needed.

###### Acknowledgements.

We thank the anonymous referee for the positive and constructive comments.
S.E. and L.L. acknowledge financial contribution from the contracts 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), INAF mainstream project 1.05.01.86.10, and from the
European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant
agreement n. 871158). The i(cm)z model described in Sect. 2 is encoded as an
IDL function and available upon request from the first author.

## References

  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Andrade-Santos et al. (2017) Andrade-Santos, F., Jones, C., Forman, W. R., et al. 2017, ApJ, 843, 76
  * Angelinelli et al. (2022) Angelinelli, M., Ettori, S., Dolag, K., Vazza, F., & Ragagnin, A. 2022, A&A, 663, L6
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Bahar et al. (2022) Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7
  * Bhattacharya et al. (2013) Bhattacharya, S., Habib, S., Heitmann, K., & Vikhlinin, A. 2013, ApJ, 766, 32
  * Campitiello et al. (2022) Campitiello, M. G., Ettori, S., Lovisari, L., et al. 2022, A&A, 665, A117
  * CHEX-MATE Collaboration et al. (2021) CHEX-MATE Collaboration, Arnaud, M., Ettori, S., et al. 2021, A&A, 650, A104
  * Diemer & Kravtsov (2015) Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108
  * Dutton & Macciò (2014) Dutton, A. A. & Macciò, A. V. 2014, MNRAS, 441, 3359
  * Eckert et al. (2017) Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astronomische Nachrichten, 338, 293
  * Eckert et al. (2022) Eckert, D., Ettori, S., Pointecouteau, E., van der Burg, R. F. J., & Loubser, S. I. 2022, A&A, 662, A123
  * Eckert et al. (2021) Eckert, D., Gaspari, M., Gastaldello, F., Le Brun, A. M. C., & O’Sullivan, E. 2021, Universe, 7, 142
  * Eckert et al. (2019) Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, Astronomy and Astrophysics, 621, A40
  * Eckert et al. (2011) Eckert, D., Molendi, S., & Paltani, S. 2011, A&A, 526, A79
  * Ettori (2015) Ettori, S. 2015, MNRAS, 446, 2629
  * Ettori & Eckert (2022) Ettori, S. & Eckert, D. 2022, A&A, 657, L1
  * Ettori et al. (2019) Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, Astronomy and Astrophysics, 621, A39
  * Ettori et al. (2020) Ettori, S., Lovisari, L., & Sereno, M. 2020, A&A, 644, A111
  * Ettori et al. (2009) Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61
  * Ghirardini et al. (2022) Ghirardini, V., Bahar, Y. E., Bulbul, E., et al. 2022, A&A, 661, A12
  * Ghirardini et al. (2019a) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019a, Astronomy and Astrophysics, 621, A41
  * Ghirardini et al. (2019b) Ghirardini, V., Ettori, S., Eckert, D., & Molendi, S. 2019b, Astronomy and Astrophysics, 627, A19
  * Ghirardini et al. (2018) Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7
  * Lovisari et al. (2021) Lovisari, L., Ettori, S., Gaspari, M., & Giles, P. A. 2021, Universe, 7, 139
  * Lovisari et al. (2020) Lovisari, L., Schellenberger, G., Sereno, M., et al. 2020, ApJ, 892, 102
  * Ludlow et al. (2016) Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214
  * Mantz et al. (2022) Mantz, A. B., Morris, R. G., Allen, S. W., et al. 2022, MNRAS, 510, 131
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Nurgaliev et al. (2017) Nurgaliev, D., McDonald, M., Benson, B. A., et al. 2017, ApJ, 841, 5
  * Oppenheimer et al. (2021) Oppenheimer, B. D., Babul, A., Bahé, Y., Butsky, I. S., & McCarthy, I. G. 2021, Universe, 7, 209
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 550, A131
  * Pratt et al. (2022) Pratt, G. W., Arnaud, M., Maughan, B. J., & Melin, J. B. 2022, A&A, 665, A24
  * Pratt et al. (2010) Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85
  * Rossetti et al. (2016) Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515
  * Sanders et al. (2018) Sanders, J. S., Fabian, A. C., Russell, H. R., & Walker, S. A. 2018, MNRAS, 474, 1065
  * Sayers et al. (2022) Sayers, J., Mantz, A. B., Rasia, E., et al. 2022, arXiv e-prints, arXiv:2206.00091
  * Voit (2005) Voit, G. M. 2005, Reviews of Modern Physics, 77, 207

##  Appendix A Direct constraints on $a_{M}$ and $a_{z}$ from thermodynamic
profiles

In Fig. 9, we present the current constraints on $a_{M,P}$ and $a_{z,P}$
obtained by fitting the gas density and temperature profiles from the binned
NFW fits (model BIN-NFW) in Sanders et al. (2018) and Ghirardini et al.
(2019a) shown in Fig. 4. The gas density profiles were fitted over the radial
range $0.15-1R_{500}$, and the gas temperature profiles were fitted over the
range $0.15-0.8R_{500}$. A grid of $\chi^{2}$ values were estimated using Eq.
10, where the error budget on the model was obtained as described in Sect. 4
for a representative object of $(M_{500},z)=(4.6\times 10^{14}M_{\odot},0.6)$
and $(6\times 10^{14}M_{\odot},0.065)$ for S18 and G19, respectively.

Figure 9: Constraints on the exponents $a_{M}$ and $a_{z}$ in $P_{\Delta}\sim
M^{a_{M}}\;E_{z}^{a_{z}}$ based upon the gas density and temperature profiles
in S18 (top panel; likelihood contour at $1$ and $5\sigma$ equivalent to
$\Delta\chi^{2}$ of 6.2 and 28.8 for 2 degrees of freedom, respectively) and
G19 (bottom panel; on $a_{M}$ only, because all the 12 objects have a redshift
$<0.1$. For the sake of comparison, $5\sigma$ corresponds to $\Delta\chi^{2}$
of 25.0 for 1 degree of freedom). The dotted red lines are the predicted self-
similar slopes, and the dashed blue lines correspond to the expected slopes
from the i(cm)z model.

##  Appendix B Hydrostatic mass bias corrected universal pressure profile in
X-COP

We have obtained a new estimate of the universal pressure profile (UPP) from
the recent analysis of the X-COP data in Eckert et al. (2022) by fitting the
nonparametric deprojected points rescaled with $M_{500}$ corrected by the
nonthermal pressure contribution as described in Eckert et al. (2019). The
rescaled profiles are shown in Fig. 10.

Figure 10: Rescaled pressure profiles for the X-COP objects after the analyses
in Eckert et al. (2019, 2022). We note that $M_{500}$ and the relative
$R_{500}$ were obtained after the correction for the estimated contribution
from the nonthermal pressure (see Eckert et al. 2019, for details).

The scatter model is the log parabola from Ghirardini et al. (2019a) (see
their Eq. 6 and Fig. 7). The constraints on the free parameter of the model
are shown in Fig. 11. The profile differs from the UPP in Arnaud et al. (2010)
only slightly, with a lower gas pressure at intermediate radii
($0.1-1R_{500}$; see Fig. 12).

Figure 11: Corner plot showing the constraints on the seven free parameters of
the model in Eq. 2. Figure 12: Comparison between the universal pressure
profiles (UPP) from the revised analysis of X-COP data and the one in Arnaud
et al. (2010).

##  Appendix C Hydrostatic bias in the X-COP

Figure 13: Radial profiles of the gas density (left) and spectroscopic
temperature (center) recovered for the four X-COP objects and, overplotted,
the predictions obtained by assuming $b=0$ (red dashed line) and $b=0.6$
(purple dashed line). (Right panel) $\Delta\chi^{2}$ for the analysis using
only the $T$ profile (dashed line) and jointly with the $n$ profile (solid
line).

Figure 14: Continued.

Figure 15: Continued.