COOL CORE CLUSTERS FROM COSMOLOGICAL SIMULATIONS

During their hierarchical assembly, galaxy clusters grow in mass via diffuse accretion processes as well as via merger events. These processes lead to the shock heating of the intra-cluster medium (ICM) to a virial temperature of up to T ∼ 10⁸ K, with central electron number densities corresponding to n_e ∼ 10⁻³ cm⁻³ (e.g., Voit 2005; Kravtsov & Borgani 2012, for a review). At this density, hot baryons in the core regions have a cooling time, t_cool ∝ T^1/2/n_e, shorter than the Hubble time. Under this condition, the ICM suffers from radiative losses and becomes colder and denser. This process might initiate a run-away cooling process, which would cause an extreme accretion rate in the central galaxy, unless some heating mechanism were to balance the cooling. Substantial evidences from X-ray observations proved that the regulating mechanism is the feedback by active galactic nuclei (AGNs; e.g., McNamara & Nulsen 2007; Voit et al. 2015).

In a situation of steady accretion, low-entropy⁸ and metal-enriched gas associated with merging substructures, or funneled by filaments, sink toward the central regions of the so-called cool-core clusters (CC; named as such by Molendi & Pizzolato 2001). On the contrary, non-cool-core (NCC) objects are observed with nearly isentropic gas cores at a higher entropy level and show no evidence of a spike in their metal abundances (De Grandi et al. 2004; Maughan et al. 2008).

While idealized simulations reproduce some of the many features of the CC systems (e.g., Gaspari 2015; Li et al. 2015), producing the diverse populations of CC and NCC clusters has so far proved to be quite a formidable challenge for cosmological hydrodynamic simulations (Borgani & Kravtsov 2011, and references therein). Earlier simulation works predicted that the origin of the two classes of CC and NCC clusters is "primeval," with CC being destroyed exclusively by early mergers or preheating phenomena (Burns et al. 2008; McCarthy et al. 2008; Poole et al. 2008; Planelles & Quilis 2009). In fact, a CC cluster was claimed to maintain its high-density and metal-rich core even after a major merger, if that happens at z ≤ 0.5. However, this prediction is in confict with observations and recent cosmological simulations by Hahn et al. (2015). In most cases, CC clusters show regular X-ray morphology while the opposite is true for NCC systems. Data, thus, favor an "evolutionary" model in which a transition between the CC and NCC status of a cluster can take place over relatively short timescales in consequence of a merger event that mixes the convective stable gas and destroys the CC (Rossetti et al. 2011). Specifically, Hahn et al. (2015) advocate that the dichotomy between CC and NCC clusters originates by low-angular-momentum major mergers.

In this Letter, we present a set of simulated clusters in which, for the first time, we find (1) co-existence of CC and NCC systems, (2) agreement with the observed thermodynamical and chemical properties of the two classes of clusters, and (3) evidence of a recent transition between CC and NCC in the evolution of individual objects. As we will discuss, this result is achieved thanks to the combined action of a new AGN feedback model (Steinborn et al. 2015) and an improved smoothed particle hydrodynamics (SPH) scheme (Beck et al. 2015).

We provide a short description of the simulations analyzed here, while we refer to a forthcoming paper (S. Planelles et al. 2015, in preparation) for a more thorough description. The simulations have been generated from the same set of initial conditions originally presented by Bonafede et al. (2011). They correspond to zoomed-in initial conditions of 29 Lagrangian regions selected around 24 massive clusters with M₅₀₀ between 5 and 20 × 10¹⁴ h⁻¹ M_⊙ and 5 poorer clusters with M₅₀₀ in the range 0.7–3 × 10¹⁴ h⁻¹ M_⊙. The simulations are carried out with the GADGET-3 code (Springel 2005) and include an upgraded version of the SPH scheme, as described in Beck et al. (2015). These developments consist of (a) a higher-order (Wendland C⁴ instead of the standard B-spline) interpolating kernel to better describe discontinuities and reduce clumpiness instability, (b) a time-dependent artificial viscosity term to capture shocks and minimize viscosity away from shock regions, (c) an artificial conduction term that largely improves the SPH capability of following gas-dynamical instabilities and mixing processes. With respect to the physical thermal conduction, it similarly promotes the transport of heat, since it is applied to contact discontinuities in internal energy, but, differently, it does not depend on thermodynamical properties of the gas. To assess the performances of this improved SPH implementation, several standard hydrodynamical tests, such as weak and strong shocks, shear flows, gas mixing, and self-gravitating gas, have been performed in Beck et al. (2015).

The physical processes included in our simulations can be described as follows. Metallicity-dependent radiative cooling and the effect of a uniform time-dependent UV background are considered as in Planelles et al. (2014; see also Wiersma et al. 2009). A sub-resolution model for star formation from a multi-phase interstellar medium is implemented as in Springel & Hernquist (2003). Kinetic feedback driven by a supernova (SN) is accounted in the form of galactic winds with a velocity of v_w ∼ 350 km s⁻¹. Metal production from SN-II, SN-Ia, and asymptotic-giant-branch stars follows the original receipt by Tornatore et al. (2007; see also Planelles et al. 2014). The AGN feedback is modeled with the implementations recently presented in Steinborn et al. (2015). The released thermal energy accounts for contributions by both mechanical outflows and radiation, separately computed in the code. Eddington-limited gas accretion onto BHs is computed by multiplying the Bondi rate by a boost factor α. In this Letter, we consider only cold accretion (α_cold = 100 and α_hot = 0). We verified that the results hold when hot accretion is included with α_hot = 10. For the feedback energy, we assume variable efficiencies for the radiation and the mechanical outflow (Equations (19) and (20) of Steinborn et al. 2015) that depend on the BH mass and on the accretion rate. Finally, we fix _f = 0.05 as the efficiency of coupling the radiated BH energy to the ICM.

We emphasize that the parameters defining the modification of SPH and the AGN feedback model have been tuned to overcome the limitations of standard SPH in passing classic hydrodynamical tests and to reproduce the observational scaling relation between BH mass and stellar mass of the host galaxies, including the brightest cluster galaxies (McConnell & Ma 2013). No attempt has been pursued to reproduce any of the observational properties of the ICM. The existing match between other simulated and observed stellar and hot gas properties will be presented in future papers.

The dark matter mass resolution is equal to 8.3 × 10⁸ h⁻¹ M_⊙ and the z = 0 Plummer-equivalent softening length for the gravitational force is = 3.75 h⁻¹ kpc. In our sample, this value is always inferior than 0.01 × R₅₀₀.

Several approaches have been used in the literature to classify clusters as CC and NCC systems based, e.g., on the central gradient of temperature or entropy, on the central entropy level, on the cooling time, on the mass deposition rate inferred from X-ray spectroscopy, on the cuspiness of the central gas density profile (comparisons among these criteria are presented in Cavagnolo et al. 2009; Hudson et al. 2010; Leccardi et al. 2010; McDonald et al. 2013; Pascut & Ponman 2015). In order to assign a degree of "cool-coreness" to our simulated clusters, we jointly use the following two criteria, which measure the shape and level of the entropy profiles in the central regions and are extensively applied to observational data.

Following Rossetti et al. (2011, see also Leccardi et al. 2010), we base the first criterion on the computation of the pseudo entropy:

$\begin{eqnarray}&&\sigma =\displaystyle \frac{({T}_{{\rm{IN}}}/{T}_{{\rm{OUT}}})}{{({\mathrm{EM}}_{{\rm{IN}}}/{\mathrm{EM}}_{{\rm{OUT}}})}^{1/3}},\end{eqnarray} \tag{ 1 }$

where the spectroscopic-like temperature (Mazzotta et al. 2004), T, and the emission measures, EM, are computed in the IN region, r < 0.05 × R₁₈₀, and in the OUT region, 0.05 × R₁₈₀ < r < 0.15 × R₁₈₀. We define CC those clusters with σ < 0.55.

In addition, we request the value of the central entropy to be K₀ < 60 keV cm². We model the entropy profile as

$\begin{eqnarray}&&K={K}_{0}+{K}_{100}\times {(R/{R}_{100k})}^{\alpha },\end{eqnarray} \tag{ 2 }$

where we treat as free parameters the core entropy, K₀, the entropy at R_100k(=100 kpc), K₁₀₀, and the slope, α, (Cavagnolo et al. 2009). For each simulated cluster, the entropy profile is measured in linearly equispaced spherical shells, with the innermost radius being that of the sphere containing at least 100 gas particles. In all cases, this radius is greater than 5 × . We use R₅₀₀ as the outermost radius to which we perform the fit. Three clusters with σ < 0.55 have steadily declining entropy profiles. In these cases, K₀ is not constrained. We, thus, consider the overall profiles and keep in the CC sample the two objects with α ∼ 1 and exclude that with α < 0.5.

The distribution of σ and K₀ is continuous and no clear bimodality is present, in contrast with recent results by Hahn et al. (2015). Considering the aforementioned observational thresholds, we find that 38% (11/29) of our simulated clusters at z = 0 are classified as CC. This fraction is quite close to the 35% value reported by Eckert et al. (2011) from a similar mass distributed sample, HIGFLUGS, after accounting for the sample incompleteness and the X-ray selection bias.

In the upper panels of Figure 1, we compare the entropy profiles of our simulated clusters with observational results obtained by Pratt et al. (2010) from XMM-Newton observations of the REXCESS clusters. Our simulations now produce a population of CC clusters, whose entropy profiles decline down to the innermost regions. This result confirms that the energy extracted from the central BHs is distributed in such a way to compensate for the radiative losses of gas at low entropy and, therefore, to keep it in the hot phase despite the fact that it formally had a short cooling time. More generally, we note that observed and simulated entropy profiles agree quite well in normalization and slope, for both CC and NCC populations.

Zoom In Zoom Out Reset image size

Dubois et al. (2011) were also able to produce the correct CC entropy structure in a high-resolution simulation of a Virgo-like cluster by including the effect of AGN feedback. However, in the same paper, these authors claimed that this conclusion is spoiled as soon as the effect of metallicity is included in the cooling function. This confirms the fragility of self-regulated AGN feedback. We, therefore, compare the simulated mass-weighted iron profiles to the metallicity of nearby clusters, z < 0.2, from Ettori et al. (2015),⁹ in the bottom panels of Figure 1. In keeping with observational results, simulated CC clusters are characterized by a rather high enrichment level in the central regions while NCC objects have a shallower peak. This agreement suggests that our simulations provide an effective description of the mixing related to gas-dynamical processes acting during the cluster formation and of the rise of enriched gas caused by feedback processes: the same processes leading to the diversity of entropy profiles are also responsible for the diversity of the profiles of metal abundance.

From an object-by-object investigation of our simulated clusters, we recognize a variety of situations for the establishment of a CC and for the transformation between CC and NCC, which can take place over quite a short time and even at low redshift. As an example, in Figure 2, we show the evolutionary sequences of pseudo-entropy maps from z = 0.8 to z = 0 for two objects. Following an observational approach, we create maps of spectroscopic-like temperature, ${{\mathcal{M}}}_{{\rm{TSL}}},$ and of X-ray soft band ([0.5–2] keV) emissivity, ${{\mathcal{M}}}_{{\rm{X}}},$ which we then combine as ${{\mathcal{M}}}_{{\rm{TSL}}}/{{\mathcal{M}}}_{{\rm{X}}}^{1/3}$ to build maps of pseudo entropy (Finoguenov et al. 2010). The upper panel shows the case of an NCC cluster at z = 0.8 that slowly and continuously develops a CC structure. On the contrary, the lower panel corresponds to a system with a low-entropy core at high redshift. The core is still in place at z = 0.25 and it is eventually destroyed by a merger between the main halo and a sub-clump at z ∼ 0.1. Remarkably, this last case of CC remnant is not an isolated event: almost half of NCC clusters at z = 0 were CC at z = 1 (Rossetti et al. 2011).

Zoom In Zoom Out Reset image size

These results vary somewhat from previous claims, also based on simulations, which stated that CCs can only be destroyed at high redshift through major mergers, while being much more resilient against mergers at low redshift (e.g., Burns et al. 2008; Poole et al. 2008). A close comparison of our results with these previous analyses is not straightforward due to substantial differences in both the hydrodynamic schemes and the implementation of feedback. Poole et al. (2008) carried out non-cosmological simulations of isolated two-body mergers using a more standard implementation of SPH, while Burns et al. (2008) used a Eulerian code with a cosmological set-up. Neither one of these two analyses included the effect of AGN feedback and, for this reason, the strength of their cores is probably affected by overcooling. In order to surpass this limitation, Burns et al. (2008) truncated star formation below a certain redshift, while Poole et al. (2008) excluded the central 40 kpc in their analysis. Both strategies, however, do not directly attack the main problem of the delicate balance between cooling and heating. This balance, essential for the creation of a low-entropy core, needs to be regulated since early times and can be achieved in our simulations only by including AGN feedback.

It has long been known that ICM evolution purely driven by gravitational processes leads to self-similar entropy profiles that scale as a power-law of radius, K(r) ∝ r^α with α ∼ 1–1.1 in the cluster outskirts (Voit et al. 2005; Nagai et al. 2007). Non-radiative simulations, however, differ in their behavior in the central regions. A flat core profile is found in grid-based codes (Frenk et al. 1999) and in recent improved versions of SPH (Biffi & Valdarnini 2015; Sembolini et al. 2015, and references therein). Standard SPH simulations, instead, produce entropy profiles that are declining toward the center.

In addition to mixing, radiative processes have the counterintuitive effect of also raising entropy in core regions (e.g., Borgani & Kravtsov 2011). If not counteracted by some energy feedback process, radiative cooling selectively removes gas with low entropy and short cooling time from the hot-X-ray emitting phase, thus fuelling star formation (e.g., Li & Bryan 2012) and leaving in the ICM only gas with a relatively higher entropy.

The AGN feedback plays a key role in explaining our results. Considering the new SPH set-up but excluding AGN feedback, we find that overcooling leads to a thermal structure of central cluster regions whose entropy profiles behave as in NCC clusters (top panel of Figure 3). At the same time, the excess of star formation produces an exceedingly high level of ICM metal enrichment, with heavy elements mostly locked around galaxies. As a consequence, simulated clusters without AGNs always display spikes of central metallicity, an enrichment pattern typical of CC clusters (bottom panel of Figure 3). As observations suggest, the AGN feedback should have the threefold effect of (1) compensating radiative losses so as to prevent gas with a short cooling time from dropping out of the ICM, (2) reducing star formation in massive halos at low redshift, and (3) raising metals away from star-forming regions (e.g., Sijacki et al. 2007; Fabjan et al. 2011; McCarthy et al. 2011).

Zoom In Zoom Out Reset image size

However, attempts to create realistic thermal structures of CC in simulated clusters so far have failed, even when including AGN feedback, which is otherwise successful in reproducing X-ray scaling relations (e.g., Puchwein et al. 2008; Planelles et al. 2014) and their scatter (Le Brun et al. 2014). The encouraging success of the simulations presented here in producing thermal and chemo-dynamical structures of CC in agreement with observations is related to the combined action of the AGN feedback model and of the artificial thermal diffusion.

Figure 4 compares the entropy profiles of the cluster shown in the upper panels of Figure 2 with those of the same object simulated without artificial conduction but includes other recent implementations such as the new kernel and the artificial viscosity. By excluding artificial conduction, the entropy profiles at z = 0.8 and z = 0 are similar to each other with a large isentropic core, characteristics of a strongly NCC object. This contrast with the previous finding showing a transition from an NCC to a CC situation (Figure 2). Besides confirming the results reported in previous works based on standard SPH (Sijacki et al. 2008; Planelles et al. 2014), this outcome highlights the role of artificial conduction in providing an effective description of the processes of mixing that lead to the redistribution of AGN feedback energy in central regions.

Zoom In Zoom Out Reset image size

In summary, we showed that our simulations produce a population of clusters that have a CC reasonable structure, both from a thermo- and chemo-dynamical point of view. Furthermore, once this result is achieved, the simulations also naturally produce a mix of CC and NCC clusters, which are similar to that observed. We trace the reason for the success of our simulations to the combined action of artificial conduction, introduced to overcome numerical limitations of standard SPH, and of AGN feedback, introduced to regulate star formation in massive halos.

As a concluding word of caution, we would like to emphasize that the success of the simulations presented here does not imply that we are self-consistently describing the details of all of the physical processes that lead to the creation of cool cores: inflation of bubbles of high-entropy gas from the shocks of sub-relativistic jets launched by the AGN, gas circulation and turbulence triggered by the buoyancy of these bubbles, possible effects of magnetic fields in stabilizing them, thermal conduction, and cosmic rays. These are all effects for which we have circumstantial evidence from observations, but that are not included in our simulations. Still, our results demonstrate that our description of AGN feedback and artificial conduction provides a realistic effective description of the above physical processes, without requiring any specific fine tuning.

We are greatly indebted to Volker Springel for the access to the GADGET3 code; to D. Fabjan, V. Fiorenzo, M. Petkova, and L. Tornatore for the simulation set-up; and to the referee, D. Eckert, S. Ettori, A. Evrard, M. Gaspari, S. Molendi, P. Monaco, P. Tozzi, M. Voit for useful discussions. We acknowledge financial support from PIIF-GA- 2013-627474, NSF AST-1210973, PRIN-MIUR 201278X4FL, PRIN-INAF 2012 "The Universe in a Box: Multi-scale Simulations of Cosmic Structures," the INFN INDARK grant, "Consorzio per la Fisica" of Trieste, DFC Cluster of Excellence "Universe," DFC Research Unit 1254, CONICET-Argentina, FonCyT. Simulations are carried out using Flux HCP Cluster at the University of Michigan, Galileo at CINECA (Italy), with CPU time assigned through ISCRA proposals and an agreement with the University of Trieste, and PICO at CINECA though our expression of interest.