Heating and enriching the intracluster medium

Heating and enriching the intracluster medium

C. J. Short, P. A. Thomas and O. E. Young
Astronomy Centre, University of Sussex, Falmer, Brighton, BN1 9QH, United Kingdom
E-mail: P.A.Thomas@sussex.ac.uk
Abstract

We present numerical simulations of galaxy clusters with stochastic heating from active galactic nuclei (AGN) that are able to reproduce the observed entropy and temperature profiles of non-cool-core (NCC) clusters. Our study uses -body hydrodynamical simulations to investigate how star formation, metal production, black hole accretion, and the associated feedback from supernovae and AGN, heat and enrich diffuse gas in galaxy clusters. We assess how different implementations of these processes affect the thermal and chemical properties of the intracluster medium (ICM), using high-quality X-ray observations of local clusters to constrain our models. For the purposes of this study we have resimulated a sample of massive galaxy clusters extracted from the Millennium Simulation. Sub-grid physics is handled using a semi-analytic model of galaxy formation, thus guaranteeing that the source of feedback in our simulations is a population of galaxies with realistic properties. We find that supernova feedback has no effect on the entropy and metallicity structure of the ICM, regardless of the method used to inject energy and metals into the diffuse gas. By including AGN feedback, we are able to explain the observed entropy and metallicity profiles of clusters, as well as the X-ray luminosity-temperature scaling relation for NCC systems. A stochastic model of AGN energy injection motivated by anisotropic jet heating – presented for the first time here – is crucial for this success.

With the addition of metal-dependent radiative cooling, our model is also able to produce CC clusters, without over-cooling of gas in dense, central regions.

keywords:
hydrodynamics – methods: N-body simulations – galaxies: clusters: general – galaxies: cooling flows – X-rays: galaxies: clusters.

1 Introduction

1.1 Background

Clusters of galaxies are believed to be the largest gravitationally-bound objects in the Universe. Their deep gravitational potential well means that the largest clusters are ’closed boxes’, in the sense that baryons ejected from cluster galaxies by supernova (SN) explosions and active galactic nuclei (AGN) do not escape the cluster completely, but instead end up in the hot, diffuse plasma that fills the space between cluster galaxies – the intracluster medium (ICM).

The thermal properties of intracluster gas, which can be measured with X-ray telescopes such as Chandra, XMM-Newton and Suzaku, thus provide a unique fossil record of the physical processes important in galaxy and galaxy cluster formation and evolution, such as radiative cooling, star formation, black hole accretion and the subsequent feedback from supernovae (SNe) and AGN. In addition, measurements of the ICM chemical abundances yield information about the production of heavy elements in stars in member galaxies, providing constraints on nucleosynthesis, and the processes responsible for their transport into the ICM.

A key diagnostic of the thermal state of intracluster gas is provided by the gas entropy111We define the gas entropy as , where is Boltzmann’s constant, is the gas temperature, is the electron number density and is the ratio of specific heats for a monoatomic ideal gas.. Entropy remains unchanged under adiabatic processes, such as gravitational compression, but increases when heat energy is introduced and decreases when radiative cooling carries heat energy away, thus providing an indicator of the non-gravitational processes important in cluster formation.

In recent years, spatially-resolved observations have facilitated a detailed examination of the radial distribution of entropy in clusters. Observed entropy profiles are typically found to scale as at large cluster-centric radii, 222We define as the radius of a spherical volume within which the mean matter density is times the critical density at the redshift of interest. The mass enclosed within this sphere is denoted by . (e.g. Ponman et al. 2003; Sun et al. 2009; Cavagnolo et al. 2009; Sanderson et al. 2009; Pratt et al. 2010). This power-law scaling agrees with that predicted by simple analytical models based on spherical collapse (Tozzi & Norman, 2001) and cosmological hydrodynamical simulations that include gravitational shock heating only (e.g. Voit et al. 2005, hereafter STY10; Nagai et al. 2007, hereafter STY10; Short et al. 2010, hereafter STY10).

However, in the inner regions of clusters, observations have unveiled the presence of a mass-dependent entropy excess with respect to theoretical expectations, and a large dispersion in central entropy values (e.g. Ponman et al. 1999; Lloyd-Davies et al. 2000; Ponman et al. 2003; Pratt et al. 2006; Morandi & Ettori 2007; Cavagnolo et al. 2009; Pratt et al. 2010). The source of this entropy excess is likely to be a combination of non-gravitational heating from astrophysical sources, such as SNe and AGN, and cooling processes. There is evidence that the distribution of central entropy values is bimodal, with morphologically disturbed non-cool-core (NCC) systems having an elevated central entropy compared to dynamically relaxed cool-core (CC) systems (e.g.  Cavagnolo et al., 2009). It is thought that the association of unrelaxed morphology with a high central entropy is an indication that either cool cores are destroyed by mergers, or that cool cores have never been able to form in these objects.

The chemical state of the ICM is characterised by its metallicity, the proportion of chemical elements present heavier than H and He. X-ray spectroscopy of galaxy clusters has revealed emission features from a variety of chemical elements, O, Ne, Mg, Si, S, Ar, Ca, Fe and Ni, all of which are synthesised in stars and transported to the ICM by processes such as ram-pressure stripping, galactic winds and AGN outflows. Type Ia SNe produce a large amount of Fe, Ni, Si, S, Ar and Ca, but compared to Type II SNe, they only produce very small amounts of O, Ne and Mg. Type Ia SN products are found to dominate in cluster cores, whereas Type II SN products are more evenly distributed (Finoguenov et al. 2000; Tamura et al. 2001; see Böhringer & Werner 2010 for a recent review). This can be explained by early homogeneous enrichment by Type II SNe, which produce -elements in the protocluster phase, and a subsequent, more centrally-peaked enrichment by Type Ia SNe, which have longer delay times and continue to explode in the central cD galaxy long after the cluster is formed.

The first detailed measurements of spatial abundance distributions were made by De Grandi & Molendi (2001), using data from BeppoSAX, who measured the radial Fe abundance profiles for a sample of massive clusters. They found that CC clusters have a sharp Fe abundance peak in central regions, whereas NCC clusters have flat Fe abundance profiles. Subsequent observations with XMM-Newton and Chandra have confirmed this dichotomy between the metallicity distribution in CC and NCC clusters (Tamura et al., 2004; Vikhlinin et al., 2005; Pratt et al., 2007; Baldi et al., 2007; Maughan et al., 2008; Leccardi & Molendi, 2008; Matsushita, 2011).

It is thought that cluster mergers, and the subsequent mixing of intracluster gas, are responsible for destroying the central abundance peak found in CC clusters. However, some systems with a highly disturbed morphology are also found to have a high central metallicity. Leccardi et al. (2010, see also ) suggest that these objects correspond to relaxed CC systems that have undergone a major merger, or a significant AGN heating event, very recently, so that mixing processes have not yet had sufficient time to fully erase low-entropy gas and the central abundance peak.

Explaining the observed thermal and chemical properties of the ICM from a theoretical perspective requires a detailed understanding of the complex interplay between large-scale gravitational dynamics and the various small-scale astrophysical processes mentioned above. Numerical cosmological hydrodynamical simulations have emerged as the primary tool with which to tackle this problem. There has been considerable effort to include the processes relevant for cluster formation and evolution in simulations in a self-consistent manner; see Borgani & Kravtsov (2009) for a recent review. However, an explicit treatment is unfeasible since these processes all occur on scales much smaller than can be resolved with present computational resources.

Hydrodynamical simulations that include models of radiative cooling, star formation, metal production and galactic winds generally fail to reproduce observed ICM temperature, entropy and metallicity profiles. Simulated temperature profiles typically have a sharp spike at small cluster-centric radii, followed by a rapid drop in temperature moving further into the core (e.g. Valdarnini 2003; Tornatore et al. 2003; Borgani et al. 2004; Romeo et al. 2006; Sijacki et al. 2007; Nagai et al. 2007), in clear conflict with the smoothly declining (flat) profiles of observed CC (NCC) clusters (e.g. Sanderson et al. 2006; Vikhlinin et al. 2006; Arnaud et al. 2010). This is due to the adiabatic compression of gas flowing in from cluster outskirts to maintain pressure support, following too much gas cooling out of the hot phase. This over-cooling causes excessive star formation in cluster cores, with predicted stellar fractions being about a factor of larger than the observed value of (Balogh et al., 2001; Lin et al., 2003; Balogh et al., 2008), which in turn leads to excessive Fe production in central regions, generating steeper abundance profiles than observed (e.g. Valdarnini 2003; Tornatore et al. 2004; Romeo et al. 2006; Tornatore et al. 2007; Davé et al. 2008).

It is generally accepted that the solution to the over-cooling problem in hydrodynamical simulations is extra heat input from AGN. Simple analytical arguments convincingly show that the energy liberated by accretion onto a central super-massive black hole is sufficient to suppress gas cooling and thus quench star formation. The precise details of how this energy is transferred to the ICM are not well understood at present, but it appears that there are two major channels via which black holes interact with their surroundings (see McNamara & Nulsen 2007 for a review).

At high redshift, mergers of gas-rich galaxies occur frequently and are expected to funnel copious amounts of cold gas towards galactic centres, leading to high black hole accretion rates and radiating enough energy to support the luminosities of powerful quasars. Quasar-induced outflows have been observationally confirmed in a number of cases (e.g. Chartas et al. 2003; Crenshaw et al. 2003; Pounds et al. 2003; Ganguly & Brotherton 2008; Dunn et al. 2010).

Evidence for another mode of AGN feedback, not related to quasar activity, can be seen in nearby CC clusters, which often contain radio-loud X-ray cavities in the ICM. These bubbles are thought to be inflated by relativistic jets launched from the central super-massive black hole (Blanton et al., 2001; Bîrzan et al., 2004; McNamara et al., 2005; Fabian et al., 2006; Morita et al., 2006; Jetha et al., 2008; Gastaldello et al., 2009; Dong et al., 2010; Giacintucci et al., 2011). Bubbles may rise buoyantly, removing some of the central cool, enriched gas and allowing it to mix with hotter gas in the outer regions of groups and clusters. Together with the accompanying mechanical heating, this can constitute an efficient mechanism for suppressing cooling flows, and redistributing metals throughout the ICM. Such flows are seen in simulations of idealised clusters, performed with hydrodynamical mesh codes (e.g. Churazov et al. 2001; Quilis et al. 2001; Ruszkowski & Begelman 2002; Brüggen et al. 2002; Brüggen 2003; Dalla Vecchia et al. 2004; Roediger et al. 2007; Brüggen & Scannapieco 2009).

Various authors have implemented self-consistent models of black hole growth and AGN feedback in cosmological simulations of galaxy groups and clusters (in addition to cooling, star formation, and thermal and chemical feedback from SNe). Springel et al. (2005a) developed a model for quasar mode AGN feedback (see also Di Matteo et al. 2005), which was used in cosmological simulations of galaxy groups by Bhattacharya et al. (2008). A model for radio mode AGN feedback based on bubble injection was proposed by Sijacki & Springel (2006), which was subsequently extended by Sijacki et al. (2007) to include quasar mode AGN feedback as well. Both Sijacki & Springel (2006) and Sijacki et al. (2007) performed cosmological simulations of a few massive clusters with their respective models.

These studies demonstrated, in a qualitative manner, that AGN feedback is effective in reducing the amount of cold baryons and star formation in the central regions of groups and clusters. Furthermore, the gas density is reduced and the temperature is increased, elevating the central entropy. Sijacki et al. (2007) also showed that AGN outflows drive metals from dense, star-forming regions to large radii, flattening ICM abundance profiles relative to those predicted by a run without AGN feedback. Such trends are precisely what is required to reconcile simulations of galaxy clusters with observations.

A more quantitative assessment of the impact of AGN feedback on the ICM was conducted by Puchwein et al. (2008). They resimulated a sample of groups and clusters with the scheme of Sijacki et al. (2007), finding that the model could reproduce the observed X-ray luminosity-temperature scaling relation, at least on average. However, since their sample size is quite small, it is unclear whether the model can generate a realistic population of CC and NCC systems and thus explain the observed scatter about the mean relation. In addition, the stellar fraction within the virial radii of their simulated objects appears larger than observed.

Another detailed study was undertaken by Fabjan et al. (2010), who resimulated a sample of groups and clusters in a cosmological setting, using a model closely related to that of Sijacki et al. (2007), but with a different implementation of radio mode AGN feedback. On group scales, they found that AGN heating was able to successfully balance radiative cooling, reproducing observed stellar fractions, but the central entropy (at ) was about a factor of too high. In addition, their predicted group Fe abundance profiles are flat for , whereas observed profiles have a negative gradient out to the largest radii for which measurements are possible (e.g. Rasmussen & Ponman 2009). There is also an indication that the Fe distribution may be too sharply peaked in central regions compared to observations. The effect of AGN feedback on galaxy groups was also investigated by McCarthy et al. (2010), who implemented the AGN feedback scheme of Booth & Schaye (2009) in a cosmological simulation. With this model they were able to explain the observed entropy, temperature and Fe abundance profiles of groups, as well as observed X-ray scaling relations.

For massive clusters, Fabjan et al. (2010) showed that their model can reproduce the entropy structure of the ICM, but a factor of too many stars were formed. The cluster Fe abundance profiles they obtained have a shape consistent with that of observed profiles, although with a higher normalisation, but the central Fe abundance may be over-estimated. Dubois et al. (2011) also examined the role of AGN feedback in establishing the properties of the ICM, using a cosmological AMR simulation of a massive cluster with a prescription for jet heating by AGN. The entropy profile of their cluster agrees well with that of observed CC clusters if metal-cooling is neglected, and when metals are allowed to contribute to the radiative cooling, the resulting profile resembles that of a NCC cluster instead. However, the metallicity profile of their cluster appears steeper than observed.

1.2 This work

In this work, we pursue a different, but complementary, approach to the theoretical study of galaxy clusters. Instead of undertaking self-consistent hydrodynamical simulations, we adopt the hybrid approach of Short & Thomas (2009, hereafter SHT09) which couples a semi-analytic model (SA model) of galaxy formation to a cosmological -body/smoothed-particle hydrodynamics (SPH) simulation. In this model, the energy imparted to the ICM by SNe and AGN is computed from a SA model and injected into the baryonic component of a non-radiative hydrodynamical simulation; see SHT09 for details. The main advantage of this approach is that feedback is guaranteed to originate from a realistic population of galaxies, since SA models are tuned to reproduce the properties of observed galaxies. As a consequence, the stellar fraction in massive clusters agrees with observations (Young et al., 2011), which is not the case in self-consistent hydrodynamical simulations.

We have extended the model of SHT09 to follow the metal enrichment of the ICM. Note that Cora (2006, see also ) have already used a similar hybrid technique to study the pollution of intracluster gas by heavy elements. However, they did not include energy injection from SNe and AGN, which are likely to affect the distribution of metals in the ICM.

In the model of SHT09, the energy liberated by SN explosions and black hole accretion is assumed to be distributed uniformly throughout the diffuse gas of the host halo. With this rather ad hoc heating model they were able to reproduce observed X-ray scaling relations for NCC clusters, but ICM entropy profiles were found to be flatter than observed within 0.5 times (STY10). These simulations do not well resolve the core (), nor do they include radiative cooling that is likely to be important in this region, at least for CC clusters. However, we would expect that they should be able to provide a much better fit to X-ray observations of NCC clusters outside the core.

The primary goal of this paper is, therefore, to formulate a new feedback model that has a clear physical motivation and that is better able to explain the radial variaton of both the thermal and chemical properties of intracluster gas outside the core of the cluster. To help us do this we test a wide variety of different models for SN and AGN feedback and metal enrichment, using a selection of X-ray data (namely, entropy and metallicity profiles and the luminosity-temperature scaling relation) to identify the features that a model should possess in order to reproduce the data.

Our conclusion is that a stochastic heating model, motivated by observations of anisotropic AGN outflows, provides a better fit to the observed properties of the ICM than more commonplace models, such as heating a fixed number of neighbours or heating particles by a fixed temperature. Using entirely plausible duty cycles and opening angles for the jets, it is possible to provide an acceptable fit to all available observations with our model.

Note that the use of SA models means that the feedback is not directly coupled to the cooling of the gas – that is why our previous work and the bulk of this paper uses non-radiative simulations and restricts its attention to NCC clusters. However, towards the end of the paper we introduce radiative cooling in an attempt to reproduced CC clusters. We estimate the degree to which the SA model fails to supply the required feedback energy and show that there can be a substantial short-fall at high redshift, but that it averages to under 10 per cent over the lifetime of the cluster. We are able to qualitatively reproduce some CC profiles, but we do not provide a detailed quantitative analysis here.

In this work, we neglect many physical effects such as magnetic fields, cosmic rays, thermal conduction, turbulent mixing, etc.. Our principal reason for doing this is to keep the model simple and ease interpretation of our results. Some of these may be important in the central regions of CC clusters () but there is little evidence that they play a significant role at the larger radii that we use to constrain our models. We discuss this further at the end of the paper.

The layout of this paper is as follows. In Section 2, we present the details of our hybrid numerical model and describe our cluster simulations. We investigate the effect of SN feedback on the thermal and chemical properties of the ICM in Section 3, and assess how our results are affected by different choices of SN feedback and metal enrichment models. We show, in agreement with previous work, that SN have little impact on the entropy structure of the intracluster gas. In Section 4 we examine the impact of additional heating from AGN: these can reproduce the correct scaling relations but give entropy profiles that are too flat. Our results motivate a new, stochastic feedback model based on jet heating, which is described in Section 5. In this section, we also discuss what this model predicts for the thermal and chemical properties of the ICM, and we conduct an exhaustive comparison with observational data in Section 6. In Section 7 we demonstrate that our model is capable of producing both CC and NCC clusters with the inclusion of metal-dependent radiative cooling. Our conclusions are presented in Section 8.

For those readers who are mostly interested in the final model itself, rather than the steps used to motivate it, we recommend skipping Sections 3 and 4, at least on first reading.

2 Simulations

We make use of hydrodynamical resimulations of a sample of massive galaxy clusters extracted from the dark-matter-only Millennium Simulation (Springel et al., 2005b). Our sample consists of objects with and forms a subset of the larger sample of groups and clusters resimulated by STY10 for their so-called FO simulation, one of the Millennium Gas Simulations333The Millennium Gas Simulations are a series of hydrodynamical simulations designed to add gas to the dark matter structures found in the Millennium Simulation (Springel et al., 2005b). At present, there are simulations, each of which employs a different physical mechanism for raising the entropy of intracluster gas. The first of these is a reference model that includes gravitational heating only (the GO run). The second includes radiative cooling and uniform preheating at as a simple model for heating from astrophysical sources (the PC run). The third simulation is the FO run, where feedback from galaxies is computed from a SA model using the hybrid model of SHT09.. See STY10 for details of the cluster selection procedure. Basic properties of our clusters are listed in Table 1.

Cluster name
C1
C2
C3
C4
C5
C6
C7
C8
C9
C10
C11
C12
C13
C14
C15
C16
C17
C18
C19
C20
C21
C22
C23
C24
C25

We define dynamical temperature as , where  kg is the mean particle mass and is the mean square velocity.

Table 1: The masses, (in units of ), and dynamical temperatures, (in units of keV), of the clusters used in this study within (second and third columns, respectively), and (third and fourth columns, respectively). Cluster C1 is our fiducial cluster, used for most of the plots in this paper.

Following STY10, the feedback model we adopt in our simulations is the hybrid scheme of SHT09, where a SA model of galaxy formation is used to compute the number of stars formed and the energy transferred to the ICM by SNe and AGN. We refer the reader to STY10 for a full description of the modelling process and simulation parameters.

Briefly, we first perform dark-matter-only simulations of each region containing a cluster in our sample using the massively parallel TreePM -body/SPH code GADGET-2 (Springel, 2005). Virialised dark matter haloes are identified at each simulation output using the friends-of-friends (FOF) algorithm, with a standard linking length of of the mean inter-particle separation (Davis et al., 1985). Only groups with at least particles are kept, yielding a minimum halo mass of . Gravitationally bound substructures orbiting within these FOF haloes are then found with a parallel version of the SUBFIND algorithm (Springel et al., 2001). From the stored subhalo catalogues we construct dark matter halo merger trees by exploiting the fact that each halo will have a unique descendant in a hierarchical scenario of structure formation; see Springel et al. (2005b) for further details.

The second stage is to generate galaxy catalogues for each resimulated region by applying the Munich L-Galaxies SA model of De Lucia & Blaizot (2007) to the halo merger trees. A full description of the physical processes incorporated in L-Galaxies and model parameters is given in Croton et al. (2006) and De Lucia & Blaizot (2007). For each galaxy in these catalogues, we use its merger tree to compute the change in stellar mass, , and mass accreted by the central black hole, , between successive model outputs. Knowledge of enables us to incorporate star formation in our simulations as described below. From and we can also calculate the energy imparted to intracluster gas by Type II SNe, , and AGN, , respectively. Details are given in SHT09.

For the purpose of this work, we have extended the model of SHT09 to also follow the enrichment of the ICM by metals ejected from galaxies in winds. In the L-Galaxies SA model  of the mass of newly formed stars is instantaneously returned, and deposited in the cold gas disc of the host galaxy (Croton et al., 2006). In other words, the model assumes that metal ejection is instantaneous and does not distinguish between emission from Type II and prompt Type Ia SNe, and that from delayed Type Ia SNe and AGB stars. This will be added in future work.

In each model galaxy, metals can reside in several distinct phases: stars, cold disc gas, hot halo gas and gas ejected by winds from the halo into an external ‘reservoir’. Only the latter two are relevant for the ICM. We define the total mass in metals in diffuse gas to be

(1)

where and are the mass in metals in hot and ejected gas, respectively.

Once a galaxy falls into a FOF group, becoming a satellite of the central galaxy of the halo, all of its metals in hot and ejected gas are assumed to be associated with the central galaxy. It follows that is non-zero only for central galaxies. Given a particular halo at some output redshift , we compute the change in metal content of the ICM, , since the previous output, , by taking for the central galaxy and subtracting the sum of for every galaxy that is a progenitor of any galaxy contained in the host FOF group and also a central galaxy of a halo at :

(2)

The quantity is used to implement metal enrichment of the ICM in our simulations, as described in subsequent sections.

Finally, we couple the L-Galaxies SA model to hydrodynamical simulations of our clusters to track the effect of feedback from galaxies on the thermal and chemical properties of the ICM. The initial conditions for these resimulations are the same as for the dark-matter-only runs described above, except that we add gas particles with zero gravitational mass. This ensures that the dark matter distribution remains undisturbed by the inclusion of baryons, so that the halo merger trees used to generate the semi-analytic galaxy catalogues will be the same. Gas particles are added at a lower resolution than the dark matter, simply to ease the computational cost of our simulations. The resolution we have adopted is sufficient to obtain numerically converged estimates of bulk cluster properties for systems with keV (SHT09).

Every time an output redshift is reached in our hydrodynamical simulations, temporary ‘galaxy’ particles are introduced at positions specified by the SA model galaxy catalogue. For each galaxy, we know the increase in stellar mass since the last output, and we remove this mass from the hot phase by converting the nearest gas particles into collisionless star particles, using a stochastic method to ensure that is an integer. Once star formation is complete, we then distribute metals and the heat energy available from SNe and AGN amongst neighbouring gas particles in some way, as described in the following sections. Following the injection of metals and entropy, the galaxy particles are removed and the simulation continues until the next output time, when the process is repeated. The main purpose of this paper is to investigate different ways of heating and enriching intracluster gas, using X-ray observations of galaxy clusters to constrain our models.

Cluster catalogues are generated at from our simulations using a procedure similar to that of Muanwong et al. (2002). Full details of our cluster extraction method are given in STY10.

Following SHT09, we choose to neglect gas cooling processes in our hydrodynamical simulations throughout most of this work. Although cooling is relatively unimportant for the majority of the ICM, we cannot expect to reproduce the low central entropy and steep entropy profiles of observed CC clusters, as demonstrated by STY10. However, in Section 7 of this paper, we make a first attempt to overcome this limitation of the model of SHT09 by including metal-dependent radiative cooling in our simulations. With the addition of cooling, we show that it is indeed possible to produce both CC and NCC systems using our hybrid approach.

3 Feedback from Type II supernovae

In this section we investigate how galactic winds driven by Type II SNe shape the chemical and thermal properties of intracluster gas. We pay particular attention to how our results are affected by varying the feedback scheme in our simulations. Sections 3.13.3 describe the models; then in Sections 3.4 and 3.5 respectively, we show that SNe simply do not provide enough energy to significantly alter the entropy and metallicity profiles of the ICM. Most of the metals in the ICM originate outside the central cluster galaxy and we argue that the metallicity profile in these models is imposed by the accretion history of the ICM.

3.1 Supernova feedback models

There are two broad classes of SN feedback models deployed in numerical simulations: thermal, where the available energy is used to raise the temperature of neighbouring gas particles (Katz, 1992; Mori et al., 1997; Thacker & Couchman, 2000; Kay et al., 2002; Brook et al., 2004; Stinson et al., 2006), and kinetic, where neighbouring particles are given a velocity ‘kick’ (Navarro & White, 1993; Mihos & Hernquist, 1994; Kawata, 2001; Kay et al., 2002; Springel & Hernquist, 2003; Oppenheimer & Davé, 2006; Dubois & Teyssier, 2008).

It is well known that simple thermal feedback schemes fail in simulations with cooling since the injected energy is radiated away before it has any hydrodynamical effect. This problem is typically evaded by suppressing radiative cooling by hand. However, this is not an issue for us since cooling processes are not included in any of our simulations until Section 7. We have experimented with a variety of both thermal and kinetic models, which we now describe.

3.1.1 Thermal models

The thermal feedback models employed in our simulations can be grouped into three categories, depending on the method used to inject the SN energy, , into the ICM.

In our first scheme, we simply heat a fixed number of the gas particles closest to each galaxy, where the number of neighbours heated is , or . This is the approach typically adopted in fully self-consistent hydrodynamical simulations with radiative cooling, star formation and thermal SN feedback.

The second method we have investigated is to heat all gas particles within some sphere centred on each galaxy, where the radius of the sphere is assumed to be some fraction, , of the halo virial radius, . We have explored , and . In this model, is defined as the number of gas particles enclosed by the sphere. If no neighbours are found, the radius of the sphere is increased until a single gas particle is found.

Our third approach is to heat neighbouring gas particles by a multiple, , of the halo virial temperature, , defined by

(3)

The values of that we adopt are , and . The number of gas particles that can be heated with the available energy is then

(4)

In GADGET-2 the thermodynamic state of each fluid element is defined in terms of the entropic function

(5)

where is the thermal energy per unit mass of a particle and is its density. Supplying heat energy to a gas particle causes to increase. Note that is related to the X-ray gas entropy via , where  kg is the mean molecular mass per free electron.

In each of our three feedback schemes, we heat particles by raising their thermal energy by a fixed amount

(6)

implemented in GADGET-2 as an entropy boost of

(7)

The product of the cosmic baryon fraction, , and the virial density, , gives the mean overdensity of baryons within the virial radius. If the required neighbours are not found within a distance of a galaxy and the search radius has to be increased, the density of some of these particles may be less than . By using , rather than , in the denominator of equation (7), we are assuming that the amount of energy used to heat such particles is ; the rest of the energy is taken to be used up as the gas does work expanding adiabatically to a density .

For the first two schemes mentioned above, we have also tested an alternative heating model where gas particles are given a fixed entropy, rather than energy, boost:

(8)

Denser particles close to a galaxy are then heated to a higher temperature than more distant, lower density particles.

3.1.2 Kinetic models

Kinetic SN feedback is implemented in our simulations by assuming that gas particles closest to a model galaxy are given a velocity kick. The number of particles that receive a kick depends on the available energy:

(9)

where the wind speed, , is a free parameter. We have considered several different values for the wind speed: km s, km s, km s and km s. For any given galaxy, we impose the constraint that the wind speed cannot be less than the virial speed444We define the virial speed of a halo to be the circular velocity at the virial radius. of the host halo, , so that the case km s is equivalent to assuming that material is ejected at the virial speed. To ensure that is an integer, we draw a random number uniformly from the unit interval and compare it with the fractional part of : if is less (greater) than the fractional part of , we round up (down) to the nearest integer.

The velocity of each kicked particle is modified according to

(10)

where is a unit vector that is either oriented in a random direction on the unit sphere, or in the direction from the galaxy to the wind particle.

We have only studied kinetic feedback models where the wind speed is a constant for all galaxies. Similar models are often employed in self-consistent hydrodynamical simulations (e.g. Navarro & White 1993; Springel & Hernquist 2003; Dalla Vecchia & Schaye 2008). However, there are other possibilities, such as momentum-driven winds, where the wind speed scales with the galaxy velocity dispersion (e.g. Martin 2005; Oppenheimer & Davé 2006), and models where the outflow velocity increases with galactocentric radius (Steidel et al., 2010).

3.2 Metal enrichment models

We distribute metals amongst gas particles in our simulations as follows. For each model galaxy, all gas particles contained within a sphere centred on the galaxy are identified. As before, the radius of the sphere is chosen to be a fraction, , of the halo virial radius, where , or . The metals in diffuse gas produced by the galaxy are then shared evenly amongst these particles, so that the metal mass associated with each particle, , increases by an amount

(11)

where is the number of gas particles inside the sphere.

Given the total mass in metals for a gas particle, we could then define its metallicity simply as

(12)

which we refer to as the particle metallicity. However, in this work, we prefer to use the smoothed metallicity (Okamoto et al., 2005; Tornatore et al., 2007), defined by

(13)

where the smoothed metal mass density, , is computed in an analogous way to the standard SPH density estimate:

(14)

Here is the number of SPH smoothing neighbours and is a spherically-symmetric smoothing kernel, which depends upon the separation of particles and , , and the smoothing length of particle , . The smoothed metallicity of a gas particle is updated whenever its SPH density is calculated. Once a gas particle is converted into a star particle, its smoothed metallicity remains fixed for the rest of the simulation. Note that the metallicity of particles does not affect the gas dynamics in our non-radiative simulations.

3.3 Naming conventions

Table 2 lists all of our SN feedback models. Note that each model, including the reference gravitational heating only model (GO), follows the conversion of gas into stars as dictated by the underlying SA model.

Model name Type Energy injection method Comments
GO - - Gravitational heating only
SN_Th_N Thermal Fixed energy , ,
SN_Th_R Thermal Fixed energy , ,
SN_Th_T Thermal Fixed energy , ,
SN_En_N Thermal Fixed entropy , ,
SN_En_R Thermal Fixed entropy , ,
SN_KiR_V Kinetic Velocity kick, random /km s, 300, 600, 1000
SN_KiD_V Kinetic Velocity kick, directed /km s
SN_KiR_Z Kinetic Velocity kick, random , , ; /km s
Table 2: Supernova feedback models. Unless otherwise stated, the radius for metal injection is (). Note that all models, including the GO model, follow the conversion of gas into stars.

3.4 Entropy profiles

To test the effect of different implementations of SN feedback on the entropy structure of the ICM, we have resimulated our fiducial cluster, C1, with each of our models. Figure 1 shows the resulting entropy profiles. For comparison, we also show the observed entropy profiles of CC and NCC clusters in the REXCESS sample (Pratt et al., 2010, hereafter PAP10) . To facilitate a fair comparison, we only plot profiles of observed clusters that have a mass, , within of that of cluster C1.

Figure 1: Entropy profiles for cluster C1 resimulated with different implementations of supernova feedback (coloured lines; see legend for model details). Note that the profile obtained from a gravitational heating only model (GO) is also shown. For comparison, we display observed profiles of similar-mass CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10).

The main point to note is that all of our models yield almost identical entropy profiles that are in good agreement with the profile obtained from the reference GO run. In all cases, the profiles scale approximately as for , consistent with spherical accretion models (e.g. Tozzi & Norman 2001) and cosmological simulations that include gravitational heating only (e.g. Voit et al. 2005; Nagai et al. 2007). For , the entropy profiles flatten off significantly, exhibiting a small spread in central entropy.

Compared to the observed entropy profiles of CC clusters, the profiles predicted by our models have a steeper slope at , and the normalisation is systematically too low. In the case of NCC clusters, it is evident that none of our models can explain the shallow profiles characteristic of these systems.

We have checked that these results hold for other clusters in our sample, so we conclude that SNe have a negligible impact on the thermodynamical properties of intracluster gas and, furthermore, the manner in which the feedback energy is injected is unimportant.

3.5 Metallicity profiles

The metal enrichment model has only one free parameter, , which controls the radius of the spherical region about a galaxy in which metals are injected. We have resimulated cluster C1 with three different values of , fixing the SN feedback scheme to be the kinetic model where gas particles are kicked in a random direction with km s. Figure 2 shows the emission-weighted metallicity profiles that result, along with observed Fe abundance profiles of CC and NCC clusters from Matsushita (2011, hereafter MAT11). We plot all observed clusters with a mass above of that of cluster C1, in order to obtain a reasonable number of both CC and NCC objects (most NCC clusters in the sample of MAT11 are considerably more massive than cluster C1).

Figure 2: Emission-weighted metallicity profiles for cluster C1 resimulated with the same kinetic supernova feedback model, but varying the radius of the region within which metals are injected (solid coloured lines). See the legend for details of the metal enrichment models adopted. Observed profiles of CC (dashed grey lines) and NCC (solid black lines) clusters from the sample of MAT11 are also shown.

We note that it is difficult to directly compare the metallicity profiles of our clusters with those of observed clusters from MAT11, for several reasons. Firstly, the observed profiles are Fe abundance profiles, but our simple metal enrichment model does not include the contribution from Type 1a SNe, a major source of Fe, nor does it track the production of individual chemical elements. Secondly, we could, in principle, adjust the yield in the SA model underpinning our simulations, which would allow us to alter the normalisation of our metallicity profiles. For these reasons, we focus on the shape of metallicity profiles, instead of their normalisation, when assessing the impact of different feedback and enrichment models on the ICM enrichment pattern.

It is apparent from Figure 2 that varying the metal injection radius has a large impact on ICM metallicity profiles in core regions, . If metals are injected in a concentrated fashion (), then we see a sharp peak in the metal distribution within that region that is not reflected in the observational data. As the injection radius is increased, the gradient of the profile becomes progressively shallower until, when , the slope provides a good match to that of the profiles of observed NCC clusters, except the few systems that have a high central abundance more typically found in CC clusters. It is possible that these systems are CC remnants. Recall that radiative cooling is not included in our simulations so we do not expect to be able to reproduce the abundance peaks seen in the core regions of CC clusters.

A metal enrichment model where metals are distributed throughout the halo can be justified if the bulk of the metals found in intracluster gas were brought in by infalling material, rather than being produced by star formation in the central galaxy of the halo. To check whether this is the case, we have modified the L-Galaxies SA model to follow what fraction of metals in diffuse halo gas are produced by the central galaxy of the halo. Figure 3 shows this fraction as a function of halo virial mass for all clusters in our sample, and for a selection of halos with taken from the Millennium Simulation galaxy catalogues. For all of our clusters, the fraction of metals in hot halo gas produced by the central galaxy is less than , implying that nearly all of the metals in diffuse gas are indeed accreted. Note that this may change somewhat when we extend L-Galaxies to track the time-dependence of metals returned by Type Ia SNe and AGB stars, as some of the metal production will be delayed until after the formation of the central cluster galaxy.

Figure 3: Fraction of metals in hot halo gas produced by the central galaxy of the halo, , as a function of virial mass. The black points are for a subsample of halos with extracted from the Millennium Simulation galaxy catalogues. The red circles correspond to our resimulated clusters. For these massive systems, under of metals in the hot gas are produced by the central galaxy, implying that nearly all of the metals are accreted.

Ideally one would like to inject metals locally about satellite galaxies falling into a halo, rather than distributing them uniformly throughout the halo. However, this is not possible with the De Lucia & Blaizot (2007) version of L-Galaxies since all the metals in the diffuse gas associated with a galaxy are assumed to be instantaneously stripped once it becomes a satellite galaxy. In future work we plan to switch to a different treatment of satellite galaxies whereby hot gas is gradually removed from infalling galaxies by tidal and ram-pressure stripping (e.g.  Henriques & Thomas, 2010).

We now examine whether changing the SN feedback scheme affects ICM metallicity profiles. To do this, we have resimulated our fiducial cluster with all of the models described in Section 3.1, fixing in each case. Figure 4 compares the emission-weighted metallicity profiles obtained from our cluster simulations with those of the same observed clusters from the MAT11 sample.

Figure 4: Emission-weighted metallicity profiles for cluster C1 resimulated with different supernova feedback schemes, assuming that metals are distributed uniformly throughout the halo (coloured lines; see the legend for feedback model details). We also show profiles of CC (dashed grey lines) and NCC (solid black lines) clusters from the observational sample of MAT11.

The main point to note is that the metallicity profiles obtained from all our various simulations are essentially the same and we have checked that this conclusion remains valid when metals are injected in a concentrated manner (). It follows that SN feedback has no impact on the metal distribution in clusters, and the precise way in which the energy available from SNe is used to heat intracluster gas is irrelevant.

The metallicity profiles of clusters in the presence of supernova feedback were investigated by Tornatore et al. (2007). They also found that changing the SNe feedback rate makes little difference to the slope of the profiles (although it does change the normalisation). We show in Figure 11 and Section 5.3 that the stronger AGN jet feedback can have a larger effect.

3.6 Summary

Our study so far has revealed that feedback from SNe has a negligible effect on both ICM entropy and metallicity profiles, regardless of the manner in which the energy is assumed to be transferred to the gas.

In light of this freedom, we choose our fiducial SN feedback scheme to be the kinetic model SN_KiR_V600, where gas particles are given a kick in a random direction with km s. This is similar to the model of Dalla Vecchia & Schaye (2008). A wind speed of km s is consistent with observations of local (e.g. Veilleux et al. 2005) and –3 (e.g. Steidel et al. 2010) starburst galaxies.

The metallicity profiles reflect the manner in which metals are injected into the diffuse gas. For our fiducial metal enrichment model we assume that the metals ejected from galaxies are distributed uniformly throughout the entire halo, since this gives a good match to the slope of the metallicity profiles of observed NCC clusters. This model is justified by the fact that nearly all of the metals in intracluster gas are accreted, rather than being produced by the central galaxy of the halo.

The model that forms the basis for the rest of the work presented in this paper is thus SN_KiR_Z1.

4 Feedback from active galactic nuclei

As shown in the previous section, the heating of intracluster gas by stellar feedback alone clearly cannot account for the excess entropy observed in cluster cores, indicating that an additional feedback mechanism must be at play. The favoured candidate is the energy liberated by the accretion of gas onto central supermassive black holes at the centres of galaxies. Our goal in this section is to assess how the properties of the ICM are altered by the inclusion of this extra heating from AGN, and how our results are affected by different numerical implementations of AGN feedback.

In Section 4.1 we describe simple AGN heating models, similar to those found in the literature, then in Section 4.2 we use comparisons with observed entropy profiles to conclude that none of these models are entirely satisfactory. One model with extreme wind speeds does provide an adequate fit to the data and that motivates the new stochastic heating model developed in Section 5.

4.1 AGN feedback models

The amount of energy available from AGN heating, , is not arbitrary but is set by the model described in Secton 3.1.2 of SHT09. Although that paper considered only a single heating model, we find that the global X-ray lumnosity-temperature relation is dependent mainly upon the normalisation of the heating, and is relatively unaffected by the particular manner in which the heat is injected. That can have a large effect on the entropy profiles, however, as we show below.

4.1.1 Thermal models

The first set of thermal AGN feedback models that we have tested are identical to the thermal SN feedback models described previously in Section 3.1.1, except with replaced by . Similar prescriptions for AGN feedback have been employed in numerous other works (e.g. Springel et al. 2005a; Di Matteo et al. 2008; Booth & Schaye 2009; Fabjan et al. 2010).

4.1.2 Kinetic models

We implement kinetic AGN feedback in our simulations in the same way as kinetic SN feedback; see Section 3.1.2. The only differences are that the number of particles kicked (equation 9) now depends on the energy available from black hole accretion, , rather than that available from SN explosions, , and we have adopted larger wind speed values, km s, km s and km s, in line with measured AGN outflow velocities (Pounds et al., 2003; Chartas et al., 2003; Crenshaw et al., 2003; Ganguly & Brotherton, 2008; Dunn et al., 2010).

4.1.3 Naming conventions

Table 3 lists all of our various AGN feedback models. In each case the SN feedback and metal injection schemes are the same as for model SN_KiR_Z1.

Model name Type Energy injection method Comments
AGN_Th_N Thermal Fixed energy , ,
AGN_Th_R Thermal Fixed energy , ,
AGN_Th_T Thermal Fixed energy , ,
AGN_En_N Thermal Fixed entropy , ,
AGN_En_R Thermal Fixed entropy , ,
AGN_KiR_V Kinetic Velocity kick, random /km s 1000, 4 500, 20 000
AGN_KiD_V Kinetic Velocity kick, directed /km s 20 000
Table 3: AGN feedback models. In each case, supernova feedback is implemented using a kinetic model where particles neighbouring a galaxy are given a kick in a random direction with velocity km s, and ejected metals are assumed to be distributed uniformly throughout the entire host halo.

4.2 Entropy profiles

In order to assess how sensitive the thermodynamical properties of the ICM are to different implementations of AGN feedback, we have resimulated our fiducial cluster with each of our thermal and kinetic AGN feedback models.

4.2.1 Thermal models

The entropy profiles obtained from our thermal models are displayed in Figure 5. For comparison, we also show the profile predicted by our fiducial SN feedback model (SN_KiR_Z1), and entropy profiles of observed clusters of a similar mass.

As for the SN feedback, all of the thermal models give very similar results, even though there are large differences in the way in which the available energy is shared amongst gas particles in the various schemes. Essentially, the profiles follow the predicted scaling at large radii, but as we move in towards the core they begin to flatten off at . At radii interior to this, the slope is shallower than seen in either CC or NCC clusters, leading to an over-estimate of the central entropy. Simple preheating models predict similarly large isentropic cores at (e.g. STY10).

Figure 5: Entropy profiles for cluster C1 resimulated with different implementations of thermal AGN feedback (coloured lines; see the legend for model details). For comparsion, we also show the profile obtained from a run with kinetic supernova feedback only, model SN_KiR_Z1, and observed profiles of CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10).

4.2.2 Kinetic models

Figure 6 compares the entropy profiles predicted by our kinetic AGN feedback models with the same set of observed cluster profiles as in Figure 5.

In cluster outskirts, , the models give similar results, but there are clear difference at radii less than this. For the lowest wind speed, km s, we see a very flat entropy profile, with a hint of an entropy inversion in the core. As the wind speed is increased, the entropy profile steadily steepens, providing a good match to observed NCC cluster profiles when km s. Almost identical results are obtained when kicks are imposed in the direction from the galaxy to the wind particle, rather than in a random direction.

To understand this behaviour, we have checked how the trajectories of kicked particles are affected by variations in the wind speed. This is done by identifying the main progenitor of our cluster at high-redshift using the halo merger trees, selecting all gas particles within of this object that have just received a kick, then tracking the cluster-centric positions of these particles to . For a high wind speed, km s, the available AGN energy is only sufficient to kick a small number of particles and we find that their large momentum carries them beyond . As we reduce the wind speed to km s, the number of particles kicked increases but their momentum gain is smaller, so they do not escape from the cluster core before their kinetic energy is converted to thermal energy. This leads to an increase of the gas entropy in the central regions, establishing a flat entropy profile as seen in Figure 6.

Figure 6: Entropy profiles for cluster C1 resimulated with different kinetic AGN feedback models (coloured lines; see legend for model details). For comparsion, we also show the profile obtained from a run with kinetic supernova feedback only, model SN_KiR_Z1, and observed profiles of CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10).

4.3 Summary

Several interesting results have emerged from our study of the effect of different AGN feedback models on the thermal properties of the ICM.

We have found that simple thermal feedback schemes, based on heating a fixed number of particles, heating particles within a fixed fraction of the virial radius, or heating particles by a fixed fraction of the virial temperature, all heat the gas in cluster central regions excessively, leading to a higher core entropy than observed. All of the thermal models we have tested give very similar results.

When AGN feedback is implemented in a kinetic manner, ICM entropy profiles are found to be sensitive to the wind speed adopted. For low wind speeds, the resulting entropy profiles are too flat, as in the thermal case. This is because the available energy is shared amongst a large number of particles, and kicked particles do not have sufficient momentum to escape central cluster regions before their kinetic energy is thermalised.

As the wind speed is increased, the number of particles kicked decreases as and kicked particles are able to reach larger cluster-centric radii before thermalisation of their kinetic energy. Consequently, more low-entropy material remains in core regions, so entropy profiles become progressively steeper, approaching observed ones. For km s, the predicted profiles agree well with observed profiles of NCC clusters. However, such a high wind speed is perhaps physically unrealistic.

From our discussion, it seems that the key ingredient of a successful AGN feedback model must be to ensure that only a small fraction of particles in central cluster regions are heated/kicked, so that these particles have sufficient entropy/momentum to reach cluster outskirts, leaving low-entropy metal-rich gas behind in the core. In the next section we formulate a new AGN feedback prescription that has this desired feature, and is motivated by the observed interaction of AGN with their environment.

5 A new model for feedback from active galactic nuclei

There is a growing body of observational evidence that AGN feedback may be mostly related to radio-loud AGN. In the local universe, observations of galaxy groups and clusters often show X-ray cavities coincident with lobes of radio emission linked to the central galaxy by radio jets (Blanton et al., 2001; Bîrzan et al., 2004; McNamara et al., 2005; Fabian et al., 2006; Morita et al., 2006; Jetha et al., 2008; Gastaldello et al., 2009; Dong et al., 2010; Giacintucci et al., 2011). It is thought that these bubbles are inflated by the central AGN, and may provide an efficient means of removing cool, enriched gas from cluster cores as they rise buoyantly through the cluster atmosphere, thus quenching star formation.

At high redshift, , emission-line kinematics of radio galaxies based on rest-frame optical integral-field spectroscopy have revealed powerful bipolar outflows with kinetic energies equivalent to of the rest mass of the central supermassive black hole (e.g. Nesvadba et al. 2006, 2008). These AGN-driven winds are energetic enough to remove copious amounts of gas from the host galaxy, preventing further accretion onto the black hole and suppressing star formation. Large-scale energetic outflows have also been observed in ultra-luminous infrared galaxies (Alexander et al., 2010), a galaxy population potentially an order of magnitude more common than distant radio galaxies.

Although it is not yet fully understood how the energy released by black hole accretion is transferred to the surrounding gas, the observational data suggests that the energy is input in a directional manner, via jets or collimated outflows, rather than isotropically. To reflect this, we have developed an anisotropic, stochastic heating model where only some of the gas particles neighbouring a galaxy are heated per duty cycle of the AGN. We note that higher-resolution models of feedback from AGN in cluster cores also favour anisotropic heating (e.g.  Gaspari et al., 2012).

In Section 5.1 we describe this heating model in detail, then in Section 5.2 we use observed entropy profiles and X-ray scaling relations to determine optimal model parameters.

5.1 Stochastic AGN feedback model

The basis of our new model for AGN heating is as follows. For each galaxy, we first identify all gas particles contained within a sphere centred on the galaxy, where the radius of the sphere is some fraction, , of the halo virial radius. We then assume that the probability that any of these particles has been heated by AGN feedback during the time elapsed, , since the previous SA model output is

(15)

where is a parameter controlling the fraction of particles heated over the AGN duty cycle, . Based on observational data, we take yrs (e.g. Bîrzan et al. 2004; Fabian et al. 2006; Jetha et al. 2008). With our choice of SA model output times we then have for .

For each gas particle neighbour, we draw a random number uniformly from the unit interval and compare it with : if the particle is given an entropy boost

(16)

and if the particle is not heated. By including the heating probability in the denominator of equation (16), we ensure that the total amount of energy injected into the gas is (approximately) the same for different choices of .

We have also experimented with supplying the AGN heat energy to particles as a fixed energy boost. However, this makes virtually no difference to our results so we do not discuss these models hereafter.

There are two free parameters in our model: and . The values of these parameters we have tested in this work are , and , and , , and . Note that in the case our model reduces to AGN_Th_R.

It is interesting to link the parameter to the opening angle of AGN jets. If we make the simple approximation that large-scale AGN outflows can be treated as biconical jets, each with opening angle , then it follows that . For the range of values of tested here, this corresponds to .

Table 3 lists all of our stochastic AGN feedback models. To distinguish these from the AGN heating models of the previous section, we have given them the label JET.

Model name Type Energy injection method Comments
JET_R_D Stochastic Fixed entropy , ,
, , ,
JET_Z Stochastic Fixed entropy , ,
Table 4: Stochastic AGN feedback models. In each case, supernova feedback is implemented using a kinetic model where particles neighbouring a galaxy are given a kick in a random direction with velocity km s. Unless otherwise stated, the radius for energy and metal injection is ( and , respectively), and the fraction of particles heated per AGN duty cycle is .

5.2 Entropy profiles and scaling relations

We now investigate whether our new physically-motivated stochastic AGN feedback scheme yields a better match to observed cluster profiles than the simple thermal and kinetic models discussed in the previous section.

Recall that our stochastic model has two free parameters: , which governs the radius about a galaxy in which energy is injected, and the fraction, , of neighbouring gas particles that are heated per AGN duty cycle. The first issue to address is how varying these parameters affects cluster properties. We then identify an optimal choice for these parameters by using a selection of observational data to constrain our model.

5.2.1 The effect of changing

Figure 7 shows the effect of varying on the entropy profile of cluster C1. In each case, is kept fixed at . It is apparent that the entropy structure of the ICM is relatively insensitive to the choice of , with only small differences between the three different runs. The best match to observed NCC cluster profiles arises when , in which case we find excellent agreement with the observational data.

Figure 7: Entropy profiles for cluster C1 resimulated with our stochastic AGN feedback model for different values of the parameter (solid coloured lines; see legend for model details). The other parameter in the model, , is fixed at . The profiles of observed CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10) are also displayed for comparison.

Figure 8 shows the - relation for our full -cluster sample for the three different values of tested. Here also, we can see that the three different models predict a very similar - relation. Note that the trend in - with is not monotonic: excessively large and excessively small values of will both leave behind low-entropy core particles. However, the variation is small and each of the chosen values yields an adequate match to the observed relation for NCC clusters, with providing the best match of the three.

Figure 8: The X-ray luminosity-temperature scaling relations predicted by our stochastic AGN feedback model with (asterisks), (triangles) and (circles), keeping the other model parameter, , fixed at . See the legend for model names. X-ray properties are computed within . For comparative purposes, we also plot observational data for CC (diamonds) and NCC (squares) clusters in the REXCESS sample (PCA09).

5.2.2 The effect of changing

We now turn our attention to the effect of the parameter . We have done four runs, with , , and , respectively, keeping fixed at unity. The entropy profile of cluster C1 in each case is displayed in Figure 9. It is immediately clear that varying has a much larger effect on the entropy of intracluster gas than . As is increased from to , the slope of the entropy profile at radii becomes progressively shallower. For the slope is too steep compared to that of observed NCC cluster profiles, whereas it is too flat for . The values and both give a good match to the observational data for NCC clusters.

To explain the variation in cluster entropy profiles with , we have again examined what happens to particles that are heated by AGN feedback in each of our runs. In the case where , the probability of a particle being heated is low, but any particle that is heated receives a large entropy boost since appears in the denominator of equation (16). The high entropy of heated particles causes them to rise buoyantly to large cluster-centric radii, , leaving the entropy profile in the core relatively undisturbed compared to a run with SN feedback only. When is increased to , many more particles in central regions are heated by AGN feedback since is larger, and the entropy boost they are given is smaller. Accordingly, the distance they move outwards from the core is less, resulting in a higher central entropy and a flatter profile.

Figure 9: Entropy profiles for cluster C1 resimulated with our stochastic AGN feedback model for several choices of the parameter (solid coloured lines; see legend for model details). The other model parameter, , is set to unity. For comparison, we also show the profiles of observed CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10).

The large impact of on cluster entropy profiles is reflected in the - scaling relation. This is demonstrated in Figure 10 where we show the - relation for our full cluster sample predicted by each of our four models. For , the low central entropy causes an enhanced X-ray luminosity, so all of our simulated clusters lie well above the mean observed relation for NCC clusters. In fact, the predicted relation in this case resembles the observed relation for CC clusters, although this is artificial since we have not included cooling processes. As is increased, the normalisation of the - relation decreases and the slope becomes steeper. A good match to the observed NCC cluster - relation is obtained when . For larger values of , the relation is too steep, so that low-temperature systems have too low a luminosity for their mass. This is because the AGN heating has raised the core entropy in these systems to an excessive level.

Figure 10: The X-ray luminosity-temperature scaling relations predicted by our stochastic AGN feedback model with (asterisks), (triangles), (crosses) and (circles), keeping the other model parameter, , set to unity. See the legend for model names. X-ray properties are computed within . Observed CC (diamonds) and NCC (squares) clusters from the REXCESS (PCA09) are also displayed.

5.2.3 Identifying optimal parameter values using observations

The next issue to address is whether observational data can help us to constrain the two free parameters of our stochastic AGN feedback model. For this analysis, we have resimulated all clusters in our sample (C1–C25) using different combinations of the pair of parameters . The values adopted are , and , and , , and , giving a grid of models in total.

We assess the suitability of each model by testing how well it reproduces the observed scaling of three fundamental ICM observables with spectroscopic temperature: (i) the entropy profile normalisation, which we take to be the entropy at (typically about ), (ii) the entropy profile shape, defined as the ratio of the entropy at to the entropy measured at (i.e. the central entropy), and (iii) the X-ray luminosity. Again, the source of the observational data is the REXCESS. Note that we are only aiming to match the observed scaling relations for NCC clusters since cooling is not included in our simulations at this stage.

In the following analysis, we neglect any clusters in our simulated sample that have large amounts of substructure. To differentiate between dynamically relaxed and disturbed systems, we use the substructure statistic

(17)

where is the location of the dark matter potential minimum and is the centre of mass of the cluster within . Following Kay et al. (2007), we say that a cluster is disturbed if , and relaxed otherwise.

For each scaling relation, the criterion we use to test how well our models reproduce the mean observed relation is the statistic

where is the number of (relaxed) simulated clusters and , or , depending on the relation being considered. The quantities and are the normalisation and slope of a power-law fit to the corresponding observed relation,

(19)

obtained by using the BCES orthogonal linear regression method (Akritas & Bershady, 1996) in log-log space, taking into account the errors in both any . The normalisation has units of cm and for and , respectively, and is dimensionless for . The factor is included to remove the predicted self-similar evolution, where the index , and for the -, - and - relations, respectively.

The scatter expected from statistical uncertainties, , is

(20)

where

(21)

and and are the errors in any , respectively.

We estimate the raw scatter, , using error-weighted distances to the regression line:

where and is the number of observed NCC clusters in the sample.

Finally, the intrinsic scatter, , about each observed mean relation is estimated as

(23)

We examine how well our models reproduce the observed scatter about each mean relation, by using the Kolmogorov-Smirnov (K-S) test to determine if the residuals for the observed and simulated samples are drawn from the same distribution in each case.

Tables 5 and 6 show the probabilities (-values) for each of our models obtained from the (equation 5.2.3) and K-S tests, respectively, for the case of the - relation. The models we deem to be acceptable are highlighted in bold. It is evident that all models with are ruled out, at least for this particular relation, and only three models provide an acceptable match to both the mean observed relation and the associated scatter.

\backslashbox 0.1 0.32 1.0
0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.20 0.96
0.69 0.00 0.01
Table 5: test probability values for the - scaling relation.
\backslashbox 0.1 0.32 1.0
0.00 0.01 0.00
0.00 0.02 0.00
0.41 0.12 0.24
0.73 0.00 0.02
Table 6: K-S test probability values for the - scaling relation.

For the sake of brevity, we do not present the corresponding tables for the - and - relations. We simply note that the test for the - relation rules out all models expect the two with and , whereas the K-S test rules out all models with or . In the case of the - relation, both the and K-S tests rule out all models with , tending to favour models that populate the upper-right corner of the table.

We have combined the results of all tests ( for each of the scaling relations) using Fisher’s method for combining -values. The overall probabilities for our models are summarised in Table 7. It is evident that only one model is now acceptable, the model with . Therefore, we choose our fiducial AGN feedback scheme to be the stochastic model JET_R1_D.

\backslashbox 0.1 0.32 1.0
0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.34
0.00 0.00 0.00
Table 7: Combined probability values for all three scaling relations.

We emphasise that the purpose of this section has not been to conduct a rigorous statistical analysis, but merely to provide us with an indication of the region of our model parameter space that is favoured by the observational data. There are several possible caveats to our analysis. For example, we are assuming Gaussian errors and that each test is strictly independent, both of which may not be the case.

5.3 Metallicity profiles

Powerful AGN outflows are an obvious candidate for transporting metal-rich material away from the central regions of halos. Indeed, there is observational evidence that enriched gas is entrained by bubbles inflated by central AGN and removed from cluster cores (e.g. Forman et al. 2005; Million et al. 2010; Kirkpatrick et al. 2011). Cosmological hydrodynamical simulations have also demonstrated that AGN are important for the metal enrichment of intracluster gas (e.g. Sijacki et al. 2007; Moll et al. 2007; Fabjan et al. 2010; Wiersma et al. 2011).

Recall from Section 3.5 that SN feedback alone has a negligible impact on the distribution of metals in the ICM, regardless of the way in which the available energy is injected into the diffuse gas. We now investigate whether the inclusion of extra energy input from AGN affects ICM metallicity profiles. To do this, we have resimulated cluster C1 with our fiducial AGN feedback model (JET_R1_D), varying the parameter in our metal enrichment scheme (recall that this parameter sets the radius of the spherical region about a galaxy within which metals are injected).

Figure 11 compares the resulting emission-weighted metallicity profiles with those of observed clusters from the sample of MAT11. As before, the best match to the gradient of the profiles of observed NCC clusters arises when metals are distributed uniformly throughout the entire halo (). As is decreased, we see the development of a sharp central abundance peak that is in conflict with the observational data. Recall that a distributed metal enrichment model is justified since the SA model underpinning our simulations predicts that almost all metals in the ICM are accreted, rather than being produced by BCGs.

Comparing the metallicity profile predicted by model JET_Z1 to that obtained from the SN_KiR_Z1 run, we can see that addition of AGN feedback has indeed displaced some metals from central cluster regions, leading to a flatter profile, but the effect is small. This is because only a small fraction of the particles in the core are heated by AGN in our model, and they receive a sufficiently large entropy boost to escape the central regions of clusters, leaving the majority of metal-rich material behind in the core.

Figure 11: Emission-weighted metallicity profiles for cluster C1 resimulated with our fiducial stochastic AGN feedback model, but with different metal enrichment schemes (solid coloured lines; see legend for model details). The profile obtained from a run with supernova feedback only is shown for comparison, as well as the profiles of observed CC (dashed grey lines) and NCC (solid black lines) clusters from the sample of MAT11.

5.4 Summary

We have developed a new AGN feedback model where heat energy is injected into intracluster gas in a stochastic manner. Our model is physically-motivated and has just two free parameters: , which sets the radius of the spherical region within which energy is injected, and , which is the fraction of particles in this region that are heated per AGN duty cycle.

We have found that ICM entropy profiles and the - scaling relation are fairly insensitive to variations in , but depend strongly on . For small values of , the resulting entropy profiles are close to those obtained from a run with SN feedback only, implying that AGN heating has little effect in this case. This is because only a few particles are heated, and they receive large entropy boosts, which causes them to rise buoyantly to large distances from the cluster centre, leaving the majority of low-entropy gas behind.

As is increased, the probability of a particle being heated also increases: hence more particles in the core are heated, they are given a smaller entropy injection, and so they do not escape central cluster regions. As the heated gas expands, the gas density drops, causing the gradient of the resulting entropy profiles to become shallower.

We have used three observed scaling relations to identify an optimal choice for the two free parameters in our model: . Setting roughly corresponds to a jet opening angle of . With these parameter choices, our model, named JET_R1_D, can explain the observed entropy profiles and - relations for NCC clusters.

Using JET_R1_D as our fiducial AGN feedback model, we have demonstrated that AGN heating has little impact on the distribution of metals in the ICM by comparing to a model with kinetic supernova feedback only. When the metals produced by stars in galaxies are distributed uniformly throughout the entire host halo (model JET_Z1), the resulting abundance gradients provide a good match to those observed in NCC clusters.

For the remainder of this paper we thus adopt model JET_Z1 as our fiducial model for star formation, metal production, black hole growth and associated stellar and AGN feedback.

6 Comparison with observations

In this section we conduct a detailed assessment of how well our fiducial model (JET_Z1) can reproduce key observed thermal and chemical properties of intracluster gas.

6.1 Thermal properties of the ICM

Figure 12 compares the predicted entropy profiles of all clusters (C1–C25) in our sample with the profiles of observed systems in the same mass range. The profiles of relaxed (disturbed) simulated clusters are shown as solid (dotted) red lines.

Figure 12: Entropy profiles for clusters resimulated with our fiducial stochastic AGN feedback model. The profiles of relaxed (disturbed) systems are shown by solid (dotted) red lines. Observed profiles of CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10) are also displayed for comparative purposes. The observed clusters span the same mass range as our simulated ones.

First impressions are that our fiducial model generates clusters whose entropy profiles agree well with those of observed NCC systems, both in terms of normalisation and gradient. The central entropy is too high in three of our objects, but two of these are classified as disturbed systems. To assess our model more quantitatively, we now examine how the entropy profile normalisation and shape scale with system temperature.

Figure 13 shows the entropy profile normalisation (defined as the entropy at ) as a function of spectroscopic-like temperature. Filled (open) circles correspond to relaxed (disturbed) objects (both here and in all subsequent figures). Observational data for NCC clusters in the REXCESS are also shown. The parameters of the accompanying predicted and observed best-fit relations are summarised in Table 8. Note that we only consider relaxed systems in our sample when performing the fit and, as before, we adopt the BCES orthogonal fitting method. It is evident that the - relation predicted by our model is a good match to the observed relation: the slope and normalisation are both within of the observed values, and the scatter about the mean relation is comparable.

Figure 13: The scaling of entropy profile normalisation, , with temperature predicted by our fiducial stochastic AGN feedback model. We compute the spectroscopic-like temperature within . Filled (open) circles correspond to relaxed (disturbed) systems in our simulated cluster sample, and the solid red line is the best-fit relation considering relaxed systems only. The observed relation for NCC clusters in the REXCESS sample (PAP10) is also shown by open squares and a solid black line.
Relation Predicted Observed
-
-
-
-

and are the best-fitting normalisation and slope of the relations, respectively (see equation 19), and is the intrinsic scatter about the mean relation (equation 23).

Table 8: Best-fit parameters (with errors) for scaling relations obtained from our full 25-cluster simulated sample, and from the REXCESS observations of PAP10. Note that we only consider relaxed clusters in our sample when deriving predicted relations. All fits were performed using the BCES orthogonal regression method.

In Figure 14 we display the variation of the ratio (a measure of the entropy profile shape) with temperature. The slope of our predicted relation is consistent with that of the observed relation and the scatter is identical to the observed value; see Table 8. However, the normalisation is slightly lower. There are two reasons for this offset. First, one of our relaxed clusters has an anomalously low value of , which lowers the normalisation of the predicted - relation. This objects correspond to the relaxed system with an excessive central entropy in Figure 12. Second, three of the observed clusters lie considerably above any of our simulated clusters on the - plane. This acts to raise the normalisation of the observed relation relative to the predicted one. Although classified as NCC systems in REXCESS, these objects actually have a low central entropy, reminiscent of CC clusters; see Figure 12. Without these outliers, there is good overall agreement between the predicted and observed - relations.

Figure 14: The scaling of entropy profile shape, , with temperature predicted by our fiducial stochastic AGN feedback model. The spectroscopic-like temperature is computed within . Relaxed (disturbed) systems in our simulated cluster sample are shown as filled (open) circles. The solid red line is the best-fit relation obtained using our relaxed clusters only. For comparison, we also display the observed relation for NCC clusters from the REXCESS (PAP10; open squares and black line).

Finally, we contrast our predicted - scaling relation with the REXCESS NCC cluster relation in Figure 15. We predict slightly less scatter about the mean relation than observed, but we recover the normalisation and slope of the observed relation to within , as summarised in Table 8. We conclude that our fiducial model yields an - relation that is a good match to the observed relation for NCC clusters.

Figure 15: The X-ray luminosity-temperature scaling relation predicted by our fiducial stochastic AGN feedback model. X-ray properties are computed within . Filled (open) circles correspond to relaxed (disturbed) systems in our simulated cluster sample, and the solid red line is the best-fit relation for relaxed objects only. For comparative purposes, we also plot the observed relation for NCC clusters in the REXCESS sample (PCA09; squares and solid black line).

6.2 Chemical properties of the ICM

The predicted emission-weighted metallicity profiles of all clusters in our sample are displayed in Figure 16, along with Fe abundance profiles of observed CC and NCC clusters from the sample of MAT11. To ensure a fair comparison, we only plot the profiles of observed clusters that lie in the mass range spanned by our simulated sample. Again, solid (dotted) red lines correspond to relaxed (disturbed) systems. We remind the reader of the limitations of our metallicity model: we assume only prompt enrichment and cannot discriminate between ejecta from core collapse and Type 1a SN. Nevertheless, it is evident that the profiles of our simulated clusters are in reasonable agreement with those of observed NCC clusters, both in terms of normalisation and slope, with the exception of the few observed NCC objects that have a sharp central abundance peak, which could be CC remnants (Leccardi et al., 2010; Rossetti & Molendi, 2010). To demonstrate this more rigorously, we now investigate how the metallicity profile shape scales with temperature.

Figure 16: Emission-weighted metallicity profiles for clusters resimulated with our fiducial stochastic AGN feedback model. The profiles of relaxed (disturbed) systems are shown by solid (dotted) red lines. For comparison, we also show observed profiles of CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10). We only show observed clusters with a mass in the same range as our simulated objects.

The measure of metallicity profile shape we adopt is the ratio of the metallicity at a radius of to that at a radius of . We chose these particular radii since nearly all of the clusters in the sample of MAT11 have a metallicity profile defined over this radial range, thereby maximising the number of observed clusters we can compare our predictions to.

Figure 17 shows the predicted scaling of the ratio with spectroscopic-like temperature. The corresponding observed relation for NCC clusters is also shown for comparison, and the parameters of both best-fit relations are presented in Table 8. We recover the observed gradient to within , but the predicted normalisation is lower than observed, and the scatter is larger. However, the observational errors are large and there is one observed cluster that lies considerably above all the others, which will act to increase the normalisation of the observed relation relative to that of ours.

Figure 17: The scaling of emission-weighted metallicity profile shape, , with temperature predicted by our fiducial stochastic AGN feedback model. The filled (open) circles correspond to relaxed (disturbed) systems in our simulated cluster sample. The solid red line is the best-fit relation obtained when considering just the relaxed systems in our sample. The observed relation for NCC clusters in the sample of MAT11 is also displayed (squares and solid black line).

7 Including radiative cooling: A first attempt

None of the simulations presented in this paper thus far incorporate cooling processes. In this section, we make a first attempt to extend our hybrid feedback scheme by allowing gas to cool radiatively. Our aim is to formulate a feedback model which can produce both CC and NCC clusters, whilst avoiding catastrophic over-cooling of gas in central cluster regions. We emphasise that this work is exploratory, intended merely to demonstrate that such a model is possible with our approach.

The addition of gas cooling is likely to lead to differences between the predictions of kinetic and thermal feedback schemes that have not been apparent in our previous non-radiative runs. However, it is not our intention here to conduct an exhaustive comparison of different feedback models when cooling processes are included; we save this for future work.

We implement AGN feedback using the -parameter stochastic heating model developed in the Section 5 (however, as we shall see, the optimal parameter choices change with the addition of cooling), and we adopt our fiducial model for SN feedback (gas particles neighbouring a galaxy are given a kick in a random direction with a speed of km s) and metal enrichment (metals produced by stars in galaxies are uniformly distributed throughout the entire host halo).

Metal-dependent radiative cooling is included in our simulations as follows. For each gas particle, we know its (smoothed) metallicity, (see equation 13), and we can compute its temperature from its entropy, , and density, . With this information we then calculate the cooling rate using the cooling function of Sutherland & Dopita (1993), and reduce the entropy of the gas particle accordingly.

Figure 18 compares the entropy profile of cluster C1 obtained from runs with our fiducial feedback model (JET_Z1) with and without cooling. It is apparent that the heating from SN and AGN has not been sufficient to prevent over-cooling in central cluster regions: there is a sharp drop in gas temperature at , leading to a steep decline in the entropy profile, and there is an entropy increase at larger radii due to hotter, lower-density gas flowing inwards from cluster outskirts to maintain pressure support in the core.

Figure 18: Entropy profiles for cluster C1 resimulated with our fiducial stochastic AGN feedback model without cooling (solid red line), and with metal-dependent radiative cooling (solid blue line). Observed profiles of similar-mass CC (dashed grey lines) and NCC (solid black lines) clusters in the REXCESS sample (PAP10) are also displayed for comparative purposes.

A priori, there is no reason to expect that the amount of SN and AGN heating provided by the underlying SA model would be sufficient to precisely balance radiative cooling in the simulation. This is because L-Galaxies employs a simple cooling recipe based on the assumption that haloes have a spherically-symmetric isothermal gas distribution, which is typically not the case in hydrodynamical simulations, so the predicted cooling rate of gas in haloes will be different to that in the simulation. Since the amount of gas that can cool to form stars and accrete onto central black holes governs the level of subsequent feedback, such differences in gas cooling rates imply that it is unlikely a self-regulating feedback loop would be established in the simulation.

To address this problem, we have developed an ad hoc extension of our stochastic AGN feedback model where we inject extra energy into cluster cores as a crude representation of additional AGN heating that would have arisen from enhanced black hole accretion due to more efficient cooling. Such a scheme is justifiable, provided the extra energy input required to balance radiative cooling is a small fraction of that originally available from the SA model.

The details of our model are as follows. At each SA model output, we identify all gas particles in the simulation residing in the central regions of haloes (). At each subsequent timestep, we test if any of these particles have cooled below a threshold temperature of K; if they have, we raise their temperature to some multiple, , of the virial temperature, (equation 3), of their host halo at the previous output time. We continue in this fashion until the next model output is reached, at which point the list of particles contained in halo cores is reset and the process is repeated.

In what follows we keep the radius of energy injection in our AGN feedback model fixed at unity (). We then have two free parameters: the fraction of particles heated per AGN duty cycle, , and , which controls the temperature cold particles in cluster cores are heated to. We now explore the effect of varying these parameters on the entropy distribution in clusters. All of our models are summarised in Table 9, where we have given them the label ZCOOL to emphasise that they include metal-dependent radiative cooling.

Table 9: Stochastic AGN feedback models with metal-dependent radiative cooling and a prescription for additional heating of cold gas in cluster cores. In each case, supernova feedback is implemented using a kinetic model where particles neighbouring a galaxy are given a kick in a random direction with velocity km s. Energy and metals are both injected within a radius of ( and , respectively). Unless otherwise stated, the fraction of particles heated per AGN duty cycle is , and cold particles in cluster cores are heated to a temperature of times the halo virial temperature.
Model name Type Energy injection method Comments
ZCOOL_D Stochastic Fixed entropy