Rhapsody-G simulations I: the cool cores, hot gas and stellar content of massive galaxy clusters

We present the Rhapsody-G suite of cosmological hydrodynamic AMR zoom simulations of ten massive galaxy clusters at the $M_{\rm vir}\sim10^{15}\,{\rm M}_\odot$ scale. These simulations include cooling and sub-resolution models for star formation and stellar and supermassive black hole feedback. The sample is selected to capture the whole gamut of assembly histories that produce clusters of similar final mass. We present an overview of the successes and shortcomings of such simulations in reproducing both the stellar properties of galaxies as well as properties of the hot plasma in clusters. In our simulations, a long-lived cool-core/non-cool core dichotomy arises naturally, and the emergence of non-cool cores is related to low angular momentum major mergers. Nevertheless, the cool-core clusters exhibit a low central entropy compared to observations, which cannot be alleviated by thermal AGN feedback. For cluster scaling relations we find that the simulations match well the $M_{500}-Y_{500}$ scaling of Planck SZ clusters but deviate somewhat from the observed X-ray luminosity and temperature scaling relations in the sense of being slightly too bright and too cool at fixed mass, respectively. Stars are produced at an efficiency consistent with abundance matching constraints and central galaxies have star formation rates consistent with recent observations. While our simulations thus match various key properties remarkably well, we conclude that the shortcomings strongly suggest an important role for non-thermal processes (through feedback or otherwise) or thermal conduction in shaping the intra-cluster medium.


INTRODUCTION
While simulations of galaxy formation in Milky Way-sized haloes (e.g. Guedes et al. 2011;Christensen et al. 2014) or small cosmic volumes (e.g. Vogelsberger et al. 2014;Dubois et al. 2014;Schaye et al. 2015) are making substantial progress, realizing the population of galaxies that reside in the most massive cosmic haloes, those hosting rich clusters of galaxies, remains a formidable challenge (Kravtsov & Borgani 2012). Compared to a 10 12 M Milky Way halo, simulating a 10 15 M cluster requires at least an extra order of magnitude in spatial Email: oliver.hahn@oca.eu resolution and three orders of magnitude in mass, thus requiring larger and longer computations. In addition, the Gaussian random field nature of the initial conditions biases the progenitors of clusters toward extreme systems at all redshifts, meaning cluster evolution is tied to that of the first stars, the earliest proto-galaxies, and the most massive quasars and their supermassive black hole interiors at high redshifts. Finally, the deep nature of the gravitational potential well retains most of the cosmic baryonic content associated with the dark matter, implying that a complex mix of coupled hydrodynamic, magnetohydrodynamic, chemical, and radiative processes must be solved in order to gain insights into the cycling and transport of mass, energy/entropy, momentum, and metals across ∼ 10 Mpc regions over most of the 13.8 Gyr history of the universe.
Early hydrodynamic simulations with gas cooling successfully formed multiple galaxies within a cluster (Katz & White 1993;Evrard et al. 1994), but the central galaxies were too massive compared to observations. Radio and X-ray observations of cavities near the central galaxy highlighted the need for strong feedback to curtail central cooling (c.f., McNamara & Nulsen 2007. Inclusion of supermassive black hole (SMBH) feedback into semi-analytic methods applied to halo ensembles from N-body simulations led to a reduction in central galaxy stellar masses (Croton et al. 2006;De Lucia & Blaizot 2007). Application within hydrodynamic simulations required estimating accretion rates onto the SMBH of active galactic nuclei (AGN), a task that in principle requires very high spatial resolution. Schemes such as that developed by Booth & Schaye (2009) were promoted to estimate Bondi-Hoyle accretion rates from simulations with roughly kiloparsec spatial resolution. Gas dynamic simulations using this approach to AGN feedback show improvements to central cluster galaxy morphology (Martizzi et al. 2012a) as well as improved scaling of hot intracluster medium (ICM) properties with halo mass (e.g. Planelles et al. 2014;Le Brun et al. 2014, and references therein).
The detailed nature of AGN feedback and its relationship to the phenomenology of the core plasma in clusters remain subjects of active investigation. X-ray observations show that the massive cluster population can be divided into cool-core (CC) and non-cool core (NCC) categories, with the distinction determined by the strength of a surface brightness cusp within the inner ∼ 100 kpc (Allen et al. 2011). The physics that controls this CC/NCC dichotomy has remained under debate. Major mergers may be capable of driving a CC to NCC transition through shock heating and/or ablation of core gas, but early cosmological gas dynamic simulations with cooling and supernova feedback found major mergers to be ineffective at heating the core gas  unless the merger had sufficiently low angular momentum (Poole et al. 2008). The observational finding that morphologically disturbed clusters rarely host cool cores offers hoewever strong empirical evidence that mergers play a role in determining core morphology ).
Still, it remains possible that internal processes, such as AGN jet-driven turbulent heating, may drive the transition from a CC to a NCC state, at least temporarily (Guo & Oh 2009;Parrish et al. 2009;Guo & Mathews 2010;Parrish et al. 2010;Ruszkowski & Oh 2010) or may play an important role in pressurising the cores. Idealized, high-resolution AMR simulations find cold, chaotic accretion onto the SMBH at rates many times the Bondi expectation (Gaspari et al. 2013b), and this has motivated models in which precipitation of turbulent core gas serves to self-regulate AGN accretion and feedback (Gaspari et al. 2013a;Li & Bryan 2014;Li et al. 2015). Voit et al. (2015) provide empirical support that such a model, coupled with moderate conductive heating from the external reservoir of ICM plasma, may help explain CC/NCC dichotomy. Ultraviolet imaging of brightest cluster galaxies in the CLASH cluster sample ) reveals knots and extended filamentary structures suggestive of the bi-polar streams that emerge in the simulations of Li et al. (2015).
In this paper, we present results from a suite of gas dynamic simulations of ten high mass haloes. These simulations extend our previous suite of N -body simulations of massive clusters falling inside a very narrow mass range at z = 0 (Wu et al. 2013a,b) to multi-physics adaptive mesh simulations. They include cooling, star formation and supermassive black holes as well as their respective feedback. These subgrid models have been shown before to reproduce realistic BCG masses (Martizzi et al. 2012a). Here, we investigate how well our simulations reproduce the observed X-ray properties of galaxy clusters, most notably their density, temperature and entropy profiles and find that, much like in observed systems, our simulations produce a clear CC/NCC dichotomy most clearly seen in the entropy profiles for which we are able to give a physical explanation in terms of major cluster mergers. Reproducing the properties of the hot cluster plasma however has to go hand-in-hand with achieving also compatible results for the full galaxy population in and around clusters. While the field is still far from predictive simulations across the full multiwavelength range covered by current and future observational studies of clusters, we are able to highlight several successes as well as important shortcomings of such state-of-the-art cluster simulations in a full cosmological context.
The structure of this article is as follows. In Section 2, we introduce the sample of zoom simulations that we study as well as the numerical methods and models that we employ. We then compare the ICM profiles of the simulated clusters with various X-ray observations in Section 3. In Section 4, we investigate the origin of the cool-core / non-cool-core dichotomy that we find in our sample. We then compare the stellar properties of the galaxies formed in the simulations with observational data in Section 5. Finally, in Section 6, we study the evolution of the simulated clusters along several mass-observable scaling relations important for cosmology. We discuss the impact of specific modelling choices and the influence of numerical effects on our results in Section 7, before we summarise our results and conclude in Section 8.

DESCRIPTION OF THE SIMULATIONS AND THE NUMERICAL APPROACH
In this section, we describe how we selected a representative sample of massive clusters at fixed mass at z = 0 from a simulation of a large cosmological volume, combining an average and a tail-biased set of objects selected from a larger sample. We also discuss the details of our numerical and algorithmic approaches, in particular the sub-grid models we have employed to account for sub-resolution physics due to cooling and energy injection by stars and massive black holes.

Initial conditions and general approach
The current Rhapsody-G simulation suite includes ten hydrodynamical zoom-in simulations selected from the original Rhapsody sample of massive galaxy clusters . Nine are chosen to have a similar final mass of Mvir ≈ 6×10 14 h −1 M and the tenth has Mvir ≈ 1.3 × 10 15 h −1 M . The original Rhapsody clusters were identified from one of the Carmen simulations from the LArge Suite of DArk MAtter Simulations (LasDamas 1 ). The Carmen simulation on which we base our re-simulations is a cosmological volume of 1 h −1 Gpc. All simulations are based on a ΛCDM cosmology with density parameters Ω b = 0.045 for baryons and Ωm = 0.25 for total matter, ΩΛ = 0.75 for the cosmological constant. The longwave spectral index is ns = 1, the amplitude normalization is σ8 = 0.8, and the Hubble parameter is h = 0.7. The original Rhapsody haloes were identified based on a Gadget-2 (Springel 2005) simulation of the full box using 1120 3 particles. The subsequent N -body only Rhapsody simulations were also carried out using Gadget-2 but using "zoom" initial conditions (using MUSIC; Hahn & Abel 2011) for a sphere of 8 h −1 Mpc centred on each selected cluster at z = 0 with an effective resolution of 4096 3 (4K) and 8192 3 (8K) particles (see Wu et al. 2013b, for details). All initial conditions (both for the original box and all subsequent zooms) were performed using second-order Lagrangian perturbation theory at z = 50.
In this paper, we re-generate the respective initial conditions using Music for ten clusters (see Section 2.2 below for details on the selection and properties of these haloes) but now including baryon perturbations. We assume here that baryons are fully tracing the dark matter perturbations already at z = 50, which is accurate enough for our purposes here (and common procedure), but strictly speaking not correct in detail (see e.g. Angulo et al. 2013, for a detailed discussion). In addition, we use also Lagrangian perturbation theory for the baryons using the local Lagrangian approximation (see Hahn & Abel 2011, for details on baryon ICs for grid and particle codes).

The Rhapsody-G subset of Rhapsody clusters
From the sample of ∼100 clusters in the Mvir = 10 14.8±0.05 h −1 M mass bin of the original Rhapsody sample ), we selected a subset of 10 clusters which we investigate using multi-physics simulations in this paper. To sample both extreme cases of formation history and more average clusters, we considered the two-dimensional ordering shown in Figure 2 of Wu et al. (2013b), where clusters are ranked first by halo concentration c, and then at similar concentration ranked a second time by the numbers of subhaloes N sub above vmax > 100 km/s. The extreme corners of this space are occupied by the systems with IDs 337, 377, 572 and 653. We additionally included system 545 which has similar properties to 337 (high concentration, high substructure fraction). We note that this naturally includes the fossil system 572 which was discussed in more detail already in Wu et al. (2013a) as being a particularly early forming system and occupying the tail in many halo properties. We complement this subset of "extreme" clusters with four more clusters taken from the central region of the c-N sub space. These are the ones with IDs 211, 348, 361, 448. Additionally we include cluster 474, which has a mass that is twice as high as the rest of the sample. These clusters, along with various properties discussed in Section 3.1 are listed in Table 1.

Hydrodynamics, N-body and gravity
In order to follow the non-linear multi-physics evolution of our initial conditions, we use here the Eulerian adaptive mesh refinement code Ramses (Teyssier 2002) and include radiative cooling, as well as sub-grid models for star formation and AGN feedback which we discuss below. Ramses is based on a MUSCL scheme with an approximate Riemann solver and gravity is solved using the multi-grid method directly on the AMR hierarchy. Dark matter as well as stellar and sink particles are evolved using standard N -body techniques.
We employ an overdensity-based (i.e. "Lagrangian") refinement strategy that splits cells if they reach an overdensity of eight, i.e. refinement of the base grid by n additional levels requires a density of 8 nρ . The maximum refinement level for the 4K runs (which constitute the main part of our analysis here) is determined by maintaining a smallest cell size of physical 3.8 h −1 kpc. The 4K dark matter N -body particle mass is 8.27 × 10 8 h −1 M , and the initial mass per hydro cell is ∼ 1.82 × 10 8 h −1 M respectively. The 8K run of RG 653, which we consider here for convergence purposes, has twice the mass and force resolution, i.e. 1.9 h −1 kpc, and an eight times smaller N -body particle mass than the 4K runs. The high-resolution Lagrangian patch from which the 8 h −1 Mpc sphere centred on each cluster will form is tagged using a passive scalar colour field that is advected with the gas. Dynamic refinement is restricted to the regions where this colour field is non-zero. We thus focus all computational resources on the forming cluster and its immediate environment.

Modelling cooling, star formation, stellar feedback and chemical evolution
The simulations include optically thin radiative cooling using the cooling rates of Sutherland & Dopita (1993) for H, He and metal line cooling. Different metals are not evolved separately, but the total gas metallicity is advected with the hydrodynamical equations as a passive scalar and is sourced by the supernova feedback model. A UV background is included assuming the parameterisation of Haardt & Madau (1996) and an instantaneous reionisation at z = 10 taking into account an earlier reionisation in the particularly overdense proto-cluster environment.
The unresolved cold and dense gas that will constitute the interstellar medium (ISM) of galaxies is approximated using a polytropic temperature floor where n * = 0.1 cm −3 is the threshold for star formation (see below), T * = 10 4 K is a characteristic temperature mimicking thermal and turbulent motions in the ISM, and γ * = 5/3 is the effective polytropic exponent. In practice, gas can be heated above the temperature floor, but cannot cool below it. We assume star formation to occur in a cell when the gas density exceeds n * . In this case a star particle is spawned carrying 20 per cent of the mean baryon mass of the cell. We set the local star formation rate per cell aṡ with the free-fall time t ff = (3π/32Gρ) 1/2 of the cell and a star formation efficiency per free-fall time * = 0.02. Stellar particles are seeded locally from a Poisson process. Supernova feedback is implemented based on the model of Dubois & Teyssier (2008), in which each newly formed star particle releases a fraction η = 0.1 of its mass and metals with a yield of y = 0.1 into its surrounding cells through supernovae (SNe) after 10 Myr (this implies that 1 per cent of the time integrated global star formation rate is returned as metals to the ISM). In addition, each supernova injects an energy of 10 51 erg into the surrounding ISM which regulates the star formation efficiency at galaxy scale halo masses. The free parameters of the supernova feedback have been calibrated to reproduce stellar masses consistent with abundance matching results (e.g. Behroozi et al. 2013a) at halo masses M halo 10 12 M for haloes resolved with 1000 particles.

Modelling AGN feedback
The deep potential wells of galaxy clusters require a stronger heating source than supernova feedback to prevent a central cooling catastrophe. Active galactic nuclei (AGNs) provide a natural source of energy in massive galaxies to offset these extreme cooling flows (see e.g. McNamara & Nulsen 2007;Fabian 2012, for reviews). In the simulations discussed in this paper, we include a purely thermal AGN feedback, based on the subgrid models of Springel et al. (2005) with the additional energy injection thresholding of Booth & Schaye (2009), which is commonly employed in SPH simulations (e.g. Le Brun et al. 2014;Schaye et al. 2015, for recent examples). We give a brief summary of our AGN feedback model below but refer the reader to Paper 2, Martizzi et al. 2015, in prep., for a more detailed exposition.

Thermal AGN feedback model
Sink particles, representing supermassive black holes, are created when contiguous regions of high density gas exceed 10 −29 g/cm 3 (identified on the fly by the clustering method of Bleuler et al. 2014) that (1) do not already contain a sink particle, (2) are gravitationally bound, and (3) have an accretion rateṀ clump = M 4 t ff > 30M /yr, where M4 is the mass contained within a sphere of 4 computational cells. This guarantees that sink particles are spawned only in the most massive haloes at high redshift. The initial black hole mass is chosen proportional to the clump accretion rate times the Salpeter time.
Sink particles then accrete at a boosted Bondi-Hoyle accretion rate (e.g. Booth & Schaye 2009) where u is the fluid velocity, cs the sound speed,rB = min(rB, 4∆x) the free-fall limited Bondi radius, and if nH > n * = 0.1 H/cc, 1 otherwise (4) the density dependent boost factor. We additionally impose the Eddington limit ontoṀBH assuming that r = 0.1 of the accreted rest mass energy is converted into radiation. Numerically, the amountṀBH∆t up to a maximum of half of the gas mass contained in each cell is removed from cells inside the sink radius. The thermal energy ∆E = c rṀBHc 2 ∆t would in principle be released into the gas in each time step, but following Booth & Schaye (2009), we accumulate EAGN = ∆E and inject the accumulated EAGN into the sink radius (i.e. the sphere of radius 4 cells) only once EAGN > 3 2 mgaskB∆T , where mgas is the gas mass enclosed by the sink radius. The coupling efficiency c 0.15 has been calibrated to reproduce the MBH − σ relation . In our fiducial model, we set ∆T = 10 7 K, but vary this by one order of magnitude up and down when we investigate the effect of trading fewer violent AGN events against more frequent less violent AGN events in Section 7.1. We emphasise however that the total energy injected by the AGN obviously remains more or less the same. The parameter ∆T has been shown to allow a tuning of the bulk properties of the ICM (e.g. Le Brun et al. 2014) in SPH simulations.

Phenomenological AGN feedback model
Since the AGN energy is thus injected close to the resolution limit into a region containing the densest cells of the simulation, we also consider a more phenomenological model, inspired by the AGN model of Battaglia et al. (2010), in which we distribute EAGN over a resolved sphere of a radius determined by the black hole mass. Battaglia et al. (2010) used a formula for the injection radius that depends on the halo mass. Since we do not track halo masses on the fly with our algorithm and prefer a local criterion, we decided to express the injection radius used by Battaglia et al. (2010) in terms of the black hole mass by using the relation MBH ∝ V 4 c , where Vc is the circular velocity of the halo (Volonteri et al. 2011, noting that for an NFW profile the circular velocity is related to halo mass as Vc ∝ M 1/3 200 ). The radius of the injection sphere is then given by where E(z) is the ratio of the Hubble constant at redshift z and z = 0, i.e. H0, RAGN is larger than 4 cell sizes at most times. This additionally reduces the thermal energy inserted into each individual cell, reducing immediate loss of energy through efficient cooling at the highest densities. Again, we only consider this alternative model, which we term the "phenomenological model", when investigating the impact of changes in the thermal feedback model parameters in Section 7.1.

Halo/galaxy finding and their properties
In order to identify haloes, including their stellar, black hole and gaseous content in our simulation, we use a heavily modified version of Rockstar-galaxies, which is a special version of the original version of Rockstar ) adapted to work also for simulations including gas and stars (see e.g. Knebe et al. 2013, where this version was used). In order to interface Ramses with Rockstar-galaxies, we convert all leaf cells of the AMR tree to pseudo-particles of variable mass, which is however typically not varying over more than an order of magnitude due to the Lagrangian refinement strategy. During the (sub-)halo finding, we then calculate all relevant halo and galaxy properties directly inside Rockstar-galaxies. Unbinding of gravitationally unbound particles and cells is performed, and we obtain the masses, radii, centres, and bulk velocities of the dark matter, stellar, gaseous, and black hole content of each (sub-)halo. We use the Consistent-Trees code (Behroozi et al. 2013c) to link the halo/galaxy catalogues across redshifts based on only the dark matter particle IDs.
We additionally compute several galaxy-related quantities, such as the star formation rate, mean stellar age, surface brightness, and magnitude in various photometric bands. For the star formation rate, we simply sum the mass of all star particles younger than 100 Myr, divided by 100 Myr. For the surface brightness and magnitude calculations, we determine the luminosity of each star particle in a given filter applied to the spectral energy distribution obtained from the simple stellar population models from Starburst99 (Leitherer et al. 1999), assuming the metallicity and age of the star particle as the mean age and metallicity of the stellar population. In order to derive the surface brightness-limited galaxy masses, we perform a one-dimensional projection and integrate the enclosed stellar mass out to a given surface brightness limit. We note that we do not take dust absorption into account, which might affect our comparison with observations. Finally, we also calculate several properties related to the hot plasma, such as the X-ray luminosity and an X-ray emissivity per cell, in order to allow us to weight quantities by their observability in X-ray observations, based on an APEC emission model, similar to Biffi et al. (2012). We neglect however additional sophistications such as actual photon sampling, PSF effects or spectral fitting in order to determine the X-ray temperature, which are clearly beyond the scope of this paper. Finally, we also compute the SZ flux decrement for each cluster. We detail how we perform X-ray and SZ measurements separately in Sections 6.1 and 6.2, respectively, in the context of the evolution of the clusters along scaling relations.

COMPARISON OF ICM PROFILES WITH OBSERVATIONS
We perform a detailed comparison of density, temperature, entropy, and gas depletion profiles with X-ray observations in this section. To show both the mean and deviations over cosmic time of the simulated profiles, we consider individually stacked profiles over time for each simulated cluster for all snapshots for which the cluster mass falls into a narrow range. We discuss this stacking procedure in detail below.

Global cluster properties: masses and cool cores
We first discuss some global properties of our sample of simulated galaxy clusters. To give a visual impression of the diversity of the sample, Figure 1 shows density-weighted temperature and density maps for all clusters at z = 0. In Table 1, we list all the simulated clusters including their z = 0 masses and their half-mass redshift z 1/2 , i.e. the highest redshift at which the most massive progenitor exceeds half the final mass. We quote two z = 0 mass measurements: the virial mass (calculated using the overdensity criterion of Bryan & Norman 1998) and M500c, the mass enclosed inside 500 times the critical density. For the comparison of the cluster ICM profiles with observational data from the ACCEPT sample (Cavagnolo et al. 2009), in Section 3, we decided to select a very narrow mass bin 4 < M500c(z)/10 14 M < 6 in both the ACCEPT sample and the simulated clusters. For the simulated clusters, we take every fifth snapshot out of the ∼ 300 snapshots that we have stored for each cluster for which it falls into this mass bin. We then perform a stack analysis for these snapshots separately for each cluster in order to quantify the mean profiles and the variance around these mean profiles (see Section 3.2 for details). In Table 1, we thus also give the number of snapshots N stack used in the stack analysis of each cluster and the mean redshiftz stack of the stack. N stack is most of the time, but not always, a larger number if the cluster has a high formation redshift. Finally, we indicate whether the cluster is classified as a cool-core (CC) or non-cool-core (NCC) cluster based on whether the mean central entropy of the stack entropy profile is below 40 keV cm 2 at r = 10 kpc. We find that RG 361, RG 377, RG 448 and RG 545 are classified as cool core systems according this classification. for the RHAPSODY-G simulations listed along with the z = 0 masses and the formation redshift z 1/2 . We also list the number of snapshots that enter our stacking analysis and the mean redshift of the stack, as well as whether the cluster is classified as cool-core in the stack analysis (see Section 3.1 for details). The system RG 572 is a fossil cluster, a significant outlier in almost all its halo properties, but is close to a cool-core system.
Halo RG 572 is close to a CC in principle as well, but we will treat it separately as its core properties are even more extreme. The other systems are NCC clusters. We note that we see no immediate connection between either z 1/2 or the time the cluster spends in the 4 < M500c(z)/10 14 M < 6 mass bin and the CC/NCC distinction. Two CC clusters are from the "extreme" sample and two from the "average". The ratio of 4 out of 9 (excluding the fossil system RG 572 as an extreme outlier) is, within the errorbars of small-number statistics, consistent with the number of cool core systems in the ACCEPT sample (6 out of 28) and with other observational estimates at these mass scales (e.g. Chen et al. 2007) of ∼ 40 per cent.

Comparison with ICM profiles from X-ray data
We next compare the gas temperature, entropy, and electron density profiles to observational data. We specifically confront the ICM X-ray profile data from the ACCEPT sample (Cavagnolo et al. 2009). As mentioned above, in order to enable a precise comparison with our Rhapsody-G clusters, we perform a stringent mass cut of the ACCEPT clusters. We combine X-ray and SZ masses to this end to reduce the influence of hydrostatic mass bias at least to some degree. Specifically, for each cluster in the ACCEPT sample, we find its M500 masses from the MCXC catalogue ) and its MSZ mass from the Planck 2015 catalogue (Planck Collaboration et al. 2015b). Cross-matching the three catalogues and selecting all clusters with 4 × 10 14 < M/M < 6 × 10 14 in both M500 and MSZ, we finally obtain a sample of 29 ACCEPT clusters (see Table 1, all but 5 from the Abell catalogue). We additionally excluded Abell 2069 which makes the mass cuts but is a strong outlier in both its entropy and electron density profile as it is a merging cluster. As for the RG clusters, we define cool-core as having a decreasing profile and a central entropy of at most 40 keV cm 2 at r = 10 kpc. We note that the overall trends of the ACCEPT CC profiles are in good enough agreement with the more recent analysis of Mantz et al. (2015) of cool-core systems, most notably the slope of the entropy profiles, for our purposes here.
In Figure 2, we show the comparison between the mass-  Figure 1. Mass weighted temperature and gas density maps at z = 0 of the simulated clusters used in this paper at the 4K resolution, i.e. ρ = ρ 2 dl/ ρ dl and T = T ρ dl/ ρ dl. Notable extreme cases are: RG 337 is a merging system and RG 572 is a fossil cluster, RG 474 has twice the mass of the other systems; RG 337 and RG 377 are richer in substructure than average and RG 653 is poorer than average. The images are 8 h −1 Mpc wide and the projection depth is 4 h −1 Mpc centred on each cluster. Comparison of ICM profiles from Rhapsody-G (red and blue ribbons) with observational data from our mass-matched subset of ACCEPT clusters (black solid and dashed lines). Top: electron density profiles, 2nd from top: entropy profiles, 3rd from top: temperature profiles, and bottom: total mass profiles. The coolcore clusters in the ACCEPT subset are indicated by dashed lines, while the cool-core Rhapsody-G clusters are shown in blue and the non-cool-core in red. The simulated haloes were stacked individually for each cluster while occupying the same mass range as the ACCEPT subset. The shaded ribbons indicate the 1σ scatter in each stack, reflecting time variations in the profiles. We excluded the fossil cluster RG 572 (see text).
selected ACCEPT sub-sample and stacked profiles from the Rhapsody sample. For each simulated cluster, we show the mean profile as well as the standard deviation around that mean profile, obtaining one ribbon for each Rhapsody cluster in Figure 2. For this analysis, we excluded the fossil cluster RG 572 since it shows an extremely fluctuating entropy and core density, caused by a very active and rapid growing central hypermassive black hole of ∼ 3 × 10 11 M at z = 0. Note that in the full Rhapsody sample of more than 100 clusters this was still a completely abnormal case. This cluster is so atypical in its core properties already in pure N -body simulations that we only note here that it is also an outlier in its baryonic core properties.
The number of profiles stacked for each cluster is different and given in Table 1. In all cases, we find that the dispersion in each stack around the mean profile is smaller than the difference between the cool and non-cool-core profiles, indicating a stable bimodality. This can be clearly seen from the small extent of the ribbons in Figure 2, which indicate a very small amount of scatter around the mean in the stack of each cluster over time. The cool cores are thus a much more significant feature than short term fluctuations in the profiles. The non-cool-core profiles from ACCEPT and Rhapsody agree well within their respective scatter, for both the electron density profiles and the entropy profiles. The ACCEPT cool-core clusters however show a much weaker cool core than the Rhapsody CC clusters: the observed clusters show only a moderate increase in core electron density inside the innermost ∼ 50 kpc, while the Rhapsody cool cores show a strong drop in entropy and increase in electron density already at scales of 150 kpc. It is clear that the cool core systems are most likely still undergoing overcooling to some degree despite the central AGN that is efficiently fuelled during the CC phase. The CC/NCC dichotomy is thus a long lived property of our clusters, consistent with the observational constraints of e.g. McDonald et al. (2013). The dichotomy arose naturally in a larger sample of cosmological cluster simulations and can be explained by differences in the assembly history and nature of major mergers of the clusters (see our analysis in Section 4).
Interestingly, outside the core, the temperature profiles show a somewhat discrepant temperature at large radii, specifically a ∼ 30 per cent difference in X-ray temperature at 200 kpc, with much better agreement at smaller radii. Upon closer inspection one notices that a similar but slightly weaker offset also exists in the power-law part of the entropy profile, where the simulated profiles are systematically offset to lower entropy. A similar discrepancy can be seen also, e.g., in Figure 7 of Dubois et al. (2011), independently of the much larger range of subgrid models employed in that reference. This may point towards additional physics missing from these simulations.
In order to investigate better the nature of the cool cores, in Figure 3, we plot cooling and free-fall time profiles (c.f. also Fig. 9 of Li et al. 2015). We see that for the non-cool core systems, the cooling time is at least a factor of 10 above the freefall time in the core, while for the cool core systems, the two time scales are much closer to each other. Normally the central AGN should increase the central cooling time dramatically, but it appears to not prevent the formation of cool cores in four of our systems. This is despite a dramatic growth that the central black hole undergoes in the CC cases. For one of them, RG 545, we demonstrate in Section 7.1 that the existence of the cool core does not depend on the details of the AGN feedback parameters. We are thus led to believe that in our simulations  . Ratio of the hydrostatic and the gravitating mass as a function of radius. Outside of the core and well inside the virial radius, the clusters are very close to hydrostatic equilibrium, with large deviations at radii 1 Mpc. There is no obvious difference between the CC systems (blue) and the NCC systems (red). The grey shaded band indicates ±20 per cent deviations from hydrostatic equilibrium.
the CC/NCC dichotomy has a cosmological origin, and will investigate this aspect further in Section 4.

Hydrostatic mass and deviations
The cluster is stabilised against gravitational forces by thermal pressure as well as bulk and random flow motions (e.g. Nelson et al. 2014). We next investigate to what degree the thermal pressure alone counteracts gravity in our simulations. In adaptive mesh simulations with a Lagrangian refinement strategy (which is however indispensable if galaxy formation is included), the numerical Reynolds number is relatively high so that no turbulent cascade develops and the energy is very quickly dissipated into thermal energy. In contrast to uniform high resolution simulations with fully developed turbulence (Miniati & Beresnyak 2015), we can focus on thermal effects here inside the outer infall region. We follow Nelson et al. (2014) and define the thermal (or hydrostatic) mass as where the angle brackets denote mass-weighted averages over spherical shells at a given radius r. In Figure 4, we show the ratio M therm (r)/Mgrav(r), where Mgrav(r) = Mtot(< r) is the total mass enclosed in radius r. All clusters are closest to hydrostatic equilibrium at radii of ∼ 300 kpc with a mean deviation of at most 20 per cent inside of ∼ 600 kpc. At large radii, the variance is however substantial, most likely due to incomplete thermalisation of the region which is perturbed by accreting gas. Most importantly for our analysis, we find no significant difference in deviations from hydrostatic equilibrium over radii 100 kpc r 1 Mpc between CC and NCC systems.
There is however an indication of systematic differences at 1 Mpc, as well as, obviously, at scales of the core.

Gas depletion profiles
A crucial observable reflecting the degree to which collisional matter follows the total mass distribution of galaxy clusters is given by the gas depletion profiles. Especially in the cluster cores, cooling and AGN feedback will affect both the gas and the dark matter distribution. In particular, if the central black hole accretes gas quasi-periodically, AGN feedback can transform the cuspy dark matter profiles into cored profiles (Martizzi et al. 2013, consistent with what we find in the RHAPSODY-G simulations). In Figure 5, we show the ratio of enclosed gas mass to total mass in units of the universal baryon fraction. We split our sample of clusters into non-coolcore (top) and cool-core (bottom) systems. We compare our results to the observational relations of Mantz et al. (2014), as indicated by the blue ribbons. The NCC systems follow the observational relation reasonably well. The fossil cluster and the CC systems however show a very large central gas fraction, in many cases even above the universal gas fraction (similar to the results of Burns et al. 2008). This result is clearly in tension with the results of Mantz et al. (2014), and reflects the high electron densities and very low central entropies we have seen for the CC systems above. The conclusion must be that in the case of the CC systems, the AGN feedback is not able to stabilise the core at levels consistent with X-ray observations. The exciting possibility is that other forms of non-thermal feedback (or processes) must be plausibly involved in order to bring these results in line with observations. We note that our results for NCC systems are consistent with published results from SPH simulations (e.g. Battaglia et al. 2013;Sembolini et al. 2013;Planelles et al. 2013) at scales of 0.1 Rvir, but are somewhat higher at larger radii.
Less visible but as significant is the discrepancy of the gas fraction at ∼ Rvir/2, where the profiles are slightly but systematically high with very little scatter. As we see below in Section 17, the thermal AGN feedback does not reach these large radii in our simulations, and it is not possible to tune the energy injection parameter ∆T to deplete the cool cores efficiently.
Despite these discrepancies, the use of the gas fraction at R2500 as a robust cosmological observable (Allen et al. 2008;Mantz et al. 2014) is strongly supported by our simulations since R2500 appears to be outside the reach of the AGN and shows a virtually unbiased thermal mass (c.f. Section 3.3).

THE ORIGIN OF THE COOL-CORE / NON-COOL-CORE DICHOTOMY
We next inspect the origin of the CC/NCC dichotomy in our simulations. We find that at early times, z ∼ 0.6 and z ∼ 0.8, respectively, both RG 348 and RG 545 are cool-core clusters and have a very similar assembly history, which we show in Figure 6. Both undergo a major merger with a similar mass ratio at similar times. However, the cooling time of the core (top panel of Figure 6) rises dramatically more in the case of RG 348 to a value of about 5 Gyr after an initial higher spike. Cluster RG 545 only experiences an increase in cooling time to 2.5 Gyr. We define the core here as the gas enclosed in an overdensity of 8 × 10 4 and evaluate the cooling time at that radius (we find that the free-fall time reduces only slightly over the same time scale so that roughly t cool /t ff ∝ t cool ). After the merger, both clusters have a more quiescent merger history Figure 6. Assembly history of RG 348 and RG 545. The first is a non-cool-core cluster at z = 0, the latter is a cool core cluster. The top panel shows the cooling time evolution in the core for the two systems, the bottom panel the mass accretion history. Both undergo a major merger between a ∼ 0.5 and a ∼ 0.6. While the core of RG 348 is substantially heated during the merger, the effect on RG 545 is not as dramatic. The vertical lines indicate the times for which we show images illustrating the evolution of RG 545 and RG 348 in Figure 7 and 8, respectively.    Figure 8. Cool core destroyed in a major merger (RG 348). We show the minimum entropy along the line-of-sight (left column) and the density (middle left) of gas with a temperature of at least 0.1 keV, as well as the density weighted temperature in linear scale (middle right) and the dark matter density (right column). The first row is at z = 0.61, the next rows are 171, 330, and 448 Myr later. The top left panel clearly shows the two cool cores of the systems before the merger. During the merger, the cores undergo a strong shock (clearly visible in 2nd temperature panel from top). After the collision, the cool cores are dissipated and undergo mixing with higher entropy gas. Solid and dashed circles indicate the cores of the main cluster and the merging cluster, respectively. The many low entropy blobs are satellite galaxies that do not carry much mass with them. The depth of the images is comoving 1 h −1 Mpc centred on the more massive cluster. and both cores start to cool again. The cooling of RG 545 is much more rapid, plausibly since it is quickly experiencing runaway cooling, while RG 348 cools at a slower rate -the time to z = 0 is still 7 Gyr, shorter than its cooling timemost likely reflecting AGN heating that is able to keep the core hot.
After inspecting the time evolution of the systems, the difference in how the mergers proceed is obvious. In the case of RG 545, see Figure 7, the cores do not collide since the merger has large angular momentum, so that the cool core is perturbed but not destroyed and is just sloshing in the centre of the halo. Quite different is the merger that RG 348 experiences: here the two cool cores of the progenitor clusters collide and dissipate within ∼ 100 Myr. We show their time evolution over ∼ 450 Myr during the merger in Figure 8: the left column shows the minimum entropy along the comoving 1 h −1 Mpc image depth of X-ray emitting gas with a temperature above 0.1 keV. We also show the density of gas above 0.1 keV, as well as the temperature (in linear scale) and the dark matter density. For RG 545, we only show entropy and density since the merger is less spectacular. From the sequence of RG 348, one sees that the cores pass right through each other (in the moment shown in the second row); after the collision, the cool core of the main cluster halo is destroyed, while the one of the less massive progenitor survives weakened. In the last time shown and subsequently, the surviving core mixes quickly with the hot cluster gas due to ram pressure and Kelvin-Helmholtz instabilities. We also note that, as expected, we detect a weak displacement between the BCG as well as the dark matter halo centre and the gas centre after the collision, consistent also with the observational findings of e.g. Semler et al. (2012).
The finding of a connection between CC/NCC and mergers in early simulations put forward by Burns et al. (2008) and Planelles & Quilis (2009) thus should be complemented by this new finding. Not only slowly accreting systems can have cool cores; the amount of angular momentum in major mergers (i.e. whether the collision is head-on in such a way that the cores do pass through each other) appears to be of high importance as well. This is also consistent with the findings from idealised non-cosmological cluster merger simulations by Poole et al. (2008). Major mergers are thus a necessary but not sufficient condition for the destruction of cool cores. The survival of RG 545 as a CC system clearly shows that the cores are not easily disrupted, in some cases, even by a major merger. For all clusters in our sample it is true that, very similar to the two examples in Figure 6, the cool core clusters never had their cooling time increased beyond at most 20 times the core free fall time. This is roughly consistent with Meece et al. (2015) who find that a ratio of at most 10 is needed in order to undergo a significant thermal instability. Since, absent other heat sources, every cluster will slowly cool again after a major merger event (and basically every cluster will experience one or more major mergers throughout its history), the sharp distinction between the CC and NCC systems in terms of their profiles then arises simply due to how much a merger perturbs the core -determined by the merger's angular momentumand whether or not rapid cooling can set in subsequently. In the case of RG 348, the cool core has been dissipated in no more than a few 100 Myr with a subsequent cooling time that is significantly higher. The disparity in these two time-scales makes this a quick transition, allowing the bimodality apparent in the cluster population to be sustained.

PROPERTIES OF THE GALAXIES IN AND AROUND THE SIMULATED CLUSTERS
We next discuss the properties of the galaxy population in and around our simulated clusters. Since the high-resolution region is an 8 h −1 Mpc sphere around each cluster at z = 0, we have a sizeable sample of lower mass haloes ranging from lower mass clusters to group to galaxy scales in the outskirts of our massive clusters. This facilitates a statistical comparison with the more field-dominated observational samples. Nevertheless all galaxies in our simulations do live in a high-density environment, which should be kept in mind when comparing these results to observations.

Comparison with abundance matching results
Apart from the hot intracluster plasma, galaxy clusters harbour a large population of satellite galaxies. A general short-coming of galaxy cluster simulations is that in these simulations it is computationally not affordable to resolve the length scales relevant to important aspects of galaxy formation (on say 100 pc and smaller). The hope is that at least the most massive galaxies are resolved well enough. A bona-fide test of the realism of our simulations is thus to investigate whether the central cluster galaxies indeed have stellar masses and other properties that are consistent with observations. Remarkably, despite the fairly poor resolution at galaxy mass halo scales, we find that the star formation efficiency and supernova feedback as dictated by the subgrid models is in fact able to reproduce accurate stellar masses across all halo masses that we resolve with 1000 particles in our simulations. In Wu et al. (2015), based on the same simulations discussed here, we find stellar mass fractions of ∼ 10% of the baryon fraction, consistent with observational constraints.
In Figure 9, we show the comparison of the stellar masses of our galaxies (satellite and central) as a function of their halo mass with results obtained using the abundance matching technique (c.f. e.g. Conroy et al. 2006;Moster et al. 2010;Behroozi et al. 2013a). The very massive systems that we consider here occupy however the very tail of the stellar mass functions, and in fact, the older abundance matching relations (e.g. Moster et al. 2010;Behroozi et al. 2013a) do not capture very well the scaling at the highest masses due to the way stellar mass was counted in the stellar mass functions used. More accurate results have been obtained by Kravtsov et al. (2014) taking explicit care of the intracluster light component; these are fully consistent with the the relations of Behroozi et al. (2013a) when updated with the stellar mass functions of Bernardi et al. (2013). We show a comparison of these relations with the stellar masses from our simulations measured using a range of techniques. In particular, we have produced mock V-band images and count the stellar mass inside of isophotes where the surface brightness is above 25 and 24 mag arcsec 2 , as well as using a simple spherical density threshold of 2.5 × 10 6 M /kpc 3 (panels from left to right in Figure 9). We see that, while we largely overproduce stellar mass com- Figure 10. As Figure 9, but a comparison between the M * -M 500c relation for galaxies from a higher-resolution simulation of cluster RG 653 (simulated at 8 times better mass and twice better spatial resolution than the full sample) at z = 0 and abundance matching constraints.
pared the original Behroozi et al. (2013a) result, our simulations are consistent with more recent estimates for massive systems. The stellar masses in massive haloes above M500c ∼ 2×10 12 M are still somewhat high, but not dramatically so, and details depend on the exact definition of how stellar mass of central group and cluster galaxies should be measured. Notably, the stellar masses at galactic halo scales seem to undershoot the observational relations. However, this is a resolution effect. In Figure 10, we show the respective plot for the galaxies in and around cluster RG 653 simulated at twice better spatial and eight times better mass resolution. With higher resolution, the full range of halo masses from 10 11 − 10 15 M is in good agreement, with the slight overproduction of stars in the most massive haloes remaining.

Star formation rates
We next want to investigate the reason for the somewhat high stellar masses in the most massive haloes in our simulation. To this end, we show each of our simulated galaxies (for all clusters) in the stellar-mass vs. star formation rate plane. We estimate our star formation rates as the mass of star particles in each simulated galaxy that have an age of less than 100 Myr divided by the 100 Myr. Galaxies which do not have any star particle fulfilling this condition get assigned a random star formation rate (SFR) between 10 −2 and 10 −3 M yr −1 . We perform this step in order to allow a visual estimate of the fraction of quenched galaxies in Figure 11.
From Figure 11, we see that our galaxies with significant star formation follow a narrow star forming sequence that is consistent in normalisation and slope with the SDSS star forming sequence as measured by Woo et al. (2013) (solid black line). In addition, we show the star formation rates determined by Liu et al. (2012) for a sample of central cluster galaxies from SDSS, and those of McDonald et al. (2015) for BCGs of SPT clusters. When comparing with the limit below Figure 11. Star formation rate as a function of stellar mass for central (orange) and satellite (blue) galaxies from all 10 Rhapsody-G haloes at z = 0. Galaxies with no detectable star formation in the previous 500 Myr have been assigned random values between 10 −3 and 10 −2 M /yr. For comparison, the main sequence of star forming galaxies in SDSS from Woo et al. (2013) (solid black line) and the division between star-forming and quenched galaxies (black dotted) along with the OII-derived star formation rates for the BCG samples of Liu et al. (2012) and the UV estimates of McDonald et al. (2015) are shown (note that for the latter, data points with only long lower error bars are upper limits).
which Woo et al. (2013) would classify galaxies as "quenched" (black dashed line), we see that many more massive galaxies in this simulation are not quenched -the SDSS quenched fraction exceeds 50 per cent above M * ∼ 10 10.5 in Woo et al. (2013) (private communication from J. Woo). While the SDSS results are averaged over all field galaxies, it may be that some biases exist in these most dense environments around the most massive clusters. However, it seems more plausible that these simulations either miss an additional physical mechanism or lack the resolution to reproduce a realistic quenched fraction at the intermediate masses. While the highest star formation rates we observe are roughly consistent with Liu et al. (2012) and perfectly consistent with McDonald et al. (2015) (note that the data points with one-sided lower error bar are upper limits only), it is clear that quenching of intermediate (and maybe even high) mass galaxies is inefficient in our simulations. This is consistent with the somewhat high stellar masses compared to the abundance matching results. Finding the mechanism(s) necessary to suppress star formation in massive galaxies is among the most pressing tasks of galaxy formation theory today and it is not surprising that our simulations are not performing significantly better than better resolved simulations dedicated to study galaxy formation.

IMPLICATIONS FOR CLUSTER COSMOLOGY: EVOLUTION ALONG SCALING RELATIONS
The use of cluster counts as a cosmological probe exploits population scaling relations that link total system mass to observable properties of the galaxies or hot plasma. In this section, we compare the evolution of the simulated clusters along a range of scaling relations with observational results from X-ray and CMB data and some other published results from simulations. We focus on the Sunyaev-Zeldovich effect as well as on X-ray luminosity and temperature. In the absence of a large statistical sample of clusters, tracking a few clusters over time allows us to relate the role of evolutionary processes to their impact on observables.
Since systems in the simulated sample have roughly the same mass at z = 0, their evolutionary tracks can be used to assess the contribution of their assembly history to scatter in the mass-proxies. In order to study the evolution of the scatter or additional biases relevant to the full population, however, larger samples would be required. While we compare our scalings to observations, we caution that this type of exercise is non-ideal in that: i) we employ true, three-dimensional spherical masses in simulations while observational estimates may be biased with respect to these values; ii) similarly, the intrinsic properties of the simulations are measured directly in the simulations rather than derived from models applied to mock observations, and; iii) the statistics of the simulated sample are not necessarily representative of a broader mass-complete sample. The purpose of this exercise is to identify areas of agreement and discrepancy, which in turn should provide insights into future adjustments to the astrophysical feedback model.

Sunyaev-Zeldovich masses
We first compare the Compton-y parameters integrated over our simulated clusters, i.e. calculated as with SZ masses and Y5R500 measurements from the Planck HFI union catalog 2015 (R2.08), see Planck Collaboration et al. (2015a). We have converted the Planck Y5R500 values to Y500 using the ratio of 1.814 between the measurement at 5 and at 1R500 quoted in Melin et al. (2011) and assumed the Planck cosmological parameters when converting Y5R500 to units of kpc 2 . We show the results of this comparison in Figure 12, along with the unbiased (i.e. setting the hydrostatic mass bias to b = 0) Planck2013 mean baseline relation E −2/3 (z) Y500 kpc 2 = 10 1.81±0.02

X-ray observables
Our simulations show good agreement with the Planck2013 baseline for all systems over the entire range from ∼ 5 × 10 13 to 10 15 M . Despite the large variance between the clusters' assembly histories, we find very little scatter, ∼ 0.2 dex around the Planck2013 mean relation (consistent with other simulations, e.g. Sembolini et al. 2013), with only the fossil cluster The Planck data is taken from the 2015 union catalog with rescaled Y 500 = Y 5R500 /1.814 according to Melin et al. (2011). Figure 13. Comparison of the M 500 − L X scaling relation with observations and other simulations. The RG clusters are again shown as evolutionary tracks, the luminosity is the integrated 0.1-2 keV luminosity estimated from APEC spectra. For the clusters from the MCXC sample we have used the 0.1-2 keV luminosities from MCXC but the respective Planck SZ masses here. The dashed-dotted line represents the fit from Arnaud et al. (2010), which is almost identical to the scaling relation used in MCXC, and the dashed line is the best fit self-similar solution (i.e. the ∝ M 4/3 500 fit to the MCXC/Planck joint data). The dotted line shows the best fit from Biffi et al. (2014), who however use the wider Chandra band, while the blue hatched region is from Le Brun et al. (2014), and uses also the 0.1-2 keV band. Figure 14. Comparison of the ACCEPT/MCXC clusters M 500 − T X relation with the evolutionary tracks of the RG clusters. From left to right, we show comparisons with different estimates of the gas temperature in the simulations: (Left) the spectroscopic-like temperature, with core excluded, (Center) volume-weighted temperature and (Right) the mass weighted temperature. The gray dots are a combination of X-ray masses from the MCXC catalog with temperature estimates taken from the ACCEPT catalog, and the black long-dashed line is the best fit power law to this data, while the short-dashed line indicates a power-law fit to RG 474, which covers the largest range in mass evolution.
RG 572 being somewhat of a significant outlier (though still insignificant given the scatter in the Planck M − Y data). In fact, it is remarkable that RG 572 does not stand out. In particular, the scatter is substantially smaller than the scatter around the M − Y scaling relation in the Planck2015 SZ catalog. The large observed scatter is mostly due to complications from projection and low angular resolution, which requires joint fitting of Y5R500 and cluster size θ500. The large and markedly asymmetric scatter in the Planck M − Y data is an artefact of our simplistic rescaling from Y5R500 to Y500 with a constant factor and appears in rather stark contrast to the simulation results. Investigating the simulation prediction for Y5R500 would require a matched filter approach to reduce interloper effects at larger radii, as employed by the Planck team, and possibly full lightcone projections (e.g., White et al. 2002), which are beyond the scope of this work.
In Figure 13, we show a comparison of the X-ray luminosity of the RG clusters during their evolution as a function of their mass with various scaling relations from the literature as well a direct comparison with clusters from the MCXC sample ) using the Planck SZ masses and the MCXC luminosity (since otherwise one simply reproduces the scaling relation used to estimate the masses, which is very close to Arnaud et al. 2010). We measure the X-ray luminosity by computing the emissivity per cell using the APEC tabulated emission model and exclude the central core region inside 0.15 R500, i.e.
where ne(x) and nH(x) are the electron and hydrogen number density, T (x) is the electron temperature, and (T, Z) is the (temperature-and metallicity-dependent) emissivity taken from the APEC tables and integrated over the same energy range 0.1 − 2 keV as the MCXC clusters. The exclusion of the core is necessary as it introduces a strong variability of the X-ray luminosity due to AGN events where the central density fluctuates and unrealistically high luminosities can be obtained. The samples that we compare against however do include the core, so that the comparisons have to be taken with a grain of salt. For comparison, we also show the results of Le Brun et al. (2014), who use the same AGN model in SPH simulations and whose X-ray luminosity appears to be somewhat low compared to the MCXC sample. For completeness, we also show the best-fit relation from Arnaud et al. (2010), which is basically identical to the scaling relation used for the M500 measurements in MCXC, as well as the best fit from Biffi et al. (2014). The latter authors use the Chandra bolometric band (0.3 − 12 keV), which explains the higher Xray luminosity compared to the MCXC data. In comparison, our clusters evolve along a scaling relation that is somewhat shallower than both Arnaud et al. (2010) and Le Brun et al. (2014), consistent with Biffi et al. (2014) (who however used the larger spectral window) but slightly steeper than a simple self-similar scaling ∝ M 4/3 (shown as a green dashed line in the figure and whose normalisation has been obtained by fitting to the MCXC data). In fact, the best-fit slope in our case is closer to 1.2, but we do not want to make a more quantitative statement at this stage since the comparison of evolving small (autocorrelated) samples should be taken with some caution for making rigorous predictions. It suffices to observe here that the simulation X-ray luminosities are consistently higher (by about a factor of two around 10 14 M and at most ∼ 20 per cent at 10 15 M ) than the MCXC luminosities, with a scaling relation that is slightly steeper than the self-similar expectation of 4/3. There is a weak indication that the normalisation of the evolutionary track of clusters is persistent, i.e. RG 211 tends to be lower than RG 361, which is lower than RG 348 for a range of masses. This is plausibly indicative of a connection between X-ray luminosity and assembly history at fixed mass. We will investigate the origin of scatter in future work. Once again it is remarkable that RG 572 does not stand out dramatically, its core-excised X-ray luminosity is only slightly lower than the rest of the sample. Since the use of Planck M500 shifts the MCXC clusters away from the Arnaud et al. (2010) relation, we expect that hydrostatic mass bias will play an important role when comparing the RG clusters to observations. To conclude our comparison with X-ray scaling relations, in Figure 14, we show the mass-temperature relation for a subset of the MCXC clusters for which we could take the Xray temperatures from the ACCEPT cluster catalog. Due to the known biases in estimating temperatures from simulated clusters that reflect temperatures measured from X-ray spectroscopy (see e.g. Biffi et al. 2014, for in-depth comparisons), we show a range of differently weighted temperature estimates in the different panels of the figure. First, we compare with the spectroscopic-like temperature as introduced by Mazzotta et al. (2004) as as a fit to spectroscopic temperature estimates with Chandra or XMM-Newton. We additionally do not include cells with a temperature below 0.5 keV in this average, which is necessary in the case of cooling to avoid the inclusion of non-X-ray emitting gas. Furthermore, we again exclude the core inside 0.15 R500 in order to avoid even larger fluctuations in the estimated cluster temperature due to the influence of the AGN model on this region. This is also common practice in observations to reduce scatter in the scaling relations (e.g. Mantz et al. 2010). We found that a temperature estimated by computing the mean energy of all photons emitted using an APEC model (see discussion above) above 0.1 keV leads to basically identical results within the scope of this first analysis, so that we do not consider these 'mock' X-ray observations in this paper further. The other two temperature estimates we consider are a simple mass and a volume-weighted average of all cells with no core excluded and no cut in minimum gas temperature, so that they reflect the true mean temperatures, not including possible observational biases.
We note that the spectroscopic-like (as an APEC emissivity weighted temperature not shown here) temperature has a much larger variance compared to the volume-and mass-weighted averages. The variance substantially increases if the core is also included due to the n 2 -dependence of the emissivity. In all cases, and thus plausibly not affected by potential biases compared to observational results, the M − T scaling that we observe is slightly steeper than the one exhibited by the MCXC/ACCEPT data. In particular, we find that the cluster temperature tends to be somewhat low for masses above 10 14 M . The discrepancy, while systematic, is not dramatic. It is consistent with our comparison of the temperature profiles with the ACCEPT data, in which we found a similar temperature discrepancy at large radii.
To summarise, in no case did we find a very strong dependence of any cluster scaling relation on the assembly history. A weak dependence might be possible for the M − LX relation in our data. This result is consistent with earlier analysis in this direction, by e.g. Jeltema et al. (2008), who however found a dependence on the dynamical state once the scaling relations are not compared against the true mass but against, e.g., the hydrostatic mass. This is very plausibly the case for our simulations as well, since the hydrostatic mass bias can be large close to R500c (see our analysis in Section 3.3). It is thus not entirely clear if and how the ACCEPT comparison may be biased by hydrostatic mass errors which could alleviate also the discrepancies we observed between the profiles in Section 3.2. Figure 15. Impact of the AGN model on X-ray scaling relations: As Figure 13, but only for the cluster RG545 using two different implementations and a range of parameters of the thermal blastwave AGN model. Specifically, we compare the phenomenological injection model ('pheno') that has a well-resolved injection bubble with the standard Booth & Schaye (2009) model with injection close to the resolution scale. In both cases we vary the thermal energy accumulation threshold ∆T . No significant impact is seen on the bolometric core-excised X-ray luminosity.

DISCUSSION OF THE NUMERICAL MODELLING APPROACHES
We have established already in the previous sections that our results are numerically converged by comparing the ten runs at 4K resolution to the higher resolution 8K run of RG 653. Next, we investigate the dependence of our results on particular choices of the AGN feedback model parameters and discuss discrepancies with published results based on Lagrangian methods.

AGN model dependence: X-ray properties
In this section, we investigate the robustness of our results to changes in the thermal AGN feedback model. In particular, we want to see whether the cold cores we find are stable to such changes. We furthermore will investigate the degree of change in the scaling relations that can result from such astrophysical changes beyond the resolution of the simulation in the subgrid model. While it is clear that the very short cooling times in the centres of massive clusters necessitate a mechanism to prevent efficient cooling in their cores in order to reproduce both the observed gas fractions and realistic galaxy masses, what this mechanism is, or how it should be modelled on the scales accessible to three-dimensional cosmological simulations is less clear. Virtually all current state-of-the-art cosmological simulations that tap into this mass regime invoke some form of AGN feedback to resolve this problem. As shown in Figure 5, the AGN model we adopted fails however to reduce the gas fraction at radii 0.5 Rvir, where our simulations predict somewhat high baryon fractions compared to e.g. Mantz et al. (2014). If AGNs however have such a dramatic effect on the entire ICM, most of cluster cosmology would have to rely on a very tight relation between cluster mass and AGN feedback energy. The analysis of Le Brun et al. (2014) e.g. shows that the entire baryon gas profile can be varied by tuning the energy accumulation threshold ∆T of the feedback model of Booth & Schaye (2009). While this is possible, it should be expected that AGN feedback could be a significant source of scatter, relating processes at the pc scale to the Mpc-scale of the ICM. In this section we thus repeat this analysis and investigate how sensitive our results are to variations of the energy accumulation parameter ∆T .
In Figure 15, we show again the evolution of the X-ray luminosity of a cluster as a function of its mass over time. For this analysis, we focus on the cluster RG 545, which we found to be a cool-core system during the time it is in the mass-range between 4 < M500/10 14 M < 6. The selection of a CC cluster for this analysis will allow us to investigate later how robust the existence of the cool core in the simulation is when parameters of the AGN model are varied. We have not found any large differences between the CC and non-CC systems in neither the M − LX or M − T scaling relations once there cores are excluded.
When changing ∆T over two orders of magnitude from 10 6 K to 10 8 K, we see no clear systematic effect on the Xray luminosity. This appears in tension with what Le Brun et al. (2014) find (their Figure 1b), where the X-ray luminosity decreases by 0.7 − 0.8 dex when ∆T is increased by a similar amount. We caution that we only studied the response of a single system here, while Le Brun et al. (2014) study the impact on the whole sample of clusters in the Cosmo-OWLS volume of 400 h −1 Mpc. As expected, the phenomenological injection model (labeled 'pheno' in the figures, cf. Section 2.5) is less bursty in X-ray luminosity compared to the fiducial injection method, since the energy is distributed over a larger volume. In conclusion, it does not appear possible in this model to tune the normalisation of the LX − M relation in any significant way using AGN feedback parameters in our simulation.
The excess X-ray luminosity is however consistent with the findings of e.g Choi et al. (2015), who find that in their simulations of group-scale haloes with thermal AGN feedback, the X-ray luminosity is a factor of ∼ 50 − 100 higher than observed. This excess luminosity is reduced, and becomes consistent once they employ a kinetic feedback model.

AGN model dependence: stability of cool cores and gas depletion
Much like in our sample, Burns et al. (2008) (and Planelles & Quilis 2009 found the CC/NCC dichotomy to naturally arise in samples of cosmological cluster simulations. Unlike in RHAPSODY-G, their simulations were of much coarser resolution (AMR min cell size of 15.6 h −1 kpc and N -body particle mass of 9 × 10 9 h −1 M ) and did not include AGN feedback which has left room for speculations about different origins of the CC/NCC dichotomy. Of course, this difference w.r.t. our simulations is crucial, since increasing the spatial resolution increases the severe central cooling catastrophe and leads to unrealistic BCG stellar masses providing the main argument for the necessity to include AGN feedback as a central energy source to compensate cooling losses and bring BCG masses into realistic ranges (c.f. Martizzi et al. 2012a,b). It is thus a non-trivial result that the cool cores survive this energy Figure 16. Impact of AGN model on ICM profiles: As the top three panels of Figure 2, but only for the CC cluster RG545 using two different implementations and a range of parameters of the thermal AGN injection threshold. The phenomenological model has the most significant impact by raising the entropy also at small radii, changing the electron density to a power law in the core, and flattening the temperature profile. The lowest injection threshold ∆T = 10 6 K produces frequent outward travelling blast waves that show up in the entropy and temperature profile but only have a modest effect on the density. In no case radii outside ∼ 100 kpc are affected however.
injection, and one may wonder whether their survival is only due to particular choices of the feedback model parameters and whether more rare but violent, or more frequent but less violent AGN events might destroy them. This question is of particular importance not only because one may wonder about the model parameter dependence but also since the dominant role of AGNs in CC/NCC transitions has been advocated in the literature (Guo & Oh 2009;Guo & Mathews 2010). We next investigate the robustness of the cool core clusters to changes in the AGN energy injection threshold. We note once again that the injection threshold ∆T does not control the total energy injected into the ICM, but only its portioning. In Figure 16, we show the dependence of the electron density, temperature and entropy profiles on ∆T -to be compared with Figure 2. We find a non-monotonic dependence of the central core slope of the electron density profile on ∆T . In particular, the lowest threshold we considered, ∆T = 10 6 K yields a cored profile with a central density of ∼ 0.1cm −1 , Figure 17. Impact of AGN model on gas fraction: As Figure 5, but only for the CC cluster RG 545 using two different implementations and a range of parameters of the thermal blastwave AGN model. The phenomenological model has the most significant impact to keep the gas fraction below the universal baryon fraction also at small radii, but in all cases the gas fraction at 0.1R vir is high compared to the constraints from Mantz et al. (2014). In contrast to the non-cool-core case, the thermal blast wave here is insufficient to substantially reduce the gas fraction in the inner parts of the halo. more consistent with the ACCEPT observational constraints. At the same time, the frequent AGN bursts show up as outward travelling shock waves in the entropy and temperature profiles inside 150 kpc. In all cases considered there is no effect outside that radius on either entropy or density, nor is the entropy profile changed to that of the average NCC systems. This result lets us conclude that the formation of the cool core cannot be prevented by the central thermal AGN feedback model, regardless of the injection threshold and region we considered. The additional overcooling of the cool core can be somewhat tuned but not alleviated by this feedback model. Whether kinetic feedback can resolve this problem (c.f. Li et al. 2015) is an interesting possibility to be investigated in future work.
A robust property of our cool core systems that we identified above was the high central gas fraction of order the universal baryon fraction inside the cool cores -inconsistent with the observations of e.g. Mantz et al. (2014). We thus also want to investigate the impact of the AGN model parameters on the gas depletion profiles. In Figure 17, we show essentially the same plot as in Figure 5, varying again the injection threshold parameter ∆T . While we saw that the lowest injection threshold ∆T = 10 6 K was able to impact the entropy and density profiles, it also reduces the central gas fraction, but only very slightly so, and only at radii Rvir/10. Similarly, the larger injection regions of the phenomenological model reduce the central gas fraction below the universal value. Once again, in no case could we affect regions outside the core leaving the gas fraction inconsistent with observations at ∼ 0.2 Rvir. The energy available from the thermal blastwave, even after accumulating energy to heat to higher temperatures is completely insufficient to reduce the gas fraction outside the core to the fgas ∼ 0.5Ω b /Ωm observed by Mantz et al. (2014). Plausibly, thermal conduction might play a role here to distribute energy better towards the outskirts.

A note on fundamental differences between AMR and SPH
Various authors using SPH simulations have argued that the parameter ∆T has a strong impact on the physical properties of both the galaxies and the intracluster gas. As we have discussed above, this is in contrast to our own findings. Le Brun et al. (2014) find changes in the gas density, Pike et al. (2014) in the gas pressure out to R500. Furthermore, Le Brun et al. (2014) found that ∆T = 10 8 K allowed them to make cluster properties compatible with a range of observables. We discussed in more detail the impact of the two flavours of energy injection as well as the energy accumulation threshold on the various observables and various physical properties in Section 7.1 above. As we demonstrated there, in our Eulerian AMR simulations, the particular choice of ∆T and even the size of the injection region does not affect larger radii. It is thus plausible that thermal feedback has a different impact in Lagrangian and Eulerian simulations. Unfortunately, to the best of our knowledge, no direct comparison exists beyond the indirect estimate of Chaudhuri et al. (2013) who find that the amount of energy per unit mass that the AGN has to provide is larger for SPH than for AMR. While the lack of entropy cores in purely adiabatic SPH simulations has been known for a long time (see Frenk et al. 1999;Power et al. 2014, for the original result and a recent explanation of its origin), more modern formulations of SPH are able to alleviate this discrepancy with Eulerian methods (see e.g. Rasia et al. 2014;Sembolini et al. 2015, for recent method comparisons). At the same time, the effect of feedback models in controlled experiments is not well documented. We hope that a comparison of the numerical discrepancies, as suggested by our results, will be undertaken by the community in the near future.

A note on metal enrichment
We find (see paper 2, Martizzi et al. 2015, in prep. for more details) that both our intracluster medium as well as the galaxy population have a metallicity which is slightly lower than what is of observed in the gas (e.g. De Grandi et al. 2004;Matsushita 2011;Werner et al. 2013) and stellar metallicity (e.g. Gallazzi et al. 2005). A metallicity distribution more consistent with observations has been reproduced in various SPH simulations (e.g. Sijacki et al. 2007;Planelles et al. 2014). Also in the Illustris simulation (not SPH) realistic metallicities were produced (Vogelsberger et al. 2014), albeit at the price of baryon fractions that are completely inconsistent with observational constraints (Haider et al. 2015, their Figure 1). On the other hand, our results are consistent with the low metallicities found by other AMR simulations, e.g. Dubois et al. (2014). This aspect is another possible systematic discrepancy between the methods and should be investigated in more detail. We note that standard SPH does not include any mixing of metals (or any tracers, although explicit metal diffusion can be added, see e.g. Shen et al. 2010). Planelles et al. (2014) included a smoothed metallicity estimate when calculating cooling times and state that their unsmoothed metallicity field is very noisy. We investigate the aspect of metal enrichment in our simulations further in paper 2, Martizzi et al. 2015, in prep.

SUMMARY AND CONCLUSIONS
We present simulations of nine clusters of mass Mvir ∼ 10 15 M and one of twice that mass, including cooling, star formation, as well as supernova and AGN feedback with a physical resolution of 3.8 h −1 kpc. The simulations include the environment of the clusters inside of spheres of 8 h −1 Mpc around each cluster at z = 0. In this paper, we compare in detail the ICM density, temperature, entropy, and gas depletion profiles with X-ray data by performing a time-ensemble analysis for each cluster over a narrow mass bin and comparing with observed clusters in the same mass range. Next, we investigate the evolution of our simulated clusters over cosmic time with a range of cosmological observables that serve as mass proxies. In this paper, we focus on the mass vs Compton-y, X-ray luminosity and X-ray temperature scaling relations that are of particular importance for cluster cosmology. We also establish the numerical convergence of our results with resolution and their robustness against changes in the AGN feedback parameters.
We summarise our findings as follows: (i) We find a persistent cool-core/non-cool core dichotomy in our clusters. The cool cores are insensitive to changes in the thermal AGN feedback model parameters.
(ii) We link the disruption of cool cores to low-angular momentum major mergers. Major mergers with enough angular momentum leave the cool cores intact. Core disruption occurs on time scales of at most a few 100 Myr, with much increased core cooling times after the disruption, leading to a quick transition and thus a stable bimodality.
(iii) Our simulations agree with the Planck M500 − Y500 scaling relation with very little scatter. We do not identify a strong dependence of the scatter on the accretion history or the AGN model parameters.
(iv) The RG clusters are more X-ray luminous than a comparison sample from the MCXC catalog. The clusters evolve along scaling relations in the M − LX plane that are consistent with self-similar scaling. There is a slight indication that the scatter in this relation correlates with details of the assembly history.
(v) The non-cool core clusters reproduce density, entropy and mass profiles of an ACCEPT comparison sample well, and are consistent with the observed gas depletion profiles. The cool core systems have excess central gas and a too low central entropy compared to the ACCEPT clusters. In addition, there is a general indication that at large radii, the simulated clusters have a slightly too low entropy and temperature, and a slightly too high density compared to observations. (vi) The galaxies forming in our simulations have realistic masses and are consistent with abundance matching results across three decades in halo mass. At higher masses, the simulated galaxies are slightly more massive than observed at a given halo mass, although we caution that observational issues complicate a detailed comparison. The star formation rates at the high mass end are consistent with recent observational constraints for BCGs.
The discrepancies we observe and listed above are plausibly related to shortcomings of the simplistic central thermal blast wave model: (i) In our AMR simulations, we find that thermal AGN feedback does not affect the ICM at significantly large radii. We see no effect on gas at scales of ∼ Rvir/2, nor is the AGN able to mildly stabilise the cool core systems. Once a cool core forms, it cools below observed core entropies and leads to a core baryon fraction inconsistent with observational constraints.
(ii) The above finding is discrepant with SPH results in the literature. In our simulations, details of the AGN energy injection are irrelevant for global ICM properties. In particular, the X-ray and SZ scaling relations are unaffected by details of the thermal AGN model. This is quite in contrast to the findings of, e.g., Le Brun et al. (2014) and points to a possible discrepancy between SPH/Lagrangian and Eulerian methods and how feedback couples to gas in such simulations.
The inability of the thermal AGN model to shape larger scales in our simulations plausibly points to other forms of energy injection (e.g. through kinetic feedback, see Li et al. 2015;Choi et al. 2015) or additional processes in shaping the intracluster medium. Several published results suggest that thermal conduction might play a central role both in stabilising cool cores and at larger scales (e.g. Guo et al. 2008;Parrish et al. 2009;Ruszkowski & Oh 2010;Arth et al. 2014). We will investigate these aspects in future research.