# The scatter and evolution of the global hot gas properties of simulated galaxy cluster populations

Amandine Le Brun, Ian G. Mccarthy, Joop Schaye, Trevor J. Ponman

<figure><img src="image_1.png" /></figure>

## To cite this version:

Amandine Le Brun, Ian G. Mccarthy, Joop Schaye, Trevor J. Ponman. The scatter and evolution of the global hot gas properties of simulated galaxy cluster populations. Monthly Notices of the Royal Astronomical Society, 2017, 466 (4), pp.4442-4469. 10.1093/mnras/stw3361. hal-01554004

HAL Id: hal-01554004

https://hal.science/hal-01554004v1

Submitted on 24 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The scatter and evolution of the global hot gas properties of simulated galaxy cluster populations

Amandine M. C. Le Brun,<sup>1,2*</sup> Ian G. McCarthy,<sup>2</sup> Joop Schaye<sup>3</sup> and Trevor J. Ponman<sup>4</sup>

Accepted 2016 December 22. Received 2016 December 16; in original form 2016 June 14

## ABSTRACT

We use the cosmo-OverWhelmingly Large Simulation (cosmo-OWLS) suite of cosmological hydrodynamical simulations to investigate the scatter and evolution of the global hot gas properties of large simulated populations of galaxy groups and clusters. Our aim is to compare the predictions of different physical models and to explore the extent to which commonly adopted assumptions in observational analyses (e.g. self-similar evolution) are violated. We examine the relations between (true) halo mass and the X-ray temperature, X-ray luminosity, gas mass, Sunyaev–Zel’dovich (SZ) flux, the X-ray analogue of the SZ flux ($Y_X$) and the hydrostatic mass. For the most realistic models, which include active galactic nuclei (AGN) feedback, the slopes of the various mass–observable relations deviate substantially from the self-similar ones, particularly at late times and for low-mass clusters. The amplitude of the mass–temperature relation shows negative evolution with respect to the self-similar prediction (i.e. slower than the prediction) for all models, driven by an increase in non-thermal pressure support at higher redshifts. The AGN models predict strong positive evolution of the gas mass fractions at low halo masses. The SZ flux and $Y_X$ show positive evolution with respect to self-similarity at low mass but negative evolution at high mass. The scatter about the relations is well approximated by log-normal distributions, with widths that depend mildly on halo mass. The scatter decreases significantly with increasing redshift. The exception is the hydrostatic mass–halo mass relation, for which the scatter *increases* with redshift. Finally, we discuss the relative merits of various hot gas-based mass proxies.

**Key words:** galaxies: clusters: general – galaxies: clusters: intracluster medium – galaxies: evolution – galaxies: formation – galaxies: groups: general – intergalactic medium.

## 1 INTRODUCTION

Galaxy clusters are potentially powerful systems for measuring fundamental cosmological parameters, such as the overall matter density of the Universe, the amplitude of the matter power spectrum, as well as the evolution of dark energy (for recent reviews, see Voit 2005; Allen, Evrard & Mantz 2011; Kravtsov & Borgani 2012).

The classical test is to compare the theoretically predicted and observed evolutions of the abundance of galaxy clusters. Since the abundance of dark matter haloes is a strong function of mass, to exploit clusters for cosmological purposes, one usually requires an accurate method for estimating the masses of individual clusters using (generally speaking) quite limited observational information

(but see Caldwell et al. 2016 and Ntampaka et al. 2016 for mass-independent tests). Furthermore, the scatter and covariance of the adopted mass–observable relations must be properly included in the cosmological modelling, and a detailed knowledge of the selection function of the survey is also necessary.

The use of theoretical models/simulations has become commonplace to calibrate mass–observable relations and to assess their scatter and biases (e.g. Kravtsov, Vikhlinin & Nagai 2006). Aside from assisting in the calibration of absolute mass measurements, theoretical models are required to predict the abundance of clusters as a function of their mass for a given set of cosmological parameters (e.g. Jenkins et al. 2001; Tinker et al. 2008). Here, we note that a number of recent studies have shown that the predicted abundance is sensitive to the details of feedback processes associated with galaxy formation (e.g. Cui, Borgani & Murante 2014; Cusworth et al. 2014; Martizzi et al. 2014; Velliscig et al. 2014), as

<footnote n="1"><sup>1</sup>Laboratoire AIM, IRFU/Service d'Astrophysique – CEA/DRF – CNRS – Université Paris Diderot, Bât. 709, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="2"><sup>2</sup>Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="3"><sup>3</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="4"><sup>4</sup>Astrophysics and Space Research Group, School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote>* E-mail: amandine.le-brun@cea.fr</footnote>
energetic feedback can eject baryons from collapsed structures (e.g. McCarthy et al. 2011) and lower their total mass, thereby reducing the number of haloes above a given mass threshold.

Ongoing and upcoming galaxy cluster surveys, such as the Dark Energy Survey (The Dark Energy Survey Collaboration 2005), *eRosita* (Merloni et al. 2012), *Euclid* (Laureijs et al. 2011), SPT-3G (Benson et al. 2014) and Advanced ACTpol (Henderson et al. 2016), will deliver samples of tens of thousands of galaxy clusters. With such large samples becoming available, the limiting uncertainties in the cosmological analyses will be due to systematic (as opposed to statistical) errors, the largest of which are likely to be associated with absolute mass calibration on the observational side and our incomplete knowledge of the effects of galaxy formation physics on the total masses and observable properties [e.g. X-ray luminosity, Sunyaev–Zel’dovich (SZ) flux, etc.] of clusters on the theoretical side.

Cosmological hydrodynamical simulations can help to address both of these issues, by calibrating the mass–observable relations (or, alternatively, by providing a testbed for direct observational mass reconstruction algorithms, such as synthetic weak lensing maps, which observers can use to calibrate the mass–observable relations empirically) and by providing a self-consistent framework for capturing the effects of galaxy formation physics on the predicted halo mass distribution. However, an important caveat is that predictions of current simulations are often sensitive to the details of the ‘sub-grid’ modelling of important feedback processes (Le Brun et al. 2014; McCarthy et al. 2016; Sembolini et al. 2016), as we will demonstrate here as well. Therefore, continual confrontation of the simulations with the observations (via production of realistic synthetic observations of the simulations) is also needed to test the realism of the former.

Encouragingly, the realism of simulations of galaxy clusters, in terms of their ability to reproduce the observed properties of local clusters, has been improving rapidly in recent years and is due in large part to the inclusion of energetic active galactic nuclei (AGN) feedback (e.g. Sijacki et al. 2007; Bhattacharya, Di Matteo & Kosowsky 2008; Puchwein, Sijacki & Springel 2008; Fabjan et al. 2010; McCarthy et al. 2010; Planelles et al. 2013; Le Brun et al. 2014; Planelles et al. 2014; McCarthy et al. 2016). However, much less is known about the realism of such simulations at higher redshifts, where there are fewer high-quality observations with which to compare the simulations, and there are significantly greater uncertainties in the role of important selection effects for the observed systems. The sparseness of high-quality observations of high-redshift systems means that cosmological tests using distant clusters will have to rely much more heavily on simulations, both to help calibrate the mass–observable relations and self-consistently predict the abundance of haloes in the presence of baryons. It is therefore crucial to examine what current simulations predict in terms of the evolution of the hot gas properties of clusters.

In this study, we use the cosmo-OverWhelmingly Large Simulation (cosmo-OWLS; Le Brun et al. 2014; McCarthy et al. 2014) suite of large-volume cosmological hydrodynamical simulations to conduct a study of the scatter and evolution of the global hot gas properties of large populations of galaxy groups and clusters as a function of the important non-gravitational physics of galaxy formation. We will examine to what extent the predicted scaling relations evolve self-similarly (both in terms of amplitude and slope), can be characterized by simple power laws with log-normal scatter and whether or not the scatter itself depends on mass and redshift. We will examine a large number of commonly used scaling relations, including the dependences on the true total halo mass of the gas

mass, (soft) X-ray luminosity, SZ flux, $Y_X$ (the X-ray analog of the SZ flux), mass-weighted and spectroscopic temperatures and X-ray hydrostatic mass. We also present an analysis of the evolution and scatter of the X-ray luminosity–temperature relation (in Appendix C).

This work extends upon an already relatively large body of previous studies of the impact of the non-gravitational physics of galaxy formation on the evolution of scaling relations such as those of Short et al. (2010), Stanek et al. (2010), Battaglia et al. (2010, 2012), Fabjan et al. (2011), Planelles et al. (2013) and Pike et al. (2014), which will be examined in the context of the results of this study in Section 9.

This paper is organized as follows. We briefly introduce the cosmo-OWLS simulation suite used here in Section 2 and summarize self-similar theory in Section 3. We then describe how we fit the mass–observable scaling relations and the scatter about them in Section 4. In Sections 5, 6 and 7, we examine the evolution of the mass slopes, amplitude and scatter, respectively, in the various mass–observable relations. We then examine the scatter and evolution of the hydrostatic bias in Section 8. Finally, we conduct a short comparison to previous studies in Section 9 and summarize and discuss our findings in Section 10.

Throughout the paper, masses, luminosities, temperatures, integrated SZ signal and X-ray analogues of the integrated SZ signal are quoted in physical $M_{\odot}$, erg s$^{-1}$, keV, Mpc$^2$ and $M_{\odot}$ keV, respectively; In denotes natural logarithm, while $\log_{10}$ corresponds to decimal logarithm.

## 2 SIMULATIONS

We take advantage of the cosmo-OWLS suite of cosmological simulations described in detail in Le Brun et al. (2014, hereafter L14; see also McCarthy et al. 2014; van Daalen et al. 2014; Velliscig et al. 2014; Le Brun, McCarthy & Melin 2015).

The cosmo-OWLS simulations constitute an extension to the OWLS project (Schaye et al. 2010). cosmo-OWLS was conceived with cluster cosmology in mind and is composed of large volume ($400\ h^{-1}\ Mpc$ on a side) periodic box simulations with $2 \times 10^{24}$ particles using updated initial conditions derived from the Planck data<fnref n="1" /> (Planck Collaboration XVI 2013) {$\Omega_m, \Omega_b, \Omega_\Lambda, \sigma_8, n_s, h$} = {$0.3175, 0.0490, 0.6825, 0.834, 0.9624, 0.6711$}. This results in dark matter (DM) and (initial) baryon particle masses of $\approx 4.44 \times 10^9\ h^{-1}\ M_{\odot}$ and $\approx 8.12 \times 10^8\ h^{-1}\ M_{\odot}$, respectively. The gravitational softening of the runs presented here is fixed to $4\ h^{-1}\ kpc$ (in physical coordinates below $z=3$ and in comoving coordinates at higher redshifts). We use $N_{ngb} = 48$ neighbours for the smoothed particle hydrodynamics (SPH) interpolation and the minimum SPH smoothing length is fixed to one tenth of the gravitational softening.

The simulations were conducted with a version of the Lagrangian TreePM-SPH code GADGET3 (last described in Springel 2005), which was significantly modified to incorporate new ‘sub-grid’ physics as part of the OWLS project. All the runs used here were started from identical initial conditions and only the included non-gravitational physics and some of its key parameters were methodically altered. We use four of the five physical models presented in L14 here:

<footnote n="1">¹ We also ran simulations using initial conditions derived from the 7-year Wilkinson Microwave Anisotropy Probe (WMAP) data. We will only present results from the Planck cosmology ones here, but will comment on any notable differences with the equivalent runs in the WMAP7 cosmology.</footnote>
<table><caption>**Table 1.** cosmo-OWLS runs presented here and their included sub-grid physics. Each model has been run in both the WMAP7 and Planck cosmologies.</caption><thead><tr><th>Simulation</th><th>UV/X-ray background</th><th>Cooling</th><th>Star formation</th><th>SN feedback</th><th>AGN feedback</th><th>ΔT<sub>heat</sub></th></tr></thead><tbody><tr><td>NOCOOL</td><td>Yes</td><td>No</td><td>No</td><td>No</td><td>No</td><td>—</td></tr><tr><td>REF</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>No</td><td>—</td></tr><tr><td>AGN 8.0</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>10<sup>8.0</sup> K</td></tr><tr><td>AGN 8.5</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td><td>10<sup>8.5</sup> K</td></tr></tbody></table>

**NOCOOL:** a non-radiative model. It includes net heating from the Haardt & Madau (2001) X-ray and ultraviolet photoionizing background, whose effects on the intracluster medium (ICM) are however negligible.

**REF:** this model also incorporates prescriptions for metal-dependent radiative cooling (Wiersma, Schaye & Smith 2009a), star formation (Schaye & Dalla Vecchia 2008), stellar evolution, mass-loss and chemical enrichment (Wiersma et al. 2009b) from Type II and Ia supernovae and asymptotic giant branch stars, and kinetic stellar feedback (Dalla Vecchia & Schaye 2008).

**AGN 8.0 and AGN 8.5:** in addition to the physics included in the REF model, these models include a prescription for supermassive black hole growth (through both Eddington-limited Bondi–Hoyle–Lyttleton accretion and mergers with other black holes) and AGN feedback (Booth & Schaye 2009, which is a modified version of the model originally developed by Springel, Di Matteo & Hernquist 2005). The black holes accumulate the feedback energy until they can heat neighbouring gas particles by a pre-determined amount $\Delta T_{\text{heat}}$. As in Booth & Schaye (2009), 1.5 per cent of the rest-mass energy of the gas that is accreted on to the supermassive black holes is employed for the feedback. This yields a satisfactory match to the normalization of the black hole scaling relations (Booth & Schaye 2009; see also L14), which is insensitive to the exact value of $\Delta T_{\text{heat}}$. The two AGN models used here only differ by their value of $\Delta T_{\text{heat}}$, which is the most important parameter of the AGN feedback model in terms of the gas-phase properties of the resulting simulated population of groups and clusters (McCarthy et al. 2011, L14). It is fixed to $\Delta T_{\text{heat}} = 10^{8.0}$ K for AGN 8.0 and $\Delta T_{\text{heat}} = 10^{8.5}$ K for AGN 8.5. Note that since the same quantity of gas is being heated in these models, more time is required for the black holes to accrete a sufficient amount of gas for heating the adjacent gas to a higher temperature. Hence, increased heating temperatures result into more episodic and more violent feedback episodes.

Table 1 provides a list of the runs used here and the sub-grid physics that they include.

Haloes are identified with a standard friends-of-friends (FoF) separation algorithm on the dark matter particles with a linking length of 0.2 in units of the mean interparticle separation. The baryonic content of the haloes is identified by locating the nearest DM particle to each baryonic (i.e. gas or star) particle and connecting it with the FoF group of the DM particle. Artificial haloes are removed by carrying out an unbinding calculation with the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009): any FoF halo that does not have at least one self-bound substructure (called subhalo) is removed from the FoF groups list.

Spherical overdensity masses $M_{\Delta}$ (where $M_{\Delta}$ is the total mass within a radius $r_{\Delta}$ that encloses a mean internal overdensity of $\Delta$ times the critical density of the Universe) with $\Delta = 200, 500$ and $2500$ have been computed for all the FoF haloes. The spheres are centred on the position of the minimum of the gravitational potential of the main subhalo (the most massive subhalo of the FoF

halo). Then, all galaxy groups and clusters with $M_{500} \ge 10^{13} M_{\odot}$ are extracted from each snapshot for analysis. There are roughly 25 000 such systems at $z=0$ in the NOCOOL run with the Planck cosmology, for example.

The SZ signal is characterized by the value of its spherically integrated Compton parameter $d_{\Delta}^{2} Y_{500} = (\sigma_{T}/m_{e}c^{2}) \int P dV$ where $d_{\Delta}$ is the angular diameter distance, $\sigma_{T}$ the Thomson cross-section, $c$ the speed of light, $m_{e}$ the electron rest mass, $P = n_{e}k_{B}T_{e}$ the electron pressure and the integration is done over the sphere of radius $r_{500}$. The X-ray equivalent of the SZ signal is characterized by $Y_{x, 500} = M_{\text{gas, 500}} T_{\text{spec, cor}}$ where $M_{\text{gas, 500}}$ is the gas mass enclosed within $r_{500}$ and $T_{\text{spec, cor}}$ is the core-excised X-ray spectroscopic temperature (computed within the $[0.15-1]r_{500}$ annulus).

Note that contrary to what was done in L14, in this study, we use true halo masses (as opposed to halo masses computed under the assumption of hydrostatic equilibrium) and that the X-ray luminosities, spectral temperatures, gas masses and integrated SZ signals were computed within the true $r_{500}$ aperture and in the case of core-excised quantities within the annulus $[0.15-1]r_{500}$ (as opposed to within $r_{500, hse}$ and $[0.15-1]r_{500, hse}$ where $r_{500, hse}$ is the value of $r_{500}$ obtained when the halo masses are computed assuming hydrostatic equilibrium). The spectral temperatures and X-ray luminosities were computed using the synthetic X-ray methodology presented in L14 though. The rationale behind these choices is that we aim to elucidate the relations between the hot gas observables and true halo mass, since those are useful for: (i) calibrating the mass–observable relations whose use is of paramount importance when doing (precision) cosmology with galaxy clusters; and (ii) making large synthetic surveys by applying template methods to large dark matter only simulations (e.g. Bode et al. 2007; Sehgal et al. 2007, 2010).

The various cosmo-OWLS models have been compared to a wide range of observational data in L14 and McCarthy et al. (2014, see also Hojjati et al. 2015; Le Brun et al. 2015). In L14, we concentrated on the comparison to low-redshift properties such as X-ray luminosities and temperatures, gas mass fractions, entropy and density profiles, integrated SZ signal, I-band mass-to-light ratio, dominance of the brightest cluster galaxy and central black hole masses and concluded that the fiducial AGN model (AGN 8.0) produces a realistic population of galaxy groups and clusters, broadly reproducing both the median trend and, for the first time, the scatter in physical properties over approximately two decades in mass ($10^{13} M_{\odot} \lesssim M_{500} \lesssim 10^{15} M_{\odot}$) and 1.5 decades in radius ($0.05 \lesssim r/r_{500} \lesssim 1.5$). In McCarthy et al. (2014), we investigated the sensitivity of the thermal SZ power spectrum to important non-gravitational physics and found that while the signal on small and intermediate scales is highly sensitive to the included galaxy formation physics, it is only mildly affected on large scales.

We note that no explicit attempt was made in cosmo-OWLS to calibrate the simulations to reproduce the hot gas properties of groups and clusters. As a consequence, some of the models perform better than others in terms of their comparison with
observational data. McCarthy et al. (2016) have recently calibrated the same simulation code to better reproduce the stellar masses of galaxies, while also reproducing the observed trend between hot gas mass fraction and halo mass as inferred from X-ray observations. The calibrated model is referred to as BAHAMAS. In terms of the local gas-phase properties, the BAHAMAS model is similar to the cosmo-OWLS AGN 8.0 model. We will briefly comment below on any significant differences in the predicted evolution of the mass–observable relations of the BAHAMAS and AGN 8.0 models, but defer a detailed study of the BAHAMAS model to future work (Barnes et al. 2017).

## 3 SELF-SIMILAR SCALINGS

The dominant force in the formation and evolution of galaxy clusters is gravity. Since gravity is scale free, galaxy clusters are, to ‘zeroth order’, expected to obey self-similarity, such that the properties only depend upon the cluster mass, and more massive galaxy clusters are just scaled up versions of less massive ones with a scaling factor that depends only upon the mass ratios (e.g. White & Rees 1978; Kaiser 1986; Voit 2005). While not expected to be strictly valid in a universe with a significant baryonic component<fnref n="2" /> (since baryons can radiate, which leads to star and black hole formation and then feedback), the self-similar model is still quite useful as a baseline for the interpretation of simulations and observations alike.

If one defines the total cluster mass (denoted as $M_\Delta$) as that contained within a region that encloses a mean overdensity $\Delta\rho_{crit}(z)^3$, then, under the assumption of self-similarity, one can predict both the redshift evolution (which is in this case only due to the increasing mean density of the Universe) and the slope of a given mass–observable relation. The redshift dependence comes from the evolution of the critical density for closure:

$$ \rho_{\text{crit}}(z) \equiv \frac{3H(z)^2}{8\pi G} = E(z)^2 \frac{3H_0^2}{8\pi G} = E(z)^2 \rho_{\text{crit}}(z=0), \quad (1) $$

where

$$ E(z) \equiv \frac{H(z)}{H_0} = \sqrt{\Omega_m(1+z)^3 + \Omega_\Delta} \quad (2) $$

gives the evolution of the Hubble parameter, $H(z)$, in a flat Lambda cold dark matter Universe. For instance, since

$$ M_{\Delta} \propto \Delta\rho_{\text{crit}}(z) r_{\Delta}^{3} \quad (3) $$

by definition, the cluster size scales as

$$ r_{\Delta} \propto M_{\Delta}^{1/3} E(z)^{-2/3}. \quad (4) $$

Gas falling into a cluster potential well is heated via shocks and will eventually settle and achieve approximate virial equilibrium

within that potential. The gas is then expected to have a temperature that is close to the virial temperature:

$$ k_B T_{\Delta} \propto -\frac{1}{2} \Phi = \frac{G M_{\Delta} \mu m_p}{2 r_{\Delta}} \quad (5) $$

where $k_B$ is the Boltzmann constant and $\mu$ is the mean molecular weight. Thus, the self-similar temperature–total mass relation can be obtained by combining equations (4) and (5) and is as follows:

$$ T_{\Delta} \propto M_{\Delta}^{2/3} E(z)^{2/3}. \quad (6) $$

The bolometric X-ray emission of massive clusters is dominated by thermal bremsstrahlung, implying that it scales as $L_X^{bol} \propto \rho^2 \Lambda(T) r_\Delta^3 \propto \rho^2 T^{1/2} r_\Delta^3$ as the cooling function $\Lambda(T) \propto T^{1/2}$, which combined with equations (4) and (6) gives the self-similar bolometric X-ray luminosity–total mass and bolometric X-ray luminosity–temperature relations (we examine the luminosity–temperature relation in Appendix C):

$$ L_{X,\Delta}^{\text{bol}} \propto M_{\Delta}^{4/3} E(z)^{7/3} \quad (7) $$

and

$$ L_{X,\Delta}^{\text{bol}} \propto T_{\Delta}^{2} E(z). \quad (8) $$

Hereafter, the bolometric X-ray luminosity $L_X^{bol}$ will be simply denoted as $L_X$.

Finally, $Y_{SZ,\Delta} \propto Y_{X,\Delta} \equiv M_{gas,\Delta} T_{\Delta} \propto M_{\Delta} T_{\Delta}$ assuming a constant gas fraction. Thus, the self-similar integrated SZ signal–total mass and $Y_X$–total mass relations follow from equation (6):

$$ Y_{X/SZ,\Delta} \propto M_{\Delta}^{5/3} E(z)^{2/3}. \quad (9) $$

With the launch of the first X-ray telescopes, such as the *Einstein Observatory*, *EXOSAT* and *ROSAT*, in the 1980s–1990s, it was quickly realized that the self-similar model was incompatible with the observations of the X-ray luminosity–temperature relation (e.g. Mushotzky 1984; Edge & Stewart 1991; Markevitch 1998; Arnaud & Evrard 1999; Lumb et al. 2004; Osmond & Ponman 2004; Pratt et al. 2009; Hilton et al. 2012), which was found to be significantly steeper than the self-similar expectation ($L_X \propto T^\alpha$ with $\alpha \simeq 2.5–3$ for clusters and likely even steeper for groups). This led to the conclusion that some non-gravitational processes, most likely connected to galaxy formation, must be breaking the self-similarity by introducing some physical scales (see for instance Evrard & Henry 1991; Kaiser 1991, for the first proposed solutions to this puzzle).

However, it should be noted that the X-ray luminosity is likely to be particularly sensitive to non-gravitational physics, since the luminosity is dominated by dense, centrally concentrated gas with short cooling times and is therefore likely to be significantly affected by radiative cooling and feedback. Other quantities, such as the mean temperature or integrated SZ flux, are not expected to deviate as strongly from the self-similar prediction, owing to their stronger contribution from gas at large radii, which is less affected by non-gravitational processes. Furthermore, even if (some of) the local relations do not strictly obey self-similarity with mass, the *evolution* can still be close to the self-similar prediction (and self-similarity with mass may be a better approximation to the truth at higher redshift than in the local Universe). At present, there are precious few constraints on the evolution of the hot gas properties of clusters that are robust to uncertainties in selection effects. For these reasons, self-similar evolution is still commonly adopted in cosmological analyses (e.g. Allen et al. 2008; Vikhlinin et al. 2009b; Planck Collaboration XX 2013; Mantz et al. 2014, 2015; Planck Collaboration XXIV 2015).

<footnote n="2">² Note that even if one completely neglects baryons and their associated non-gravitational physics (as in the case of a dark matter only simulation), haloes still do not strictly obey self-similarity. That is because in cold dark matter models the smallest objects collapse first, while the most massive (galaxy clusters) are still collapsing today. The internal structure of haloes is sensitive to the time of collapse via the evolution of the background density (e.g. Wechsler et al. 2002). Thus, while gravity may be scale free, the finite age of the Universe imprints a scale in structure formation.</footnote>

<footnote n="3">³ as is often the convention. One could also use a multiple of the mean matter density $\rho_m(z)$, which leads to predictions for self-similar evolution in terms of powers of $1+z$ instead of $E(z)$. As discussed later in the text, we have tried fitting for evolution using powers of $1+z$ and found that it decreases the quality of the fits.</footnote>
## 4 FITTING OF RELATIONS

Our characterization of the scaling relations is a two-stage process, in which we first derive the median mass–observable relations (i.e. T–M, L$_X$–M, M$_{gas}$–M, Y$_X$–M, Y$_{SZ}$–M and M$_{hse}$–M) as a function of redshift and then measure the scatter about these relations. Specifically, we first compute the median values of the observable Y (where hereafter Y denotes one of T, L$_X$, M$_{gas}$, Y$_X$, Y$_{SZ}$ and M$_{500, hse}$) in 10 equal-width logarithmic mass bins over the range 13.0 ≤ log$_{10}$[M$_{500}$(M$_\odot$)] ≤ 15.5 at nine different redshifts (z = 0, 0.125, 0.25, 0.375, 0.5, 0.75, 1.0, 1.25 and 1.5). To characterize the scatter about the median relations, we simply divide each cluster’s observable value by that expected from (a cubic spline interpolation of) the median mass–observable relation, Y$_{spline}$. We then fit the residuals (i.e. measure the scatter) in four bins of M$_{500}$ (chosen to be log$_{10}$[M$_{500}$(M$_\odot$)] = 13.0–13.5, 13.5–14.0, 14.0–14.5 and 14.5–15.5) with a log-normal distribution of the form:

$$P(X) = \frac{\text{binsize}}{X\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(\ln X - \mu)^2}{2\sigma^2}\right), \quad (10)$$

where $X = Y/Y_{spline}$.

When fitting the log-normal distribution, we fix $\mu = 0$ (i.e. the central value of the histogram is imposed to be the median observable in the corresponding mass bin i.e. $Y = Y_{spline}$) and use the MPFIT least-square minimization package in IDL (Markwardt 2009) to perform the fitting. Note that the binsize prefactor is needed to convert the log-normal Probability Distribution Function into a histogram whose bin width is equal to binsize. The binsize is set to the width (in linear units) of the distribution enclosing the central 90 per cent of the simulated systems divided by a factor of 10 (i.e. the width of the distribution is resolved with 10 bins).

For comparison, we have also computed the root mean square (RMS) dispersion about the median scaling relations as a function of mass and redshift as follows:

$$\sigma_{j,rms}(z) = \sqrt{\frac{\sum_{i=1}^{N_j} (\ln Y_i(z) - \ln Y_{spline,i}(z))^2}{N_j(z)}}, \quad (11)$$

where $N_j(z)$ is the number of system in mass bin j at redshift z and $Y_{spline,i}(z)$ is the value of Y(z) obtained for M$_{500,i}$ (z) using the best-fitting spline. Since the trends obtained using the RMS dispersion are nearly identical to the ones obtained with the log-normal scatter (the RMS scatter tends to be only slightly larger than the log-normal estimate), in the remainder of the paper, we only present the results obtained for the log-normal scatter described above.

From the above procedure, we obtain a median mass–observable scaling relation and the log-normal scatter about it (and its dependence on mass) as a function of redshift for each of the simulation models. We then fit an evolving power law of the form

$$Y = 10^A E(z)^\alpha \left( \frac{M_{500}}{10^{14} M_{\odot}} \right)^\beta, \quad (12)$$

to the median relations and the log-normal scatter about them.

Note that self-similar theory predicts that the mass–observable relations should be single power laws in both mass and E(z). However, as we will demonstrate below, a single power law in mass does not describe the mass–observable relations particularly well. For this reason, we also consider an evolving broken power law of the form

$$Y = 10^{A'} E(z)^{\alpha'} \left( \frac{M_{500}}{10^{14} M_{\odot}} \right)^{\epsilon'} , \quad (13)$$

where

$$\epsilon' = \begin{cases} \beta' & \text{if } M_{500} \le 10^{14} M_{\odot} \\ \gamma' & \text{if } M_{500} > 10^{14} M_{\odot} \end{cases} \quad (14)$$

Finally, we also try an evolving broken power law with a redshift-dependent low-mass power-law index of the form

$$Y = 10^{A''} E(z)^{\alpha''} \left( \frac{M_{500}}{10^{14} M_{\odot}} \right)^{\epsilon''} , \quad (15)$$

where

$$\epsilon'' = \begin{cases} \beta'' + \delta'' E(z) & \text{if } M_{500} \le 10^{14} M_{\odot} \\ \gamma'' & \text{if } M_{500} > 10^{14} M_{\odot} \end{cases} \quad (16)$$

The motivation for a redshift-dependent low-mass power-law index is that groups are more sensitive to non-gravitational physics than more massive clusters. Consequently, their evolution deviates more strongly from the self-similar prediction, which we demonstrate in detail below. The functional form of the redshift dependence of the low-mass power-law index was empirically inferred by looking at the redshift dependence of $\beta'$ (and $\gamma'$; see Fig. 3 and corresponding text).

Note that in all three cases (equations 12, 13 and 15), we held the mass pivot point fixed to $10^{14} M_{\odot}$. There are two main reasons for this: virtually all of the scaling relations appear to break at $M_{500} \sim 10^{14} M_{\odot}$ for the radiative models, and having a fixed pivot point makes the comparison between different physical models and fitting formulae straightforward.

Note that here $\chi^2$ is defined as

$$\chi^2 \equiv \sum_{i=1}^{N_{\text{bin}}} (Y_i - Y_{\text{bf},i})^2, \quad (17)$$

where $Y_{bf,i}$ is given by one of the following equations: (12) or (13) or (15); as no errors can straightforwardly be assigned to the variables (this is especially true when we fit the median scaling relations since, in this case $Y_i$ is held fixed to $Y_{i,spline}$ as is explained at the beginning of this section).<fnref n="4" /> For the same reason, the quoted errors for the best-fitting parameters should be taken ‘with a pinch of salt’ as they have been computed by rescaling the diagonal values of the covariance matrix computed by MPFIT so that the $\chi^2$ per degree of freedom is equal to 1 (according to the suggestion in the documentation of MPFIT). It should nevertheless be noted that they are, most of the time, much smaller than the differences between the different physical models. Additionally, error bars are not relevant when comparing different physical models, which were started from the same initial conditions as is the case here.

The results of fitting the evolving power law and broken power law with an evolving low-mass power-law index for the AGN 8.0 simulation (our most realistic model, see L14) are presented in Tables 2 and 3, respectively.<fnref n="5" /> Hereafter, $\beta$, $\epsilon'$ and $\epsilon''$ are called the mass slope.

It is worth mentioning that we also experimented with characterizing the evolution of median relations and the scatter about them using powers of $1+z$ instead of E(z), as is sometimes done in the literature for simulations and observations alike (e.g. Ettori et al. 2004; Maughan et al. 2006; Kay et al. 2007; Sehgal et al. 2011;

<footnote n="4">⁴ Therefore, the values of $\chi^2$ should only be used to compare the respective quality of the fits for the evolving power law and broken power law models at fixed scaling relation and physical model, as their differences and/or ratios are otherwise meaningless.</footnote>

<footnote n="5">⁵ The results for the other simulations are given in Appendix B.</footnote>
<table><caption>**Table 2.** Results of fitting the evolving power law (equation 12) to both the median relation and the log-normal scatter about it for the AGN 8.0 simulation. The scatter is the log-normal scatter in the natural logarithm of the *Y* variable. The results for the other simulations are given in Appendix B. Note that *α*<sub>SS</sub> and *β*<sub>SS</sub> represent the self-similar predictions for the evolution and mass slope, respectively.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th><th>α<sub>SS</sub></th><th>β<sub>SS</sub></th></tr></thead><tbody><tr><td>T<sub>spec, cor</sub>–M<sub>500</sub></td><td>Median</td><td>0.280 ± 0.004</td><td>0.356 ± 0.024</td><td>0.577 ± 0.006</td><td>0.044</td><td>30</td><td>2/3</td><td>2/3</td></tr><tr><td>L<sub>bol</sub>–M<sub>500</sub></td><td>Median</td><td>43.440 ± 0.015</td><td>2.920 ± 0.083</td><td>1.812 ± 0.019</td><td>0.496</td><td>30</td><td>7/3</td><td>4/3</td></tr><tr><td>M<sub>gas, 500–M<sub>500</sub></sub></td><td>Median</td><td>12.851 ± 0.012</td><td>0.576 ± 0.066</td><td>1.317 ± 0.015</td><td>0.322</td><td>30</td><td>0</td><td>1</td></tr><tr><td>Y<sub>X, 500–M<sub>500</sub></sub></td><td>Median</td><td>13.137 ± 0.014</td><td>0.909 ± 0.077</td><td>1.888 ± 0.018</td><td>0.436</td><td>30</td><td>2/3</td><td>5/3</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>–M<sub>500</sub></td><td>Median</td><td>-5.754 ± 0.014</td><td>0.981 ± 0.077</td><td>1.948 ± 0.018</td><td>0.431</td><td>30</td><td>2/3</td><td>5/3</td></tr><tr><td>M<sub>500, hse, spec</sub>–M<sub>500</sub></td><td>Median</td><td>13.902 ± 0.009</td><td>-0.027 ± 0.047</td><td>0.952 ± 0.011</td><td>0.161</td><td>30</td><td>0</td><td>1</td></tr><tr><td>T<sub>spec, cor</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.967 ± 0.026</td><td>-0.492 ± 0.142</td><td>-0.174 ± 0.033</td><td>1.473</td><td>30</td><td>–</td><td>–</td></tr><tr><td>L<sub>bol</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.430 ± 0.028</td><td>-0.464 ± 0.155</td><td>-0.133 ± 0.036</td><td>1.740</td><td>30</td><td>–</td><td>–</td></tr><tr><td>M<sub>gas, 500–M<sub>500</sub></sub></td><td>Scatter</td><td>-0.976 ± 0.013</td><td>-0.422 ± 0.069</td><td>-0.412 ± 0.016</td><td>0.349</td><td>30</td><td>–</td><td>–</td></tr><tr><td>Y<sub>X, 500–M<sub>500</sub></sub></td><td>Scatter</td><td>-0.750 ± 0.027</td><td>-0.441 ± 0.148</td><td>-0.250 ± 0.034</td><td>1.599</td><td>30</td><td>–</td><td>–</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.862 ± 0.025</td><td>-0.163 ± 0.136</td><td>-0.269 ± 0.032</td><td>1.344</td><td>30</td><td>–</td><td>–</td></tr><tr><td>M<sub>500, hse, spec</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.646 ± 0.018</td><td>0.352 ± 0.096</td><td>-0.118 ± 0.022</td><td>0.678</td><td>30</td><td>–</td><td>–</td></tr></tbody></table>

<table><caption>**Table 3.** Results of fitting the evolving broken power-law (equation 15) to both the median relation and the log-normal scatter about it for the AGN 8.0 simulation. The scatter is the log-normal scatter in the natural logarithm of the *Y* variable. The results for the other simulations are given in Appendix B.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A''</th><th>α''</th><th>β''</th><th>γ''</th><th>δ''</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec, cor</sub>–M<sub>500</sub></td><td>Median</td><td>0.314 ± 0.005</td><td>0.257 ± 0.020</td><td>0.703 ± 0.018</td><td>0.514 ± 0.009</td><td>-0.063 ± 0.011</td><td>0.013</td><td>28</td></tr><tr><td>L<sub>bol</sub>–M<sub>500</sub></td><td>Median</td><td>43.469 ± 0.018</td><td>2.590 ± 0.078</td><td>2.163 ± 0.069</td><td>1.846 ± 0.036</td><td>-0.259 ± 0.042</td><td>0.189</td><td>28</td></tr><tr><td>M<sub>gas, 500–M<sub>500</sub></sub></td><td>Median</td><td>12.940 ± 0.010</td><td>0.227 ± 0.042</td><td>1.738 ± 0.037</td><td>1.181 ± 0.020</td><td>-0.240 ± 0.023</td><td>0.056</td><td>28</td></tr><tr><td>Y<sub>X, 500–M<sub>500</sub></sub></td><td>Median</td><td>13.246 ± 0.010</td><td>0.504 ± 0.044</td><td>2.383 ± 0.039</td><td>1.712 ± 0.020</td><td>-0.275 ± 0.024</td><td>0.060</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>–M<sub>500</sub></td><td>Median</td><td>-5.637 ± 0.011</td><td>0.689 ± 0.045</td><td>2.336 ± 0.040</td><td>1.710 ± 0.021</td><td>-0.175 ± 0.024</td><td>0.062</td><td>28</td></tr><tr><td>M<sub>500, hse, spec</sub>–M<sub>500</sub></td><td>Median</td><td>13.895 ± 0.017</td><td>-0.011 ± 0.072</td><td>0.930 ± 0.063</td><td>0.967 ± 0.033</td><td>0.009 ± 0.039</td><td>0.160</td><td>28</td></tr><tr><td>T<sub>spec, cor</sub>–M<sub>500</sub></td><td>Scatter</td><td>-1.171 ± 0.024</td><td>-0.027 ± 0.103</td><td>-0.807 ± 0.091</td><td>0.256 ± 0.048</td><td>0.269 ± 0.055</td><td>0.327</td><td>28</td></tr><tr><td>L<sub>bol</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.623 ± 0.014</td><td>-0.313 ± 0.060</td><td>-0.448 ± 0.053</td><td>0.381 ± 0.028</td><td>0.013 ± 0.032</td><td>0.112</td><td>28</td></tr><tr><td>M<sub>gas, 500–M<sub>500</sub></sub></td><td>Scatter</td><td>-1.025 ± 0.021</td><td>-0.390 ± 0.088</td><td>-0.485 ± 0.077</td><td>-0.280 ± 0.041</td><td>-0.002 ± 0.047</td><td>0.240</td><td>28</td></tr><tr><td>Y<sub>X, 500–M<sub>500</sub></sub></td><td>Scatter</td><td>-0.936 ± 0.032</td><td>-0.084 ± 0.135</td><td>-0.762 ± 0.119</td><td>0.171 ± 0.062</td><td>0.189 ± 0.072</td><td>0.562</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>–M<sub>500</sub></td><td>Scatter</td><td>-1.037 ± 0.031</td><td>0.273 ± 0.130</td><td>-0.849 ± 0.115</td><td>0.085 ± 0.060</td><td>0.263 ± 0.070</td><td>0.525</td><td>28</td></tr><tr><td>M<sub>500, hse, spec</sub>–M<sub>500</sub></td><td>Scatter</td><td>-0.756 ± 0.024</td><td>0.566 ± 0.102</td><td>-0.424 ± 0.090</td><td>0.127 ± 0.047</td><td>0.115 ± 0.055</td><td>0.324</td><td>28</td></tr></tbody></table>

Lin et al. 2012). However, we found that this generally leads to poorer fits and we will thus not discuss the results obtained using powers of 1 + *z* any further.

In Fig. 1, we show examples of our fits at various redshifts for one of the variables studied here: the gas mass, which is a representative variable. In each panel, the grey dots correspond to the individual simulated groups and clusters with $\log_{10}[M_{500}(M_{\odot})] \ge 13.0$ taken from the AGN 8.0 model, the solid blue, green and red lines, respectively, correspond to the best-fitting evolving power law (equation 12), broken power law (equation 13) and broken power law with a redshift-dependent low-mass power-law index (equation 15) to the median gas mass in bins of mass and redshift and the dashed red lines correspond to the best-fitting evolving broken power law with a redshift-dependent low-mass power-law index to the log-normal scatter in bins of mass and redshift. For $z \le 1.5$, the median relations and the scatter about them are reasonably well modelled by evolving broken power laws with redshift-dependent low-mass power-law indices of the form given by equation (15), whereas power laws and broken power laws of the form given by equations (12) and (13) fail to reproduce the median relations, especially at the low-mass end.

## 5 EVOLUTION OF THE MASS SLOPE

We start by examining the evolution of the logarithmic slope (the mass slope) of the total mass-observable relations. Note that, by

definition, no evolution of the mass slope is expected in the context of the self-similar model. The slope of a particular relation is fixed and can be predicted assuming only Newtonian gravity and that the gas is in virial equilibrium (see Section 3). Any evolution or deviation at any redshift from the predicted mass slope signals that either some non-gravitational physics is at play, or that the gas is not virialized (or both).

In Figs 2 and 3, we show the evolution of the mass slopes from $z=0$ to $1.5$ for the scaling relations between total mass and core-excised<fnref n="6" /> temperature (for both mass-weighted and X-ray spectroscopic temperature), bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal for each of the four physical models (different coloured curves). Fig. 2 shows the evolution of the mass slope obtained when we fit a single power law (equation 12) at each individual redshift [and so omitting the $E(z)$ factor for the moment], while Fig. 3 shows the evolutions of the low-mass (left-hand panel) and high-mass (right-hand panel) mass slopes resulting from the fitting of the broken power law in equation (13) to the median scaling relations at each individual redshift. In each panel, the solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation.

Starting first with the mass-X-ray temperature relation, the fitted mass slope is slightly shallower than the self-similar expectation of 2/3 for all the models. This result is mostly independent of redshift

<footnote n="6">⁶ The results for non-core excised temperatures are presented in Appendix A.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 1.** Gas fraction (in units of the universal baryon fraction $\Omega_b/\Omega_m$)–$M_{500}$ relation at six different redshifts ($z = 0, 0.25, 0.5, 0.75, 1.0$ and $1.5$ from top left to bottom right). In each panel, the grey dots correspond to the individual simulated groups and clusters with $\log_{10}[M_{500}(M_{\odot})] \ge 13.0$, the solid blue, green and red lines, respectively, correspond to the best-fitting evolving power law (equation 12), broken power law (equation 13) and broken power law (equation 15) with a redshift-dependent low-mass power-law index to the median gas mass. The dashed red lines correspond to the best-fitting evolving broken power law with a redshift-dependent low-mass power-law index to the log-normal scatter in bins of mass and redshift. The median relations and the scatter about them are reasonably well modelled by evolving broken power laws with redshift-dependent low-mass power-law indices of the form given by equation (15), whereas the other functional forms fail to reproduce the median relations, especially at the low-mass end.

and mass, but does depend somewhat upon the included sub-grid physics. The sensitivity to sub-grid physics is stronger when using the observable spectroscopic temperature, as opposed to the mass-weighted temperature.

We note that changing the sub-grid physics can affect the mean temperature in several ways. First, the mean temperature profile can be altered, because the mean entropy of the gas can be raised or lowered by including feedback and radiative cooling. The degree of scatter about the mean temperature profile (i.e. ‘multiphase’ structure) will also be affected. Energetic feedback processes can drive outflows and introduce turbulence, so that the temperature of the gas is no longer just determined by the entropy configuration of the gas and the potential well depth. In addition, the degree of clumpiness of the gas, which is affected by feedback, will influence the observable mean spectroscopic temperature, since denser gas contributes more to the X-ray emissivity. (The same is true for the gas-phase

<figure><img src="image_5.png" /></figure>

**Figure 2.** Evolution of the mass slope of the scaling relations between total mass and core-excised temperature (for both mass-weighted and X-ray spectroscopic temperature), bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from top left to bottom right). In each panel, we plot the evolution of the best-fitting power-law indices obtained by fitting the power law given by equation (12) at each individual redshift. The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The AGN models show significant deviations from self-similarity for all the scaling relations, except for the mass–temperature relation for which only a mild deviation is predicted (independent of the included sub-grid physics). The deviations from self-similarity increase with decreasing redshift for the AGN models. The models that lack efficient feedback (i.e. NOCOOL and REF) show approximately self-similar behaviour for most scaling relations.

metallicity.) However, in spite of all of these possible effects on the mean temperature, fundamentally, the mean temperature of the gas cannot deviate very strongly from the virial temperature, otherwise it will not be close to equilibrium with the total potential well, in which case the gas will readjust itself in a few sound-crossing or dynamical times (in practice, the minimum of the two) in order to achieve equilibrium.

The picture is not as clear-cut in the case of the bolometric X-ray luminosity–total mass relation, for which the mass slope displays strong simultaneous redshift and non-gravitational physics dependencies. When the self-similar prediction of 4/3 (see equation 7) is considered, even the non-radiative simulation (NOCOOL) and the simulation that neglects AGN feedback altogether (REF) have bolometric X-ray luminosity–mass scaling relations that are slightly steeper than self-similar; this is nearly independent of redshift for both models but only at high mass for REF. The deviations from self-similarity are probably due to the fact that the gas does not trace
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Evolution of the mass slopes of the scaling relations between total mass and core-excised temperature (for both mass-weighted and X-ray spectroscopic temperature), bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from top left to bottom right). In each subpanel, we plot the evolution of the low-mass (left-hand panel) and high-mass (right-hand panel) best-fitting power-law indices obtained by fitting the broken power law given by equation (13) at each redshift independently. The solid curves (red, orange, blue and green) correspond to the different models and the horizontal dashed lines to the self-similar expectation, respectively. With the exception of the X-ray temperature and the soft X-ray luminosity, for both the low-mass and the high-mass slopes, the non-radiative (NOCOOL) model and the REF model (which neglects AGN feedback) are approximately consistent with the self-similar expectation, whereas the models that include AGN feedback deviate significantly from self-similarity. The main differences between the two power-law indices at fixed scaling relation and physical model is that the low-mass one displays a stronger redshift dependence and tends to deviate more strongly from self-similarity.</figcaption></figure>

the dark matter (it has e.g. a different mass–concentration relation compared to the dark matter) which affects both the density and the temperature, and hence the X-ray luminosity. The temperature and density are also potentially affected by non-thermal pressure support. It is worth mentioning that this scaling relation is steeper in the calibrated BAHAMAS model than for both of the AGN models discussed here and that this is only the case at the low-mass end.

The total mass–gas mass relation is steeper than the self-similar slope of 1 (which corresponds to gas tracing total mass, with a constant gas fraction) for all of the radiative models, whereas it is consistent with self-similarity for the non-radiative (NOCOOL) model. These results are approximately independent of mass and redshift for the NOCOOL and REF models. For the models that include AGN feedback, the mass slope steepens with decreasing redshift and mass, moving away from the self-similar expectation at low redshifts. The low-mass regime shows a particularly large deviation with respect to the self-similar result (see Fig. 3), with the slope approaching 3/2 at late times. The fact that the mass slope deviates most strongly from the self-similar result at low masses makes sense, since the AGN feedback is more efficient at ejecting gas from the high-redshift progenitors of groups than those of massive clusters (McCarthy et al. 2011). That the mass slope

moves further away from the self-similar result with decreasing redshift is also understandable, since haloes of fixed mass (as we consider here) have shallower potential wells at low redshift, as we discuss further below. It is noteworthy that the gas mass–total mass relation is the only other scaling (besides the soft X-ray luminosity–total mass one) whose slope has been affected by the calibrations undertaken as part of the BAHAMAS project: it is slightly steeper at the low-mass end than both AGN 8.0 and 8.5.

In the bottom panels of Figs 2 and 3, we examine the slopes of the relations between total mass and the integrated SZ signal and its X-ray analogue, $Y_X$. The models that neglect feedback (NOCOOL and REF) exhibit approximately self-similar behaviour over the full range of redshifts we explore. The inclusion of AGN feedback, however, results in mass slopes that are significantly steeper than the self-similar expectation of 5/3, but move closer to the self-similar result with increasing mass and redshift, in analogy to the gas mass trends discussed above. Taken at face value, the results obtained here for the SZ signal in the AGN models conflict with the observational results of Planck Collaboration XI (2013) and Greco et al. (2015). But, as discussed in detail in Le Brun et al. (2015) (see especially their fig. 5 and also fig. 20 of McCarthy et al. 2016), due to its limited resolution (~5 arcmin beam), Planck does not measure the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of the normalizations of the scaling relations between total halo mass and core-excised temperature (for both mass-weighted and X-ray spectroscopic temperature), bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from the top left subpanel to the bottom right subpanel). The amplitude of each scaling relation in the four $\log_{10}[M_{500}(M_{\odot})]$ bins (denoted by solid lines of different colours) has been normalized by the self-similar expectation for the redshift evolution at fixed mass (shown as a horizontal dashed line). The different panels (continued over the page) correspond to the different physical models. The mass-temperature relation evolves in a negative fashion with respect to the self-similar model (i.e. slower), independently of the included ICM physics. The gas mass evolves approximately self-similarly for the non-radiative simulation but shows a positive evolution (i.e. faster than self-similar) when radiative cooling and particularly AGN feedback is included, a result which is strongly mass dependent. The SZ flux and $Y_X$ evolve negatively with respect to the self-similar expectation for models that neglect efficient feedback, driven by the negative evolution of the temperature. When AGN feedback is included, the *sign* of evolution of the SZ flux and $Y_X$ with respect to self-similar depends on halo mass, driven by the strong halo mass dependence of the gas mass evolution combined with the negative evolution of the temperature.</figcaption></figure>

SZ flux within $r_{500}$ but within a much larger aperture of size $5r_{500}$. This effectively diminishes the sensitivity of the SZ signal to the impact of the non-gravitational physics of galaxy formation.

Finally, it is worth pointing out that all the above results are independent of the choice of cosmology.

In short, we find that our most realistic models that include efficient feedback from AGN predict significant deviations in the mass slope from the self-similar expectation for all of the scaling relations we have examined. The one exception to this is the mass-temperature relation (especially in the mass-weighted case), where only a small deviation from the self-similar expectation is found (this is generally true, independent of the details of the included ICM physics). For the other scaling relations, all of which depend directly on the gas density/mass, the deviations from the self-similar prediction are strongest at low redshifts and low halo masses. The models that neglect efficient feedback (NOCOOL and REF), on the other hand, have mass slopes that are approximately consistent with self-similar expectations.

# 6 EVOLUTION OF THE NORMALIZATION

In Fig. 4, we show the evolution of the normalization of the scaling relations between total mass and core-excised temperature (for both mass-weighted and X-ray spectroscopic temperature), bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from top left to bottom right). The amplitude of each scaling relation in the four $\log_{10}[M_{500}(M_{\odot})]$ bins (denoted by solid lines of different colours) has been normalized by the self-similar expectation for the redshift evolution at fixed mass (shown as an horizontal dashed line). In the remainder of the paper, a scaling relation whose $E(z)$ exponent is smaller (larger) than the self-similar expectation listed in Section 3 will be referred to as having negative (positive) evolution.

We first consider the evolution of the mass-temperature relation. Interestingly, the normalization of this relation evolves in a negative sense with respect to the self-similar expectation, such that the spectroscopic (mass-weighted) temperature is predicted to be $\approx 30$ per cent ($\approx 15$ per cent) lower than predicted using self-similar arguments by $z=1$. This negative evolution is approximately
<figure><img src="image_8.png" /><figcaption>Figure 4 - continued</figcaption></figure>

independent of halo mass and only weakly dependent on the included gas physics. That the negative evolution is approximately independent of the included physics implies that it is mainly a consequence of the merger history of clusters. Two possible causes are: a change in the mass structure of haloes with redshift (i.e. evolution of the mass-concentration relation) and/or a change in the degree of virialization of the gas with redshift. We expect the former effect to be quite weak, as it only affects the mean temperature through a slight change in the weighting of the particles when calculating the mean.

We have examined the degree of virialization of the gas as a function of redshift, by tracking the evolution of the ratio of kinetic to thermal energy of the hot gas in different fixed mass bins (see Fig. 5). We find that, independently of the included gas physics, haloes of fixed mass have strongly increasing kinetic-to-thermal energy ratios with increasing redshift, evolving from a typical value of 10–15 per cent at $z=0$ up to 20–30 per cent by $z=1$. These kinetic motions (in both bulk flows and turbulence) act as a source of non-thermal pressure support that increases with increasing redshift, implying that a lower temperature (with respect to low-redshift clusters of the same mass) is required to achieve equilibrium within the overall potential well at high redshift and therefore likely drive the negative evolution (with respect to self-similarity) of the mass-temperature relation. We therefore warn against the ‘simplistic’ interpretation of deviations of any scaling relation from self-similar evolution as indication of the effects of non-gravitational physics.

The amplitude of the bolometric X-ray luminosity–total mass relation evolves positively for all the models that include efficient feedback (i.e. the AGN models). The amplitude of the evolution

is strongly mass dependent, slightly redshift dependent (it flattens out as redshift increases) and is strongly sensitive to the non-gravitational physics of galaxy formation (it becomes more positive as the feedback intensity is increased with a reversal from mostly negative to mostly positive when AGN feedback is included). The slight negative evolution in the models without AGN feedback is most likely due to the negative evolution of the mass–temperature relation whereas the positive evolution is linked to ‘ease’ of gas execution (see the discussion about the evolution of the total mass–gas mass relation below).

The total mass–gas mass relation is approximately consistent with self-similar evolution for the non-radiative model, but exhibits positive evolution when non-gravitational physics is included, particularly when the feedback ‘intensity’ is increased (going from REF to AGN 8.0 to AGN 8.5). For our most realistic models (i.e. the AGN models), the evolution is strongly mass dependent and somewhat redshift dependent.

A likely explanation for the strong positive evolution of the gas mass (and X-ray luminosity) is that, since haloes of fixed mass are denser at higher redshifts, more energy is required to eject gas from these higher redshift haloes. More precisely, the binding energy can be approximated by

$$E_{\text{bind}} \propto \frac{GM_{\Delta}^2}{r_{\Delta}}, \qquad (18)$$

which combined with equation (4) gives

$$E_{\text{bind}}(z) \propto M_{\Delta}^{5/3} E(z)^{2/3}. \qquad (19)$$
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Evolution of the kinetic-to-thermal ratio from $z=0$ to $z=1.5$ in each of the four $\log_{10}[M_{500}(M_{\odot})]$ bins (denoted by solid lines of different colours). The kinetic-to-thermal ratio increases with increasing redshift and that independently of mass and physical models. The introduction of AGN feedback notably increases the importance of the kinetic motions in the galaxy groups regime. These results are independent of cosmology.</figcaption></figure>

Hence, the binding energy increases with redshift making expulsion of gas due to outflows more difficult. Note that the supermassive black holes that power the AGN could ‘know’ about the evolution of the binding energy of their host dark matter halo in the sense that the black hole masses are determined by their halo binding energy through their self-regulation (see Booth & Schaye 2010, 2011). However, the growth of black holes in massive galaxies may be governed by black hole mergers rather than by self-regulated gas accretion (e.g. Peng 2007).

The integrated SZ flux and $Y_X$ show perhaps the most interesting behaviour, in terms of the evolution of the amplitude of their relations with halo mass. For models with inefficient feedback, there is a mild negative evolution with respect to the self-similar expecta-

tion, which is driven by the negative evolution of the temperature combined with the nearly self-similar evolution of the gas mass (note that $Y_X$, $SZ$, $\Delta \equiv M_{\text{gas}}, \Delta T_{\Delta}$). Things get more interesting when AGN feedback is included. In particular, low-mass haloes display a positive evolution with respect to the self-similar result, whereas high-mass haloes show negative evolution. This change in the sign of the effect with respect to the self-similar model is driven by the strong halo mass dependence of the gas mass evolution. In particular, low-mass haloes show a strong positive evolution in the gas mass that more than compensates for the negative evolution in the temperature, leading to a positive evolution in the SZ flux and $Y_X$. High-mass haloes, however, show little evolution in the gas mass (their gas fractions are already near the universal fraction $\Omega_b/\Omega_m$)
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Evolution of the log-normal scatter about the scaling relations between total mass and core-excised X-ray spectroscopic temperature, bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from the top subpanel to the bottom subpanel). For each simulation and each scaling relation, we plot the log-normal scatter as a function of $M_{500}$ and denote the redshift using lines of different colours. The different columns (continued over the page) correspond to the different physical models. For most scaling relations, the log-normal scatter varies only mildly with mass, is somewhat sensitive to non-gravitational physics, but displays a moderately strong redshift dependence (it tends to decrease with increasing redshift). With the exception of the X-ray luminosity, all the examined hot gas mass proxies have a similar scatter at fixed total mass of about 10 per cent.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6 - continued</figcaption></figure>

and, when combined with the negative evolution in the temperature, this leads to a negative evolution of the SZ flux and $Y_X$ with respect to the self-similar expectation.

We note that all the results described in this section are approximately independent of the choice of cosmology, in the sense that

the general trends are preserved but the exact values of e.g. the $E(z)$ exponents are slightly different. The only noteworthy difference is that the highest mass bin evolves somewhat faster for all the scalings considered in the WMAP7 runs compared to the Planck runs. The calibrated BAHAMAS model has an evolution that is
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Log-normal scatter of the scaling relations between total mass and core-excised X-ray spectroscopic temperature, X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from top to bottom) at $z = 0$ (left-hand panel) and $z = 1$ (right-hand panel). For each scaling relation and each redshift, we plot the log-normal scatter as a function of $M_{500}$ and denote the physical model using lines of different colours. The amplitude of the scatter tends to increase with increasing feedback (be it from SN or AGN) intensity at both redshifts.</figcaption></figure>

bracketed by that of AGN 8.0 and AGN 8.5 for all the scalings but the mass–temperature one (for both mass-weighted and X-ray temperature).xref ref-type="fig" rid="fig8"/>

To summarize, the mass–temperature relation evolves in a negative fashion with respect to the self-similar model, independently of the included ICM physics. The gas mass evolves approximately self-similarly for the non-radiative simulation but shows a positive evolution when radiative cooling and particularly AGN feedback is included, a result that is strongly mass dependent. The SZ flux and $Y_X$ evolve negatively with respect to the self-similar expectation for

models without efficient feedback, driven by the negative evolution of the temperature. When AGN feedback is included, however, the sign of evolution of the SZ flux and $Y_X$ with respect to self-similar depends on halo mass (positive evolution for low-mass haloes, negative evolution for high-mass haloes), driven by the strong halo mass dependence of the gas mass evolution combined with the negative evolution of the temperature.

## 7 EVOLUTION OF THE SCATTER

Having quantified the evolution of the median relations between total mass and the various hot gas observables, we now move on to an examination of the scatter about these relations. Quantifying the degree of scatter (and its possible dependence on halo mass, redshift and ICM physics) is extremely important for cosmological studies using cluster abundances, since the associated Eddington bias can be significant due to the steepness of the halo mass function.

In Fig. 6, we show the evolution of the log-normal scatter (at fixed total mass) about the scaling relations between total halo mass and core-excised X-ray spectroscopic temperature, bolometric X-ray luminosity, gas mass, $Y_X$ and the integrated SZ signal (from top to bottom). For each simulation and each scaling relation, we plot the log-normal scatter as a function of $M_{500}$ and denote the redshift using lines of different colours.

For most scaling relations, the log-normal scatter varies only mildly with mass, usually changing by less than a factor of ≈2–3 over 1.5 decades in $M_{500}$. Nevertheless, it is clearly mass dependent in all of the models, whereas virtually all current observational studies adopt mass-invariant scatter.

The amplitude of the scatter at fixed redshift is sensitive to the included non-gravitational physics, with the amplitude tending to increase with increasing complexity (and realism) of the included galaxy formation physics and with increasing AGN feedback intensity, as can be seen more easily in Fig. 7 (for both $z = 0$ and 1). This is perhaps not so surprising, as we know that AGN feedback is having a large effect on the mean properties.

On the positive side (for observational studies), the amplitude of the scatter tends to decrease with increasing redshift, which is likely due to the lessening importance of non-gravitational physics at higher redshifts when examining the evolution of haloes of fixed mass. Typically, the scatter is reduced by ~50 per cent from $z = 0$ to 1.

We note that the results above are generally robust to changing the cosmology from Planck to WMAP7.

Table 4 summarizes the scatter about the median scaling relations at $z = 0$. The scatter in the observable at fixed mass, $\sigma_{\ln Y|M}$, is

<table><caption>**Table 4.** Scatter about the median scaling relations at $z = 0$. $\sigma_{\ln Y|M}$ was computed using the best-fitting $A''$ values for the log-normal scatter of the AGN 8.0 model listed in Table 3. $\sigma_{\ln M|Y}$ was computed as a sigma-clipped RMS deviation from the (cubic spline interpolation of) median $M-Y$ relation (analogous to equation 11) of the AGN 8.0 model. The zero-point uncertainty of $Y$ was estimated as the RMS dispersion of the normalization $A$ of the best-fitting evolving power law obtained for the median relations of the three radiative models (REF, AGN 8.0 and 8.5) and by multiplying it by $\ln 10$. The observational constraints on $\sigma_{\ln Y|M}$ are summarized in the final column.</caption><thead><tr><th>Scaling relation</th><th>σ<sub>ln Y|M</sub></th><th>σ<sub>ln M|Y</sub></th><th>Zero-point uncertainty in Y</th><th>Observational constraints on σ<sub>ln Y|M</sub></th></tr></thead><tbody><tr><td>T<sub>spec, cor</sub>–M<sub>500</sub></td><td>≈5 per cent</td><td>≈20 per cent</td><td>≈5 per cent</td><td>≈15 per cent</td></tr><tr><td>L<sub>bol</sub>–M<sub>500</sub></td><td>≈25 per cent</td><td>≈25 per cent</td><td>≈40 per cent</td><td>≈40 per cent</td></tr><tr><td>M<sub>gas, 500</sub>–M<sub>500</sub></td><td>≈10 per cent</td><td>≈10 per cent</td><td>≈25 per cent</td><td>≈15 per cent</td></tr><tr><td>Y<sub>X, 500</sub>–M<sub>500</sub></td><td>≈10 per cent</td><td>≈15 per cent</td><td>≈25 per cent</td><td>≈20 per cent</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>–M<sub>500</sub></td><td>≈10 per cent</td><td>≈10 per cent</td><td>≈20 per cent</td><td>≈20 per cent</td></tr><tr><td>M<sub>500, hse, spec</sub>–M<sub>500</sub></td><td>≈15 per cent</td><td>≈25 per cent</td><td>≈5 per cent</td><td>—</td></tr></tbody></table>
computed using the best-fitting $A''$ values for the log-normal scatter of the AGN 8.0 model listed in Table 3. $\sigma_{\ln M^Y}$ was computed as a sigma-clipped (using a threshold value of $2\sigma$ for the clipping) RMS deviation from the (cubic spline interpolation of) median $M-Y$ relation (analogous to equation 11) of the AGN 8.0 model. The zero-point uncertainty of $Y$ was estimated by computing the RMS dispersion of the normalization $A$ of the best-fitting evolving power law obtained for the median relations of the three radiative models (REF, AGN 8.0 and 8.5) and by multiplying it by $\ln 10$. The unphysical non-radiative simulation (NOCOOL) was excluded from the computation.

All but one of the hot gas proxies examined here have a similar scatter at fixed total mass of order 10 per cent. Note, however, that the mass dependence is significant (see Fig. 7). The bolometric X-ray luminosity, which is the odd one out, has a significantly larger scatter at fixed total mass (it is about three times higher). The theoretical relations exhibit values of $\sigma_{\ln Y|M}$ that are slightly smaller than some of the recent observational constraints, which are summarized in the rightmost column of Table 4 (see e.g. Pratt et al. 2009; Vikhlinin et al. 2009a; Mantz et al. 2010b; Andersson et al. 2011; Planck Collaboration XI 2011; Lin et al. 2012; Planck Collaboration III 2013; Planck Collaboration XX 2013; Allen et al. 2011 for a review). We caution that this difference is at least partly due to the fact that we are not using full synthetic X-ray observations of the ‘depth’ typical of these observational data sets. A careful comparison would also require us to fold in the full observational selection functions of these data sets, which is beyond the scope of this work. We note, however, that the values obtained for $\sigma_{\ln Y|M}$ are similar to those found in several previous simulation studies (see for instance Kravtsov et al. 2006; Nagai 2006; Stanek et al. 2010; Fabjan et al. 2011; Kay et al. 2012).

## 8 SCATTER AND EVOLUTION OF THE HYDROSTATIC BIAS

High-quality X-ray observations that allow for the measurement of spatially resolved temperature and density profiles can be used to derive a direct estimate of the total mass via the assumption of hydrostatic equilibrium. However, clusters are not expected to be perfectly in equilibrium, and simulations are often relied upon to calibrate the degree of ‘hydrostatic bias’ that is expected. While such hydrostatic analysis is generally applicable to local systems only, owing to the rapidly declining X-ray surface brightness with increasing redshift, it is nevertheless interesting to examine the evolution and scatter in the hydrostatic bias and its sensitivity to the included gas physics, as future X-ray missions such as eRosita and eventually Athena will be able to extend this kind of analysis to higher redshifts. (Stacking of large numbers of high-redshift systems in current data, to derive mean hydrostatic halo mass estimates, is also possible.) We note that there has already been considerable attention devoted to quantifying the hydrostatic bias in simulated clusters (e.g. Evrard, Metzler & Navarro 1996; Kay et al. 2007, 2012; Nagai, Vikhlinin & Kravtsov 2007a; Piffaretti & Valdarnini 2008; Meneghetti et al. 2010; Battaglia et al. 2012; Rasia et al. 2012, 2014). However, most of these studies examined only relatively small samples of clusters in zoomed simulations and could not explore the system-to-system scatter in the hydrostatic bias and its evolution, which we do below.

The hydrostatic masses were computed as in L14. Briefly, the parametric models of Vikhlinin et al. (2006) were fitted to either the temperature and density profiles coming from the results of the spectral fitting or the mass-weighted profiles over the radial range

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Evolution of the mass slope of the total mass–hydrostatic mass scaling relations obtained when using either the density and temperature profiles derived from spectral fitting (right-hand panel) or the mass-weighted profiles (left-hand panel) to compute the hydrostatic masses. In each panel, we plot the redshift evolution of the best-fitting power-law indices obtained by fitting the broken power law given by equation (12). The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The mass-weighted hydrostatic mass tracks the absolute mass well (i.e. $M_{500, hse, mw} \propto M_{500}$), whereas the spectroscopically derived hydrostatic mass tracks the true halo mass slightly less well, with $M_{500, hse, spec} \propto M_{500}^{0.9}$ typically.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Evolution of the mass slope of the total mass–hydrostatic mass scaling relations obtained when using either the density and temperature profiles derived from spectral fitting (right-hand subpanels) or the mass-weighted profiles (left-hand subpanels) to compute the hydrostatic masses. In each subpanel, we plot the redshift evolution of the low-mass (top panel) and high-mass (bottom panel) best-fitting power-law indices obtained by fitting the broken power law given by equation (13). The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The mass-weighted hydrostatic mass tracks the absolute mass well (i.e. $M_{500, hse, mw} \propto M_{500}$), whereas the spectroscopically derived hydrostatic mass tracks the true halo mass less well, with the exact scaling depending on the included gas physics and halo mass range.</figcaption></figure>

$[0.1–1.1]r_{500}$ and hydrostatic equilibrium was assumed to compute the mass profiles.

In Figs 8 and 9, we show the evolution of the mass slopes from $z=0$ to $1.5$ for the total mass–hydrostatic mass scaling relations for each of the four physical models. In Fig. 8, we plot the redshift evolution of the best-fitting power-law indices obtained by fitting
the power law given by equation (12), whereas in Fig. 9, we plot the redshift evolutions of the low-mass (top panel) and high-mass (bottom panel) best-fitting power-law indices obtained by fitting the broken power law given by equation (13). The solid curves (red, orange, blue and green) represent the various simulations. The horizontal dashed lines correspond to the self-similar expectation (i.e. hydrostatic mass tracks the true mass). We consider two varieties of hydrostatic mass: one where we use the mass-weighted density and temperature profiles directly from the simulations (left-hand panels) and one where we use the spectroscopically derived density and temperature profiles in the hydrostatic equation (right-hand panels). The former is a better indicator of the true degree of hydrostatic bias, whereas the latter is also sensitive to observational biases in the derived gas density and temperature profiles (e.g. due to gas clumping and multiphase structure).

When the hydrostatic masses are computed using the mass-weighted profiles, the mass slopes are close to the self-similar expectation of 1, independent of mass, redshift and the nature of the included ICM physics. The use of the spectroscopic profiles, however, breaks the self-similarity, by making the relation shallower than the self-similar expectation of 1 in the vast majority of cases (the only exceptions being the high-mass end slope of the AGN 8.0 model for $z < 0.6$). This is generally true for all redshifts, masses and physical models. It is worth mentioning that the model that displays the largest deviation from self-similarity is the unphysical non-radiative model NOCOOL, likely due to the excessive gas clumping and high degree of multiphase structure present in this simulation. Because the gas is unable to radiate, large quantities of gas can exist at very high densities/short cooling times in the NOCOOL model. Such gas is thermally unstable when radiative cooling is switched on and it quickly cools down to very low (non-X-ray-emitting) temperatures, where it is either converted into stars or heated by the central AGN. Thus, inclusion of radiative cooling and feedback generally leads to a less clumpy X-ray emitting ICM. Additionally, the mass slope is shallower for clusters than for groups and, in the clusters' case, the slope gets steeper with increasing feedback intensity. This is probably due to the fact that clusters, being less relaxed, tend to have more mass in substructures (be clumpier; see e.g. Neto et al. 2007) and to the degree of clumpiness being dependent on feedback intensity (see e.g. Battaglia et al. 2015). Note that both are only apparent dependences that are due to the use of the spectroscopic profiles since they disappear in the left-hand panels of Fig. 9 when the mass-weighted profiles are used.

In Fig. 10, we show the evolution of the normalization of the total mass-hydrostatic equilibrium bias scaling relations obtained when using both the spectroscopic (right-hand subpanels) and mass-weighted (left-hand subpanels) density and temperature profiles to compute the hydrostatic masses. The results are shown for the four $\log_{10}[M_{500}(M_{\odot})]$ bins (denoted by solid lines of different colours).

When the hydrostatic masses are computed using the true mass-weighted profiles, the bias $1-b \equiv \frac{M_{500, hse, mw}}{M_{500}}$ is nearly independent of redshift and only weakly dependent on mass and included gas physics, with a typical value of $\approx 0.8$; i.e. $M_{500, hse, mw} \approx 0.8 M_{500}$. The use of spectroscopic profiles changes the situation, however, introducing a strong mass dependence in the bias for the NOCOOL and REF models. As already noted, this is likely due to the effects of gas clumping and the high degree of multiphase temperature structure present in these models, particularly the non-radiative one. Inclusion of energetic AGN feedback reduces the degree of gas clumping, which removes some (but not all) of the mass dependence of the spectroscopic 'hydrostatic' bias parameter. Typically, the spectroscopic bias parameter ranges between 0.7 to 1.0 de-

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Evolution of the normalization of the total mass-hydrostatic equilibrium bias scaling relations obtained when using either the density and temperature profiles derived from spectral fitting (right-hand subpanels) or the mass-weighted profiles (left-hand subpanels) to compute the hydrostatic masses. The four $\log_{10}[M_{500}(M_{\odot})]$ bins are denoted by solid lines of different colours and the different panels correspond to the different physical models. A mass and non-gravitational physics dependence is introduced when the profiles coming from the spectral fitting are used. Galaxy groups and clusters have different sensitivity to galaxy formation physics: their hydrostatic equilibrium biases vary in opposite directions with increased feedback intensity.</figcaption></figure>

pending on the exact mass range under consideration for the AGN feedback models, and is approximately independent of redshift. This is similar to the range ($1-b \approx 0.7-0.95$) suggested by recent measurements of the scaling relation between weak lensing (see e.g. von der Linden et al. 2014; Hoekstra et al. 2015; Penna-Lima et al. 2016; Smith et al. 2016) or CMB lensing (see e.g. Melin & Bartlett 2015) and hydrostatic masses.
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Evolution of the log-normal scatter about the total mass–hydrostatic mass scaling relations obtained when using either the density and temperature profiles coming from the results of the spectral fitting (bottom subpanels) or the mass-weighted profiles (top subpanels) to compute the hydrostatic masses. For each simulation and each scaling relation, we plot the log-normal scatter as a function of $M_{500}$ and denote the redshift using lines of different colours. For both types of hydrostatic masses, the log-normal scatter varies only mildly with mass, but displays slightly stronger redshift (it tends to increase with increasing redshift contrarily to most of the hot gas scalings; see Fig. 6) and gravitational physics dependences. About 0.04 dex additional scatter is introduced in the relation by the use of spectral fitting to obtain the density and temperature profiles.</figcaption></figure>
Finally, in Fig. 11, we show the evolution of the log-normal scatter (at fixed total mass) about the total mass–hydrostatic mass scaling relations obtained when using either the spectroscopic (bottom sub-panels) or mass-weighted (top subpanels) density and temperature profiles to compute the hydrostatic masses. For each simulation and each scaling relation, we plot the log-normal scatter as a function of $M_{500}$ and denote the redshift using lines of different colours.

For both types of hydrostatic masses, the log-normal scatter varies only mildly with mass, but displays a slightly stronger redshift dependence, tending to increase with increasing redshift. This trend is contrary to that of the scaling relations we investigated in Section 7, which show decreasing scatter with increasing redshift. Only a modest amount of additional scatter ($\approx0.04$ dex, typically) is introduced when we use the spectroscopically derived bias as opposed to the mass-weighted hydrostatic bias.

It is noteworthy that, on average, the hydrostatic bias has a larger overall scatter than the other mass proxies (apart from the X-ray luminosity) we have considered (see Table 4) and becomes increasingly less competitive with the other proxies with increasing redshift. However, together with the core-excised X-ray temperature, it has the smallest zero-point uncertainty due to the uncertain non-gravitational physics (see Table 4 and discussion in Section 10).

All these results are robust to changing the cosmology from Planck to WMAP7.

## 9 COMPARISON WITH PREVIOUS WORK

Hot gas scaling relations and their (potential) deviations from self-similar expectations have been explored in a number of previous studies using cosmological hydrodynamical simulations, some of which also included AGN feedback (or preheating; Short et al. 2010; Stanek et al. 2010; Battaglia et al. 2010, 2012; Fabjan et al. 2011; Planelles et al. 2013; Pike et al. 2014). Here we compare our findings with previous studies.

In terms of the mass slopes, Stanek et al. (2010), Fabjan et al. (2011), Biffi et al. (2014) and Pike et al. (2014) all found that the slope of the mass–temperature scaling is shallower than the self-similar expectation, in agreement with the results of this study (but see Kravtsov et al. 2006; Nagai, Kravtsov & Vikhlinin 2007b who found no deviation from self-similarity). It is worth noting though that Pike et al. found that the slope decreases with increasing redshift whereas we found a steepening and that Biffi et al. also concluded that some of the deviation is due to X-ray spectral fitting. Kay et al. (2012) and Battaglia et al. (2012) similarly saw a steepening of the SZ flux–mass relation when AGN feedback is included (but see Pike et al. 2014), whereas Nagai (2006) concluded that galaxy formation physics has a small impact by comparing a non-radiative zoom simulation with one that included cooling and star formation but no AGN feedback. Our results on the SZ flux–mass relation are consistent with these studies, in that we find that inclusion of cooling, star formation and stellar feedback (as in the REF model) does not significantly affect the slope of the relation, but the further inclusion of energetic AGN feedback, which lowers the gas fractions of groups, steepens the relation.

Short et al. (2010), Stanek et al. (2010) and Fabjan et al. (2011), in agreement with the results discussed here, concluded that the mass–temperature relation evolves negatively (but see Nagai et al. 2007b). Stanek et al. similarly attributed it to the fact that the kinetic-to-thermal ratio increases with redshift. In terms of the gas mass–total mass scaling, Planelles et al. (2013) found that it evolves self-similarly even when AGN feedback is included and Fabjan et al. (2011) found a mild positive evolution even in the adiabatic case and

in a model similar in spirit to REF (see also Kravtsov et al. 2006), both in clear disagreement with this study. Note that for Planelles et al., this is due to the fact that the baryon fraction nearly tracks the universal baryon fraction $\Omega_b/\Omega_m$ even in their simulations with AGN feedback (i.e. the feedback was not sufficiently energetic to lower the gas fractions of groups). Short et al. (2010) and Battaglia et al. (2012) reported noticeable deviations from self-similarity in the evolution of $Y_X$ and SZ flux, respectively, in agreement with the results presented here, whereas Nagai (2006), Kay et al. (2012), Pike et al. (2014), Kravtsov et al. (2006), Nagai et al. (2007b) and Fabjan et al. (2011) saw no deviations from the self-similar expectation for the normalization of the SZ flux (for the first three references) and of $Y_X$ (for the last three) even when non-gravitational physics is included. We again attribute these differences to inefficient feedback in these studies, that we argue, in McCarthy et al. (2011) and L14, is necessary to reproduce the observed low baryon fractions of groups.

Nagai et al. (2007a), in reasonable agreement with this study, found a 5–20 per cent hydrostatic bias (see also e.g. Meneghetti et al. 2010; Battaglia et al. 2012; Rasia et al. 2012). Additionally, Kay et al. (2012) reported a 20–30 per cent hydrostatic bias for the non-radiative runs and an $\sim10$ per cent decrease when non-gravitational physics is included, in agreement with the results presented here.

## 10 SUMMARY AND DISCUSSION

We have used the cosmo-OWLS suite of large-volume hydrodynamical simulations (Le Brun et al. 2014; McCarthy et al. 2014) to investigate the scatter and evolution of the global hot gas properties of large populations of simulated galaxy groups and clusters. The simulations use $2 \times 10^{24}$ particles in a 400 $h^{-1}$ Mpc on a side box, and assume either a WMAP7-year or a Planck 2013 cosmology. We have studied all the haloes with $M_{500} \ge 10^{13} M_{\odot}$ over the redshift range $z = 0–1.5$. For instance, there are over 25 000 such systems in the Planck cosmology version of the NOCOOL model. The most realistic models, which include efficient feedback from AGN, have been shown by Le Brun et al. to reproduce a wide range of observable properties of local systems. Note that cosmo-OWLS forms an extension of the OWLS project (Schaye et al. 2010) and was designed to help quantify the importance of uncertainties in ‘sub-grid’ physics for cluster cosmology efforts.

From the study presented here, we conclude the following.

(i) The median relations between (true) halo mass and the (core-excised) mass-weighted and X-ray spectroscopic temperature, bolo-metric X-ray luminosity, gas mass, SZ flux, its X-ray analogue and the hydrostatic mass and the scatter about them are reasonably well modelled by evolving broken power laws (equation 15), whereas single power laws of the form given by equation (12) perform less well (Fig. 1). In the AGN models, the physical origin of this result is tied to the break of the self-similarity in the relation between gas mass fraction and halo mass (which is also observed): massive clusters have approximately constant gas fractions near the universal value of $\Omega_b/\Omega_m$, while for lower mass systems ($M_{500} \lesssim 10^{14} M_{\odot}$), the gas fraction is a steadily declining function of halo mass.

(ii) In terms of the mass slopes (i.e. the logarithmic slopes of the various mass–observable relations), we find that our most realistic (AGN) models predict large deviations from the self-similar expectation for all of the scaling relations we have examined. The one exception is the mass–temperature relation, where only a weak, though still significant, deviation from self-similarity is found (independently of the details of the included ICM physics). For the
other scaling relations, all of which depend directly on the gas density/mass, the deviation from the self-similar prediction is strongest at low redshifts and low halo masses. Models without efficient feedback (NOCOOL and REF), on the other hand, have mass slopes that are approximately consistent with self-similar expectations (Figs 2–3).

(iii) In terms of the evolution of the normalization of the scaling relations, we find that the mass–temperature relation evolves in a negative fashion with respect to the self-similar model (i.e. slower than predicted), independently of the included ICM physics. This is likely due to the increasing contribution of non-thermal pressure support with increasing redshift (the kinetic-to-thermal energy ratio of the ICM increases from a typical value of 10–15 per cent at $z=0$ up to 20–30 per cent by $z=1$; see Fig. 5). The gas mass evolves approximately self-similarly (i.e. remains constant) for the non-radiative simulation, but shows a positive evolution when radiative cooling and particularly AGN feedback is included, a result which is strongly mass dependent. We hypothesize that, in the case of the AGN models, the positive evolution is due to the increasing binding energy of haloes of fixed total mass with increasing redshift (making gas expulsion more difficult at high redshift). The SZ flux and $Y_X$ evolve negatively with respect to the self-similar expectation for models without efficient feedback, driven by the negative evolution of the temperature. When AGN feedback is included, however, the sign of the evolution of the SZ flux and $Y_X$ with respect to self-similar depends on halo mass (positive evolution for low-mass haloes, negative evolution for high-mass haloes), driven by the strong halo mass dependence of the gas mass evolution combined with the negative evolution of the temperature (Fig. 4).

(iv) The scatter about the various mass–observable relations is mass dependent, varying by a factor of 2–3 over 1.5 decades in $M_{500}$ (Fig. 7). The overall amplitude of the scatter at fixed redshift is somewhat sensitive to the included non-gravitational physics, with the amplitude tending to increase with increasing complexity (and realism) of the included galaxy formation physics and with increasing AGN feedback intensity. Encouragingly, the overall amplitude of the scatter tends to decrease with increasing redshift, which is likely due to the lessening importance of non-gravitational physics at higher redshifts when examining the evolution of haloes of fixed mass. Typically, the scatter is reduced by ~50 per cent from $z=0$ to 1 (Fig. 6). At $z=0$, the X-ray temperature, gas mass, SZ signal and $Y_X$ all have similar amounts of scatter at fixed total mass (typically 10 per cent), whereas the bolometric X-ray luminosity has considerably larger scatter (Table 4).

(v) We have also separately analysed the scatter and evolution in the hydrostatic mass–true mass relation (see Section 8 and Figs 8–11). We find that, overall, the hydrostatic mass tracks the true mass reasonably well (such that $M_{500, hse} \propto M_{500}$). The typical value of the true bias is $1-b \equiv \frac{M_{500, hse, mw}}{M_{500}} \approx 0.8$ and is insensitive to redshift, true halo mass, and included ICM physics. When we derive hydrostatic masses using synthetic X-ray spectroscopy to infer the density and temperature profiles, however, the bias becomes more dependent on the true halo mass range and the included ICM physics, with typical median values ranging between $1-b \equiv \frac{M_{500, hse, spec}}{M_{500}} \approx 0.7$ and 1.0 with little evolution. At $z=0$, the typical level of scatter in the hydrostatic mass is 15 per cent (i.e. slightly larger than for most of the other proxies we have considered, barring the bolometric X-ray luminosity). However, contrary to what was found for the other proxies, the scatter in the hydrostatic mass *increases* with increasing redshift, making it an increasingly less competitive proxy at higher redshifts.

(vi) The bolometric X-ray luminosity–X-ray temperature relation qualitatively behaves in approximately the same way as the bolometric X-ray luminosity–mass relation (compare Figs 2–6 to Figs C1–C4). The differences can be explained by the behaviour of the mass–temperature relation.

(vii) The vast majority of these results are robust to changing the cosmology from Planck to WMAP7. The only noteworthy exceptions are the high-mass ends of all the studied mass–observable relations and the bolometric X-ray luminosity–temperature relation, which evolve slightly faster in the WMAP7 cosmology.

Our results have highlighted that, while some of the trends are sensitive to the sub-grid modelling, there are many robust trends that are directly relevant for ongoing and upcoming surveys. For example, for samples including both galaxy groups and clusters, none of the mass–observable relations are well-fit by a single power law for any of the physical models (apart from the unphysical non-radiative model). We advocate the use of broken power laws as an alternative. In terms of amplitude evolution, with the exception of the gas mass–halo mass relation restricted to high halo masses, none of the scaling relations we have examined evolves self-similarly in any of the models: we therefore warn against the ‘simplistic’ interpretation of deviations of any scaling relation from self-similar evolution as an indication of the effects of feedback processes. The precise evolution generally depends on the included ICM physics, for which we advocate using the predictions of the AGN models (which at least reproduce the $z=0$ cluster population well). In addition, we find that the scatter about the relations generally decreases with increasing redshift, independently of the included ICM physics. The one important exception is the hydrostatic mass, whose performance as a mass proxy worsens (scatter increases) with increasing redshift.

Our detailed analysis of the various mass–observable relations also allows us to comment on what is the ‘best’ total mass proxy. We argue that an ideal mass proxy should: (i) be easy to measure (including with relatively shallow data); (ii) be strongly correlated with mass; (iii) have a small intrinsic scatter; (iv) be insensitive to the choice of cosmology; (v) be insensitive to the cluster dynamical state; (vi) be insensitive to uncertain baryon physics (i.e. the zero-point is well-known); and (vii) evolve in a manner that is easy to characterize. On the basis of the above criteria, we argue that the X-ray temperature performs the best. Its main strength relative to the other mass proxies is that it is least sensitive to uncertain baryon physics (in terms of its zero-point, evolution and intrinsic scatter) while having amongst the smallest intrinsic scatter in mass at fixed proxy (see Table 4).

If one's sole criterion for judging the best mass proxy is instead based on the intrinsic scatter (for example, because one plans to use a high-quality sub-sample with mass proxy information in combination with 'self-calibration' to remove any zero-point uncertainties; e.g. Levine, Schulz & White 2002; Hu 2003; Majumdar & Mohr 2003, 2004; Mantz et al. 2010a; Wu, Rozo & Wechsler 2010), one might instead conclude from inspection of Table 4 that the integrated SZ signal, $Y_X$ and the gas mass have a slight advantage over X-ray temperature (see Allen et al. 2011, and references therein). We caution, however, that our predictions for the intrinsic scatter of these mass proxies (and their dependences on halo mass and redshift) are sensitive to unresolved feedback physics, whereas the scatter in the mass–temperature relation is generally insensitive to baryonic physics (see Fig. 7).

In this study, we have examined the scaling relations between various global hot gas properties and halo mass. In a future study,
we plan to examine the evolution and scatter in the structure of the hot gas (i.e. profiles).

## ACKNOWLEDGEMENTS

The authors would like to thank the members of the OWLS team for their contributions to the development of the simulation code used here and the anonymous referee for a constructive report. AMCLB is grateful to Monique Arnaud, Andrey Kravtsov, Daisuke Nagai and Gabriel Pratt for constructive discussions and to Remco van der Burg for useful comments on an earlier version of the manuscript. AMCLB acknowledges support from an internally funded PhD studentship at the Astrophysics Research Institute of Liverpool John Moores University, from the French Agence Nationale de la Recherche under grant ANR-11-BS56-015 and from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 340519. IGM is supported by an STFC Advanced Fellowship. JS is sponsored by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 278594-GasAroundGalaxies. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

## REFERENCES

* Allen S. W., Rapetti D. A., Schmidt R. W., Ebeling H., Morris R. G., Fabian A. C., 2008, MNRAS, 383, 879

* Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409

* Andersson K. et al., 2011, ApJ, 738, 48

* Arnaud M., Evrard A. E., 1999, MNRAS, 305, 631

* Barnes D. J., Kay S. T., Henson M. A., McCarthy I. G., Schaye J., Jenkins A., 2017, MNRAS, 465, 213

* Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., Sijacki D., 2010, ApJ, 725, 91

* Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 74

* Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2015, ApJ, 806, 43

* Benson B. A. et al., 2014, in Holland W. S., Zmuidzinas J., eds, Proc. SPIE Conf. Ser. Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII. SPIE, Bellingham, p. 91531P

* Bhattacharya S., Di Matteo T., Kosowsky A., 2008, MNRAS, 389, 34

* Biffi V., Sembolini F., De Petris M., Valdarnini R., Yepes G., Gottlöber S., 2014, MNRAS, 439, 588

* Bode P., Ostriker J. P., Weller J., Shaw L., 2007, ApJ, 663, 139

* Booth C. M., Schaye J., 2009, MNRAS, 398, 53

* Booth C. M., Schaye J., 2010, MNRAS, 405, L1

* Booth C. M., Schaye J., 2011, MNRAS, 413, 1158

* Caldwell C. E., McCarthy I. G., Baldry I. K., Collins C. A., Schaye J., Bird S., 2016, MNRAS, 462, 4117

* Cui W., Borgani S., Murante G., 2014, MNRAS, 441, 1769

* Cusworth S. J., Kay S. T., Battye R. A., Thomas P. A., 2014, MNRAS, 439, 2485

* Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431

* Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497

* Edge A. C., Stewart G. C., 1991, MNRAS, 252, 414

* Ettori S. et al., 2004, MNRAS, 354, 111

* Evrard A. E., Henry J. P., 1991, ApJ, 383, 95

* Evrard A. E., Metzler C. A., Navarro J. F., 1996, ApJ, 469, 494

* Fabjan D., Borgani S., Tornatore L., Saro A., Murante G., Dolag K., 2010, MNRAS, 401, 1670

* Fabjan D., Borgani S., Rasia E., Bonafede A., Dolag K., Murante G., Tornatore L., 2011, MNRAS, 416, 801

* Greco J. P., Hill J. C., Spergel D. N., Battaglia N., 2015, ApJ, 808, 151

* Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays (Available at http://moriond.in2p3.fr)

* Henderson S. W. et al., 2016, J. Low Temp. Phys., 184, 772

* Hilton M. et al., 2012, MNRAS, 424, 2086

* Hoekstra H., Herbonnet R., Muzzin A., Babul A., Mahdavi A., Viola M., Cacciato M., 2015, MNRAS, 449, 685

* Hojjati A., McCarthy I. G., Harnois-Deraps J., Ma Y.-Z., Van Waerbeke L., Hinshaw G., Le Brun A. M. C., 2015, JCAP, 10, 047

* Hu W., 2003, Phys. Rev. D, 67, 081304

* Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS, 321, 372

* Kaiser N., 1986, MNRAS, 222, 323

* Kaiser N., 1991, ApJ, 383, 104

* Kay S. T., da Silva A. C., Aghanim N., Blanchard A., Liddle A. R., Puget J.-L., Sadat R., Thomas P. A., 2007, MNRAS, 377, 317

* Kay S. T., Peel M. W., Short C. J., Thomas P. A., Young O. E., Battye R. A., Liddle A. R., Pearce F. R., 2012, MNRAS, 422, 1999

* Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353

* Kravtsov A. V., Vikhlinin A., Nagai D., 2006, ApJ, 650, 128

* Laureijs R. et al., 2011, ESA report (arXiv:1110.3193)

* Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270 (L14)

* Le Brun A. M. C., McCarthy I. G., Melin J.-B., 2015, MNRAS, 451, 3868

* Levine E. S., Schulz A. E., White M., 2002, ApJ, 577, 569

* Lin Y.-T., Stanford S. A., Eisenhardt P. R. M., Vikhlinin A., Maughan B. J., Kravtsov A., 2012, ApJ, 745, L3

* Lumb D. H. et al., 2004, A&A, 420, 853

* McCarthy I. G. et al., 2010, MNRAS, 406, 822

* McCarthy I. G., Schaye J., Bower R. G., Ponman T. J., Booth C. M., Dalla Vecchia C., Springel V., 2011, MNRAS, 412, 1965

* McCarthy I. G., Le Brun A. M. C., Schaye J., Holder G. P., 2014, MNRAS, 440, 3645

* McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2016, MNRAS, 465, 2936

* Majumdar S., Mohr J. J., 2003, ApJ, 585, 603

* Majumdar S., Mohr J. J., 2004, ApJ, 613, 41

* Mantz A., Allen S. W., Rapetti D., Ebeling H., 2010a, MNRAS, 406, 1759

* Mantz A., Allen S. W., Ebeling H., Rapetti D., Drlica-Wagner A., 2010b, MNRAS, 406, 1773

* Mantz A. B., Allen S. W., Morris R. G., Rapetti D. A., Applegate D. E., Kelly P. L., von der Linden A., Schmidt R. W., 2014, MNRAS, 440, 2077

* Mantz A. B. et al., 2015, MNRAS, 446, 2205

* Markevitch M., 1998, ApJ, 504, 27

* Markwardt C. B., 2009, in Bohlender D. A., Durand D., Dowler P., eds, ASP Conf. Ser. Vol. 411, Astronomical Data Analysis Software and Systems XVIII. Astron. Soc. Pac., San Francisco p. 251

* Martizzi D., Mohammed I., Teyssier R., Moore B., 2014, MNRAS, 440, 2290

* Maughan B. J., Jones L. R., Ebeling H., Scharf C., 2006, MNRAS, 365, 509

* Melin J.-B., Bartlett J. G., 2015, A&A, 578, A21

* Meneghetti M., Rasia E., Merten J., Bellagamba F., Ettori S., Mazzotta P., Dolag K., Marri S., 2010, A&A, 514, A93

* Merloni A. et al., 2012, preprint (arXiv:1209.3114)

* Mushotzky R. F., 1984, Phys. Scr. T, 7, 157

* Nagai D., 2006, ApJ, 650, 538

* Nagai D., Vikhlinin A., Kravtsov A. V., 2007a, ApJ, 655, 98

* Nagai D., Kravtsov A. V., Vikhlinin A., 2007b, ApJ, 668, 1

* Neto A. F. et al., 2007, MNRAS, 381, 1450

* Ntampaka M., Trac H., Cisewski J., Price L. C., 2016, ApJ, preprint (arXiv:1602.01837)

* Osmond J. P. F., Ponman T. J., 2004, MNRAS, 350, 1511
Peng C. Y., 2007, ApJ, 671, 1098

Penna-Lima M., Bartlett J. G., Rozo E., Melin J.-B., Merten J., Evrard
A. E., Postman M., Rykoff E., 2016, preprint (arXiv: 1608.05356)

Piffaretti R., Valdarnini R., 2008, A&A, 491, 71

Pike S. R., Kay S. T., Newton R. D. A., Thomas P. A., Jenkins A., 2014, MNRAS, 445, 1774

Planck Collaboration XI, 2011, A&A, 536, A11

Planck Collaboration III, 2013, A&A, 550, A129

Planck Collaboration XI, 2013, A&A, 557, A52

Planck Collaboration XVI, 2013, A&A, 571, A16

Planck Collaboration XX, 2013, A&A, 571, A20

Planck Collaboration XXIV, 2015, A&A, 594, A24

Panelles S., Borgani S., Dolag K., Ettori S., Fabjan D., Murante G., Tornatore L., 2013, MNRAS, 431, 1487

Panelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., Ragone-Figueroa C., Dolag K., 2014, MNRAS, 438, 195

Pratt G. W., Croston J. H., Arnaud M., Böhringer H., 2009, A&A, 498, 361

Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53

Rasia E. et al., 2012, New J. Phys., 14, 055018

Rasia E. et al., 2014, ApJ, 791, 96

Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210

Schaye J. et al., 2010, MNRAS, 402, 1536

Sehgal N., Trac H., Huffenberger K., Bode P., 2007, ApJ, 664, 149

Sehgal N., Bode P., Das S., Hernandez-Monteagudo C., Huffenberger K., Lin Y.-T., Ostriker J. P., Trac H., 2010, ApJ, 709, 920

Sehgal N. et al., 2011, ApJ, 732, 44

Sembolini F. et al., 2016, MNRAS, 459, 2973

Short C. J., Thomas P. A., Young O. E., Pearce F. R., Jenkins A., Muanwong
O., 2010, MNRAS, 408, 2213

Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877

Smith G. P. et al., 2016, MNRAS, 456, L74

Springel V., 2005, MNRAS, 364, 1105

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS,
328, 726

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Stanek R., Rasia E., Evrard A. E., Pearce F., Gazzola L., 2010, ApJ, 715, 1508

The Dark Energy Survey Collaboration, 2005, preprint (astro-ph/0510346)

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

van Daalen M. P., Schaye J., McCarthy I. G., Booth C. M., Dalla Vecchia
C., 2014, MNRAS, 440, 2997

Velliscig M., van Daalen M. P., Schaye J., McCarthy I. G., Cacciato M., Le Brun A. M. C., Dalla Vecchia C., 2014, MNRAS, 442, 2641

Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray
S. S., Van Speybroeck L., 2006, ApJ, 640, 691

Vikhlinin A. et al., 2009a, ApJ, 692, 1033

Vikhlinin A. et al., 2009b, ApJ, 692, 1060

Voit G. M., 2005, Rev. Mod. Phys., 77, 207

von der Linden A. et al., 2014, MNRAS, 443, 1973

Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 393, 99

Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b, MNRAS, 399, 574

Wu H.-Y., Rozo E., Wechsler R. H., 2010, ApJ, 713, 1207

APPENDIX A: RESULTS FOR THE
MASS-NON-CORE-EXCISED TEMPERATURES
SCALING RELATIONS

As it is sometimes impossible to compute core-excised temperatures for observed groups and clusters (especially for instance at high redshift or for short exposure times), we also present the results for the mass-temperature relation in the case of non-core excised tempera-

<figure><img src="image_17.png" /></figure>

**Figure A1.** Evolution of the mass slope from *z* = 0 to 1.5 for the total mass–temperature scaling relations (for both non-core-excised mass-weighted and X-ray spectroscopic temperatures). In each panel, we plot the redshift evolution of the best-fitting power-law indices obtained by fitting the broken power law given by equation (12) at each redshift independently. The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. Whether one excises the core or not when computing the temperature does not noticeably affect the results presented in the main part of the paper (compare with Fig. 2).

<figure><img src="image_18.png" /></figure>

**Figure A2.** Evolution of the mass slope from *z* = 0 to 1.5 for the total mass–temperature scaling relations (for both non-core-excised mass-weighted and X-ray spectroscopic temperatures). In each subpanel, we plot the redshift evolution of the low-mass (top panel) and high-mass (bottom panel) best-fitting power-law indices obtained by fitting the broken power law given by equation (13) at each individual redshift. The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The high-mass end of the mass–X-ray temperature relation of the AGN 8.5 model is steeper in the absence of core excision (compare with Fig. 3).

tures (be it mass-weighted or X-ray) and conclude that whether one
excises the core or not when computing the temperature does not
noticeably affect the results presented in the main part of the pa-
per, with a few noteworthy exceptions: (i) the high-mass end of the
mass–X-ray temperature relation of the AGN 8.5 model is steeper in
the absence of core excision and it even goes from being shallower
than the self-similar expectation when the core is excised to being
steeper (compare the left-hand panel of Fig. 3 to the bottom panel
of Fig. A2); (ii) the evolution of the normalization being slightly
more negative for the X-ray temperature when the core is excised
(compare the top right subpanels of Fig. 4 to the panels of Fig. A3)
<figure><img src="image_19.png" /><figcaption>**Figure A3.** Evolution of the normalization from $z=0$ to 1.5 for the total mass–temperature scaling relations (for both non-core-excised mass-weighted and X-ray spectroscopic temperatures) for each simulation. The normalizations of each scaling relation in the four $\log_{10}[M_{500}(M_{\odot})]$ bins (denoted by solid lines of different colours) have been normalized by the self-similar expectation for the redshift evolution at fixed mass (shown as an horizontal dashed line). The evolution of the normalization is slightly more negative for the X-ray temperature when the core is excised (compare with Fig. 4). The fact that, for the highest mass bin, the mass–temperature scalings have a slightly faster evolution in the *WMAP7* cosmology than in the *Planck* one is also unaffected by core excision.</figcaption></figure>

and feedback is included (be it from SNe or AGN); and (iii) as might be expected, a global increase (at all masses and redshift) of the scatter for all the physical models (compare the top panels of Fig. 6 to the panels of Fig. A4) in the absence of core excision. Finally, it is worth noting that the fact that, for the highest-mass bin, the finding that the mass–temperature scalings have a slightly faster evolution in the *WMAP7* cosmology than in the *Planck* one is also unaffected by core excision.

<figure><img src="image_20.png" /><figcaption>**Figure A4.** Evolution of the log-normal scatter from $z=0$ to 1.5 for the total mass–non-core-excised X-ray spectroscopic temperature scaling relation. For each simulation, we plot the log-normal scatter as a function of $M_{500}$ and denote the redshift using lines of different colours. There is a global increase (at all masses and redshift) of the scatter for all the physical models (compare to the top panels of Fig. 6) in the absence of core excision.</figcaption></figure>

## APPENDIX B: RESULTS FOR THE OTHER PHYSICAL MODELS

In Tables B1–B6, we present the results of the fitting procedures described in Section 4 for three of the other cosmo-OWLS physical models in *Planck* cosmology.
<table><caption>**Table B1.** Results of fitting the evolving power law (equation 12) to both the median relation and the log-normal scatter about it for the NOCOOL simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Median</td><td>0.173 ± 0.006</td><td>0.377 ± 0.032</td><td>0.547 ± 0.007</td><td>0.074</td><td>30</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Median</td><td>44.129 ± 0.009</td><td>2.267 ± 0.051</td><td>1.411 ± 0.012</td><td>0.190</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Median</td><td>13.136 ± 0.002</td><td>0.066 ± 0.011</td><td>1.020 ± 0.002</td><td>0.008</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Median</td><td>13.311 ± 0.004</td><td>0.471 ± 0.021</td><td>1.569 ± 0.005</td><td>0.033</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Median</td><td>−5.573 ± 0.003</td><td>0.647 ± 0.017</td><td>1.645 ± 0.004</td><td>0.021</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Median</td><td>13.904 ± 0.008</td><td>−0.006 ± 0.041</td><td>0.849 ± 0.010</td><td>0.124</td><td>30</td></tr><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Scatter</td><td>−0.983 ± 0.008</td><td>0.303 ± 0.045</td><td>0.280 ± 0.010</td><td>0.147</td><td>30</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Scatter</td><td>−0.381 ± 0.011</td><td>−0.325 ± 0.060</td><td>0.025 ± 0.014</td><td>0.259</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−1.381 ± 0.026</td><td>−0.360 ± 0.139</td><td>−0.124 ± 0.032</td><td>1.400</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−0.863 ± 0.012</td><td>0.111 ± 0.064</td><td>0.192 ± 0.015</td><td>0.301</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Scatter</td><td>−0.930 ± 0.016</td><td>0.010 ± 0.089</td><td>0.035 ± 0.021</td><td>0.581</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Scatter</td><td>−0.657 ± 0.012</td><td>0.607 ± 0.065</td><td>0.165 ± 0.015</td><td>0.307</td><td>30</td></tr></tbody></table>

<table><caption>**Table B2.** Results of fitting the evolving broken power law (equation 15) to both the median relation and the log-normal scatter about it for the NOCOOL simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A''</th><th>α''</th><th>β''</th><th>γ''</th><th>δ''</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Median</td><td>0.200 ± 0.010</td><td>0.290 ± 0.042</td><td>0.656 ± 0.037</td><td>0.499 ± 0.019</td><td>−0.056 ± 0.023</td><td>0.054</td><td>28</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Median</td><td>44.178 ± 0.009</td><td>2.260 ± 0.039</td><td>1.457 ± 0.035</td><td>1.269 ± 0.018</td><td>0.026 ± 0.021</td><td>0.047</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Median</td><td>13.143 ± 0.003</td><td>0.038 ± 0.015</td><td>1.054 ± 0.013</td><td>1.008 ± 0.007</td><td>−0.019 ± 0.008</td><td>0.007</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Median</td><td>13.328 ± 0.007</td><td>0.406 ± 0.028</td><td>1.648 ± 0.025</td><td>1.543 ± 0.013</td><td>−0.044 ± 0.015</td><td>0.024</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Median</td><td>−5.577 ± 0.006</td><td>0.633 ± 0.024</td><td>1.654 ± 0.021</td><td>1.663 ± 0.011</td><td>−0.014 ± 0.013</td><td>0.017</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Median</td><td>13.951 ± 0.025</td><td>−0.156 ± 0.046</td><td>1.035 ± 0.041</td><td>0.765 ± 0.021</td><td>−0.097 ± 0.025</td><td>0.064</td><td>28</td></tr><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Scatter</td><td>−1.011 ± 0.015</td><td>0.393 ± 0.064</td><td>0.168 ± 0.057</td><td>0.330 ± 0.029</td><td>0.058 ± 0.034</td><td>0.125</td><td>28</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Scatter</td><td>−0.374 ± 0.019</td><td>−0.231 ± 0.080</td><td>−0.062 ± 0.071</td><td>−0.030 ± 0.036</td><td>0.084 ± 0.043</td><td>0.196</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−1.378 ± 0.050</td><td>−0.427 ± 0.214</td><td>−0.055 ± 0.188</td><td>−0.108 ± 0.097</td><td>−0.055 ± 0.115</td><td>1.383</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−0.943 ± 0.015</td><td>0.348 ± 0.066</td><td>−0.107 ± 0.058</td><td>0.341 ± 0.030</td><td>0.149 ± 0.035</td><td>0.131</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Scatter</td><td>−1.018 ± 0.026</td><td>0.207 ± 0.109</td><td>−0.231 ± 0.096</td><td>0.225 ± 0.049</td><td>0.109 ± 0.058</td><td>0.358</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Scatter</td><td>−0.705 ± 0.019</td><td>0.629 ± 0.079</td><td>0.105 ± 0.070</td><td>0.297 ± 0.036</td><td>−0.012 ± 0.042</td><td>0.189</td><td>28</td></tr></tbody></table>

<table><caption>**Table B3.** Results of fitting the evolving power law (equation 12) to both the median relation and the log-normal scatter about it for the REF simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Median</td><td>0.314 ± 0.003</td><td>0.248 ± 0.019</td><td>0.530 ± 0.004</td><td>0.027</td><td>30</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Median</td><td>43.546 ± 0.007</td><td>2.349 ± 0.038</td><td>1.365 ± 0.009</td><td>0.106</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Median</td><td>12.917 ± 0.004</td><td>0.352 ± 0.021</td><td>1.091 ± 0.005</td><td>0.032</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Median</td><td>13.229 ± 0.006</td><td>0.588 ± 0.034</td><td>1.609 ± 0.008</td><td>0.085</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Median</td><td>−5.720 ± 0.006</td><td>0.723 ± 0.034</td><td>1.745 ± 0.008</td><td>0.082</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Median</td><td>13.912 ± 0.006</td><td>0.041 ± 0.031</td><td>0.907 ± 0.007</td><td>0.069</td><td>30</td></tr><tr><td>T<sub>spec</sub>, cor−M<sub>500</sub></td><td>Scatter</td><td>−0.950 ± 0.022</td><td>−0.285 ± 0.118</td><td>−0.025 ± 0.028</td><td>1.014</td><td>30</td></tr><tr><td>L<sub>bol</sub>−M<sub>500</sub></td><td>Scatter</td><td>−0.423 ± 0.014</td><td>−0.654 ± 0.076</td><td>−0.191 ± 0.018</td><td>0.415</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−1.069 ± 0.009</td><td>−0.805 ± 0.051</td><td>−0.137 ± 0.012</td><td>0.185</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500−M<sub>500</sub></td><td>Scatter</td><td>−0.866 ± 0.017</td><td>−0.538 ± 0.092</td><td>−0.109 ± 0.022</td><td>0.619</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>−M<sub>500</sub></td><td>Scatter</td><td>−0.937 ± 0.018</td><td>−0.290 ± 0.100</td><td>−0.150 ± 0.023</td><td>0.728</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec−M<sub>500</sub></td><td>Scatter</td><td>−0.690 ± 0.015</td><td>0.426 ± 0.082</td><td>0.080 ± 0.019</td><td>0.485</td><td>30</td></tr></tbody></table>
<table><caption>**Table B4.** Results of fitting the evolving broken power law (equation 15) to both the median relation and the log-normal scatter about it for the REF simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A''</th><th>α''</th><th>β''</th><th>γ''</th><th>δ''</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Median</td><td>0.317 ± 0.007</td><td>0.221 ± 0.029</td><td>0.559 ± 0.025</td><td>0.531 ± 0.013</td><td>-0.021 ± 0.015</td><td>0.025</td><td>28</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Median</td><td>43.505 ± 0.010</td><td>2.495 ± 0.043</td><td>1.187 ± 0.038</td><td>1.435 ± 0.019</td><td>0.097 ± 0.023</td><td>0.056</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Median</td><td>12.934 ± 0.005</td><td>0.256 ± 0.021</td><td>1.199 ± 0.019</td><td>1.077 ± 0.010</td><td>-0.070 ± 0.011</td><td>0.014</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Median</td><td>13.254 ± 0.009</td><td>0.454 ± 0.040</td><td>1.762 ± 0.035</td><td>1.585 ± 0.018</td><td>-0.097 ± 0.021</td><td>0.049</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Median</td><td>-5.684 ± 0.009</td><td>0.643 ± 0.039</td><td>1.853 ± 0.034</td><td>1.669 ± 0.018</td><td>-0.045 ± 0.021</td><td>0.047</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Median</td><td>13.935 ± 0.010</td><td>-0.052 ± 0.041</td><td>1.018 ± 0.036</td><td>0.871 ± 0.019</td><td>-0.063 ± 0.022</td><td>0.050</td><td>28</td></tr><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Scatter</td><td>-1.111 ± 0.024</td><td>0.233 ± 0.101</td><td>-0.665 ± 0.088</td><td>0.264 ± 0.046</td><td>0.332 ± 0.054</td><td>0.305</td><td>28</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Scatter</td><td>-0.467 ± 0.025</td><td>-0.482 ± 0.108</td><td>-0.397 ± 0.095</td><td>-0.120 ± 0.049</td><td>0.116 ± 0.058</td><td>0.352</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-1.048 ± 0.014</td><td>-0.987 ± 0.058</td><td>0.060 ± 0.051</td><td>-0.133 ± 0.026</td><td>-0.140 ± 0.031</td><td>0.103</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-0.936 ± 0.030</td><td>-0.331 ± 0.127</td><td>-0.369 ± 0.112</td><td>0.024 ± 0.058</td><td>0.129 ± 0.068</td><td>0.489</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Scatter</td><td>-1.011 ± 0.030</td><td>-0.232 ± 0.126</td><td>-0.266 ± 0.111</td><td>0.046 ± 0.057</td><td>0.001 ± 0.067</td><td>0.480</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Scatter</td><td>-0.763 ± 0.024</td><td>0.560 ± 0.104</td><td>-0.109 ± 0.091</td><td>0.245 ± 0.047</td><td>0.066 ± 0.055</td><td>0.324</td><td>28</td></tr></tbody></table>

<table><caption>**Table B5.** Results of fitting the evolving power law (equation 12) to both the median relation and the log-normal scatter about it for the AGN 8.5 simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Median</td><td>0.265 ± 0.008</td><td>0.390 ± 0.042</td><td>0.654 ± 0.010</td><td>0.127</td><td>30</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Median</td><td>43.130 ± 0.019</td><td>3.191 ± 0.103</td><td>1.767 ± 0.024</td><td>0.770</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Median</td><td>12.684 ± 0.012</td><td>0.753 ± 0.064</td><td>1.364 ± 0.015</td><td>0.295</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Median</td><td>12.951 ± 0.018</td><td>1.132 ± 0.096</td><td>2.010 ± 0.022</td><td>0.667</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Median</td><td>-5.905 ± 0.011</td><td>1.166 ± 0.061</td><td>2.055 ± 0.014</td><td>0.269</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Median</td><td>13.874 ± 0.013</td><td>-0.043 ± 0.068</td><td>0.964 ± 0.016</td><td>0.340</td><td>30</td></tr><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Scatter</td><td>-0.951 ± 0.032</td><td>-0.323 ± 0.175</td><td>-0.310 ± 0.041</td><td>2.232</td><td>30</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Scatter</td><td>-0.470 ± 0.025</td><td>-0.321 ± 0.136</td><td>-0.268 ± 0.032</td><td>1.351</td><td>30</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-0.875 ± 0.027</td><td>-0.489 ± 0.146</td><td>-0.457 ± 0.034</td><td>1.558</td><td>30</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-0.661 ± 0.020</td><td>-0.309 ± 0.109</td><td>-0.289 ± 0.025</td><td>0.862</td><td>30</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Scatter</td><td>-0.764 ± 0.020</td><td>-0.145 ± 0.107</td><td>-0.328 ± 0.025</td><td>0.830</td><td>30</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Scatter</td><td>-0.672 ± 0.056</td><td>0.469 ± 0.302</td><td>-0.201 ± 0.070</td><td>6.646</td><td>30</td></tr></tbody></table>

<table><caption>**Table B6.** Results of fitting the evolving broken power law (equation 15) to both the median relation and the log-normal scatter about it for the AGN 8.5 simulation. The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Median or scatter</th><th>A''</th><th>α''</th><th>β''</th><th>γ''</th><th>δ''</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Median</td><td>0.298 ± 0.013</td><td>0.254 ± 0.053</td><td>0.816 ± 0.047</td><td>0.604 ± 0.025</td><td>-0.094 ± 0.029</td><td>0.088</td><td>28</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Median</td><td>43.062 ± 0.022</td><td>3.084 ± 0.092</td><td>1.816 ± 0.080</td><td>2.005 ± 0.042</td><td>-0.130 ± 0.049</td><td>0.258</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Median</td><td>12.713 ± 0.015</td><td>0.489 ± 0.062</td><td>1.649 ± 0.054</td><td>1.374 ± 0.028</td><td>-0.204 ± 0.033</td><td>0.117</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Median</td><td>13.004 ± 0.025</td><td>0.752 ± 0.106</td><td>2.429 ± 0.092</td><td>1.991 ± 0.048</td><td>-0.286 ± 0.056</td><td>0.342</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Median</td><td>-5.827 ± 0.012</td><td>0.887 ± 0.053</td><td>2.395 ± 0.047</td><td>1.927 ± 0.024</td><td>-0.187 ± 0.028</td><td>0.087</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Median</td><td>13.833 ± 0.023</td><td>0.013 ± 0.096</td><td>0.875 ± 0.084</td><td>1.064 ± 0.044</td><td>0.022 ± 0.051</td><td>0.281</td><td>28</td></tr><tr><td>T<sub>spec</sub>, cor-M<sub>500</sub></td><td>Scatter</td><td>-1.024 ± 0.061</td><td>-0.016 ± 0.258</td><td>-0.673 ± 0.226</td><td>-0.205 ± 0.118</td><td>0.213 ± 0.138</td><td>2.036</td><td>28</td></tr><tr><td>L<sub>bol</sub>-M<sub>500</sub></td><td>Scatter</td><td>-0.589 ± 0.039</td><td>-0.156 ± 0.165</td><td>-0.532 ± 0.145</td><td>0.027 ± 0.076</td><td>0.066 ± 0.088</td><td>0.835</td><td>28</td></tr><tr><td>M<sub>gas</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-0.848 ± 0.053</td><td>-0.600 ± 0.223</td><td>-0.325 ± 0.196</td><td>-0.496 ± 0.103</td><td>-0.076 ± 0.119</td><td>1.533</td><td>28</td></tr><tr><td>Y<sub>X</sub>, 500-M<sub>500</sub></td><td>Scatter</td><td>-0.773 ± 0.031</td><td>0.031 ± 0.130</td><td>-0.717 ± 0.114</td><td>-0.081 ± 0.060</td><td>0.216 ± 0.069</td><td>0.518</td><td>28</td></tr><tr><td>d<sup>2</sup><sub>A</sub>Y<sub>500</sub>-M<sub>500</sub></td><td>Scatter</td><td>-0.868 ± 0.030</td><td>0.237 ± 0.127</td><td>-0.791 ± 0.112</td><td>-0.160 ± 0.058</td><td>0.256 ± 0.068</td><td>0.497</td><td>28</td></tr><tr><td>M<sub>500</sub>, hse, spec-M<sub>500</sub></td><td>Scatter</td><td>-0.733 ± 0.109</td><td>0.576 ± 0.461</td><td>-0.356 ± 0.404</td><td>-0.062 ± 0.212</td><td>0.053 ± 0.246</td><td>6.532</td><td>28</td></tr></tbody></table>

## APPENDIX C: SCATTER AND EVOLUTION OF THE BOLOMETRIC X-RAY LUMINOSITY-TEMPERATURE RELATION

In order to do an analogous analysis for the bolometric X-ray luminosity-temperature scaling relation (for both core-excised and non-core-excised X-ray spectroscopic temperatures), we apply a method similar to the one described in Section 4. In brief, we obtain a value for the median observable and the scatter at fixed temperature about the median relation as a function of temperature and

redshift as well as simple functional forms for their temperature and redshift dependences by fitting log-normal distributions in a few temperature bins at several redshifts. Specifically, the median scaling relation at the same nine different redshifts as previously was obtained by fitting a spline to the median bolometric X-ray luminosity-temperature relation computed in 10 equally logarithmically spaced temperature bins in the range $0.5 \le k_B T \le 5.0$ keV. The best-fitting spline is used to factor out the median relation when fitting a log-normal distribution to the distribution of $L_X$ in four temperature bins (chosen to be $0.5 \le k_B T \le 1.3$ keV, $1.3 \le k_B T \le 5.0$ keV).
<table><caption>**Table C1.** Results for all the physical models of fitting the evolving power law (equation C1) to both the median bolometric X-ray luminosity–temperature relation and the log-normal scatter about it (for both core-excised and non-core-excised X-ray spectroscopic temperatures). The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Simulation</th><th>Median or scatter</th><th>A</th><th>α</th><th>β</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>NOCOOL</td><td>Median</td><td>44.413 ± 0.010</td><td>1.206 ± 0.049</td><td>2.590 ± 0.024</td><td>0.242</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>REF</td><td>Median</td><td>43.519 ± 0.008</td><td>1.474 ± 0.041</td><td>2.442 ± 0.023</td><td>0.166</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.0</td><td>Median</td><td>43.478 ± 0.019</td><td>1.909 ± 0.092</td><td>3.211 ± 0.049</td><td>0.832</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.5</td><td>Median</td><td>43.255 ± 0.018</td><td>2.239 ± 0.087</td><td>2.834 ± 0.043</td><td>0.755</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>NOCOOL</td><td>Scatter</td><td>-0.569 ± 0.011</td><td>-0.647 ± 0.055</td><td>-0.053 ± 0.027</td><td>0.298</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>REF</td><td>Scatter</td><td>-0.377 ± 0.010</td><td>-0.561 ± 0.051</td><td>-0.494 ± 0.028</td><td>0.256</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.0</td><td>Scatter</td><td>-0.374 ± 0.018</td><td>-0.233 ± 0.091</td><td>-0.644 ± 0.048</td><td>0.814</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.5</td><td>Scatter</td><td>-0.300 ± 0.015</td><td>-0.429 ± 0.075</td><td>-0.584 ± 0.037</td><td>0.552</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>NOCOOL</td><td>Median</td><td>44.470 ± 0.011</td><td>1.227 ± 0.056</td><td>2.557 ± 0.026</td><td>0.253</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>REF</td><td>Median</td><td>43.491 ± 0.010</td><td>1.755 ± 0.048</td><td>2.563 ± 0.027</td><td>0.225</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.0</td><td>Median</td><td>43.416 ± 0.017</td><td>2.021 ± 0.084</td><td>3.246 ± 0.044</td><td>0.696</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.5</td><td>Median</td><td>43.079 ± 0.018</td><td>2.417 ± 0.089</td><td>2.730 ± 0.044</td><td>0.781</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>NOCOOL</td><td>Scatter</td><td>-0.561 ± 0.018</td><td>-0.592 ± 0.096</td><td>-0.156 ± 0.044</td><td>0.732</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>REF</td><td>Scatter</td><td>-0.398 ± 0.011</td><td>-0.637 ± 0.056</td><td>-0.202 ± 0.032</td><td>0.315</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.0</td><td>Scatter</td><td>-0.428 ± 0.015</td><td>-0.373 ± 0.074</td><td>-0.515 ± 0.039</td><td>0.542</td><td>33</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.5</td><td>Scatter</td><td>-0.446 ± 0.038</td><td>-0.092 ± 0.186</td><td>-0.544 ± 0.091</td><td>3.398</td><td>33</td></tr></tbody></table>

<table><caption>**Table C2.** Results for all the physical models of fitting the evolving broken power law (equation C4) to both the median bolometric X-ray luminosity–temperature relation and the log-normal scatter about it (for both core-excised and non-core-excised X-ray spectroscopic temperatures). The scatter is the log-normal scatter in the natural logarithm of the Y variable.</caption><thead><tr><th>Scaling relation</th><th>Simulation</th><th>Median or scatter</th><th>A''</th><th>α''</th><th>β''</th><th>γ''</th><th>δ''</th><th>χ<sup>2</sup></th><th>d.o.f.</th></tr></thead><tbody><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>NOCOOL</td><td>Median</td><td>44.426 ± 0.013</td><td>1.305 ± 0.048</td><td>2.353 ± 0.096</td><td>2.386 ± 0.066</td><td>0.223 ± 0.061</td><td>0.138</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>REF</td><td>Median</td><td>43.495 ± 0.013</td><td>1.488 ± 0.049</td><td>2.318 ± 0.115</td><td>2.586 ± 0.066</td><td>0.027 ± 0.073</td><td>0.141</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.0</td><td>Median</td><td>43.490 ± 0.019</td><td>1.606 ± 0.070</td><td>4.122 ± 0.153</td><td>3.450 ± 0.095</td><td>-0.700 ± 0.097</td><td>0.289</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.5</td><td>Median</td><td>43.276 ± 0.023</td><td>2.011 ± 0.087</td><td>3.510 ± 0.175</td><td>2.933 ± 0.111</td><td>-0.483 ± 0.108</td><td>0.451</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>NOCOOL</td><td>Scatter</td><td>-0.591 ± 0.018</td><td>-0.663 ± 0.067</td><td>-0.064 ± 0.132</td><td>0.121 ± 0.091</td><td>-0.047 ± 0.084</td><td>0.262</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>REF</td><td>Scatter</td><td>-0.391 ± 0.016</td><td>-0.618 ± 0.060</td><td>-0.366 ± 0.141</td><td>-0.340 ± 0.081</td><td>-0.147 ± 0.089</td><td>0.211</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.0</td><td>Scatter</td><td>-0.398 ± 0.030</td><td>-0.253 ± 0.115</td><td>-0.663 ± 0.249</td><td>-0.456 ± 0.154</td><td>-0.054 ± 0.158</td><td>0.771</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub></td><td>AGN 8.5</td><td>Scatter</td><td>-0.274 ± 0.025</td><td>-0.458 ± 0.093</td><td>-0.433 ± 0.188</td><td>-0.725 ± 0.119</td><td>-0.054 ± 0.117</td><td>0.522</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>NOCOOL</td><td>Median</td><td>44.478 ± 0.014</td><td>1.353 ± 0.056</td><td>2.277 ± 0.096</td><td>2.370 ± 0.070</td><td>0.246 ± 0.061</td><td>0.132</td><td>29</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>REF</td><td>Median</td><td>43.470 ± 0.016</td><td>1.839 ± 0.057</td><td>2.233 ± 0.139</td><td>2.614 ± 0.078</td><td>0.201 ± 0.086</td><td>0.189</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.0</td><td>Median</td><td>43.420 ± 0.018</td><td>1.760 ± 0.069</td><td>4.016 ± 0.150</td><td>3.489 ± 0.092</td><td>-0.610 ± 0.096</td><td>0.279</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.5</td><td>Median</td><td>43.035 ± 0.023</td><td>2.294 ± 0.088</td><td>2.945 ± 0.178</td><td>3.148 ± 0.112</td><td>-0.285 ± 0.113</td><td>0.463</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>NOCOOL</td><td>Scatter</td><td>-0.591 ± 0.032</td><td>-0.547 ± 0.130</td><td>-0.293 ± 0.221</td><td>0.024 ± 0.161</td><td>0.046 ± 0.141</td><td>0.699</td><td>29</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>REF</td><td>Scatter</td><td>-0.427 ± 0.019</td><td>-0.633 ± 0.069</td><td>-0.309 ± 0.167</td><td>-0.010 ± 0.094</td><td>-0.005 ± 0.104</td><td>0.273</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.0</td><td>Scatter</td><td>-0.488 ± 0.018</td><td>-0.420 ± 0.070</td><td>-0.570 ± 0.152</td><td>-0.059 ± 0.094</td><td>-0.128 ± 0.098</td><td>0.287</td><td>31</td></tr><tr><td>L<sub>bol</sub>-T<sub>spec</sub>, cor</td><td>AGN 8.5</td><td>Scatter</td><td>-0.497 ± 0.061</td><td>0.052 ± 0.233</td><td>-1.089 ± 0.472</td><td>-0.352 ± 0.297</td><td>0.304 ± 0.300</td><td>3.239</td><td>31</td></tr></tbody></table>

≤ 2.17 keV, 2.17 ≤ $k_B T$ ≤ 3.0 keV and 3.0 ≤ $k_B T$ ≤ 8.0 keV). The rest of the procedure is identical to the one described in Section 4. We can then fit an evolving power law of the form

$$L_X = 10^A E(z)^\alpha \left( \frac{k_B T}{2 \text{keV}} \right)^\beta, \qquad (C1)$$

an evolving broken power law of the form

$$L_X = 10^{A'} E(z)^{\alpha'} \left( \frac{k_B T}{2\text{keV}} \right)^{\epsilon'}, \qquad (C2)$$

where

$$\epsilon' = \begin{cases} \beta' & \text{if } k_B T \le 2 \text{ keV} \\ \gamma' & \text{if } k_B T > 2 \text{ keV}, \end{cases} \qquad (C3)$$

and finally an evolving broken power law with a redshift-dependent low-mass power-law index of the form

$$L_X = 10^{A''} E(z)^{\alpha''} \left( \frac{k_B T}{2 \text{keV}} \right)^{\epsilon''}, \qquad (C4)$$

where

$$\epsilon'' = \begin{cases} \beta'' + \delta'' E(z) & \text{if } k_B T \le 2 \text{ keV} \\ \gamma'' & \text{if } k_B T > 2 \text{ keV} \end{cases} \qquad (C5)$$

to the median relation and the log-normal scatter about it as a function of both temperature and redshift. The results of the evolving power law and broken power law with an evolving low-mass power-law index fitting for all the physical models are presented in Tables C1 and C2, respectively.

In Figs C1 and C2, we show the evolution of the mass slopes from $z=0$ to $1.5$ for the bolometric X-ray luminosity–X-ray
<figure><img src="image_21.png" /><figcaption>**Figure C1.** Evolution of the temperature slope from $z = 0$ to 1.5 for the temperature–bolometric X-ray luminosity scaling relations (for both core-excised and non-core-excised X-ray spectroscopic temperatures). In each panel, we plot the redshift evolution of the best-fitting power-law indices obtained by fitting the broken power law given by equation (C1) at each individual redshift. The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The bolometric X-ray luminosity–X-ray spectroscopic temperature (be it core excised or not) relation is steeper than self-similar for all the models considered and this independently of redshift.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure C2.** Evolution of the temperature slope from $z = 0$ to 1.5 for the temperature–bolometric X-ray luminosity scaling relations (for both core-excised and non-core-excised X-ray spectroscopic temperatures). In each subpanel, we plot the redshift evolution of the low-temperature (top panel) and high-temperature (bottom panel) best-fitting power-law indices obtained by fitting the broken power law given by equation (C2) at each redshift independently. The solid curves (red, orange, blue and green) correspond to the different simulations and the horizontal dashed lines to the self-similar expectation, respectively. The bolometric X-ray luminosity–X-ray spectroscopic temperature (be it core excised or not) relation is steeper than self-similar for all the models considered and this independently of temperature and redshift.</figcaption></figure>

s spectroscopic temperature (be it core excised or not) for each of the four physical models. Fig. C1 shows the redshift evolution of the temperature slope obtained by fitting the power law given by equation (C1) at each individual redshift whereas Fig. C2 shows the redshift evolutions of the low-temperature (top panel) and high-temperature (bottom panel) temperature slopes resulting from the fitting of the broken power law given by equation (C2) at each redshift independently. The solid curves (red, orange, blue and green) corresponding to the different simulations and the horizontal

<figure><img src="image_23.png" /><figcaption>**Figure C3.** Evolution of the normalization from $z = 0$ to 1.5 for the temperature–bolometric X-ray luminosity scaling relations (for both core-excised and non-core-excised X-ray spectroscopic temperatures) for each of the simulations. The normalizations of each scaling relation in the four temperature bins (denoted by solid lines of different colours) have been normalized by the self-similar expectation for the redshift evolution at fixed temperature (shown as an horizontal dashed line). The amplitude of the bolometric X-ray luminosity–X-ray spectroscopic temperature relation evolves positively for all the physical models with an amplitude which is strongly temperature dependent, slightly redshift dependent, strongly sensitive to the non-gravitational physics of galaxy formation (it becomes more positive with increasing feedback intensity) and slightly sensitive to cosmology (slightly faster evolution in the WMAP7 cosmology, especially for the highest temperature bin).</figcaption></figure>

dashed lines to the self-similar expectation, respectively. The bolometric X-ray luminosity–X-ray temperature (be it core excised or not) relation is steeper than the self-similar expectation of 2 (see equation 8) for bolometric X-ray luminosity for all the models considered here and this independently of temperature and redshift. The bolometric X-ray luminosity–X-ray temperature relation therefore qualitatively behaves in exactly the same way as the bolometric X-ray luminosity–total mass (see Section 5 and Figs 2 and 3). Hence,
<figure><img src="image_24.png" /><figcaption>**Figure C4.** Evolution of the log-normal scatter from $z=0$ to 1.5 for the temperature–bolometric X-ray luminosity scaling relations (for both core-excised and non-core-excised X-ray spectroscopic temperatures). For each simulation and each scaling relation, we plot the log-normal scatter as a function of $T_{\text{spec, (cor)}}$ and denote the redshift using lines of different colours. The log-normal scatter of the bolometric X-ray luminosity–X-ray spectroscopic temperature relation varies only mildly with temperature (typically changing by a factor of ≈2–3 over ~1 decade in temperature), is somewhat sensitive to non-gravitational physics (the amplitude tends to increase with increasing feedback intensity), but displays a moderately strong redshift dependence (it tends to decrease with increasing redshift).</figcaption></figure>
the deviations from self-similarity are most likely due to the same reasons, namely that the gas does not trace the dark matter affecting both density and temperature, which are also potentially further affected by non-thermal pressure support.

In Fig. C3, we show the evolution of the normalization from $z = 0$ to 1.5 for the bolometric X-ray luminosity–X-ray spectroscopic temperature (be it core excised or not) relation for the four physical models considered here. The normalizations of each scaling relation in the four temperature bins (denoted by solid lines of different colours) have been normalized by the self-similar expectation for the redshift evolution at fixed mass (shown as an horizontal dashed line). As was the case for the bolometric X-ray luminosity–total mass relation and the models that include AGN feedback (see Section 6 and Fig. 4), the amplitude of the bolometric X-ray luminosity–X-ray temperature relation evolves positively for all the physical models with an amplitude which is strongly temperature dependent, slightly redshift dependent (it slightly flattens out with increasing redshift) and strongly sensitive to the non-gravitational physics of galaxy formation (it becomes more positive with increasing feedback intensity). The reversal of the direction of the deviation from self-similarity for the models without AGN feedback (from negative for the bolometric X-ray luminosity–total mass relation to positive for the bolometric X-ray luminosity–X-ray spectroscopic temperature relation is straightforwardly explained by the negative evolution of the mass–temperature relation. Short et al. (2010) also found positive evolution for zoom simulations which include a semi-analytic prescription for AGN feedback.

In Fig. C4, we show the evolution of the log-normal scatter (at fixed temperature) from $z = 0$ to 1.5 for the bolometric X-ray

luminosity–X-ray spectroscopic temperature [be it core excised (bottom subpanels) or not (top subpanels)] relation for the four physical models considered here. For each simulation and each scaling relation, we plot the log-normal scatter as a function of $T_{\text{spec. (cor)}}$ and denote the redshift using lines of different colours. As was the case for the bolometric X-ray luminosity–total mass relation (see Section 7 and Fig. 6), the log-normal scatter of the bolometric X-ray luminosity–X-ray spectroscopic temperature relation varies only mildly with temperature (typically changing by a factor of ≈2–3 over ~1 decade in temperature), is somewhat sensitive to non-gravitational physics (the amplitude tends to increase with increasing feedback intensity), but displays a moderately strong redshift dependence (it tends to decrease with increasing redshift).

All the differences between the behaviours of the bolometric X-ray luminosity–mass and bolometric X-ray luminosity–X-ray spectroscopic temperature can be easily explained by the behaviour of the mass–temperature relation.

Finally, it is worth noting that the calibrated BAHAMAS model has a much steeper bolometric X-ray luminosity–X-ray spectroscopic temperature relation than the cosmo-OWLS AGN models apart at the high-temperature end and that the bolometric X-ray luminosity–X-ray spectroscopic temperature relation evolves slightly faster in the WMAP7 cosmology, especially for the highest temperature bin.

This paper has been typeset from a TeX/LA TeX file prepared by the author.