---
abstract: |
    We analyze hydrodynamical and cosmological simulations of galaxy clusters to study scaling relations between the cluster total masses and observable quantities such as gas luminosity, gas mass, temperature, and $Y_X$, i.e., the product of the last two properties. Our simulations are performed with the Smoothed-Particle-Hydrodynamic GADGET-3 code and include different physical processes. The twofold aim of our study is to compare our simulated scaling relations with observations at low ($z\thickapprox0$) and intermediate ($z\thickapprox0.5$) redshifts and to explore their evolution over the redshift range $z=0-2$.

    The result of the comparative study shows a good agreement between our numerical models and real data. We find that AGN feedback significantly affects low-mass haloes at the highest redshifts resulting in a reduction of the slope of the mass $-$ gas mass relation $(\sim13\%)$ and the mass $- Y_X$ relation $(\sim10\%)$ at $z=2$ in comparison to $z=0$. The drop of the slope of the mass $-$ temperature relation at $z=2$ $(\sim14\%)$ is, instead, caused by early mergers. We investigate the impact of the slope variation on the study of the evolution of the normalization.

    We conclude that the observed scaling relations should be limited to the redshift range $z=0-1$ for cosmological studies because in that redshift range the slope, the scatter, and the covariance matrix of the relations do not exhibit significant evolution. The mass$-Y_X$ relation continues to be the most suitable relation for this goal. Extending the analysis to the redshift range between 1 and 2 will be crucial to evaluate the impact generated by the AGN feedback.
author:
- |
    N. Truong$^{1}$, E. Rasia$^{2,3}$, P. Mazzotta$^{1,4}$, S. Planelles$^{5,6}$, V. Biffi$^{3,5}$, D. Fabjan$^{3,7}$, A. M. Beck$^{8}$, S. Borgani$^{3,5,9}$, K. Dolag$^{8}$, G. L. Granato$^{3}$, G. Murante$^{3}$, C. Ragone-Figueroa$^{3,10}$, L. K. Steinborn$^{8}$\
     \
    $^1$ Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca Scientifica, I-00133, Roma, Italy\
    $^2$ Department of Physics, University of Michigan, 450 Church St., Ann Arbor, MI 48109\
    $^3$ INAF, Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131, Trieste, Italy\
    $^4$ Harvard-Smithsonian Center for Astrophysics, 60 Gardner Street, Cambridge, MA 02138, USA\
    $^5$ Dipartimento di Fisica dell’ Università di Trieste, Sezione di Astronomia, via Tiepolo 11, I-34131 Trieste, Italy\
    $^{6}$ Departamento de Astronomía y Astrofísica, Universidad de Valencia, c/Dr. Moliner, 50, 46100 Burjassot, Valencia, Spain\
    $^{7}$ Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia\
    $^{8}$ University Observatory Munich, LMU, Scheinerstr. 1, D-81679, Munich, Germany\
    $^{9}$ INFN, Instituto Nazionale di Fisica Nucleare, Trieste, Italy\
    $^{10}$ Instituto di Astronomía Teórica y Experimental (IATE), Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina (CONICET),\
         Observatorio Astronómico, Universidad Nacional de Córboda, Laprida 854, X50000BGR, Córdoba, Argentina
bibliography:
- 'ref.bib'
title: 'Hydrodynamic Simulations of Galaxy Clusters: Scaling Relations and Evolution'
---

[galaxies: clusters: general — galaxies: clusters: intracluster medium — X-ray: galaxy: clusters — methods: numerical]{}

Introduction {#sec:intro}
============

Clusters of galaxies are the latest gravitationally bound structures to form. They begin their collapse at around redshift $z \sim 2$ when the dark-energy component starts to be relevant (even if not dominant yet) and just after the baryonic Universe has experienced an extremely energetic phase due to the high star formation level and the intense AGN (Active Galactic Nuclei) activity [@aird.etal.2015]. For the reason that clusters’ formation happens in such peculiar time, these massive objects are an invaluable laboratory for the study of astrophysical processes and for the derivation of cosmological parameters.

High-redshift investigations and the comparison with local observations are particularly advantageous. One of the most powerful cosmological measure is, indeed, the [*evolution*]{} of the mass function [@borgani_guzzo; @voit2005; @2009ApJ...692.1033V; @allen.etal.2011; @2011ASL.....4..204B; @2015SSRv..188...93P]. As other example, two direct probes of the AGN-activity history are $(i)$ the entropy profile of local and [*high-z*]{} objects [@dubois.etal.2011; @chaudhuri.etal.2012; @ettori.etal.2013; @pointecouteau.etal.2013] and $(ii)$ the relation between the clusters X-ray luminosity and their temperature and its evolution [@2012MNRAS.424.2086H; @takey.etal.2013; @2015arXiv151203833G]. To decode such important measurements and to extend high-redshift studies of galaxy clusters, the X-ray community is planning and working on future missions such as e-ROSITA[^1][@2012arXiv1209.3114M] and Athena[^2] [@nandra.etal.2013]. These will be complementary to millimetric surveys, e.g. SPT-G3 [@benson.etal.2014], or optical ones, e.g. EUCLID[^3] [@laureijs.etal.2011] or LSST[^4][@Ivezic.etal.2008]. Both X-ray missions will build copious samples of clusters and groups of galaxies at $z \sim 1$ or above to advance high-z investigations that are currently limited mostly to $z \leq 0.7$ objects with sparse data from most distant times (). The simplest approach to measure the mass of a large sample of objects is to access to [*scaling relations*]{} between the total mass and some observable quantities that are easier to derive. In X-rays studies, the most common mass-proxies are the [*gas mass*]{}, $M_{\rm gas}$, which can be extracted from the surface brightness profile, the [*temperature*]{}, $T$, which is solidly estimated from X-ray spectra with, at least, one thousand of counts, and their combination, $Y_X={M_{\rm g}}\times T$ [e.g., @2014MNRAS.437.1171M; @mantz.etal.2016 for recent and detailed studies on X-ray scaling relations]. The $Y_X$ parameter was firstly introduced as proxy of pressure by [@2006ApJ...650..128K] who proved its advantages through the analysis of hydrodynamical simulations. Since then, it has been diffusely adopted in X-ray cosmological studies. The biggest benefit of this parameter relies on the different reactions by the gas mass and temperature in non-equilibrium situations such as a powerful AGN burst activity [@2011MNRAS.416..801F] that can be triggered by gravitational mergers. The AGN reduces the gas mass by expelling some gas from the core and, at the same time, it heats the ICM. Due to these behaviors the $Y_X$ parameter does not substantially depend on the dynamical state or on the central AGN activity.

One of the X-ray measurements most represented in scaling-relation studies is the [*X-ray luminosity*]{} because it can easily and immediately be obtained from shallow observations. Scaling relations involving the luminosity, in particular the $L-M$ relation, are used to determine the connection between the flux limits of a planned survey with the minimum mass that can be observed. This relation, therefore, plays a crucial role in establishing the selection function of X-ray samples [e.g. @nord.etal.2008]. The history of the $L-T$ relation has been equally important because from its first determination it was clear that this relation provides information on the physics of the cluster core [e.g., @fabian.etal.1994] and the phenomena of feedback by stars or AGN [e.g., @maxim98; @2012MNRAS.421.1583M].

Over the last decade, a significant number of realistic hydrodynamical simulations enabled the study of scaling relations (@2011ASL.....4..204B for a review). A special attention was dedicated to the effects of feedback from stars [@nagai.etal.2007] and AGNs [@short_thomas_2009; @puchwein.etal.2008; @2010MNRAS.408.2213S; @2014MNRAS.438..195P; @lebrun.etal.2014; @pike.etal.2014], to the evolution of the relations up to $z\sim1$ [@2012ApJ...758...74B; @2011MNRAS.416..801F; @planelles.etal.2016; @lebrun.etal.2016], and to develop a theoretical framework to exploit the simultaneous analysis of multiple signals [@2010ApJ...715.1508S; @evrard.etal.2014]. This paper, based on realistic models of ICM quantities obtained through a new BH accretion model, AGN feedback and hydro scheme, extends the analysis to $z=2$ in vision of up-coming observational results. Most importantly, this work analyzes a simulated sample that recently proved to naturally form in a cosmological context cool-core (CC) and non-cool-core (NCC) clusters [@2015ApJ...813L..17R]. The diversity of the core entropy level in the centermost regions is mostly ascribed to the dynamical history of the system and to its overall AGN activity. We expect that in our simulated clusters the balance between radiative cooling, which forms stars, and AGN feedback, which heats the gas, is realistically reached as the systems evolve in and interact with the embedded cosmological environment. Therefore, it is compelling to study the relation between the ICM properties in an epoch right after the maximum level of the AGN activity.

The paper is arranged in the following way: in Section 2 we provide a short description of the simulated sample and refer to companion papers for a more detailed explanation of the code. Section 3 presents the computation of ICM structural quantities, the mass-proxy relations, the luminosity-based relations, and fitting methods. In Section 4, we examine the validity of our simulated data by comparing to observations at low ($z\leq0.25$) and intermediate ($0.25\leq z \leq 0.6$) redshifts. Section 5 is dedicated to explore the theoretical evolution of scaling relations from $z=0$ to $z=2$. Finally, a summary of the obtained results and conclusions is given in Section 6. All the quantities for the scaling relations are evaluated at $R_{500}$ defined as the radius of the sphere whose density is 500 times the critical density of the Universe at the considered redshift. Throughout the paper, the symbol $\log$ indicates the [*decimal*]{} logarithm and the uncertainty at 1 $\sigma$ on the best-fit parameters represents the 68.4 per cent confidence maximum-probability interval.

\[sec:cmrr\]

Simulations {#Simulation}
===========

Our analysis is based on three runs of simulated clusters. These are selected from a parent DM-only cosmological volume of 1 $h^{-3}$ Gpc$^3$ [@yoshida.etal.2001; @2011MNRAS.418.2234B] and re-simulated at higher resolution and with the inclusion of baryons.

The re-simulated Lagrangian regions are chosen around 24 massive clusters and five isolated groups and extend a radius equal to 5 times their virial radius. Their particles, identified at redshift $z=0$, are traced back to high redshifts [@2011MNRAS.418.2234B] and their number is increased to achieve a better spatial and mass resolution. The resimulations are carried out with an improved version of the GADGET-3 code [@springel05] where we activated the full set of SPH modifications as described in [@2015arXiv150207358B]. In short, these allow the SPH method to perform better in hydrodynamical standard tests including weak and strong shocks, gas mixing, and self-gravitating clouds.

We analyze three different runs of the same cluster sample. They differ on the level of sophistication on the description of the ICM physics. The joint analysis of the three samples allow us to qualitatively conjecture whether a possible evolution seen on the scaling relations is caused by the environment and the interactions with the large-scale structure distribution, in case of sample-independent results, or by the astrophysical phenomena affecting the ICM, in the other situation of different outcomes depending by the sample.

[$\bullet$  [[[nr]{}]{}]{} (Non-radiative)]{}. These simulations are carried out by the exact code used in [@2015arXiv150207358B]. The main variations with respect to the standard GADGET code include: the choice of a higher-order Wendland $C^4$ kernel function; a time-dependent artificial viscosity scheme; a thermal diffusion term (or artificial conduction) that improves the treatment of contact discontinuities and promotes fluid mixing. The performance of this code with respect to other particle- or grid-based codes is described in [@sembolini.etal.2016b] and [@sembolini.etal.2016a]. These works show how the cluster simulated with the improved GADGET version (G3-XArt) performs in a similar way to grid codes, AREPO, and the most modern SPH schemes.

[$\bullet$  [[[csf]{}]{}]{} (Cooling, Star formation and stellar Feedback)]{}. The radiative runs consider metal-dependent radiative cooling rates accordingly to [@2009MNRAS.399..574W] where 15 different elements (H, He, C, Ca, O, N, Ne, Mg, S, Si, Fe, Na, Al, Ar, Ni) are followed; the effect of a uniform UV/X-ray background radiation [@2001cghr.confE..64H]; the feedback by supernovae, as originally prescribed by [@2003MNRAS.339..289S]; the chemical model by [@2007MNRAS.382.1050T] to account for the metal enrichment (for further details see @2014MNRAS.438..195P). The wind velocity associated to the kinetic supernovae feedback is 350 km s$^{-1}$.

[$\bullet$  [[[[agn]{}]{}]{}]{}.]{} This last run is built on the [[[csf]{}]{}]{} simulations with the addition of the AGN feedback, modeled following [@steinborn.etal.2015] who improved the original model by [@springel05] by introducing a detailed and phenomenologically-driven description. The thermodynamical and metal properties and the dynamical state of the main clusters of this sample are presented also in [@2015ApJ...813L..17R], [@planelles.etal.2016], and [@biffi.etal.2016].

The cosmological setting is a $\Lambda CDM$ model with the matter density $\Omega_m=0.24$, the dark energy density $\Omega_{\Lambda}=0.76$, $n_s=0.96$, the amplitude of the density fluctuation power spectrum $\sigma_8=0.8$, and the Hubble parameter $H=72$ km s$^{-1}$ Mpc$^{-1}$. The observational data, presented in Section 4, are rescaled to the simulated cosmology.

The Plummer gravitational softening of the DM particles is equivalent to 3.75 $h^{-1}\ \rm{kpc}$ in physical units up to $z=2$ and in comoving units at higher redshifts. The gas softening is, instead, fixed to the same value in comoving coordinates at all times. The mass of the DM particle is $8.47 \times 10^{8} \ h^{-1} M_{\odot}$ and the initial mass of the gas particle is $1.53 \times 10^{8} \ h^{-1} M_{\odot}$.

The sample includes all the objects in the high-resolution Lagrangian regions with $M_{500} > 7 \times 10^{13}  E(z)^{-1} h^{-1} M_{\odot}$ with $E(z)=H(z)/H_0=(\Omega_M\times(1+z)^3+\Omega_{\Lambda})^{1/2}$. The number of [[[agn]{}]{}]{} clusters selected at any redshift of interest and the corresponding mass range are presented in Fig. \[fig:a2\]. Notice that at the highest redshifts, $z=1, z=1.5$ and $z=2$, our sample includes 74, 60, and 36 objects. The smallest systems at $z=2$ contains more than $3.5 \times 10^4$ particles providing good estimates of global quantities. We do not extend the sample to smaller systems, even if we would have obtained numerically robust global measures, because none of the planned missions will reach such small masses (see @borm.etal.2014 for e-ROSITA, @pointecouteau.etal.2013 for Athena and @benson.etal.2014 for SPT-3G).

Method of Analysis {#Method of Analysis}
==================

Computing ICM Quantities
------------------------

In the following, we briefly describe how we compute the relevant quantities from our simulated datasets. [**Masses.**]{} The total mass, $M$, is calculated by summing the contribution of all the species of particles (dark-matter, gas, and stars) within $R_{500}$. For the gas mass, ${M_{\rm g}}$, we sum up only the hot gas component[^5] of all multi-phase particles not containing more than 10% in cold mass and, therefore, no star-forming particles.

[**Temperatures.**]{} We consider both the mass-weighted temperature and the spectroscopic-like temperature [@2004MNRAS.354...10M]. The mass-weighted temperature is provided by: $${T_{\rm mw}}=\frac{\sum_i m_i T_i}{\sum_i m_i},
      \label{eq:tmw}$$ where $m_i$ and $T_i$ are the hot-gas mass and temperature of the $i^{th}$ gas particle. The spectroscopic-like temperature is introduced to ease the comparison with X-ray observations and the formula was derived considering the non-flat response of the instruments on board of Chandra and XMM-Newton: $${T_{\rm sl}}= \frac{\sum_i \rho_i m_i T_i^{0.25}}{\sum_i\rho_i m_i T_i^{-0.75}}, 
\label{eq:tsl}$$ where $\rho_i$ is the particle gas density [@2004MNRAS.354...10M see also @vikh2006]. For this computation, we used only particles emitting in the X-ray band with $T_i> 0.5$ keV.

[**The parameter $Y_{X}$.**]{} As previously said, this parameter is equivalent to the product of gas mass and the global temperature within $R_{500}$ and it is a powerful proxy for the total thermal content of the ICM because it is almost insensitive to the physical processes included in simulations [@2010ApJ...715.1508S; @2011MNRAS.416..801F; @2012ApJ...758...74B; @sembolini.etal.2014] and to the dynamical status of the clusters [@rasia.etal.2011; @kay.etal.2012]. $Y_X$ is derived from X-ray observations, therefore we adopt the spectroscopic-like temperature in its expression: $$Y_{X} = {M_{\rm g}}\times {T_{\rm sl}}. 
    \label{eq:yx}$$ We stress, however, that in this new set of simulations, which includes the artificial diffusion term, [*the difference between the mass-weighted temperature and the spectroscopic-like one is largely reduced*]{} [see also @biffi.etal.2016] with respect to previous results [@rasia.etal.2014]. In the current $z=0$ AGN sample, the mean and standard deviation values for the ratio ${T_{\rm mw}}/{T_{\rm sl}}$ computed in the sphere of radius $R_{500}$ are equal to 1.02 and 0.10, respectively.

[**X-ray Luminosity.**]{} The bolometric luminosity is computed by summing up the contribution of the emissivity, $\epsilon_i$, of all gas particles within the sphere of radius $R_{500}$: $$L = \sum_{i}\epsilon_i= \sum_i n_{e,i}n_{H,i}\Lambda(T_i,Z_i)\Delta V_i, 
 \label{eq:lx}$$ where $n_{e,i}$, $n_{H,i}$ are number densities of electrons and hydrogen atoms, respectively, $\Delta V_i = m_i/\rho_i$ is the particle’s volume and $\Lambda_i$ is the interpolation of the cooling function pre-calculated in a fine grid of temperatures and metallicities. The cooling-function tables are created by adopting the APEC model [@2001ApJ...556L..91S] in XSPEC and by referring to the \[0.01-100\] keV energy band.

The Scaling Relations and the Self-Similar Prediction
-----------------------------------------------------

The total mass of a cluster can be related to the various ICM quantities presented in the previous sections through simple power-law models. [@kaiser86] analytically derived the functional shapes of the expected scaling relations under the assumption of virial equilibrium between the kinetic and thermal energy of a galaxy cluster and its gravitational potential [ a recent extension of the formalism in dependance on the ICM physics is provided in @ettori2015]. According to this model, called self-similar (SS), the cluster total mass is the only parameter that defines both thermal and dynamical properties of the ICM (e.g., @2013SSRv..177..247G and references therein).

[**The $M-{M_{\rm g}}$ relation.**]{} The total mass and the gas mass are linearly related: $$M = C_{{M_{\rm g}}} E(z)^{0} {M_{\rm g}}, 
   \label{eq:mmg}$$ without any dependence with the redshift. Here and in the following equations, $C$ is a constant.

[**The $M-T$ relation.**]{} The total mass is related to the temperature according to: $$M =  C_{T} E(z)^{-1} T^{3/2}. 
 \label{eq:mt}$$

[**The $M-{Y_{\rm X}}$ relation.**]{} From Eqs. \[eq:yx\], \[eq:mmg\], and \[eq:mt\], the self-similar scaling of ${Y_{\rm X}}$ with cluster mass is given by $$M = C_{{Y_{\rm X}}} E(z)^{-2/5} {Y_{\rm X}}^{3/5}. 
  \label{eq:myx}$$

[**The luminosity relations: $L-T$ and $L-M$.**]{} Assuming thermal bremsstrahlung emission, the luminosity can be related to the total mass by $$L \propto E(z)^2 T^{1/2} f_{gas}^2 M, \label{8}$$ where $f_{gas}\equiv {{M_{\rm g}}}/M\equiv1/C_{M_g}$ is the gas fraction. Given the expression of the $M-{M_{\rm g}}$ and $M-T$ relations (Eqs. \[eq:mmg\] and \[eq:mt\]), the self-similar form for the $L-T$ relation is expressed as $$L = C_{LT}E(z)T^2, 
 \label{eq:lt}$$ and for the $L-M$ relation as $$L = C_{LM}E(z)^{7/3}M^{4/3}.
 \label{eq:lm}$$

Fitting Method
--------------

       Relation                         $X_0$                     $\beta$   $\gamma$
  ------------------ ------------------------------------------- --------- ----------
   $M-{M_{\rm g}}$        $2\times10^{13}\ M_\odot h^{-1}$           1         0
        $M-T$                     $3\ \textrm{keV}$                 3/2        -1
   $M-{Y_{\rm X}}$    $4\times10^{13}\ M_\odot h^{-1} \rm{keV}$     3/5       -2/5
   $L-{T_{\rm sl}}$               $3\ \textrm{keV}$                  2         1
        $L-M$             $5\times10^{14}\ M_\odot h^{-1}$          4/3       7/3

  : For each scaling relation (1st column): the pivot point, $X_0$ (2nd column), the self-similar slope, $\beta$ (3rd column), and the self-similar evolution parameter, $\gamma$ (4th column).

\[tabella1\]

To study the evolution of the parameters of the scaling relations, we fit them through a generic expression: $$\log F = \log C + \gamma \log E(z) + \beta \log (X/X_0).
   \label{11}
   $$ When $F$ represents the total mass, $M_{500}$, $X$ refers to one of the following mass-proxy: ${M_{\rm g}}$, ${T_{\rm mw}}$, ${T_{\rm sl}}$, ${Y_{\rm X}}$. Instead, when $F$ is the X-ray luminosity, $X$ stands for the spectroscopic-like temperature or the total mass. The value of each pivot point, $X_0$ is listed in Table \[tabella1\]. At first, we fix $\gamma$ to the SS expectation values (also listed in Table \[tabella1\]) and we derive $C$ and $\beta$ at each [*independent*]{} redshift (Sections 5.1 and 5.2). Later, to study the evolution of the normalization, we let $\gamma$ free to vary (Section 5.3).

We employ the routines provided by [@2007ApJ...665.1489K] to investigate the parameter space and to perform the linear regression in logarithmic space. The `linmix_err.pro` routine is applied to the single linear regression with $C$ and $\beta$ as parameters. The `mlinmix_err.pro` is adopted when we include $\gamma$ as free parameter. The independent variables are modeled as Gaussian distribution from which we extract the mean and standard deviation that we refer as the parameter best-fit value with its uncertainty. The Bayesian approach allows us to treat the intrinsic scatter, $\sigma$, estimated via the method of Monte Carlo Markov Chains (MCMC), as a free parameter. The MCMC is carried out by using Metropolis-Hastings (in `linmix_err.pro`) and Gibb Samplers (in `mlinmix_err.pro`) algorithms.

Comparison between simulated and observed scaling relations
===========================================================

Observational data sets
-----------------------

In this Section, we compare our theoretical predictions for the scaling relations with some observational results. For the mass-proxy relations, due to uncertainties associated to the amplitude of the X-ray hydrostatic mass bias derived both in observations [@vonderlinden.etal.2014; @maughan.etal.2015] and simulations [@2007ApJ...655...98N; @2012NJPh...14e5018R; @biffi.etal.2016], we prefer to refer to masses estimated via the gravitational lensing technique. More specifically, we refer to [@2013ApJ...767..116M], [@kettula.etal.2013], and [@2015arXiv151203857L]. The first sample contains 50 galaxy clusters in the redshift range $0.155<z<0.55$. The optical data are taken from the Canadian Clusters Comparison Project (CCCP, see @hoekstra.etal.2015), while the X-ray properties are based on the combined data from the $Chandra$ Observatory and the $XMM-Newton$ telescope. [@kettula.etal.2013] analyze 10 groups with redshifts between 0.1 and 0.8 taken from the COSMOS field and observed by the Hubble Space Telescope and the XMM-Newton telescope. [@2015arXiv151203857L] provides X-ray observations of 38 clusters with $z<0.6$ as part of the XXL survey. These data are associated with independent mass measurements from the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS, e.g. @2012MNRAS.427..146H; @2013MNRAS.433.2545E). The temperature measurements are restricted to the inner 300 kpc of each object.

Comparison to Observations at local and intermediate redshifts
--------------------------------------------------------------

[![image](fig2a){width="40.00000%"}]{} [![image](fig2b){width="40.00000%"}]{} [![image](fig2c){width="40.00000%"}]{} [![image](fig2d){width="40.00000%"}]{}

![image](fig3a){width="40.00000%"} ![image](fig3b){width="40.00000%"} ![image](fig3c){width="40.00000%"} ![image](fig3d){width="40.00000%"}

![The $z=0$ [[[agn]{}]{}]{} simulated luminosity-temperature relation is shown in blue and black asterisks, respectively, for cool-core and non-cool-core systems [@2015ApJ...813L..17R]. For comparison, the cool-core and non-cool-core clusters observed by @2013ApJ...767..116M are reported with filled and open magenta circles, respectively. []{data-label="fig:a1"}](fig4){width="50.00000%"}

We present the comparison between our [[[agn]{}]{}]{} simulated scaling relations at $z=0$ and $z=0.5$ and those obtained from observations at $z\leq 0.25$ and at $0.25\leq z \leq 0.6$ in Figs. \[fig:1\] and \[fig:2\], respectively. Each figure shows four scaling relations: $M-M_{g}$, $M-{T_{\rm sl}}$, $M-{Y_{\rm X}}$, and $L-M$. In this Section, all properties are measured within $R_{500}$ without any core excision to enlighten the good match that we obtain even including the region which is the most difficult to model. For the comparison, the observational measurements of $M$, ${M_{\rm g}}$, ${Y_{\rm X}}$ and $L$, which depend on $h^{-1}, h^{-5/2}, h^{-5/2}$ and $h^{-2}$, have been rescaled according to the value of the Hubble parameter of the simulation. In general, we observe a satisfactory consistency between simulations and observations for all relations, especially considering that [@2013ApJ...767..116M] measurements of the masses should be increased by $\sim20\%$ accordingly to the revised work from the group [@hoekstra.etal.2015].

The perfect accordance in the normalization of both the $M-{M_{\rm g}}$ and $L-M$ relations assures that the simulated baryon fraction is realistic over the mass range investigated. The overall $M-T$ relation is also consistent at both redshifts. It should be noted that the temperature in [@2015arXiv151203857L] (brown circles in the Figures) is measured within the aperture of $300$ kpc rather than within $R_{500}$. To compare with them, we report in the plot the simulated $M-T_{\rm sl,300kpc}$ with a red line. Since the temperature level of the ICM is mostly reached through the shocks generated by mergers and accretion during the formation of the cluster, the good match found in the two Figures implies that these processes are realistically reproduced in our models.

To investigate in a deeper manner the influence of the central regions, we compare the $L-T$ relation measured in the two categories of cool-core and non cool-core (Fig. \[fig:a1\]) clusters. The two simulated classes are taken from [@2015ApJ...813L..17R] and refer to the main halos of the 29 re-simulated regions. The distinction between the two types was established on the basis of the pseudo entropy level of the objects. We find reassurance in that the simulated clusters reproduce the same behavior of real systems: at a given temperature CC objects appear more luminous than NCC systems because of the higher central gas density that boosts their X-ray luminosity. The separation between CC and NCC is clearly reflected in our [[[agn]{}]{}]{} sample, where we also find a good reproduction of the observed $L-T$ slope. This represents a step forward with respect to earlier works. Indeed, previously, two situations emerged: either the self-regulated AGN mechanism preserved the CC structures by largely enhancing the luminosities of clusters colder than 3 keV and, therefore, failing on the $L-T$ relation [@gaspari.etal.2014; @hahn.etal.2015] or, viceversa, the $L-T$ was adequately created at the expenses of a reasonable baryon census in groups and clusters cores [@puchwein.etal.2008; @gaspari.etal.2014].

Evolution of ICM Scaling Relations
==================================

In this Section, we explore the evolution from $z=0$ up to $z=2$ of the six scaling relations: $M-{M_{\rm g}}$, $M-{T_{\rm mw}}$, $M-{T_{\rm sl}}$, $M-{Y_{\rm X}}$, $L-{T_{\rm sl}}$, and $L-M$. All the temperatures are obtained excising the central region ($<0.15$ $R_{500}$) even if excising the core produced a minimal variation. Indeed, the temperature difference is below 1% in our [[[agn]{}]{}]{} sample and below 2% in observational samples [@2012MNRAS.421.1583M] once we compare the temperature measured in the entire sphere within $R_{500}$ or excising the inner core region. Nevertheless, we decide to follow the standard choice, usually made to avoid the influence of the uncertainties related to the status of the core on the evolution study.

We apply the fitting procedure described in Section 3.4 and use the fitting function expressed by Eq. \[11\]. The pivot points are equal to the values of Table 1 and the power of evolution, $\gamma$, is here fixed to the self-similar expectations introduced in Section 3. The results, in terms of best-fitting parameters and 1 $\sigma$ uncertainties, are reported in Tables 2 and 3 for all relations, ICM physics treatments, and for different redshifts.

The Slope
---------

The evolution of the slopes of the scaling relations, $\beta$, is shown in Fig. \[fig:3\] for all six scaling relations and ICM physics.

[**$M-{M_{\rm g}}$ relation.**]{} The slope of the $M-{M_{\rm g}}$ relation is approximately constant over time for the [[[nr]{}]{}]{}  and [[[csf]{}]{}]{} runs confirming the expectations: the total mass and the gas mass grow simultaneously. Indeed, dark matter and gas increase in mass by the same fraction when the mass growth happens via slow accretion (due to constant ratio between the densities of the gas and DM components in the cluster outskirts, e.g., @rasia.etal.2004) or via major merger (due to small dependence of the gas fraction on the total mass, e.g., @2013MNRAS.431.1487P [@eckert.etal.2015], and reference therein). For the same reasons, no significant evolution is measurable in the [[[agn]{}]{}]{} runs when extending the analysis to $z\sim1$ (in agreement with previous analysis by @2011MNRAS.416..801F and @battaglia.etal.2013). However, the value of the [[[agn]{}]{}]{} slope, $\beta_{{M_{\rm g}}}$, decreases at $z=1.5$ and even more at $z=2$ to a maximum difference of $\sim13\%$ with respect to $z=0$. The cause of this variation relies on the scarse ability to retain the gas component by the lowest mass systems after the intense AGN activity that took place at $z>2$. On the other hand, systems that have a total gravitational mass of $M_{500}> 5 \times 10^{13} h^{-1} M_{\odot}$ at this high redshift had enough gravitational potential to accelerate the accretion of the necessary gas mass that brings the $z\sim 2$ most massive objects close to the $z \sim 0$ scaling relation. In other words, at redshifts $z > 1$ the gas fraction within $R_{500}$ evolves at different rates depending on the cluster mass. This causes the flattening of the $M-{M_{\rm g}}$ relation at $z=2$ with respect to $z=0$.

Finally, we confirm previous results that compared radiative and non-radiative runs and find that the slope of the radiative runs is shallower [@2010ApJ...715.1508S; @battaglia.etal.2013]. This is due to the effect of radiative cooling that converts part of the hot gas into stars and that is most efficient in low massive systems. In the local universe, the AGN balance the cooling thus the slope is closer to the gravitational-only model.

[**$M-{T_{\rm mw}}$ and $M-{T_{\rm sl}}$ relations.**]{} The slopes of both temperature relations, $M-{T_{\rm mw}}$ and $M-{T_{\rm sl}}$, drop at high redshifts. The values $\beta_{{T_{\rm sl}}}$ at $z=2$ are reduced with respect to the $z=0$ slope by $\sim23\%$, $\sim25\%$, and $\sim14\%$ in the [[[nr]{}]{}]{}, [[[csf]{}]{}]{}, and [[[agn]{}]{}]{} runs, respectively.

![image](fig5a){width="31.00000%"} ![image](fig5b){width="30.00000%"} ![image](fig5c){width="30.00000%"} ![image](fig5d){width="31.00000%"} ![image](fig5e){width="30.00000%"} ![image](fig5f){width="30.00000%"}

For either temperatures, the decrease in the $M-T$ slope does not strongly depend on the ICM physics. This indicates that the origin of the variation is due to macroscopic events, linked to the global evolution of the clusters.

To facilitate the explanation, we plot in Fig. \[fig:4\] the $M-{T_{\rm mw}}$ relation for the [[[nr]{}]{}]{} clusters at $z=0$ and $z=2$. As expected, high-redshift clusters are not thermalized yet and, therefore, they typically exhibit a lower value of temperature at fixed total mass, i.e. for $E(z) M < 2 \times 10^{14} h^{-1} M_{\odot}$ the $z=2$ temperatures (red points) are systematically on the left side of the $z=0$ ones (black points). The phenomenon is present in the entire mass range but it is less pronounced for the three most massive systems at $z=2$ that, indeed, lie extremely close to the black solid line representing the $z=0$ scaling relation. These 3 objects have recently experienced a major merger. As consequence, their mass is doubled (notice the mass separation from the rest of the sample), and their gas has been strongly heated by the induced shocks. To secure that our result was not influenced by the limited number of objects, we verified which fraction, among the highest mass systems in the $z \sim 2.3 $ MUSIC sample[^6] [@sembolini.etal.2013; @sembolini.etal.2014], has recently experienced a major merger. We found that this condition is verified for 10 systems out of the 11 objects with mass $E(z) M_{500} > 2 \times 10^{14} h^{-1} M_\odot$. We conclude that if such a massive system is already present at $z=2$ it is extremely likely (90% probability accordingly to the MUSIC sample) that it just went through a major merger phase that generated strong shocks to heat the gas. This eventuality is not surprising from a theoretical point of view since the increase steepness of the mass function at high redshift implies that it is extremely rare to grow a massive object via slow accretion. To summarize, we expect that while small clusters are not thermalized yet, the big ones can be already located in the $z=0$ scaling relation causing a shallower relation slope.

The change in slope is less prominent in the [[[agn]{}]{}]{} runs. In this circumstance, the gas of the smallest systems at $z=2$ is warmer because of the recent and intense AGN feedback activity. The phenomenon brings the smaller objects closer to the $z=0$ relation and reduces the amplitude of the slope evolution.

[**$M-{Y_{\rm X}}$ relation.**]{} The slope of the $M-{Y_{\rm X}}$ relation is constant until $z=1$ and then shows a mild decrease ($\leq 10\%$) for all the ICM runs consistent with previous results by [@sembolini.etal.2014] and by [@pike.etal.2014]. The origin of this variation can be understood by decomposing the $\beta_{{Y_{\rm X}}}$ slope into the two slopes of the $M-{M_{\rm g}}$ and $M-{T_{\rm sl}}$ [@2014MNRAS.437.1171M]: $$\beta_{{Y_{\rm X}}}=\frac{1}{1/\beta_{{M_{\rm g}}}+1/\beta_{{T_{\rm sl}}}}.$$ The complete derivation is presented in Appendix A. The mildness of the $\beta_{{Y_{\rm X}}}$ deviation is generated by the fact that none of the ICM physics runs present a coincident strong variation in both $\beta_{{M_{\rm g}}}$ and $\beta_{{T_{\rm sl}}}$. The [[[nr]{}]{}]{} and [[[csf]{}]{}]{} runs have the largest variations in $\beta_{{T_{\rm sl}}}$ but they have a constant $\beta_{{M_{\rm g}}}$, viceversa, the change in $\beta_{{M_{\rm g}}}$ for the [[[agn]{}]{}]{} run is accompanied by the mildest drop of the $\beta_{{T_{\rm sl}}}$ value. From the equation is also clear why the value of the $M-{Y_{\rm X}}$ slope is independent from the ICM prescriptions. Indeed, the [[[agn]{}]{}]{}  runs have lower $\beta_{{M_{\rm g}}}$ values than the [[[nr]{}]{}]{} simulations, because the [[[agn]{}]{}]{} smallest systems have less gas mass at fixed total mass but they also present an higher $\beta_{{T_{\rm sl}}}$ slope for the AGN extra heating that is evident in the group regime (see also discussion in @2011MNRAS.416..801F). This feature makes the $M-{Y_{\rm X}}$ relation the most suitable for cosmological studies that might include clusters presenting various astrophysical properties [see also @biffi.etal.2014].

[**$L-{T_{\rm sl}}$ and $L-M$ relations.**]{} The slopes of both luminosity-based relations are almost constant up to redshift $z=1$. At higher redshifts, $\beta_{\rm LT}$ decreases by $\sim 10\%$ for the [[[nr]{}]{}]{} and [[[csf]{}]{}]{} runs and $\beta_{\rm LM}$ grows by 20% in the [[[agn]{}]{}]{} case. These changes can be again explained by the decomposition of the luminosity-based relations’ slopes: $$\beta_{\rm LT}=\beta_{{T_{\rm sl}}}\bigg(\frac{2}{\beta_{{M_{\rm g}}}}-1\bigg)+\frac{1}{2},$$ $$\beta_{\rm LM}=\frac{2}{\beta_{{M_{\rm g}}}}+\frac{1}{2\beta_{{T_{\rm sl}}}}-1.$$ In the first case, $\beta_{{M_{\rm g}}}$ and $\beta_{{T_{\rm sl}}}$ carry a similar weight leading to a counterbalance between concurrent changes. In the second case, $\beta_{{M_{\rm g}}}$ is the dominant factor, implying a variation whenever this occurs in the $M-{M_{\rm g}}$ relation. Therefore, the steepening of the $L-M$ relation at $z>1$ is also due to the gas depletion within the potential well of small mass systems caused by the AGN ejecta at even higher redshifts.

Another notable feature in the luminosity relations is the separation of the slope values among the three physics. The slope has a self-similar behavior for the [[[nr]{}]{}]{} simulations while $\beta_{\rm LT}$ and $\beta_{\rm LM}$ increase by 50–80% and 20–50%, respectively, in the radiative runs. The change is consistent with the previous argument on the gas fractions: radiative processes reduce the amount of gas in the simulated smallest systems and, thus, strongly lower also their total luminosity similarly to what happens in real systems. The gap between the normalization is already well established at $z=2$ implying that at that epoch the stellar and AGN feedback in our simulated Universe already had a significant impact in establishing the energetics of clusters.

![Redshift evolution of the $M_{500}-{T_{\rm mw}}$ scaling relation. The result is shown for the NR simulation at redshifts $z=0$ (black) and $z=2$ (red). The solid lines are the best fit relations corresponding to each redshift.[]{data-label="fig:4"}](fig6){width="45.00000%"}

The Intrinsic Scatter
---------------------

![image](fig7a){width="32.00000%"} ![image](fig7b){width="31.00000%"} ![image](fig7c){width="31.00000%"}

We define the scatter, $\sigma$, as the mass variance in decimal logarithm at fixed signal (${M_{\rm g}}, {T_{\rm mw}}, {T_{\rm sl}}, {Y_{\rm X}}$) for the mass-proxy relations and as the luminosity variance in decimal logarithm at fixed ${T_{\rm sl}}$ and $M_{\rm tot}$ for the remaining relations. Since none of the mass$-$proxy scaling relations presents a significant trend with redshift, we limit the presentation in Fig. \[fig:5\] to the two cases with the lowest scatter: $\sigma_{{M_{\rm g}}}$ and $\sigma_{{Y_{\rm X}}}$. The luminosity-based relations, instead, present a variation of the scatter that decreases with redshift. We show in Fig. \[fig:5\] the scatter of the $L-M$ relation. The scatter values for the remaining relations can be found in the last columns of Tables 2 and 3.

The $M-{M_{\rm g}}$ relation always presents the minimum amount of scatter in line with the results from [@2010ApJ...715.1508S] and [@2010MNRAS.401.1670F]. The values are around $2-3\%$ which are 1.5 times smaller than $\sigma_{{Y_{\rm X}}}$ and a factor of $2$–$4$ smaller than the scatter of the two temperatures ($\sigma_{{T_{\rm mw}}}=0.05$ and $\sigma_{{T_{\rm sl}}} \sim 0.10$, respectively). The scatter of these four relations at $z=2$ is always consistent within 1 $\sigma$ with the $z=0$ scatter.

The two luminosity-related scaling relations $L-T$ and $L-M$ have the highest intrinsic scatter. They reach 0.2 at $z=0$ and decrease to $0.1$–$0.15$ at $z=2$. The largest variation affects the $L-M$ relation. The increase of the scatter in more recent times indicates that the luminosity is sensible to the entire merging history of the clusters and that the most significant deviations from the global scaling relation are originated by recent ($z<1$) massive mergers. A similar trend in the luminosity scatter was recently found by the Weight-the-Giant team [@mantz.etal.2016]. The variation of the $L-M$ scatter with redshift is an important factor that needs to be considered for cosmological studies. Likely, the change goes in the direction of reducing the Eddington bias caused by the different scattering of objects across the threshold of a flux limited sample [e.g., @stanek.etal.2006; @2012MNRAS.421.1583M]. Finally, radiative phenomena as stellar or AGN feedback have also effect onto diversifying the systems’ luminosity. Indeed, both the [[[csf]{}]{}]{} and [[[agn]{}]{}]{} runs show a 50% higher scatter than the [[[nr]{}]{}]{} simulations.

![Redshift evolution of the covariance matrix between the various signals of the [[[agn]{}]{}]{} run at three redshifts: $z=0$ (red asterisks), $z=1$ (orange squares), $z=2$ (brown diamonds). The Pearson correlation coefficients for $z=0,1,$ and $2$ are reported in each panel, respectively, from the left to the right.[]{data-label="fig:cov"}](fig8){width="45.00000%"}

As second step, we investigated the shape of the deviation of each signal, $\delta$, [*at fixed mass*]{} and the covariance matrix between two sets of deviations of two different signals. The scatter at fixed mass can always be accurately described by log-normal distributions [see also @lebrun.etal.2016] with a dispersion that is almost constant over time. The deviations from the scaling relations of the diverse signals are compared to each other in Fig. \[fig:cov\] where we report the results of the [[[agn]{}]{}]{} sample for $z=0,1,$ and 2. This analysis allows to identify which couples of signals have high correlation or anti-correlation ($|r| > 0.5$) and, therefore, can be used jointly to reduce the mass scatter with respect to the scatter of the individual scaling relations [@2010ApJ...715.1508S].

All the couples of X-ray quantities present positive correlation at all redshifts with the exception of no correlation found between $\delta_{{M_{\rm g}}}$ and $\delta_{{T_{\rm sl}}}$ and between $\delta_{L}$ and $\delta_{{T_{\rm sl}}}$ at $z=2$. The former behavior is consistent with the already-discussed argument about the different timescales on the variation of ${M_{\rm g}}$ and ${T_{\rm sl}}$ in reaction to mergers and accretion: since the two quantities increase at independent times we do not expect particular correlation. No correlation between $\delta_{{T_{\rm sl}}}$ and $\delta_L$ is also expected at high redshifts. Indeed, from one side, mergers produce an increase of both luminosity and temperature generating a positive correlation between the two deviations. From the other side, strong AGN bursts, common phenomena at these redshifts, are likely to cause an increase of temperature joint by a temporary [*decrease*]{} of luminosity, due to the gas lost as ejected material. This produces a negative correlation. The correlation value that we find between $\delta_L$ and $\delta_{{T_{\rm sl}}}$ at $z\leq1$ is in good agreement with the value $r=0.56\pm0.10$ found by [@mantz.etal.2016] in a $z<0.5$ sample comprehensive of both relaxed and unrelaxed clusters. In this observational paper, the authors claim a clean separation between the two classes of systems: CC relaxed systems tend to have higher $\delta_L$ and $\delta_{T}$ than the highly-disturbed NCC objects in their sample. In our simulations, we find that the six CC systems which are mostly X-ray regular have, indeed, the largest positive deviation on both quantities. However, we also have 5 CC systems that have negative variations in both quantities, particularly in temperature. Most of these objects have (colder) substructures that are entering in the cluster and are located in system outskirts. These increased the total mass of the object and lowered its overall ICM temperature.

The ${Y_{\rm X}}$ parameter shows a strong positive correlation with Pearson coefficient $r>0.5$. Clusters deviate from the $M-{Y_{\rm X}}$ mostly due to the changes in their temperature [@rasia.etal.2011] and this is confirmed by the fact that the highest correlation value is that between ${Y_{\rm X}}$ and ${T_{\rm sl}}$.

The luminosity is also firmly correlated with the other quantities. Finally, we confirm that a combination of the luminosity and temperature ($M\propto L\times T$) as proposed in [@ettori2013] is a good approach to reduce the scatter with respect to the scatters of the individual relations $M-L$ and $M-T$.

The Normalization
-----------------

![image](fig9a){width="40.00000%"} ![image](fig9b){width="40.00000%"}

The evolution of the normalization is mostly expressed as a power ${\gamma}$ of the factor $E(z)$ [@mantz.etal.2016]. In this Section, we want to enlighten how this procedure is inaccurate [*whenever the slope of the considered scaling relation varies with time*]{} [@branchesi.etal.2007]. In this case, indeed, it is not possible to quantify any change in the normalization because the relation does not move in a parallel way in the scaling-relation plane.

From observations, we still do not have any indication of evolution of $\beta$ for the paucity of observed $z\sim 1$ clusters. The situation will improve thanks to the collection of high-z data from SPT-3G [@benson.etal.2014] or e-Rosita [@2013SSRv..177..247G] and thanks also to the development of Bayesian techniques (see LIRA by @sereno_ettori2015 and @sereno2016 and also @andreon.etal.2012). To provide forecast on the precision with which these missions will measure cosmological parameters, there are studies that perform the analysis of samples including clusters at different redshifts. The usual approach is to fix the slope either to the local $\beta$ (e.g. @2014MNRAS.444.2723C) or to simultaneously fit low-, intermediate- and high-redshift objects (e.g. @2012MNRAS.424.2086H, @2015arXiv151203833G) letting $\beta$ free to vary.

To warn about the application of these methods we focus on the $M-{M_{\rm g}}$ and $L-M$ relations because they exhibit the largest variation in their slopes and, therefore, they are the best suited to evidence the possible mis-interpretation of the data. We build a sample that includes all the [[[agn]{}]{}]{} simulated clusters from $z=0$ to $z=2$ and following the observational approach we fit this combined sample with the Eqn. (\[11\]) where the evolution factor $\gamma$, the slope $\beta$ and the normalization $C$ are all free to vary and the $E(z)$ factor is appropriately considered. With this procedure, we are therefore implying that the slope value is independent from the redshift. The best-fit parameters and their 1$\sigma$ errors are reported in Table \[tb8\].

In Fig.\[fig:7\], we focus on 4 redshift bins ($z=0$, 0.5, 1, 1.5, and 2) of the $M-{M_{\rm g}}$ and $L-M$ best-fit relations of the [[[agn]{}]{}]{} runs and we report the difference between Eq. 11 expressed through the parameters reported in Table 4 and the redshift-related factor $E(z)$) and Eq. 11 expressed with the redshift-dependent values listed in Tables 2 and 3. From the Figure, it is clear that the results obtained by enforcing the same slope can significantly deviate from the results obtained in the limited redshift bin. The most severe deviations are present at higher redshifts and for the $L-M$ relation. The $M-{M_{\rm g}}$ and similarly the $M-{Y_{\rm X}}$ are mostly deviant at $z=2$ but do not vary substantially up to $z=1$ as expected since their $\beta$ value is mostly constant that redshift bin.

For the cases where the slope does not vary with time, such as $M-{T_{\rm mw}}$ and $L-T$, the evolution in normalization is a well defined measurement and the usage of a unique fitting function is well justified. In particular, we find that the predicted slope and evolution of both relations are close to the self-similar expectations with the slight tendency of less negative evolution in the $M-T$ relation ($\gamma = -0.9$) and more positive evolution in the $L-T$ relation ($\gamma=1.2$). In particular, the $L-T$ from our [[[agn]{}]{}]{} sample, $E(z)^{1.28 \pm 0.12}$ and $\beta=2.88$, is in line with recent results from [@2015arXiv151203833G] who also reported a positive evolution $E(z)^{1.64\pm0.77}$ associated with the slope $\beta =3.08 \pm 0.15$. This result is, however, in contrast with the measurements of the bolometric luminosity, $L_{\rm bol}\sim5$–$10 \times 10^{44}$ erg s$^{-1}$ of ISCS J1438+3414 at $z=1.4$ and JKCS041 at $z=2.2$ by [@andreon.etal.2011], of XDCP J0044.0-2033 at $z=1.579$ by [@tozzi.etal.2015], and of IDCS J1426.5+3508 at $z=1.75$ by [@brodwin.etal.2016]. If local relations and self-similar evolution are assumed, all these objects appeared under-luminous for their measured temperature ($4.9^{+3.4}_{-1.6}$ keV, $7.3^{+6.7}_{-2.6}$ keV, $6.7^{+1.3}_{0.9}$ keV, and $7.6^{+8.7}_{-1.9}$ keV, respectively) or viceversa hotter for their luminosity. It is difficult to make a significant comparison for low statistics of objects, the large error bars, and the impossibility to evaluate any selection effect. It is, despite, interesting to notice that all these cases are [*exceptionally*]{} massive and almost twice as hot as our most massive clusters at $z\ge 1.5$.

                               $\log C$             $\beta$            $\gamma$            $\sigma$
  ------------------ ------------------ ------------------- ------------------- -------------------
   $M-{M_{\rm g}}$     $14.208\pm0.002$   $ 0.916\pm 0.002$   $-0.047\pm 0.011$   $ 0.025\pm 0.001$
   $M-{T_{\rm mw}}$    $14.297\pm0.003$   $ 1.623\pm 0.003$   $-0.902\pm 0.017$   $ 0.045\pm 0.001$
   $M-{T_{\rm sl}}$    $14.297\pm0.005$   $ 1.661\pm 0.005$   $-0.847\pm 0.027$   $ 0.069\pm 0.002$
      $M-Y_{X}$        $14.133\pm0.002$   $ 0.597\pm 0.002$   $-0.314\pm 0.013$   $ 0.032\pm 0.001$
   $L-{T_{\rm sl}}$    $ 0.551\pm0.012$   $ 2.877\pm 0.012$   $ 1.161\pm 0.063$   $ 0.162\pm 0.005$
        $L-M$          $ 1.228\pm0.015$   $ 1.673\pm 0.015$   $ 2.520\pm 0.077$   $ 0.176\pm 0.005$

Summary and Discussion
======================

In this paper, we address the evolution of scaling relations in simulated galaxy clusters. Our study is based on a large sample of clusters, simulated in a cosmological context. The simulations are carried out with an up-graded implementation of SPH and with an improved version of the AGN physics. We first examine the reliability of our newly performed [[[agn]{}]{}]{} runs by comparing our predictions to those derived from local ($z<0.25$) and intermediate ($0.25<z< 0.6$) observations. Subsequently, we investigate the evolution of six scaling relations: $M-{M_{\rm g}}$, $M-{T_{\rm mw}}$, $M-{T_{\rm sl}}$, $M-Y_X$, $L-{T_{\rm sl}}$, and $L-M$ by comparing the [[[agn]{}]{}]{} model with two other parallel runs performed with non-radiative physics or radiative ICM but without AGN feedback. We characterize how the scaling relation features (namely the slope, intrinsic scatter, and normalization) and the covariance matrix between couples of signals change as a function of redshift. We summarize in the following our main results.

-   The simulated scaling relations at low and intermediate redshifts reproduce well the observations. Distinctly, the $L-T$ relation is well reproduced by the [[[agn]{}]{}]{} runs that are able to reproduce its slope in a mass range from small groups to clusters and the observed separation between CC and NCC clusters.

-   From $z=0$ and $z=1$, we do not detect any appreciable change on the slopes of the relations. At higher redshifts, however, all the relations exhibit some degree of evolution with the only exception of the luminosity-temperature slope which remains unchanged. In the [[[agn]{}]{}]{} runs, the gas slope $\beta_{{M_{\rm g}}}$ at $z=2$ is reduced by $\sim13\%$ with respect to the present-time value. This is caused by the effect of intense high-z AGN activity that has more impact on the lowest mass systems. At $z=2$, a shallower slope is found also for the $M-T$ relation that declines by $\sim15\%$ with respect to $z=0$. At these high redshifts, the smallest groups have a systematic low temperature due to their incomplete thermalization. The decrease of $Y_X$ ($\sim10\%$) and increase of $L-M$ ($\sim20\%$) slopes can be explained by analytically decomposing their slopes in the two contributors: $\beta_{{M_{\rm g}}}$ and $\beta_{{T_{\rm sl}}}$.

-   We do not find any significant trend in redshift for the scatter of all the mass-proxy relations. Consistent with previous theoretical studies, the $M-{M_{\rm g}}$ relation has the smallest scatter around $2-3\%$. Instead, the scatter of the two luminosity relations $L-{T_{\rm sl}}$ and $L-M$ is the largest and around $15-20\%$, over the redshift range $[0-2]$. The $L-M$ scatter increases with the decrease of redshift enlightening the significant impact on the X-ray luminosity by recent major mergers [@torri.etal.2004]. When a merger occurs the luminosity registers a permanent increase [@rowley.etal.2004]. The scatter of all relations can be well-described by a log-normal distribution whose widths are mostly constant over the redshift up to $z\sim1.5$.

-   No correlation is evident between the pair of deviation in ${M_{\rm g}}$ and ${T_{\rm sl}}$ at fixed mass and at all redshift. No correlation is registered also between $L$ and ${T_{\rm sl}}$ at $z=2$. In all other cases, positive correlations are found with Pearson coefficients always greater than 0.4. Interestingly, the correlation coefficients are almost constant from $z=0$ to $z=1$ assuring that the reasons causing coupled deviations from the scaling relations are already in place at high redshift.

-   Regarding the study of the evolution of the normalization, we stress that in situations where the slopes vary with redshift the evolution on the normalization cannot be established. On the other hand, we find that the $M-{T_{\rm mw}}$ and $L-{T_{\rm sl}}$ relations, whose slopes do not have any evolution, exhibit a negative and positive evolution of the normalization for the redshift range $[0-2]$, with values for the evolution parameters in line with recent observational studies and close to the self-similar predictions.

-   Overall, we confirm that the $M-{Y_{\rm X}}$ relation evaluated from $z=0$ to $z=1$ is the best suited for cosmological studies for the combination of its properties: in that redshift range, the slope does not vary, the evolution of the normalization can be robustly determined, the scatter is small and constant over time and most importantly the relation is solid and independent from the source of feedback either from star, supernovae, or AGN.

On the basis of our analysis of the intrinsic variations of simulated clusters, we also conclude that pushing the study of the scaling relations at higher redshifts does not seem to be an advantage because the intense AGN activity, peaking at $z \sim 2$, could have a significant impact and produce deviations from the self-similar behavior at redshifts $z>1$. This cosmic epoch is still almost an unexplored territory where the predictions of higher-resolution simulations can help in designing observational strategies for future missions such as e-Rosita and Athena. From an observational perspective, the scaling relations are expected to be calibrated by measuring the mass from weak-lensing analysis [@marrone.etal.2012; @hoekstra.etal.2015; @sereno_ettori2015; @mantz.etal.2016]. Even if this procedure is not expected to introduce mass bias [@meneghetti.etal.2010; @becker_kravtsov; @rasia.etal.2012] it will likely enlarge the scatter of the relations [@sereno_ettori.2015] and attention will need to be places for clusters close to the flux limit threshold [@nord.etal.2008]. Certainly, more efforts will need to be devoted to reducing the uncertainties on the mass calibration to few percent level.

ACKNOWLEDGEMENTS {#acknowledgements .unnumbered}
================

We are greatly indebted to Volker Springel for giving us access to the developer version of the GADGET3 code and to the entire MUSIC team to give us access to the simulated data for the comparison shown in Section 5.1. We would like to thank Gustavo Yepes, Stefano Ettori, and Mauro Sereno for providing useful discussions and comments. We acknowledge support from: PIIF-GA- 2013-627474 (E.R.); NSF AST-1210973 (E.R.); the [*Spanish Ministerio de Econom[í]{}a y Competitividad*]{}, MINECO, grant AYA2013-48226-C3-2-P (S.P.); the Generalitat Valenciana, grant GVACOMP2015-227 (S.P.); the [*Slovenian Research Agency*]{} (D.F.); the DFG Cluster of Excellence “Origin and Structure of the Universe” (A.M.B.); the DFG Research Unit 1254 “ Magnetisation of interstellar and intergalactic medium” (A.M.B.); CONICET, FonCyT and SeCyT-UNC, Argentina (C.R.F.). This work is also supported by the PRIN-MIUR 201278X4FL “Evolution of cosmic baryons” funded by the Italian Ministry of Research, by the PRIN-INAF 2012 grant “The Universe in a Box: Multi-scale Simulations of Cosmic Structures", by the INDARK INFN grant, by the “Consorzio per la Fisica di Trieste” Simulations were performed using Flux HCP Cluster at the University of Michigan and the Galileo Machine at CINECA with CPU time assigned through ISCRA proposals and an agreement with the University of Trieste. Data, post-processing, and storage has been done on the CINECA facility PICO, granted to us thanks to our expression of interest.

Decomposition of the slopes of composited signals
=================================================

In this section we derive the features of $M-{Y_{\rm X}}$, $L-{T_{\rm sl}}$, and $L-M$ based on the two fundamental relations: $M-{{M_{\rm g}}}$ and $M-{T_{\rm sl}}$:

$$M \propto E(z)^{\gamma_{{M_{\rm g}}}} \times {M_{\rm g}}^{\beta_{{M_{\rm g}}}}  \to  {M_{\rm g}}\propto E(z)^{-\frac{\gamma_{{M_{\rm g}}}}{\beta_{{M_{\rm g}}}}}\times M^{\frac{1}{\beta_{{M_{\rm g}}}}}$$

$$M \propto E(z)^{\gamma_{{T_{\rm sl}}}}\times {{T_{\rm sl}}}^{\beta_{{T_{\rm sl}}}}  \to  {T_{\rm sl}}\propto E(z)^{-\frac{\gamma_{{T_{\rm sl}}}}{\beta_{{T_{\rm sl}}}}}\times M^{\frac{1}{\beta_{{T_{\rm sl}}}}}.$$

[**The $M-{Y_{\rm X}}$ relation.**]{} By definition the ${Y_{\rm X}}$ is the product of gas mass and temperature, and thus $$Y_{X} = {M_{\rm g}}\times {T_{\rm sl}}\propto E(z)^{-\gamma_{{M_{\rm g}}}-\gamma_{{T_{\rm sl}}}}\times M^{-\beta_{{M_{\rm g}}}}.$$ From eqns. (1)-(3), one deduces $$\begin{aligned}
Y_{X}\propto M^{\frac{1}{\beta_{{M_{\rm g}}}}}\times E(z)^{-\frac{\gamma_{{M_{\rm g}}}}{\beta_{{M_{\rm g}}}}}\times M^{\frac{1}{\beta_{{T_{\rm sl}}}}}\times E(z)^{-\frac{\gamma_{{T_{\rm sl}}}}{\beta_{{T_{\rm sl}}}}}.\end{aligned}$$ Writing $M-Y_X$ in general form $$M = C_{{Y_{\rm X}}}\times E(z)^{\gamma_{{Y_{\rm X}}}}\times\bigg(\frac{Y_X}{Y_{X0}}\bigg)^{\beta_{{Y_{\rm X}}}},$$ the slope and power of evolution of the $M-Y$ are given by $$\beta_{{Y_{\rm X}}}=\frac{1}{1/\beta_{{M_{\rm g}}}+1/\beta_{{T_{\rm sl}}}},$$ $$\gamma_{{Y_{\rm X}}} = \frac{\gamma_{{M_{\rm g}}}/\beta_{{M_{\rm g}}}+\gamma_{{T_{\rm sl}}}/\beta_{{T_{\rm sl}}}}{1/\beta_{{M_{\rm g}}}+1/\beta_{{T_{\rm sl}}}}.$$ [**The $L-{T_{\rm sl}}$ relation.**]{} The general form of $L-{T_{\rm sl}}$ is given by $$L = C_{LT}\times E(z)^{\gamma_{LT}}\times \bigg(\frac{{T_{\rm sl}}}{T_0}\bigg)^{\beta_{LT}}.$$ Recalling that $L$ is related to the total mass and temperature via eqn. \[8\] $$\begin{aligned}
L &\propto& E(z)^2f_g^2\times M\times T_{sl}^{1/2}\nonumber\\
&\propto& E(z)^2\times\frac{M_g^2}{M}\times T_{sl}^{1/2}\nonumber\\
&\propto& E(z)^2\times E(z)^{-\frac{2\gamma_{M_g}}{\beta_{M_g}}}\times M^{\frac{2}{\beta_{M_g}}-1}\times T_{sl}^{1/2}.\end{aligned}$$ In the last derivation we use eqn. (1) to rewrite $M_g$ in term of $M$. In turn, by substituting $M$ in term of $T_{sl}$ via eqn. (2) one obtains $$\begin{aligned}
L&\propto& E(z)^2\times E(z)^{-\frac{2\gamma_{M_g}}{\beta_{M_g}}}\times E(z)^{\gamma_{T_{sl}}\times(\frac{2}{\beta_{Mg}}-1)}\times T_{sl}^{\beta_{T_{sl}}(\frac{2}{\beta_{M_g}}-1)}\nonumber\\
&\propto& E(z)^{2(1-\frac{\gamma_{M_g}}{\beta_{M_g}})+\gamma_{T_{sl}}(\frac{2}{\beta_{M_g}}-1)}\times T_{sl}^{\beta_{T_{sl}}(\frac{2}{\beta_{M_g}}-1)+\frac{1}{2}}.\end{aligned}$$ By comparing eqn. (10) to eqn. (8) one deduces $$\beta_{LT}=\beta_{{T_{\rm sl}}}\bigg(\frac{2}{\beta_{{M_{\rm g}}}}-1\bigg)+\frac{1}{2},$$ $$\gamma_{LT}=2\bigg(1-\frac{\gamma_{{M_{\rm g}}}}{\beta_{{M_{\rm g}}}}\bigg)+\gamma_{{T_{\rm sl}}}\bigg(\frac{2}{\beta_{{M_{\rm g}}}}-1\bigg).$$

[**The $L-M$ relation.**]{} We define the general form of $L-M$ as $$L = C_{LM}\times E(z)^{\gamma_{LM}}\times \bigg(\frac{M}{M_0}\bigg)^{\beta_{LM}}.$$ In eqn. (9), we rewrite $T_{sl}$ in term of $M$ via eqn. (2) and obtain $$\begin{aligned}
L&\propto& E(z)^{2}\times E(z)^{\frac{-2\gamma_{M_g}}{\beta_{M_g}}}\times M^{\frac{2}{\beta_{M_g}}-1}\times E(z)^{-\frac{\gamma_{T_{sl}}}{2\beta_{T_{sl}}}}\times M^{\frac{1}{2\beta_{T_{sl}}}}\nonumber\\
L&\propto& E(z)^{2(1-\frac{\gamma_{M_g}}{\beta_{M_g}})-\frac{\gamma_{T_{sl}}}{2\beta_{T_{sl}}}}\times M^{\frac{2}{\beta_{M_g}}+\frac{1}{2\beta_{T_{sl}}}-1}.\end{aligned}$$ From the last equation we can derive the slope and evolution of the $LM$ relation given as $$\beta_{LM}=\frac{2}{\beta_{{M_{\rm g}}}}+\frac{1}{2\beta_{{T_{\rm sl}}}}-1,$$ $$\gamma_{LM}=2\bigg(1-\frac{\gamma_{{M_{\rm g}}}}{\beta_{{M_{\rm g}}}}\bigg)-\frac{1}{2}\frac{\gamma_{{T_{\rm sl}}}}{\beta_{{T_{\rm sl}}}}.$$

[^1]: http://www.mpe.mpg.de/eROSITA

[^2]: http://www.the-athena-x-ray-observatory.eu

[^3]: http://www.euclid-ec.org

[^4]: http://www.lsst.org

[^5]: The gas particles can be multi-phase, carrying information on both the hot and cold gas. The cold phase, not contributing to the gas mass, provide a reservoir for stellar formation.

[^6]: http://music.ft.uam.es
