The redshift evolution of massive galaxy clusters in the MACSIS simulations
Abstract
We present the MAssive ClusterS and Intercluster Structures (MACSIS) project, a suite of 390 clusters simulated with baryonic physics that yields realistic massive galaxy clusters capable of matching a wide range of observed properties. MACSIS extends the recent BAHAMAS simulation to higher masses, enabling robust predictions for the redshift evolution of cluster properties and an assessment of the effect of selecting only the hottest systems. We study the observablemass scaling relations and the Xray luminositytemperature relation over the complete observed cluster mass range. As expected, we find the slope of these scaling relations and the evolution of their normalization with redshift departs significantly from the selfsimilar predictions. However, for a sample of hot clusters with coreexcised temperatures the normalization and slope of the observablemass relations and their evolution are significantly closer to selfsimilar. The exception is the temperaturemass relation, for which the increased importance of nonthermal pressure support and biased Xray temperatures leads to a greater departure from selfsimilarity in the hottest systems. As a consequence, these also affect the slope and evolution of the normalization in the luminositytemperature relation. The median hot gas profiles show good agreement with observational data at and , with their evolution again departing significantly from the selfsimilar prediction. However, selecting a hot sample of clusters yields profiles that evolve significantly closer to the selfsimilar prediction. In conclusion, our results show that understanding the selection function is vital for robust calibration of cluster properties with mass and redshift.
keywords:
galaxies: clusters: general  galaxies: clusters: intracluster medium  Xrays: galaxies: clusters  galaxies: evolution  methods: numerical  hydrodynamics1 Introduction
Galaxy clusters form from large primordial density fluctuations that have collapsed and virialised by the present epoch, with more massive clusters forming from larger and rarer fluctuations. This makes them especially sensitive to fundamental cosmological parameters, such as the matter density, the amplitude of the matter power spectrum and the equation of state of dark energy (see Voit, 2005; Allen et al., 2011; Kravtsov & Borgani, 2012; Weinberg et al., 2013). The observable properties of a galaxy cluster result from a nontrivial interplay between gravitational collapse and astrophysical processes. The diverse range of formation histories of the cluster population leads to scatter in the observablemass scaling relations and, as surveys select clusters based on an observable, this can lead to a biased sample of clusters, resulting in systematics when using them as a cosmological probe (e.g. Mantz et al., 2010). Many previous studies have shown that the relationship between a cluster observable, such as its temperature or Xray luminosity, and a quantity of interest for cosmology, e.g. its mass, has a smaller scatter for more massive, dynamically relaxed objects (Eke et al., 1998; Kay et al., 2004; Crain et al., 2007; Nagai et al., 2007b; Planelles et al., 2013). Therefore, the fundamental requirement when probing cosmological parameters with galaxy clusters is a sample of relaxed, massive clusters with well calibrated massobservable scaling relations.
However, galaxy clusters are rare objects, becoming increasingly rare with increasing mass, and to observe a sample large enough to be representative of the underlying population requires a survey with significant size and depth. Currently ongoing and impending observational campaigns, such as the Dark Energy Survey (The Dark Energy Survey Collaboration, 2005), eRosita (Merloni et al., 2012), Euclid (Laureijs et al., 2011), SPT3G (Benson et al., 2014) and Advanced ACTpol (Henderson et al., 2016), will be the first to have sufficient volume to yield significant samples of massive clusters. Due to their rarity, the majority of these massive clusters will be at high redshift and it is therefore critical to understand how the cluster observables and their associated scatter evolve. Additionally, the most massive clusters will be the brightest and easiest to detect objects at high redshift, making it vital to understand the selection function of the chosen cluster observable and whether the most massive clusters are representative of the underlying cluster population. Theoretical modelling of the formation of clusters and their observable properties is required to understand these issues and to further clusters as probes of cosmology. Due to the range of scales involved in cluster formation, the need to incorporate astrophysical processes and to selfconsistently predict observable properties, cosmological hydrodynamical simulations are the only viable option.
Recent progress in the modelling of largescale structure formation has been driven mainly by the inclusion of supermassive black holes and their associated Active Galactic Nucleus (AGN) feedback, which has been shown to be critical for reproducing many cluster properties (Bhattacharya et al., 2008; Puchwein et al., 2008; McCarthy et al., 2010; Fabjan et al., 2010). A number of independent simulations are now able to produce realistic clusters that simultaneously reproduce many cluster properties in good agreement with the observations (Le Brun et al., 2014; Pike et al., 2014; Planelles et al., 2014). Results from the recent BAryons and HAloes of MAssive Systems (BAHAMAS) simulations (McCarthy et al., 2016) have shown that by calibrating the subgrid model for feedback to match a small number of key observables, in this case the global galaxy stellar mass function and the gas fraction of clusters, simulations of largescale structure are now able to reproduce many observed scaling relations and their associated scatter over two decades in halo mass. However, full gas physics simulations of largescale structure formation, with sufficient resolution, are still computationally expensive. This has limited previous studies to either small samples with objects or to volumes of , all of which are too small to contain the representative sample of massive clusters that is required for cosmological studies above .
This paper introduces the Virgo consortium’s MACSIS project, a sample of massive clusters selected from a large volume dark matter simulation and resimulated with full gas physics to enable selfconsistent observable predictions. The simulations extend the BAHAMAS simulations to the most massive clusters expected to form in a cosmology. In this paper we study the cluster scaling relations and their evolution. We combine the MACSIS and BAHAMAS simulations to produce a sample that spans the complete mass range and that can be studied to high redshift, using the progenitors of the MACSIS sample. We also select the hottest clusters from the combined sample and a relaxed subset of them to examine the impact of such selections on the scaling relations and their evolution. We then study the gas profiles to further understand the differences between the samples.
This paper is organised as follows. In Section 2 we introduce the MACSIS sample and discuss the parent dark matter simulation from which the sample was selected, the selection criteria used, the model used to resimulate the haloes, how we produced the observable quantities and the three samples we use in this work. In Section 3 we investigate how the scaling relations evolve and how this evolution changes when a hot cluster sample or relaxed, hot cluster sample is selected. We then study the hot gas profiles to understand the differences in the evolution of the relations for the different samples in Section 4. Finally, in Section 5 we discuss our results and summarise our main findings.
2 Parent simulation and sample selection
In this section we describe the parent simulation, the selection of the MACSIS sample, the baryonic physics used in the resimulation of the sample and the calculation of the observable properties of the resimulated clusters. Additionally, we describe how MACSIS and BAHAMAS clusters were selected to produce the combined sample and the cuts made to yield a hot sample and its relaxed subset.
2.1 The parent simulation
To obtain a population of massive clusters we require a simulation with a very large volume . With current computational resources it is unfeasible to simulate such a volume with hydrodynamics and the required gas physics, such as radiative cooling, star formation and feedback, at a resolution high enough to accurately capture the cluster properties. An alternative option is to apply the zoomed simulation technique to a representative sample of objects from a larger volume. Therefore, we select a sample of massive haloes from a dark matter only simulation that has sufficient volume to yield a population of massive clusters and the resolution to ensure they are well characterized. We label this simulation the ‘parent’ simulation.
The parent simulation is a periodic cube with a side length of . Its cosmological parameters are taken from the Planck 2013 results combined with baryonic acoustic oscillations, WMAP polarization and high multipole moments experiments (Planck Collaboration et al., 2014a) and are , , , , , and . We note that there are minor differences between these values and the Planckonly cosmology used for the BAHAMAS simulations, but this has negligible impact on the results presented here. The simulation contained dark matter particles that were arranged in an initial glasslike configuration and then displaced according to secondorder Lagrangian perturbation theory using the ic_2lpt_gen code (Jenkins, 2010) and the public Gaussian white noise field Panphasia (Jenkins, 2013; Jenkins & Booth, 2013). ^{1}^{1}1The phase descriptor for this volume is [Panph1, L14, (2152, 5744, 757), S3, CH1814785143, EAGLE_L3200_VOL1]. The particle mass of this simulation is and the comoving gravitational softening length was set to . The simulation was evolved from redshift using a version of the Lagrangian TreePMSPH code pgadget3 (last described in Springel, 2005). Haloes were identified at using a FriendsofFriends (FoF) algorithm with a standard linking length of in units of the mean interparticle separation (Davis et al., 1985).
We plot the FoF mass function of the parent simulation at in Fig. 1. We compare it to the published relations of Jenkins et al. (2001), Angulo et al. (2012), Watson et al. (2013) and Heitmann et al. (2015). We plot the scaled differential mass function
(1) 
where is halo mass, is the mean density of the Universe at , is the number of haloes per unit volume, and is the variance of the linear density field when smoothed with a tophat filter. We plot the mass function as a function of the variable as it is insensitive to cosmology (Jenkins et al., 2001). For we find that all of the mass functions show reasonable agreement with differences of between them, with the small differences likely due to the mass function not being exactly universal (Tinker et al., 2008; Courtin et al., 2011). However, for larger values the mass functions begin to diverge, as the parent simulation has an excess of massive clusters compared to the other simulations. This is likely due to two effects. First, the MACSIS simulation is the only one to use when generating the initial conditions. It has been shown that not using results in a significant underestimation of the abundance of the rarest objects (Crocce et al., 2006; Reed et al., 2013). The second effect is simply statistics: even in a very large volume there are still low numbers of the rarest and most massive clusters, where there is likely to be significant variance between the simulation volumes.
2.2 The MACSIS sample
To select the MACSIS sample, all haloes with were grouped in logarithmically spaced bins , with . If a bin contained less than one hundred haloes then all of the objects in that bin were selected. For bins with more than one hundred objects the bin was then further subdivided into bins of 0.02 dex and ten objects from each subbin were then selected at random. The subdividing of the bins ensured that our random selection was not biased to low masses by the steep slope of the mass function. This selection procedure results in a sample of haloes that is mass limited above and randomly sampled below this limit. Table 1 shows the fraction of haloes selected from the parent simulation in each mass bin. We have compared the properties of the selected haloes with those of the underlying population and found the MACSIS sample to be representative. Additionally, in Appendix A we demonstrate that selecting by a halo’s FoF mass does not bias our results when binning clusters by their .
Mass Bin  Sample  Total  Fraction 

Size  Haloes  Selected  
Due to current computational constraints the BAHAMAS simulations are limited to periodic cubes with a side length of . There are very few clusters with a mass greater than in a volume of this size, and those that are present may be affected by the loss of power from largescale modes that are absent due to their wavelengths being greater than the box size. The zoom simulations of the MACSIS project provide an extension to the BAHAMAS periodic simulations. They provide the most massive clusters and allow the massobservable scaling relations to be studied across the complete cluster mass range.
We use the zoomed simulation technique (Katz & White, 1993; Tormen et al., 1997) to resimulate the chosen sample at increased resolution. We perform both DM only and full gas physics resimulations. The Lagrangian region for every cluster was selected so that its volume was devoid of lower resolution particles beyond a cluster centric radius of .^{2}^{2}2We define as the radius at which the enclosed average density is two hundred times the critical density of the Universe. The resolution of the Lagrangian region was increased such that the particles in the DM only simulations had a mass of and in the hydrodynamic resimulations the dark matter particles had a mass of and the gas particles had an initial mass of . In all simulations the Plummer equivalent gravitational softening length for the highresolution particles was fixed to in comoving units for and in physical coordinates thereafter. The smoothed particle hydrodynamics interpolation used neighbours and the minimum smoothing length was set to one tenth of the gravitational softening. A schematic view of the zoom approach is shown in Fig. 2.
The resolution and softening of the zoom resimulations were deliberately chosen to match the values of the periodic box simulations of the BAHAMAS project (McCarthy et al., 2016), which is a calibrated version of the OWLS code (Schaye et al., 2010), which was also used for cosmoOWLS (Le Brun et al., 2014). The subgrid models for feedback from star formation and AGN used in the BAHAMAS simulations was calibrated to obtain a good fit to the observed galaxy stellar mass function and the amplitude of the gas fractiontotal mass relation, respectively, at . Without any further tuning, the simulations then produce a population of groups and clusters that shows excellent agreement with the observations for a range of galaxyhalo, hot gashalo and galaxyhot gas relations.
2.3 Baryonic physics
The BAHAMAS simulations were run with a version of pgadget3 that has been heavily modified to include new subgrid physics as part of the OWLS project (Schaye et al., 2010). We now briefly describe the subgrid physics, but refer the reader to Schaye et al. (2010), Le Brun et al. (2014) and McCarthy et al. (2016) for greater detail, including the impact of varying the free parameters in the model and the calibration strategy. Radiative cooling is calculated on an elementbyelement basis following Wiersma et al. (2009a), interpolating the rates as a function of density, temperature and redshift from precomputed tables generated with cloudy (Ferland et al., 1998). It accounts for heating and cooling due to the primary cosmic microwave background and a Haardt & Madau (2001) ultraviolet/Xray background. The background due to reionization is assumed to switch on at .
Star formation is modelled stochastically in a way that by construction reproduces the observations, as discussed in Schaye & Dalla Vecchia (2008). Lacking the resolution and physics to correctly model the cold interstellar medium, gas particles with a density follow an imposed equation of state with . These gas particles then form stars at a pressuredependent rate that reproduces the observed KennicuttSchmidt law (Schmidt, 1959; Kennicutt, 1998). Stellar evolution and the resulting chemical enrichment are implemented using the model of Wiersma et al. (2009b), where chemical elements (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) are followed. The mass loss rates are calculated assuming Type Ia and Type II supernovae and winds from massive and asymptotic giant branch stars. Stellar feedback is implemented via the kinetic wind model of Dalla Vecchia & Schaye (2008). The BAHAMAS simulations used the calibrated massloading factor of and wind velocity . This corresponds to percent of available energy from Type II supernovae, assuming a Chabrier (2003) IMF, and yields an excellent fit to the observed galaxy mass function.
The seeding, growth and feedback from supermassive black holes (BH) is implemented using the prescription of Booth & Schaye (2009), a modified version of the method developed by Springel et al. (2005). A FoF algorithm is run onthefly and BH seed particles, with , are placed in haloes that contain at least DM particles, which corresponds to a halo mass of . BHs grow via Eddingtonlimited accretion of gas at the BondiHoyleLittleton rate, with a boost factor that is a powerlaw of the local density for gas above the star formation density threshold. They also grow by direct mergers with other BHs. A fraction, , of the rest mass energy of the accreted gas is then used to heat neighbour particles by increasing their temperature by . Changes to these parameters have a significant impact on the hot gas properties of clusters. The calibrated values of these parameters in the BAHAMAS simulations are and . The feedback efficiency , where is the radiative efficiency and is the fraction of that couples to the surrounding gas. The choice of the efficiency, assuming it is nonzero, is generally of little consequence as the feedback establishes a selfregulating scenario, but determines the black hole masses (Booth & Schaye, 2009).
2.4 Calculating observable properties
Previous studies have shown that there can be significant biases in the observable properties of clusters due to issues such as multitemperature structures and gas inhomogeneities (e.g. Nagai et al., 2007a; Khedekar et al., 2013). Therefore, when investigating cluster properties it is critical that, as far as possible, we make a likewithlike comparison with the observations. Following Le Brun et al. (2014), we do this by producing synthetic observational data for each cluster and analysing it in a manner similar to what is done for real data. Using the particle’s temperature, density and metallicity, where the metallicity is smoothed over a particle’s neighbours, we first compute a restframe Xray spectrum in the band for all gas particles, using the Astrophysical Plasma Emission Code (apec; Smith et al., 2001) via the pyatomdb module with atomic data from atomdb v3.0.2 (last described in Foster et al., 2012). A particle’s spectrum is a sum of the individual spectra for each chemical element tracked by the simulations, scaled by the particle’s elemental abundance. We ignore particles with a temperature lower than as they make a negligible contribution to the total Xray emission.
We then estimate the density, temperature and metallicity of the hot gas in logarithmically spaced radial bins by fitting a single temperature apec model, with a fixed metallicity, to the summed spectra of all particles that fall within that radial bin. We then scale the spectra by the relative abundance of the heavy elements as the fiducial spectra assume solar abundance (Anders & Grevesse, 1989). The spectra have an energy resolution of in the range and are logarithmically spaced between . To get a closer match to the observations, we multiply the spectra by the effective area of Chandra. To derive temperature and density profiles of a cluster, we fit the spectrum in the range for each radial bin with a single temperature model using a leastsquares approach.
The temperature and density profiles derived from the Xray spectra are then used to perform a hydrostatic mass analysis of the cluster. The profiles are fit with the density and temperature models proposed by Vikhlinin et al. (2006) to produce a hydrostatic mass profile. We then derive various mass and radius estimates, such as and , from the hydrostatic mass profiles. With these estimates we calculate quantities, such as or , by summing the properties of the particles that fall within the set. Coreexcised quantities are calculated in the radial range of the aperture. Luminosities are calculated by integrating the spectra of all particles within the aperture in the requisite energy band, for example, bolometric luminosities are calculated in the range . Averaged Xray temperatures are calculated by fitting a single temperature model to the sum of the spectra of all particles within the aperture. We repeat this analysis for all clusters in the combined sample at all redshifts of interest. All quantities derived in this manner are labelled with the subscript ‘spec’.
2.5 Cluster sample selection
We select clusters from MACSIS and BAHAMAS to form a ‘combined’ sample with which we can investigate the cluster scaling relations. We perform our analysis at and . We create this sample at each redshift by selecting all clusters with a mass . Additionally, we introduce a mass cut at every redshift below which we remove any MACSIS clusters. For example, at () this cut is made at (). This removes a tail of clusters with low , but have high ratios (see Appendix A). For the luminositytemperature relation, we use the temperaturemass relation of the combined sample to convert the mass cut into a temperature cut. At this results in a sample of clusters, containing clusters from BAHAMAS and MACSIS clusters, and at a sample of clusters, from BAHAMAS and from MACSIS.
The MACSIS clusters enable the investigation of the behaviour of the most massive clusters at low redshift. These clusters are commonly selected in cosmological analyses because their deep potentials are expected to reduce the impact of nongravitational processes and as the brightest clusters they require shorter exposures. We select a hot, and therefore massive, cluster sample by selecting all clusters in the combined sample with a coreexcised Xray temperature greater than . At () this yields a sample of () clusters, with () coming from the MACSIS sample. Finally, we examine the impact of selecting a relaxed subset of the hot cluster sample. Theoretically, there are many ways to define a relaxed halo (see Neto et al., 2007; Duffy et al., 2008; Klypin et al., 2011; Dutton & Macciò, 2014; Klypin et al., 2016). For this study we use the following criteria
where is distance between the cluster’s minimum gravitational potential and centre of mass, divided by its virial radius; is the mass fraction within the virial radius that is bound to substructures; and is the spin parameter for all particles inside . These criteria are not designed to select a small subset that comprises the most relaxed objects, but to simply remove those clusters that are significantly disturbed. This results in a subsample at () that contains () clusters, with () coming from the MACSIS sample.
3 The scaling relations of massive clusters
In this section we present our main results, measuring the scaling relations of our cluster samples across a range of redshifts.
3.1 Comparison to observational data
Fig. 3 shows the gas mass, , the integrated SunyaevZel’dovich (SZ) signal, , measured in a aperture as a function of estimated total mass, , (at and ) and the coreexcised bolometric Xray luminosity, , as a function of coreexcised Xray temperature, , for the combined sample. We compare the sample to the relevant observational data. At all redshifts the MACSIS sample provides a consistent extension to the BAHAMAS clusters with similar scatter. At low redshift, McCarthy et al. (2016) have shown that the BAHAMAS sample shows good agreement with the observed median relations and shows similar intrinsic scatter. The MACSIS sample continues this agreement to observed highmass clusters, though there are significantly fewer clusters to compare against. In detail, it appears that the  and  relations are slightly steeper than observed. However, we would exercise caution as we have not applied the same selection criteria as was used for the observational Xray analyses.
At high redshift observational data becomes sparse and currently only SZ surveys have detected a reasonable number of clusters. At these clusters are all significantly more massive than any cluster in the BAHAMAS volume. However, the progenitors of the very massive MACSIS clusters provide a sample that can be compared with these observations. We find that the median relation shows good agreement with the observations and the intrinsic scatter of the clusters about the median relation is consistent with the scatter in the observations. Overall, we find that all quantities computed in a likewithlike manner show good agreement with the observations.
Scaling relation  Combined sample  Hot Clusters  Relaxed, Hot Clusters  

3.2 Modelling cluster scaling relations
As a baseline for understanding how the scaling relations evolve as a function of mass and redshift, we adopt the following selfsimilar scalings
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
where , is the chosen overdensity relative to the critical density and is the Xray analogue of the integrated SZ effect. These are derived in Appendix B. Although shown to be too simplistic by the first Xray studies of clusters (Mushotzky, 1984; Edge & Stewart, 1991; David et al., 1993), the selfsimilar relations allow us to investigate if astrophysical processes are less significant in more massive clusters or at higher redshift. To enable a comparison with the selfsimilar predictions, and previous work, we fit the scaling relations of our samples at each redshift. We derive a median relation by first binning the clusters into bins of log mass (width dex) or log temperature (width dex) and then computing the median in each bin with more than ten clusters. We also remove the evolution in normalization predicted by selfsimilar relations. The medians of the bins are then fit with a powerlaw of the form
(8) 
where and describe the normalization and slope of the best fit respectively, removes the expected selfsimilar evolution with redshift, is either the total mass or temperature and is the observable quantity (, , etc.). is the pivot point, which we set to for observablemass relations and to for observabletemperature relations. We note that we fix the pivot for all samples and all redshifts. Fitting to the medians of bins, rather than individual clusters, prevents the fit from being dominated by lowmass objects, which are significantly more abundant due to the shape of the mass function. For the hot sample and its relaxed subset there are too few bins with ten or more clusters to reliably derive a bestfit relation at . By limiting our sample to systems with we avoid any breaks in the powerlaw relations that have been seen both observationally and in previous simulation work (Le Brun et al., 2016).
We compute the scatter about the bestfit relation at each redshift by calculating the root mean squared (rms) dispersion in each bin according to
(9) 
where runs over all clusters in the bin, is the best fit relation for a cluster with a value and we note that . We obtain the uncertainties for our fit parameters by bootstrap resampling the clusters times. The bestfit values of all the scaling relations considered for the three samples (combined, hot and relaxed) at are summarized in Table 2 and other redshifts are listed in Appendix C. We now discuss each relation in turn.
3.3 Gas MassTotal Mass
We plot the hot gas masstotal mass scaling relation for the three samples in Fig. 4. The bestfit normalization for the combined sample shows significant evolution with redshift, with clusters of a fixed mass containing more hot gas at than at . With the inclusion of star formation, radiative cooling and feedback from supernovae and AGN, the departure from selfsimilarity is not unexpected. The increasing normalization with redshift is due to either the impact of AGN feedback or the conversion of gas to stars. As the normalization of the baryonic mass exhibits a similar trend, this evolution is being driven by AGN feedback. A plausible explanation is as follows. The mean density of the Universe increases with redshift and cluster potentials at a fixed mass get deeper with increasing redshift. This reduces the efficiency with which AGN expel gas from the cluster with increasing redshift, leading to a higher gas mass at higher redshift for clusters at a fixed mass. In addition, AGN have less time to act on and expel gas from clusters that form at higher redshifts. The AGN breaks the selfsimilar assumption of a constant gas fraction, resulting in the normalization of the gas masstotal mass relation increasing with increasing redshift. However, we note that this behaviour appears to be dependent on the implementation of the subgrid physics. Le Brun et al. (2016) use the same subgrid implementation, but with different parameters, and obtain similar behaviour. However, Planelles et al. (2013) see a constant baryon fraction with redshift suggesting that feedback is not expelling gas beyond .
The bottom left panel of Fig. 4 shows that the normalizations of the bestfit relations for the hot sample of clusters and for the relaxed subset of hot clusters are higher at than the normalization of the combined sample and evolve less with redshift. This is because hotter clusters are generally more massive and have deeper potential wells, reducing the amount of gas the AGN can permanently expel from the cluster during its formation. This flattens the slope of the relation leading to a higher normalization at the pivot.
The bottom right panel of Fig. 4 shows that the slope of the bestfit relation of the combined sample is significantly steeper than the selfsimilar prediction of unity. At a given redshift AGN feedback has expelled more gas from lower mass clusters, due to their shallower potentials, leading to a tilt in the relation. We find a slope of . Our slope is mildly shallower than found in previous a simulation work, where Le Brun et al. (2016) find a slope of for their AGN8.0 simulation, but consistent with observations, where Arnaud et al. (2007) found a slope of for a sample of clusters observed with XMM. We find negligible evolution in the slope of the relation for the combined sample.
The hot cluster sample and the relaxed subset have bestfit slopes that are consistent with the selfsimilar prediction. The increased depth of the potential well in massive clusters means that their gas mass is approximately a constant fraction of their total mass. Specifically, we find that most massive clusters have a median gas fraction of the universal baryon fraction at . This results in slopes of and for the hot cluster sample and the relaxed subset respectively. We find good agreement with the slope of found by Mantz et al. (2016) and the selfsimilar slope found by Vikhlinin et al. (2009) for relaxed cluster samples. The slope of the bestfit relation for both samples shows no significant evolution with redshift.
The top right panel of Fig. 4 shows that the scatter about the bestfit relation is independent of both mass and redshift. Averaged over all mass bins it has a value of at . The scatter reduces slightly for the hot cluster sample, with a value of , and further still for the relaxed subset, with a value of . The scatter is in reasonable agreement with the scatter of found by Arnaud et al. (2007) for a sample of clusters observed with XMM.
3.4 Xray TemperatureMass
The evolution of the coreexcised spectroscopic temperaturetotal mass scaling relations, and their scatter, for the three samples are shown in Fig. 5. The normalization of the bestfit relation of the combined sample shows a minor evolution with redshift, being lower at compared to (bottom left panel). In the selfsimilar model the temperature of the ICM is related to the depth of the gravitational potential of the cluster, under the assumption of hydrostatic equilibrium. Previous simulation work has shown that the nonthermal pressure in masslimited samples grows with redshift due to the increasing importance of mergers and resulting incomplete thermalisation (Stanek et al., 2010; Le Brun et al., 2016). Therefore, clusters increasingly violate the assumption of hydrostatic equilibrium with redshift and require a lower temperature at a fixed mass to balance gravitational collapse, which leads to a normlization that decreases with redshift compared to selfsimilar. The effective temperature of the nonthermal pressure can be estimated via
(10) 
where is the 1D velocity dispersion of the gas particles, is the mean molecular weight, is the mass of the proton and is the Boltzmann constant. Fig. 6 shows the evolution of the temperaturemass normalization once this effective kinetic temperature has been added to the spectral temperature. For all three samples the addition of the kinetic temperature results in a normalization that shows significantly reduced evolution with respect to selfsimilar.
The normalizations of the bestfit relations for the hot cluster and the relaxed hot samples are slightly higher than for the combined sample, but they show a similar trend with redshift that is removed when the kinetic temperature is included. The higher normalization occurs because, again, the hot sample has a flatter slope with mass. This flatter slope is driven by two processes. First, nonthermal pressure support becomes more important in higher mass clusters at a fixed redshift, as they have had less time to thermalise, and this lowers their temperatures. Second, we find that the bias between the spectroscopic and massweighted temperatures increases mildly with mass. This does not appear to be caused by cold clumps due to the SPH method, but is due to the presence of cooler gas in the outskirts of massive clusters that is hotter than the lower limit, contributing to the Xray spectrum, and biasing the measured temperature low for the most massive clusters. Fig. 7 shows the fractional difference between the spectroscopic and massweighted coreexcised temperatures as a function of mass. Similar to Biffi et al. (2014), we find that for lowmass clusters the spectroscopic temperature estimate agrees well with the massweighted estimate at . However, as cluster mass increases we find that the spectroscopic estimate is increasingly biased low compared to the massweighted estimate. This will also impact the hydrostatic mass estimate of the cluster and we refer the reader to Henson et al. (2016) for a more indepth study. Both of these effects lead to a flattening of the slope with mass and a higher normalization for the hot samples. We note that removing the most disturbed clusters produces a marginal decrease in the normalization of the relation, which is due to the steeper slope yielding a lower normalization at the pivot point.
We find the slope of the bestfit relation for the combined sample to be at . This is in good agreement with the slope found by previous simulation work, where values of (Short et al., 2010), (Stanek et al., 2010), (Planelles et al., 2014), (Biffi et al., 2014), (Pike et al., 2014) and (Le Brun et al., 2016) were found. All of these are in agreement with the observed temperaturetotal mass relation found for volumelimited samples, with values of for a sample of clusters observed with XMM (Arnaud et al., 2007) and for a sample of lowredshift clusters (Giles et al., 2015). We note that a caveat to these comparisons is the differing mass ranges will alter the slope as the relation is not a perfect power law. All of these relations are slightly flatter than the predicted selfsimilar slope of due to nonthermal pressure support and temperature bias.
Selecting only hot clusters produces a bestfit relation with a slope of , flatter than the combined relation. The bestfit slope of for the relaxed subset, is compatible with the combined sample. The slope of the relaxed subset is compatible with the slope found by Mantz et al. (2016) of and the slope of found by Vikhlinin et al. (2009) for relaxed clusters. However, we note that our relaxation criteria only remove the most disturbed objects, as opposed to the criteria of Mantz et al. (2015) which select the most relaxed objects. Therefore, we would likely recover a steeper slope with stricter relaxation criteria. Both samples are equally affected by the spectroscopic temperature being biased low. The slopes of the hot sample and the relaxed subset show no clear trend with redshift.
The temperaturemass scaling relations shows very low scatter, which is independent of both mass and redshift. The average scatter across all mass bins is , and for the combined sample, hot sample and relaxed subset, respectively, at . These values are consistent with the values found by both observations and previous simulations (Arnaud et al., 2007; Giles et al., 2015; Stanek et al., 2010; Short et al., 2010).
3.5  Mass
The powerlaw fits to the Xray analogue of the integrated SZ effecttotal mass relations for the three samples, and their scatter, are shown in Fig. 8. The Xray analogue signal, , is the product of the core excised spectral temperature and the gas mass, and the relation should reflect the combination of the two previously presented relations. We indeed find this to be the case. For the combined sample, the decreasing temperaturetotal mass normalization with increasing redshift offsets the increasing gas masstotal mass normalization, producing almost no evolution of the normalization for the total mass relation. The same trend was found by Le Brun et al. (2016). Therefore, the normalization evolves in a close to selfsimilar manner.
Selecting a sample of hot clusters or a relaxed subset of them leads to higher overall normalization of the bestfit relation. This is mainly due to the reduced impact of AGN feedback on the gas masstotal mass relation, which flattens the relation and leads to a higher normalization at the pivot. Both samples agree very well with the predicted selfsimilar evolution of the normalization of the relation, with the normalization of the relaxed subset changing by less than one percent between and .
The slope of the total mass relation is simply the sum of the slopes of the temperaturemass and gas masstotal mass relations and for the combined sample the slope is significantly steeper than the value predicted by selfsimilar theory. We find a value of at . The slope of our bestfit relation is consistent with those of previous simulations, who found values of (Short et al., 2010), (Planelles et al., 2014) and (Le Brun et al., 2016). Our result is also in agreement with the observational value found by Arnaud et al. (2007) of using the REXCESS cluster sample. The physical reason for the steeper slope is that gas is preferentially removed from lower mass clusters by feedback. In response to gas expulsion the remaining gas increases in temperature, offsetting some of the losses, but the loss of gas dominates and steepens the relation. The value of the slope for the bestfit relation is approximately constant with redshift, within the uncertainty of the fits.
Selecting a sample of hot clusters leads to a significant flattening of the slope of the relation, slightly flatter than the selfsimilar prediction of . With the gas masstotal mass relations of the hot sample and relaxed subset being very close to selfsimilar, the shallower than selfsimilar slope is due to the temperaturemass relation. The bestfit slope of both samples shows no significant trend with redshift.
The scatter about the bestfit relation is independent of both mass and redshift for all three samples, but it is noisy. We find an average value of at for the scatter for the combined sample, for the hot cluster sample and for the relaxed subset. These values are larger than those found previously for both simulations, where values of (Short et al., 2010), (Planelles et al., 2014) and (Le Brun et al., 2016) were found, and observations, where a value of was found for a sample of clusters observed with XMM (Arnaud et al., 2007).
3.6 Total Mass
The integrated SZ effecttotal mass relations for the three samples are shown in Fig. 9. Both the integrated SZ signal and its Xray analogue measure the total energy of the hot gas in the ICM, however the SZ signal depends on the massweighted temperature rather than the Xray spectral temperature. Our bestfit relation for the combined sample shows a mild evolution with redshift, with clusters at yielding an integrated signal that is higher than clusters at for a fixed mass. The evolution reflects the evolution in the gas masstotal mass relation. The increased evolution of its normalization compared to its Xray analogue suggests that the normalization of the massweighted temperature evolves more selfsimilarly then the spectroscopic Xray temperature and is indeed confirmed by the study of the massweighted temperaturetotal mass relation.
Selecting a sample of hot clusters or a relaxed subset of them significantly reduces the evolution in the normalization. The normalization of both samples, within the uncertainty of the fits, evolves in agreement with the selfsimilar prediction. Selecting a hot sample leads to a higher normalization than the combined sample at , due to the flatter slope of the gas masstotal mass relation yielding a flatter slope and a higher normalization at the pivot point.
The bestfit relation for the combined sample produces a slope of at , which is significantly steeper than the value predicted by the selfsimilar model. The value for the slope of the relation is consistent with previous values from both simulations, where values of (Stanek et al., 2010), (Battaglia et al., 2012), (Planelles et al., 2014), (Pike et al., 2014), (Yu et al., 2015) and (Le Brun et al., 2016) have been found, and observations, where was found for the Planck clusters (Planck Collaboration et al., 2014b) and was found for the clusters in the SPT survey. The steeper than selfsimilar slope is the result of the gas masstotal mass relation having a steeper slope. We find that the slope of the relation is independent of redshift.
The bestfit slopes of the hot cluster sample and the relaxed subset are consistent with the slope predicted by selfsimilar theory. The slopes of both samples are consistent with no evolution.
The scatter of the clusters about the bestfit relation shows no trend with either mass or redshift for all three samples. We find an average scatter of , and for the combined, hot and relaxed samples, respectively, at . This is larger than the scatter reported by previous simulations, where Battaglia et al. (2012), Pike et al. (2014), Planelles et al. (2014) and Le Brun et al. (2016) found values of , , and respectively, but in reasonable agreement with the values of and observed by Yu et al. (2015) and Planck Collaboration et al. (2014b) respectively.
3.7 Bolometric Xray LuminosityTotal Mass
Fig. 10 shows the coreexcised bolometric Xray luminositytotal mass scaling relations for the three samples and their evolution with redshift. The normalization of the bestfit relation for the combined sample shows significant evolution with redshift, being higher at compared to . The same physics driving the gas masstotal mass relation, increased binding energy, is driving the departure from selfsimilar. The Xray emission of a cluster is particularly sensitive to the thermal structure of the ICM, which depends on processes such as radiative cooling and feedback. Therefore, it is not surprising that the luminositymass relation shows significantly more evolution than other observablemass relations.
Selecting a sample of hot clusters significantly reduces the evolution in the normalization. Both the hot sample and the relaxed subset have a normalization that is larger at compared to the combined sample. The deeper potentials of more massive clusters reduces the impact of the AGN feedback and flattens the relation. This flattening leads to a higher luminosity at the pivot point. The normalizations of the bestfit relations for both the hot sample and its relaxed subset show very minor evolution, which is consistent with the selfsimilar prediction.
The slope of the bestfit relation for the combined sample is significantly steeper than the slope predicted by selfsimilar theory. At we find a slope of for the combined sample. This steepening is driven by AGN feedback being more effective in lower mass clusters. The slope at is in reasonable agreement with the slopes found in volumelimited observational samples, such as Pratt et al. (2009) who found a slope of for the REXCESS sample and Giles et al. (2015) who found a slope of for a sample of lowredshift clusters. Previous simulation work by Short et al. (2010), using the semianalytic feedback model of the Millennium Gas project, found a bolometric luminositytotal mass slope of and Stanek et al. (2010), using the preheating model of the Millennium Gas project, found a slope of . Biffi et al. (2014) found a slope of for the MUSIC simulations. The slope of the bestfit relation for the combined sample is approximately independent of redshift, with a very mild steepening of the slope with redshift occurring due to the reduction in fitting range with increasing redshift.
The slopes of the bestfit relation follow the same trend as the gas masstotal mass relation, with the hot sample and its relaxed subset producing shallower slopes that are in much better agreement with selfsimilar theory. Our bestfit slope is consistent with the observational result of Mantz et al. (2016), who found a selfsimilar slope for the coreexcised luminositytotal mass relation for a sample of relaxed clusters with .
The scatter about the best fit relation is approximately independent of both mass and redshift for all three samples, although it is relatively noisy. Averaging the scatter for the combined sample across all mass bins produces a value of . This is in reasonable agreement with the scatter found in lowredshift observational samples, where Pratt et al. (2009) find a value of and Giles et al. (2015) find a value of . Selecting hot clusters and a relaxed subset produces a small reduction in scatter about the bestfit relation with values of and respectively.
3.8 Xray LuminosityTemperature
Finally, we study the redshift evolution of the Xray luminosityspectroscopic temperature relation. Both quantities of the luminositytemperature scaling relation are observable, with the temperature tracing the depth of the potential of the cluster. This makes it a useful relation to study the impact of nongravitational physics. In Fig. 11 we plot the bolometric Xray luminosityspectroscopic temperature scaling relation for the three samples of clusters. The normalization of the bestfit relation for the combined sample shows significant evolution with redshift relative to selfsimilar. Clusters with a temperature of at have a luminosity greater than clusters with the same temperature at . This evolution can be thought of as being due to a combination of the evolution of the temperaturemass and luminositymass relations. The decreasing temperaturemass normalization and increasing luminositymass normalization with redshift combine to yield a significant evolution of the luminositytemperature normalization relative to selfsimilar.
Selecting a sample of hot clusters, or a relaxed subset of them, reduces the evolution, but there is still a mild evolution in the normalization. Hot clusters at a fixed temperature at are more luminous than those at . Combining equations (3) and (6), but allowing the slope of the relations to vary from their selfsimilar values yields
(11) 
where and are the slopes of the luminositymass and temperaturemass relations respectively. Hence, deviations of their slopes from selfsimilar leads to evolution of the normalization of the luminositytemperature relation that is not selfsimilar. With the luminositymass relation being selfsimilar for the hot cluster sample and its relaxed subset, the evolution of the normalization is being driven by the flatter than selfsimilar slope of the temperaturemass relation, which is due to the increased importance of nonthermal pressure support and the increasingly biased spectroscopic temperatures of more massive clusters.
We find a slope of for the bestfit relation at . This is significantly steeper than the slope of predicted by selfsimilar theory. However, this value is reasonable agreement with previous simulation work, (Short et al., 2010), and those found by observations, for the REXCESS sample (Pratt et al., 2009) and for a sample of clusters observed with Chandra (Maughan et al., 2012). It is clear from equation (11) that the slope of the relation depends on the slopes of the luminositymass and temperaturemass relations. The steeper than expected slope for the combined sample is due to the combined effects of AGN feedback on the luminosity slope and nonthermal pressure support and temperature bias on the temperature slope, both of which lead to a steepening of the relation compared to the selfsimilar prediction. We find that the bestfit relation steepens slightly with redshift, increasing to at . This evolution is due to the removal of highmass objects with redshift.
The bestfit slope of the hot cluster sample and the relaxed subset are flatter than the combined relation with slopes of and . This is still significantly steeper than the slope predicted by selfsimilar theory, but in good agreement with the slope of observed by Maughan et al. (2012) for their relaxed cool core cluster sample. With both samples exhibiting selfsimilar slopes for the luminositymass relations, the deviation from selfsimilarity is being driven by their temperaturemass relations.
The scatter about the bestfit relation demonstrates a trend with both temperature and redshift. Although somewhat noisy, the scatter appears to increase with decreasing temperature. The average scatter at for the combined sample is . This scatter is consistent with the simulations of Short et al. (2010), who found a scatter , and the intrinsic observational scatter of found by Pratt et al. (2009). However, it is significantly lower than the scatter of found by Maughan et al. (2012). The scatter reduces for the hot cluster sample and the relaxed subset to and respectively.
3.9 Summary
Overall, the scaling relations of the combined sample show good agreement with previous work, both simulations and observations. Departures from selfsimilarity are driven by the increased efficiency of gas expulsion by AGN feedback in clusters with shallower potentials, due to being less massive or forming at a lower redshift; the increased contribution of nonthermal pressure that supports the ICM against gravity in more massive clusters or those at higher redshifts; and the increase in the spectroscopic temperature bias for the most massive clusters. The MACSIS sample enabled the scaling relations to be studied to higher redshifts, as their progenitors are still clusters at high redshift, and the examination of the impact of selecting a sample of hot clusters on the evolution of the scaling relations. This demonstrated that massive clusters are more selfsimilar and evolve more selfsimilarly with redshift compared to the overall cluster population, as the efficiency of gas expulsion by AGN feedback is reduced due to their deeper potentials. However, it also highlighted that nonthermal pressure support becomes more important in these clusters and that their spectroscopic temperatures are biased low.
4 Evolution of gas profiles
Most of the scaling relations of hot, and therefore massive, clusters evolve in a way that is consistent with the predictions of the selfsimilar model. However, the combined sample showed significant deviations from the selfsimilar model due to the impact of nongravitational processes. To further understand the differences between the samples in the evolution of their scaling relations, we now examine the gas profiles of the different cluster samples. To enable a quantitative comparison with the observational data requires us to compare likewithlike. Therefore, we restrict the mass range of the combined sample to , yielding a sample with a median mass of . We compare this to the REXCESS cluster sample which has a median mass of and a sample of clusters from Giles et al. (2015) with a median mass of . Although this mass matching does not account for selection effects, it should allow for a quantitative comparison. We do not alter the hot sample or the relaxed subset. We factor out the expected selfsimilar evolution in the profiles by dividing by the appropriate quantity, e.g. , , or . We define these quantities as
(12) 
(13) 
(14) 
(15) 
where is the Hubble constant, is the gravitational constant, is the mean atomic weight per free electron and is the universal baryon fraction. Therefore, any changes in the profiles are due to nongravitational physics, such as AGN feedback or nonthermal pressure support.
4.1 Density profiles
The threedimensional dimensionless density profiles for the three cluster samples at and are shown in Fig. 12. We have scaled the profiles by to reduce the dynamic range. At , we compare the median profile of the combined sample with the observed median profiles from Croston et al. (2008) for the REXCESS sample and Giles et al. (2015) for a sample of lowredshift clusters observed with Chandra. The combined sample shows good agreement with the observed profiles and has similar intrinsic scatter. Beyond a radius of the median profiles of the hot sample and its relaxed subset have a similar shape as the combined sample, but the densities are higher as they are on average more massive clusters. Inside this radius the profiles of both samples have a shallower gradient compared to the combined sample. This is caused by the accretion of lowentropy, highdensity gas that sinks to the centre of the cluster potential, becoming increasingly important below (Power et al., 2014). This effect is not offset in massive clusters by the AGN feedback and so their density profiles have a shallower gradient in the core. We note that this effect can potentially impact the relations we presented in section 3. However, we presented coreexcised temperatures and luminosities, which should minimise any bias introduced by the accretion of poorly mixed gas.
At , we compare the median density profiles of the three samples to the observed profile from McDonald et al. (2013), which has been derived from a sample of clusters with a mean redshift of . These clusters were selected from the SPT survey catalogue and observed with Chandra. There is a reasonable agreement between the combined sample’s median profile and the observations, but the observations are higher between . The observed profile is in better agreement with the median profiles of the hot sample and its relaxed subset. This suggests that the observed clusters are more representative of more massive objects at . There is better agreement between the density profiles of the three samples at because the mass cut of causes the samples to converge with increasing redshift. Selecting relaxed hot clusters leads to a median profile that is slightly more centrally concentrated than for all hot clusters.
In the bottom panel of the plot we show the of the ratio of the median density profile at and the median profile at for each sample. For the hot cluster sample and the relaxed subset the profiles have evolved in a selfsimilar way beyond , showing very little change. Inside of this radius the impact of accreting lowentropy, highdensity gas that sinks to the centre of the cluster is apparent as an increase in the density profiles from to . For the combined sample the difference between the two profiles shows the increase of the depth of the potential with redshift. This leads to higher densities at and a negative change density profile at all radii with decreasing redshift.
4.2 Temperature profiles
Fig. 13 shows the threedimensional temperature profiles divided by the predicted selfsimilar temperature. At the profiles all have a similar shape, but the normalization of the combined sample is somewhat higher than those of the hot sample and its relaxed subset. This is due to the lower gas density of the combined sample, which requires a higher temperature to balance gravitational collapse. Also, there is likely to be a small effect due to the mass dependence of nonthermal pressure support, with more massive clusters having more nonthermal support and lower temperatures. The accretion of lowentropy, cold gas that sinks to the cluster core produces a steeper temperature gradient in the central profiles of the hot sample and its relaxed subset. Overlaid are the observed median temperature profiles from two cluster samples, the REXCESS sample (Arnaud et al., 2010) and a sample of clusters observed with Chandra (Giles et al., 2015). The median profile of the combined sample and its intrinsic scatter show good agreement with the observed temperature profiles and their scatter.
At all samples have a similar profile shape, but the hot sample has a lower normalization compared to the combined and relaxed hot sample. This is because nonthermal pressure support becomes increasingly important in clusters of a fixed mass with redshift, leading to a lower temperature in hot clusters. The relaxed sample removes the most disturbed objects with greatest level of nonthermal support, producing a higher median temperature profile. We compare to the observed median profile of McDonald et al. (2014). The median profiles of the combined sample and the relaxed hot sample slightly under predict the observations at and over predict the observations at large radii, but the observed profile is within the scatter of the combined sample.
Within the median temperature profiles show significantly less evolution between the two redshifts than the density profiles. The combined and hot samples deviate from selfsimilarity and show an increase in temperature from to at all radii, consistent with the decreasing temperaturemass normalization with increasing redshift found in Section 3.4. This is because nonthermal pressure support decreases with increasing redshift. Therefore, as clusters thermalise their temperatures must increase to balance gravitational collapse, resulting in a hotter temperature profile at compared to . Selecting a relaxed subset reduces the nonthermal pressure support and the median profile changes significantly less from to inside .
4.3 Pressure profiles
The dimensionless pressure profiles, scaled by , of the three cluster samples are shown in Fig. 14. The increased mass of the hot sample and its relaxed subset lead to median pressure profiles that are higher in the centre at due to their higher densities. We compare the median profiles to the observed median pressure profiles from Arnaud et al. (2010) and Giles et al. (2015) and the bestfit profile from Planck Collaboration et al. (2013). We note that the Planck result is based on the stacked profile of nearby systems. For Giles et al. (2015) we have combined their published density and temperature profiles to produce a pressure profile for each cluster. There is good agreement between the combined sample and the observed profiles, with a slight over prediction at large radii. For comparison to the Planck bestfit parameters we fit the mean profiles of our clusters at both redshifts with a generalised NavarroFrenkWhite pressure profile (Navarro et al., 1997; Nagai et al., 2007b) of the form
(16) 
We fit a four parameter model with fixed. The results are shown in table 3.
Sample  

Planck  
Combined  
Hot  
Relaxed Hot  
Combined  
Hot  
Relaxed Hot 
At the median profiles of the three samples are in closer agreement with each other, because the minimum mass limit of causes the samples to converge at high redshift. We compare our median pressure profiles with the observed profile of McDonald et al. (2014). They find a median pressure profile that is in good agreement with the median profiles, but it is most consistent with the relaxed hot sample of massive clusters.
The pressure profile of the relaxed subset shows very little evolution between and , except for the core where the increasing density leads to an increased pressure with decreasing redshift. The hot sample shows an increased pressure in the core with decreasing redshift, due to the increased density, but a negative change in pressure from to at larger radii. The combined sample shows a negative pressure change between and at all radii. The decreased pressure with decreasing redshift is caused by the decrease in density from to .
4.4 Entropy Profiles
The median entropy profiles are shown in the bottom right panel of Fig. 15 and they have been normalized by the predicted selfsimilar entropy. We note that we define entropy as
(17) 
where is the electron number density and is the chosen overdensity relative to the critical density of the Universe. At the the combined sample shows a higher normalization compared to the hot sample and its relaxed subset. This is due to it lower density profile and higher temperature profile. The gradients of the hot sample and the relaxed subset profiles steepen in the centre due to the accretion of low entropy gas. We compare with the observed median profiles of Pratt et al. (2010) and Giles et al. (2015), and the baseline profile of Voit et al. (2005) derived from nonradiative SPH simulations. The combined sample is in good agreement with the observations and tends to the nonradiative predictions at large radii.
At the three samples are in reasonable agreement with each other, all having a similar shape with the hot sample showing a marginally lower normalization. This change from is in agreement with the evolution in their density and temperature profiles. We compare the profiles to the observations of McDonald et al. (2014). The combined and relaxed hot sample show good agreement with the observed profile for , but over predict the entropy at larger radii. In contrast the median profile of the hot sample is consistent with the observations at large radii, but under predicts the entropy in the centre of the cluster.
The departure from selfsimilarity for the three samples is due to a combination of the evolution in their temperature and density profiles. The relaxed hot sample shows a mild increase in entropy from to at large radii, due to change in its temperature profile, and a decrease in entropy in the core due to the increase in density at . The increased normalization of the hot sample’s temperature profile at compared to leads to an increased entropy profile with decreasing redshift, except in the core. The combined sample shows an increase in entropy at all radii at compared to and is produced by the decreased density and increased temperature with decreasing redshift.
5 Summary & Discussion
In this work we have presented the MACSIS clusters, a sample of 390 zoomed simulations of the most massive and rarest clusters run with the stateoftheart, calibrated baryonic physics model from the BAHAMAS project (McCarthy et al., 2016) that yields realistic clusters. Such massive clusters are absent from the BAHAMAS simulation volumes of as the simulated volume is too small. After introducing the selection of the sample from the parent volume simulated with the Planck 2013 cosmology, and demonstrating the agreement of the properties of our massive cluster sample with the properties of observed massive clusters, we examined the evolution of the cluster scaling relations and the evolution of the cluster gas profiles.
By combining the MACSIS sample with the clusters in the BAHAMAS volume, we were able to examine the cluster scaling relations over the full observed mass range for the first time. Additionally, the MACSIS clusters enabled the study of the evolution of the cluster scaling relations to unprecedentedly high redshifts. Finally, the MACSIS sample enabled clusters to be selected in ways which mimic a cosmological study, such as selecting the hottest clusters, to examine if the scaling relations of such objects evolve differently from the underlying cluster population. Our main results are:

As shown in Fig. 3, the MACSIS simulations yield realistic massive clusters at low redshift and their progenitors are in good agreement with the limited observational data that is available at high redshift (i.e. ).

Scaling relations for the combined sample that spans the full observed cluster mass range show significant deviations from the simple selfsimilar theory (see Figs. 411). Both the slope of the relations and the redshift evolution of the normalization are significantly affected by nongravitational physics. The low redshift relations are in good agreement with observations and with most previous simulation work.

The main drivers of nonselfsimilar evolution are AGN feedback, nonthermal pressure support and a mild mass dependence of the spectroscopic temperature bias. Shallower potentials of clusters that are less massive or form at lower redshifts allows feedback from AGN to eject more gas. Nonthermal pressure lowers a cluster’s temperature for a given potential and is more important in more massive clusters that have had less time to thermalise. We found that the spectroscopic temperature bias increases for the most massive clusters.

With the exception of the luminositytemperature relation, we found the scatter about the bestfit scaling relations is insensitive to mass and redshift for all of the cluster samples.

Selecting a hot cluster sample, i.e. coreexcised spectroscopic temperatures , significantly alters the scaling relations and their evolution. Excluding the spectroscopic temperaturetotal mass relation, we find that the scaling relations of the hot cluster sample evolve in a much more selfsimilar manner. After accounting for the expected selfsimilar evolution with redshift, we find that the normalizations are consistent with no evolution. The slopes of the bestfit relations at each redshift are also broadly consistent with the slopes predicted by selfsimilar theory. However, the spectroscopic temperaturetotal mass relation of the hot sample deviates further from selfsimilarity than the combined sample. Selecting hot clusters removes the less massive clusters from the sample, so the hot sample is dynamically younger than the combined sample as more massive clusters form later in the hierarchical merger scenario. This increases the average level of nonthermal support in the hot sample, leading to a flatter spectroscopic temperaturetotal mass relation. Additionally, the spectroscopic temperature bias flattens the relation for the most massive clusters and this has a larger impact in a sample of only hot clusters.

Selecting a relaxed subset of hot clusters, where the most dynamically disturbed objects are removed, leads to a small reduction in the scatter for most scaling relations. Removing the most disturbed objects also leads to a reduction in the level of nonthermal support in the sample compared to the complete hot sample. This leads to steeper slope of the spectroscopic temperaturetotal mass relation compared to the hot sample and a value that is closer to the selfsimilar prediction.

The median hot gas profiles of the combined sample in general shows good agreement with observed radial profiles. The low redshift data is in very good agreement, while the data at shows reasonable agreement with the relaxed hot sample.

Comparison of the hot gas profiles at and show evolution different from selfsimilar prediction (see Figs. 1215). The combined sample shows a decreasing density profile with decreasing redshift, suggesting the impact of AGN feedback. Selecting a sample of hot clusters produces a median density profile that evolves in much more selfsimilar manner. The combined and hot samples have a median temperature profile that increases with decreasing redshift. This is likely driven by decreasing importance of nonthermal pressure support with decreasing redshift. Selecting relaxed hot cluster sample produces a median profile that evolves in better agreement with the selfsimilar prediction.
MACSIS enables the study of the observable properties of the most massive and rarest galaxy clusters. We have demonstrated that their progenitors provide a good match to the currently limited observational data at high redshift and that their observable properties evolve in a significantly more selfsimilar manner than for lowermass and lessrelaxed clusters. We have shown how the selection function can impact the derived scaling relations and radial profiles.
The size of the parent simulation enables the creation of synthetic lightcones with an area comparable to currently ongoing surveys. This will allow the impact of selection biases to be fully examined and the covariance of observable properties to be studied. Another route for future work is to improve our understanding of structure in the ICM, as the limited resolution and traditional SPH scheme used in this work limits our ability resolve structures and understand their impact on observable properties.
Acknowledgements
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 Einfrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National EInfrastructure. DJB and STK acknowledge support from STFC through grant ST/L000768/1. MAH is supported by an STFC quota studentship. IGM is supported by a STFC Advanced Fellowship. The research was supported in part by the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013) / ERC Grant agreement 278594GasAroundGalaxies. ARJ acknowledges support from STFC through grant ST/L00075X/1.
References
 Allen et al. (2011) Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409
 Anders & Grevesse (1989) Anders E., Grevesse N., 1989, Geochimica Cosmochimica Acta, 53, 197
 Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Arnaud et al. (2007) Arnaud M., Pointecouteau E., Pratt G. W., 2007, A&A, 474, L37
 Arnaud et al. (2010) Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, A92
 Battaglia et al. (2012) Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 74
 Benson et al. (2014) Benson B. A., et al., 2014, in Millimeter, Submillimeter, and FarInfrared Detectors and Instrumentation for Astronomy VII. p. 91531P (arXiv:1407.2973), doi:10.1117/12.2057305
 Bhattacharya et al. (2008) Bhattacharya S., Di Matteo T., Kosowsky A., 2008, MNRAS, 389, 34
 Biffi et al. (2014) Biffi V., Sembolini F., De Petris M., Valdarnini R., Yepes G., Gottlöber S., 2014, MNRAS, 439, 588
 Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
 Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
 Courtin et al. (2011) Courtin J., Rasera Y., Alimi J.M., Corasaniti P.S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911
 Crain et al. (2007) Crain R. A., Eke V. R., Frenk C. S., Jenkins A., McCarthy I. G., Navarro J. F., Pearce F. R., 2007, MNRAS, 377, 41
 Crocce et al. (2006) Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369
 Croston et al. (2008) Croston J. H., et al., 2008, A&A, 487, 431
 Dalla Vecchia & Schaye (2008) Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431
 David et al. (1993) David L. P., Slyz A., Jones C., Forman W., Vrtilek S. D., Arnaud K. A., 1993, ApJ, 412, 479
 Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
 Duffy et al. (2008) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64
 Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
 Edge & Stewart (1991) Edge A. C., Stewart G. C., 1991, MNRAS, 252, 414
 Eke et al. (1998) Eke V. R., Navarro J. F., Frenk C. S., 1998, ApJ, 503, 569
 Fabjan et al. (2010) Fabjan D., Borgani S., Tornatore L., Saro A., Murante G., Dolag K., 2010, MNRAS, 401, 1670
 Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
 Foster et al. (2012) Foster A. R., Ji L., Smith R. K., Brickhouse N. S., 2012, ApJ, 756, 128
 Giles et al. (2015) Giles P. A., et al., 2015, preprint, (arXiv:1510.04270)
 Haardt & Madau (2001) Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in Xrays. (arXiv:astroph/0106018)
 Heitmann et al. (2015) Heitmann K., et al., 2015, ApJS, 219, 34
 Henderson et al. (2016) Henderson S. W., et al., 2016, Journal of Low Temperature Physics,
 Henson et al. (2016) Henson M. A., Barnes D. J., Kay S. T., McCarthy I. G., Schaye J., 2016, preprint, (arXiv:1607.08550)
 Jenkins (2010) Jenkins A., 2010, MNRAS, 403, 1859
 Jenkins (2013) Jenkins A., 2013, MNRAS, 434, 2094
 Jenkins & Booth (2013) Jenkins A., Booth S., 2013, preprint, (arXiv:1306.5771)
 Jenkins et al. (2001) 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 (1986) Kaiser N., 1986, MNRAS, 222, 323
 Katz & White (1993) Katz N., White S. D. M., 1993, ApJ, 412, 455
 Kay et al. (2004) Kay S. T., Thomas P. A., Jenkins A., Pearce F. R., 2004, MNRAS, 355, 1091
 Kennicutt (1998) Kennicutt Jr. R. C., 1998, ARA&A, 36, 189
 Khedekar et al. (2013) Khedekar S., Churazov E., Kravtsov A., Zhuravleva I., Lau E. T., Nagai D., Sunyaev R., 2013, MNRAS, 431, 954
 Klypin et al. (2011) Klypin A. A., TrujilloGomez S., Primack J., 2011, ApJ, 740, 102
 Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
 Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353
 Laureijs et al. (2011) Laureijs R., et al., 2011, preprint, (arXiv:1110.3193)
 Le Brun et al. (2014) Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270
 Le Brun et al. (2016) Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2016, preprint, (arXiv:1606.04545)
 Lin et al. (2012) Lin Y.T., Stanford S. A., Eisenhardt P. R. M., Vikhlinin A., Maughan B. J., Kravtsov A., 2012, ApJ, 745, L3
 Mantz et al. (2010) Mantz A., Allen S. W., Rapetti D., Ebeling H., 2010, MNRAS, 406, 1759
 Mantz et al. (2015) Mantz A. B., Allen S. W., Morris R. G., Schmidt R. W., von der Linden A., Urban O., 2015, MNRAS, 449, 199
 Mantz et al. (2016) Mantz A. B., Allen S. W., Morris R. G., Schmidt R. W., 2016, MNRAS, 456, 4020
 Maughan et al. (2008) Maughan B. J., Jones C., Forman W., Van Speybroeck L., 2008, ApJS, 174, 117
 Maughan et al. (2012) Maughan B. J., Giles P. A., Randall S. W., Jones C., Forman W. R., 2012, MNRAS, 421, 1583
 McCarthy et al. (2010) McCarthy I. G., et al., 2010, MNRAS, 406, 822
 McCarthy et al. (2016) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2016, preprint, (arXiv:1603.02702)
 McDonald et al. (2013) McDonald M., et al., 2013, ApJ, 774, 23
 McDonald et al. (2014) McDonald M., et al., 2014, ApJ, 794, 67
 Merloni et al. (2012) Merloni A., et al., 2012, preprint, (arXiv:1209.3114)
 Mushotzky (1984) Mushotzky R. F., 1984, Physica Scripta Volume T, 7, 157
 Nagai et al. (2007a) Nagai D., Vikhlinin A., Kravtsov A. V., 2007a, ApJ, 655, 98
 Nagai et al. (2007b) Nagai D., Kravtsov A. V., Vikhlinin A., 2007b, ApJ, 668, 1
 Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Neto et al. (2007) Neto A. F., et al., 2007, MNRAS, 381, 1450
 Pike et al. (2014) Pike S. R., Kay S. T., Newton R. D. A., Thomas P. A., Jenkins A., 2014, MNRAS, 445, 1774
 Planck Collaboration et al. (2013) Planck Collaboration et al., 2013, A&A, 550, A131
 Planck Collaboration et al. (2014a) Planck Collaboration et al., 2014a, A&A, 571, A1
 Planck Collaboration et al. (2014b) Planck Collaboration et al., 2014b, A&A, 571, A20
 Planck Collaboration et al. (2015) Planck Collaboration et al., 2015, preprint, (arXiv:1502.01597)
 Planelles et al. (2013) Planelles S., Borgani S., Dolag K., Ettori S., Fabjan D., Murante G., Tornatore L., 2013, MNRAS, 431, 1487
 Planelles et al. (2014) Planelles S., Borgani S., Fabjan D., Killedar M., Murante G., Granato G. L., RagoneFigueroa C., Dolag K., 2014, MNRAS, 438, 195
 Power et al. (2014) Power C., Read J. I., Hobbs A., 2014, MNRAS, 440, 3243
 Pratt et al. (2009) Pratt G. W., Croston J. H., Arnaud M., Böhringer H., 2009, A&A, 498, 361
 Pratt et al. (2010) Pratt G. W., et al., 2010, A&A, 511, A85
 Puchwein et al. (2008) Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
 Reed et al. (2013) Reed D. S., Smith R. E., Potter D., Schneider A., Stadel J., Moore B., 2013, MNRAS, 431, 1866
 Sarazin (1986) Sarazin C. L., 1986, Reviews of Modern Physics, 58, 1
 Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
 Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
 Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
 Short et al. (2010) Short C. J., Thomas P. A., Young O. E., Pearce F. R., Jenkins A., Muanwong O., 2010, MNRAS, 408, 2213
 Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
 Stanek et al. (2010) Stanek R., Rasia E., Evrard A. E., Pearce F., Gazzola L., 2010, ApJ, 715, 1508
 The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration 2005, preprint, (arXiv:astroph/0510346)
 Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
 Tormen et al. (1997) Tormen G., Bouchet F. R., White S. D. M., 1997, MNRAS, 286, 865
 Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
 Vikhlinin et al. (2009) Vikhlinin A., et al., 2009, ApJ, 692, 1033
 Voit (2005) Voit G. M., 2005, Reviews of Modern Physics, 77, 207
 Voit et al. (2005) Voit G. M., Kay S. T., Bryan G. L., 2005, MNRAS, 364, 909
 Watson et al. (2013) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
 Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
 White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
 Wiersma et al. (2009a) Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 393, 99
 Wiersma et al. (2009b) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b, MNRAS, 399, 574
 Yu et al. (2015) Yu L., Nelson K., Nagai D., 2015, ApJ, 807, 12
Appendix A Selection Effects
The selection of the MACSIS sample was done using the mass of a halo in the parent simulation, but the scaling relations are presented using the cluster mass. This could potentially lead to a selection bias that impacts on the scaling relations presented in this work. Therefore, we make a mass cut to remove those clusters with large ratio. To assess the impact of our sample selection we plot the log of the ratio of four cluster observables: the coreexcised bolometric Xray luminosity, the coreexcised spectroscopic temperature, the gas mass and the integrated SZ signal, against their expected value from the bestfit relation as a function of the log of the ratio of the cluster’s and for the combined sample of clusters, with the cut included. Fig. 16 shows the result of these plots. Any differences in the correlations of these ratios between MACSIS and BAHAMAS could indicate that selection effects were impacting the scaling relations. We have calculated the Spearman’s rank correlation coefficients for both samples for all four observable ratios and find that only the gas mass shows a significant difference between the two samples. However, all of the quantities show only weak correlations. We therefore conclude that the cut to remove clusters with extremely high ratios has minimised any bias due to selection by .
Appendix B Selfsimilar relations
If galaxy clusters were to form from a purely monolithic gravitational collapse, astrophysical processes were negligible and they were virialised, then we would expect them to be selfsimilar objects. This would mean that their properties would depend only on their mass (White & Rees, 1978; Kaiser, 1986). The critical density of the Universe is defined as
(18) 
where is the Hubble constant, is the gravitational constant and
(19) 
A cluster can then be defined as an overdensity with mass, , inside a sphere of radius, , with some average density, , relative to the critical density
(20) 
As gas collapses into the potential, , of the cluster it is heated and, assuming that it is a collapsed iosthermal sphere, it will reach a temperature, , of
(21) 
where is the Boltzmann constant, is the mass of the proton and is the mean molecular weight. Therefore, the selfsimilar temperature of the cluster is proportional to its mass via
(22) 
Under the assumption that main cooling mechanism of the cluster is thermal bremsstrahlung, the cluster gas will emit Xrays and its bolometric emission is is proportional to
(23) 
where the cooling function for the bolometric case (e.g. Sarazin, 1986). Using equation (22), we can derive the selfsimilar prediction for the Xray luminositytemperature relation
(24) 
Assuming a constant gas fraction, the integrated SunyaevZel’dovich signal, , and its Xray analogue, , of the cluster can be predicted by
(25) 
and the selfsimilar relations are
(26) 
(27) 
Appendix C Fit parameters
The tables 49 below list the parameter values for the bestfit relations of the scaling relations presented in this paper. For there are too few clusters in too many bins to reliably measure a bitfit relation for the hot cluster sample and the relaxed subset and these values are not presented.
Redshift  Combined sample  Hot Clusters  Relaxed, Hot Clusters  

Redshift  Combined sample  Hot Clusters  Relaxed, Hot Clusters  

Redshift  Combined sample  Hot Clusters  Relaxed, Hot Clusters  
