Exploring the Energetics of Intracluster Gas with a Simple and Accurate Model

# Exploring the Energetics of Intracluster Gas with a Simple and Accurate Model

Paul Bode and Jeremiah P. Ostriker Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544; bode@astro.princeton.edu, ostriker@princeton.edu Alexey Vikhlinin Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138; avikhlinin@cfa.harvard.edu
###### Abstract

The state of the hot gas in clusters of galaxies is investigated with a set of model clusters, created by assuming a polytropic equation of state () and hydrostatic equilibrium inside gravitational potential wells drawn from a dark matter simulation. Star formation, energy input, and nonthermal pressure support are included. To match the gas fractions seen in non-radiative hydrodynamical simulations, roughly 5% of the binding energy of the dark matter must be transferred to the gas during cluster formation; the presence of nonthermal pressure support increases this value. In order to match X-ray observations, scale-free behavior must be broken. This can be due to either variation of the efficiency of star formation with cluster mass , or the input of additional energy proportional to the formed stellar mass . These two processes have similar effects on X-ray scalings. If 9% of the gas is converted into stars, independent of cluster mass, then feedback energy input of (or keV per particle) is required to match observed clusters. Alternatively, if the stellar mass fraction varies as then a lower feedback of is needed, and if the stellar fraction varies as steeply as then no additional feedback is necessary. The model clusters reproduce the observed trends of gas temperature and gas mass fraction with cluster mass, as well as observed entropy and pressure profiles; thus they provide a calibrated basis with which to interpret upcoming SZ surveys. One consequence of the increased gas energy is that the baryon fraction inside the virial radius is % of the cosmic mean, even for the most massive clusters.

cosmology:theory — galaxies:clusters:general — intergalactic medium — X-rays:galaxies:clusters

## 1 Introduction

In the near future, clusters of galaxies will prove to be an increasingly important probe into the abundance and properties of both dark matter and dark energy. Data on large numbers of clusters are becoming available from surveys at many wavelengths: in the optical/infrared, from galaxy overdensities (Koesterea07; YeeGGMHE07; Eisenhardtea08; LopesCKJ09) or from weak gravitational lensing of background objects (Dahle07z); in the X-ray (Bohringerea07; BureninVHEQM07; Finoguenovea07; Pacaudea07); and in the microwave, by locating temperature decrements caused by the Sunyaev-Zeldovich (SZ) effect (Dobbsea06; Kosowsky06; MalteSchaferB07; Muchovejea07; LinWUKLNHLWH08; BonamenteJLCNM08; Staniszewskiea08z).

Properly interpreting such surveys will require a sound understanding of the state of the Intra-Cluster Medium (ICM). Based on the assumption that the thermal energy in the gas comes solely from gravitational collapse, one may derive expected self-similar scalings between halo mass and temperature, X-ray luminosity, etc.(Kaiser91; BryanN98; EkeNF98). However, these scalings do not match the actual behavior of observed clusters (e.g ArnaudEvrard99; VikhlininKFJMMVS06; ArnaudPP07; PrattCAB08z; SunVDJFV09; VikhlininBEFHJKMNQV09). Additional modification of the thermal state of the gas from non-gravitational sources will bring these predictions into better agreement with observed scalings (Kaiser91). One such modification is star formation; by removing gas with the shortest cooling times, star formation leaves behind gas with a higher mean entropy (VoitBryan01; TozziNorman01; VoitBBB02). Energy injection from preheating or stellar and AGN feedback likewise affects scalings, as numerous simulations which include such feedback have shown (e.g. EttoriDBM06; BorganDMCSDMTTT06; MuanwongKT06; NagaiKV07; KaydSABLPST07; SijackiSMH07; BhattacharyaMK08; BorganiDDS08; DaveOS08; PuchweinSS08; ShortT08z, and references therein). The importance of these processes can be investigated analytically by considering hydrostatic, and usually polytropic, gas residing in various dark matter potentials (SutoSM98; BaloghBP99; WuFN00; Loewenstein00; TozziNorman01; KomatsuSeljak01; BabulBLP02; VoitBBB02; DosSantosDore02; AscasibarYMS03; ShimizuKSS04; LapiCM05; AfshordiLS05; SolanesMGS05; AscasibarD08; BowerMB08; CiottiP08; MoodleyWGT08z).

OstrikerBB05 and BodeOWS07 developed a simple method for populating the halo Dark Matter (DM) potentials from N-body simulations with gas, based on the assumptions of hydrostatic equilibrium and a polytropic equation of state. This method can be calibrated against local X-ray observations, and then be used to create realistic sky maps (SehgalTHB07). Rather than assuming a simple analytic profile, the distributions of halo concentration, substructure, and triaxiality are included, as well as any spatial correlations, since the halos are drawn from an accurate large scale dark matter simulation. This procedure has advantages over performing full hydrodynamic simulations. Because less computational expense is required, large volumes can be simulated. Also, one can vary parameters, such as the amount of star formation and feedback, without having to redo an entire simulation.

This paper presents a number of improvements to the model presented in OstrikerBB05 and BodeOWS07. A more sophisticated description of star formation is included, which accounts for gas recycling and allows for a calculation of metallicity evolution. BodeOWS07 assumed that every cluster halo has the same stellar mass fraction; here we allow this fraction to vary with halo mass. Also BodeOWS07 did not allow for transfers of energy from the dark matter to the gas. It is evident from simulations that this does occur, so we allow for this transfer here. Because of these limitations, BodeOWS07 required a large amount of feedback energy to match X-ray observations. As we shall see, the improvements presented here significantly reduce the feedback required.

With this method, we examine how star formation and different forms of energy input affect the ICM. It is possible to produce model clusters which match existing observational constraints, and thus make predictions for upcoming surveys. The method is described in Sec. 2 and Appendix A. In Sec. 3 it is then applied to a simulation reflecting the current best-fit cosmological parameters and the results are compared to X-ray data; in particular the effects of varying star formation and feedback are examined. In Sec. 4 a standard model is compared to further X-ray and SZ observation; we summarize and discuss implications in Sec. 5.

## 2 Method

To create a mock catalog, cluster-sized halos are extracted from a dark matter only N-body simulation, and then the gas distribution inside each halo potential is determined. This section describes the gas model and the construction of the halo catalog which is used.

### 2.1 The Gas Prescription

A given DM halo is enclosed in a cubic mesh, using a cell size, , twice the N-body particle spline softening length. The DM density, , and gravitational potential, , are computed for each cell using a Particle-Mesh code with non-periodic boundary conditions. The cell with the lowest potential, , is considered the center of the cluster. Particle velocities in the rest frame of the cluster (defined as the mean velocity of the 125 particles closest to the cluster center) are used to find the KE per unit volume, , in each cell. The virial radius (from spherical top-hat collapse; BryanN98) can then be computed, as well as the halo virial mass, , and energy, .

It is assumed that gas originally followed the DM inside the virial radius, with density and KE , where . A certain amount of the gas, , will have turned into stars. We parameterize this mass by the fraction at z=0; the evolution of with redshift is described in Appendix A. This stellar mass presumably formed from the portion of the gas with the lowest entropy and shortest cooling time, so cells are ranked by binding energy and then are checked off, starting with the most bound, until the sum of their gas masses equals ; the gas mass and energy in these cells are set to zero. Thus the initial mass and energy of the remaining gas are:

 Mg=∑kfcρDkl3, (1)
 Eg=∑kfc{ϕkρDk+@frac12tDk}l3, (2)

summing over all cells inside except those marked off for star formation. The surface pressure on the gas at the virial radius is estimated from the kinetic energy in a buffer region eleven cells thick:

 (3)

summing over the cells in the buffer region .

From this initial state, it is assumed the gas is rearranges itself into hydrostatic equilibrium in the DM potential, with a polytropic equation of state of index (as indicated by both observation and simulation; OstrikerBB05, ). The resulting gas pressure and density are given by and , where the polytropic variable

 θk≡1+Γ−1Γ(1+δrel)ρ0P0(ϕ0−ϕk). (4)

Here is a nonthermal component of pressure, assumed to be everywhere proportional to thermal pressure; the total is thus . The pressure and density at the cluster center are found using two constraints: requiring conservation of energy, and matching the external surface pressure. For a given and , the final radius of the gas initially inside is found by summing outwards from the cluster center until the initial mass is enclosed. Note this means that gas may expand or contract, changing the gas fraction inside , and also doing mechanical work. Assuming surface pressure changes little with radius, the change in energy is proportional to the change in volume: . The constraint equation for conservation of energy is thus (see BodeOWS07, ):

 Ef=∑rk

Here there are two additional inputs to the energy budget. The penultimate term represents energy input to the gas from dynamical processes, assumed to be proportional to ; the final term represents feedback energy from collapsed objects, taken to be proportional to the mass of formed stars, (different from once gas recycling is taken into account— see Appendix A). These additional energy inputs will be described in Sec. 3. Finally, requiring the final surface pressure to match the external pressure yields the second constraint

 (1+δrel)N−1rfNrf∑k=1P0θΓΓ−1k=Ps, (6)

summing over the bins in the radial range . (This is a change from BodeOWS07, which used a wider zone, thereby giving a lower measurement of surface pressure).

### 2.2 Creating the model cluster catalog

The gas prescription was applied to clusters selected from a large-scale N-body simulation. We will refer the reader to BodeOWS07 and SehgalTHB07 for further details, as the numerical methods used to carry out this simulation are the same as described in those papers. The cosmological parameters were chosen to be () = (0.044, 0.264, 0.736, 0.71, 0.96, 0.80), consistent with the WMAP 5-year results (KomatsuWMAP08). The simulation volume is Mpc on a side and contains N= particles, making the particle mass ; the particle spline softening length is kpc. At each simulation step particles were saved in one octant of a thin spherical shell centered on the origin, thus building up the matter distribution in a past light cone covering one eighth of the sky. All halos with a friends-of-friends mass of at least (293 particles) were selected from this light cone, out to redshift . Masses and radii at the virial and other overdensities are computed, in particular and at a mean overdensity of 500 times critical. The amount of star formation and iron production is determined based on the cluster redshift, as described in Appendix A. The gas prescription was then run on each halo; there are over 15,000 halos in this catalog with .

To test that resolution is sufficient, a smaller sample of halos was selected from a higher resolution simulation. This has the same cosmological parameters, but the volume is Mpc on a side, and the softening length is kpc; thus mass resolution is increased by a factor of 30.5, and spatial resolution by a factor of 5. No significant difference was seen when comparing results for lower mass halos () at the two resolutions; halos more massive than this are difficult to deal with at the higher resolution because of the large grid required.

Once the density, temperature, and metallicity in each cell are known, the X-ray luminosity in the 0.5–2 keV band is computed with the MEKAL model (MeweGO85). For the cluster temperature, the mixT code of Vikhlinin06 was used to compute the spectroscopic temperature in the projected annulus ; this is the definition of cluster temperature used throughout this paper.

There are a number of parameters which must be specified: , and . We will calibrate the model by comparing to X-ray data, in particular by matching the gas fraction as a function of . Two subsamples from are selected from the halo catalog to match the properties of the two data sets we use for the calibration (this data is shown as a function of in Figs. 1 and 2). The first data set is from VikhlininKFJMMVS06, which consists mostly of higher clusters; so, for comparison, we select all halos with in the light cone. There are 76 halos in this subsample, with temperatures above 4 keV for most parameter choices. For a given model, we compute a , summing for each observed cluster over all model points having a within 4- of the observed . The second observational set is of low redshift groups (keV) from SunVDJFV09; so, for the second subsample, we select all model halos with and . This subsample also contains 76 halos, most with keV. A is computed from this sample in an analogous manner as above. When searching for a best fit model we will minimize ; the factor of 0.5 ensures that the two terms contribute roughly equally near the minimum value.

## 3 The Role of Feedback and Star Formation

It is instructive to first consider the case of no star formation or energy input, i.e. ; we will refer to this as the “Zero” model. The resulting temperature and are shown as a function of mass as dotted lines in Fig. 4. VikhlininBEFHJKMNQV09 found the relation is well fit by a power law, , with and . Fig. 1 shows the fractional difference , where is the value obtained by plugging into the observed relation. This figure shows that at a given , the predicted temperature from this simple model is lower than observed, by over 10%. At the same time, it shows that the predicted gas fraction remains near the cosmic mean, and is independent of mass (quite different from what is observed). This results in a which is too high (by at higher masses) and an relation with a slope near the self-similar value , slightly shallower than observed. Thus we see the simplest set of assumptions yields an incorrect gas energy. There are three main processes which can alter this energy: transfer of energy from the DM through gravitational interactions, star formation, and feedback from astrophysical sources. In this section we will consider each of these in turn.

So far (as in BodeOWS07), it has been assumed that the gas originally had the same specific energy as the DM. However, this ignores an important process during cluster formation: when mergers occur, energy is transferred from the DM to the gas. In fact, the DM can lose up to 10% of its energy in an equal mass collision, and even a 10:1 merger can reduce the DM energy by a couple of percent (McCarthyBBVPTBLF07). Simulations show such mergers are common: CohnWhite05 calculate that a typical cluster will experience four 5:1 accretion events since , and FakhouriMa08 find that the mean rate for 10:1 or larger mergers in clusters is roughly two per halo per unit redshift. The magnitude of this effect can be estimated from hydrodynamical simulations, where it has long been found that baryon fraction in clusters is lower than the universal value (e.g. CenO94; EkeNF98; EttoriDBM06; GottloberY07; StanekRE09z). For simulations with no star formation or radiative processes, CrainEFJMNP07 find that the average baryon fraction within is 90% of the cosmic mean at =0. This result is independent of cluster mass, and at =1 the typical fraction is only 2 or 3% higher. Without star formation and energy feedback, the method of Sec. 2.1 gives a mean for halos with (at low redshift, ). Instead assuming that the energy transferred to the gas is 5% of the DM energy, i.e. using in Eqn. 5, this fraction is reduced to . As we are using a lower than CrainEFJMNP07, a slightly higher gas fraction at seems appropriate. Thus we adopt % for the rest of the paper.

Fig. 1, and also Fig. 2, display the change in inside after including this added energy, as short-dashed lines. While lower than before, it is still higher than observed, and still roughly independent of mass. The temperature at a given mass increases, but only slightly, so the relation is still too high and shallow. Including a term proportional to in this fashion does not change the self-similar nature of the model.

Star formation will obviously decrease the gas fraction, but it has a second effect as well: since it is the most bound, lowest entropy gas that is turned into stars, the remaining gas will have a higher specific energy than before (VoitBryan01). The simplest assumption is to make the efficiency of star formation independent of cluster mass; Fig. 2 shows, as dot-dashed lines, a model after also setting =9% for all halo masses. The predicted gas fraction drops further— for massive halos, near the observed value. Note the change in is more than 9%; the higher specific energy of the remaining gas is a significant effect, causing that gas to expand and further reducing the gas fraction in the inner parts of the halo. Temperatures also increase, so at higher masses the predicted values are close to those observed. However, as the method up to this point has been scale-free, the gas fraction still does not depend on mass; thus, for less massive clusters, the predicted and are still too large, since the slope is still the self-similar value.

There are two ways to change the scale-free nature of the models: vary star formation efficiency with halo mass, or include feedback energy proportional to stellar mass. There is observational support for both options. LinMS03, with a sample of 27 clusters, measured the fraction of stellar mass to be larger for lower mass clusters, following with . To apply this to our model, we will assume this relation also holds inside the virial radius, i.e. at . A model varying with mass in this manner is shown in Fig. 1 as dot-dashed lines. As one would expect, this decreases the gas fraction more in lower mass halos; it also changes the slope of the relation slightly.

The slopes of the and relations will clearly depend on how strongly varies with halo mass. Recently GonzalezZZ07 found a much steeper relation in a sample of 23 clusters: . A model where follows this relation is shown as long-dashed lines in Fig. 1. The effect is quite dramatic in low mass clusters because so much gas is converted into stars ( or more), and the resulting and gas fractions tend to overshoot the observed trends at low masses. If we start with the LinMS03 and vary the slope, searching for a best-fit model as described in Sec. 2.2, we find an exponent of yields a cluster sample with and relations quite close to those observed; this model is shown as solid lines in Fig. 1. We will refer to this choice of parameters () as our “SF-only” model.

The second effect which can break the scale-free nature of the models is additional energy injection, e.g. from SNe or AGN. We will assume such feedback is proportional to the number of stars formed, (see Appendix A). As an estimate of the amount of energy under consideration, for =9% and feedback energy , the energy added amounts to 0.8 keV per particle initially inside . Core collapse SNe could only add 0.4 keV/particle (assuming the model in Appendix A and erg per event), so in this case there must be some additional source such as AGN energy from accreting black holes. If the mass in black holes is , and 2% of the rest mass energy of the accreted material is converted into thermal energy (AllenDFTR06), then an additional 1.0 keV/particle could be added to the gas.

Fig. 2 shows the effect of adding feedback, starting with and . The long-dashed lines are for an extreme model with . For the most massive halos, the feedback energy is not significant compared to the binding energy of the cluster; hence there is little change evident. For less massive clusters, the added energy is enough to cause the gas to expand, lowering the gas fraction considerably. The relation becomes less steep as well, as lower mass halos become hotter; the relation is thus less steep also, though in this example the change is so large as to be incompatible with observations.

Varying the amount of feedback, the best-fit search yields a more moderate level of , shown as solid lines in Fig. 2. We will call this () the “Feedback” model. Like the SF-only model, this reproduces quite well the observed and relations. The efficiency is smaller than the equivalent number used in BodeOWS07 for several reasons. Part of the difference is merely bookkeeping: BodeOWS07 did not consider mass loss from stars, and expressed feedback as being proportional to , while here we instead use the larger mass . More substantially, BodeOWS07 did not include the term, instead effectively considering it as part of the feedback. Finally, here the gas surface pressure is calculated over a smaller radial buffer zone than in BodeOWS07, which leads to a slightly lower gas fraction initially.

Figs. 1 and 2 show that there is a degeneracy between steepening the star formation efficiency as a function of halo mass and increasing the amount of feedback per unit stellar mass. This is shown schematically in Fig. 3. We have found two models which both reproduce the X-ray data reasonably well: the SF-only model, which varies star formation efficiency as (steeper than the values found by LinMS03 or Giodiniea09z, but less steep than GonzalezZZ07), and the Feedback model, which instead has constant and adds feedback . A third option is of course some combination of these two. Since a stellar mass fraction independent of mass does appear to be contradicted by observations, we will adopt the star formation efficiency as given by LinMS03 and vary to find the best-fit model. This results in a lower level of feedback, . We will call this our “Standard” model (), which we will use for the rest of the paper. Fig. 4 shows the and relations for this third model. This choice of parameters yields a population with thermal energy as a function of mass quite close to that seen by VikhlininBEFHJKMNQV09 and SunVDJFV09. However, keep in mind that there will be a range of models around this one which will fit the X-ray data equally well– it would be difficult to distinguish between these three models based solely on the gas profiles. The shaded region in Fig. 3 gives a rough idea of the acceptable range of and .

There is one further physical effect to be considered, namely nonthermal pressure support, e.g. from cosmic rays or magnetic fields (e.g. PfrommerESJD07). We will conclude this section by considering how a moderate level of nonthermal support affects our model. If in fact , then the derived from X-rays are too low, because they are derived by assuming hydrostatic equilibrium. Fig. 5 shows our Standard model, after adjusting masses in the observational points by 10% to account for this; in this situation the model now has too high and at a given . Adding to this model means that the gas requires less thermal pressure to maintain hydrostatic equilibrium, so the solution to Eqns. 5 and 6 yields gas which is cooler and correspondingly denser. This is shown as dashed lines in Fig. 5; the model is now too cold, and the gas fraction even higher. Comparing to Fig. 1, it can be seen that increasing has the opposite effects of increasing . If we take the model and allow the dynamical energy input to vary, the best-fit result has , shown as solid lines in Fig. 5. This puts our predicted clusters back into agreement with the (adjusted) observations, at the same level as the Standard model. We will call this set of parameters () the “” model.

Increasing in this way assumes that non-radiative hydrodynamic simulations are missing some processes which transfer energy to gas during cluster formation. This certainly could be the case. For example, the dynamical friction experienced by the stars assembling into a cD galaxy would add significant amounts of energy to the ICM (ElZantKK04; Kim07; ConroyO08).

Having fixed on a Standard model, in this section we will examine how well it matches observations at other radii than .

Given that our DM halos have the standard density profile, the relation of our model should also agree with observations. Fitting all halos with keV, the best fit power law relation gives a normalization of at =5keV. This agrees well with observed X-ray samples, using masses derived from either weak lensing (Hoekstra07) or X-rays (ArnaudPP05; VikhlininKFJMMVS06).

The Standard model also has an X-ray luminosity in good agreement with the observations of VikhlininBEFHJKMNQV09, unlike the Zero model which is a factor of 2-3 more luminous. This is hardly surprising, since we calibrated the model in part on this data. Above , the predicted , but below this the relation is much steeper, approaching . VikhlininBEFHJKMNQV09 find the observed scatter at a given mass is , while the model scatter is roughly half of this. There are many processes occurring in the cores of clusters which could significantly alter not included in our model, which would increase the scatter and somewhat alter the mass dependence.

A fundamental property of the gas is its entropy, usually defined in cluster studies as , where is the electron number density. Based on Chandra observations of nearby relaxed clusters, NagaiKV07 determined how entropy varies with cluster temperature at several characteristic radii, including , , and . Their best-fit power-law approximations of these data are shown in Fig. 6 as dashed lines, which shows the entropy at these three radii as a function of . The Zero model gives too low entropy values at all radii; true clusters have shallower entropy profiles than the scale-free model predicts. Once star formation and feedback are included, however, there is much improved agreement with observed clusters; the main disagreement is at , where the model is roughly 10% too high. SunVDJFV09 also determined relations at these radii, using Chandra archival data and including many lower temperature groups. For keV, the median value from the model clusters is within 10% of their power-law fits at and . At the model is lower than the SunVDJFV09 relation by roughly 15%, for all but the most massive clusters. Our model does not consider the history of cluster formation, nor the fact that the nonadiabatic processes which can alter the entropy will likely be more pronounced in the core.

AfshordiLNS07 assembled a heterogeneous catalog of 193 X-ray clusters and then examined WMAP data to derive a mean pressure profile using the SZ effect. This profile is shown as points with error bars in Fig. 7; the pressure is normalized to the critical density times the global cluster . For comparison, we computed the mean profile of the 284 clusters from our Standard sample with keV. This is shown as solid line in Fig. 7; the mean profile of the same halos without star formation or energy input is shown as a dotted line. To scale the radius, is estimated from (as in Eqn. 1 of AfshordiLNS07), though using the radius measured directly from the DM distribution gives the same results. Outside of , changes to the gas energy show little effect on the mean profile. At smaller radii, star formation and feedback reduce the density more than they affect the temperature profile, thus lowering the pressure profile in the core. This gives better agreement with the AfshordiLNS07 profile, though it is still higher. One source of this disagreement is that they assumed inside when constraining their profile, while our model is shallower than this; had they assumed a similarly shallow inner profile, this point would have been higher. Note also the redshift and temperature distributions of the model and observed samples are quite different.

BonamenteJLCNM08 present a joint SZ and X-ray analysis of 38 clusters; here we will compare to their low-redshift sample of 19 clusters in the range . For our Standard sample with we calculated the integrated Compton y-parameter, (Eqn. 6 of BonamenteJLCNM08), out to projected ; this is shown in Fig. 8. Most of our sample is at lower gas masses than the clusters used by BonamenteJLCNM08, but where the two samples overlap the agreement is good. The best fit power-law found by BonamenteJLCNM08, with slope 1.60, is shown as a dashed line in Fig. 8. A power-law fit to our model halos with yields a slope of 1.61. Also shown as a dotted line is the Zero model. The effect of increasing the gas energy is to cause gas to expand and reduce the gas mass; however, this does not change the integrated pressure profile much, since the gas in the center must still hold up all the gas above it. This effect is larger for smaller mass halos, making the relation less steep.

In Sec. 3 we found several different models which, given the degeneracy among input parameters, gave similar results. Since one goal of this is to provide a means of predicting the SZ signal from clusters, it is useful to compare these models. Fig. 9 shows the fractional difference in the value from the Standard model, as a function of , for the SF-only, Feedback, and models. The median SZ signal is within 10% of the Standard model in each case, and they are all roughly a factor of 2 higher than the Zero model, which does not include star formation and energy injection. This highlights the importance of normalizing to X-ray data when modeling the SZ signal.

## 5 Summary and Discussion

In this paper we have described an improved method for calculating the gas distribution in a DM potential well, constraining the model to match the observed X-ray gas fractions; then we made further comparisons with X-ray and SZ observations in order to check, where possible, our conclusions. Extracting the greatest return from upcoming cluster surveys will require understanding the state of the intracluster gas. It is encouraging that the simple model presented here can so well reproduce observed profiles and trends with cluster mass. With this model the following main points have demonstrated:
Simulations show that during mergers the energy of the DM component will be tapped to provide shock heating to the gas. Transferring % of the binding energy of the DM to the gas results in gas fractions at large radii which are in good agreement with non-radiative hydrodynamical simulations.
Breaking the scale-free nature of the model to bring it into accord with observed scalings requires either varying the efficiency of star formation with cluster mass or increasing the specific energy of the gas by an amount that is not proportional to the cluster binding energy (e.g. by being proportional to the stellar mass), or both.
Without including feedback energy, varying the efficiency of star formation as reproduces observed and relations. This is steeper than found by LinMS03 and Giodiniea09z, but less steep than found by GonzalezZZ07. A relation as steep as , as seen by the latter, would leave too small gas fractions (and that gas at too high temperatures) in less massive clusters; however, the difference is less than two sigma. Also, keep in mind that what is measured is the stellar fraction inside , whereas in our model is the fraction inside ; the two need not be identical. We are assuming that stars have the same radial distribution as DM, but if stars are instead more centrally located, this would change our predicted relation. GonzalezZZ07 selected for clusters with a dominant central galaxy, and for the least massive clusters is in the outskirts of this galaxy. Because a central dominant galaxy takes up a larger fraction of the stellar mass in smaller clusters, one could argue that the stellar fraction inside changes more rapidly with cluster mass than does the fraction inside , i.e. should vary with halo mass at a larger rate than .
A model with independent of cluster mass but including feedback energy proportional to the mass of formed stars can also reproduce the observed and . However, this requires adding roughly one keV per particle initially inside . This might be possible if the output from accreting black holes is converted to thermal energy with high enough efficiency ().
The most likely possibility is that both star formation and feedback are at work. Assuming (as in LinMS03) and including a lower feedback level of also yields the proper scalings. Assuming the mass in black holes is , this energy is of their rest mass (or 0.33 keV per particle; note this does not take into account the feedback supplied by supernovae). This model was compared to further observations in Sec. 4, but the other two possibilities discussed here behave in quite a similar fashion.
Including nonthermal pressure results in a cooler, denser gas distribution. To bring such a model back into agreement with X-ray observations requires increased energy input. Once this is done, the pressure or entropy for a given is quite similar to models with , although of course the relation is different.

A couple caveats are in order. We start with halos drawn from a DM-only simulation; adding baryons to such a simulation may change the properties of the halo population (StanekRE09z). Also, we have neglected many pertinent physical processes, such as cooling in the core (note our definition of cluster temperature excludes the core). Because of this, the model halo population likely has too low a scatter at a given mass. This is exacerbated by the tight relation we assume between halo mass and star formation plus feedback. For example, if we allowed for a variation in at a given halo mass, clusters with higher would have higher temperatures and lower gas fractions, and thus lower X-ray luminosities. This would explain why catalogs of optically selected clusters tend to be underluminous in X-rays at the same mass (RykoffMBEJKRSW08).

For most halos in our standard model, the gas initially assumed to be inside has expanded outwards. For hotter (keV) clusters only, the mean final radius is (one sigma); for all clusters the mean is . A result of this is that for the most massive clusters the baryon fraction inside is roughly 90% of the cosmic mean, . Fig. 10 shows the the baryon fraction inside as a function of cluster mass for differing input parameters. The models shown have similar gas fractions (by design), but varying stellar mass fractions. This results in differing total baryon fractions at lower masses (though, as pointed out above, this may change if stars are more centrally concentrated than DM). In every case, however, the baryon fraction inside tops out near at the highest cluster masses. Our sample contains four halos with virial masses above ; their gas fractions within vary from to , with a mean of . This agrees well with the detailed multiwavelength analysis of four massive clusters by Umetsuea09z. GonzalezZZ07 also found baryon fractions inside near at high masses. As this work was nearing completion, Giodiniea09z published a new determination of the baryon fraction (not including intra-cluster light) as a function of cluster mass. Our Standard model is in good agreement with their result, shown as a dotted line in Fig. 10.

The “missing” gas in this model is distributed in the outskirts of clusters. There it would be difficult to observe directly, being at low temperatures and densities. However, Prokhorov08 has predicted in this case there would be emission in the extreme ultraviolet from these outer shells. Alternatively, as this gas is associated with clusters, its kinetic SZ signature should be correlated with the cluster velocity field (HoDS09z; HMH09z).

Computer simulations and analysis were supported by the National Science Foundation through TeraGrid resources provided by Pittsburgh Supercomputing Center and the National Center for Supercomputing Applications under grant AST070015; computations were also performed at the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology. PB was partially supported by NSF Grant 0707731.

## Appendix A Star formation, gas recycling, and metal formation

We will follow the “fossil” model of NagamineOFC06, which matches local stellar populations; this represents a bulge and a disk component with delayed exponentials. Let denote the mass in stars of type formed by time , and the total mass of such stars ever formed by a time corresponding to redshift zero, i.e. today. (At this point no allowance is made for mass lost from stars.) For the time evolution of this stellar mass we will adopt

 (A1)

This choice means the star formation rate is given by a delayed exponential

 ˙MF(τj,t)=Fjχ∗(τj,tH)˙χ∗(τj,t)=Fjχ∗(τj,tH)1τjtτje−t/τj. (A2)

NagamineOFC06 found a combination of a bulge population with decay time =1.5 Gyr and a disk population with =4.5 Gyr is most consistent with local stellar populations (hence the term fossil). Since stellar populations in clusters are older than average, we will assume 90% of the total stars formed by are from the bulge component, and the remaining 10% from the younger disk component. The resulting star formation history is shown in Fig. 11; there is little new star formation for 1. This figure assumes =0.1 at =0; note that if increases then the amount of gas inside will decrease, changing the final value of this ratio.

The supernova rate also depends upon the delay time distribution. For a burst of star formation at =0, suppose the fraction of SNe of type occurring by time is given by . Then the number of SNe per unit stellar mass per unit time will follow the commonly used distribution (e.g. Strolger+04; MannucciDVP06)

 Φ(τi,t)=νi˙fSN(τi)=νiτie−t/τi, (A3)

where is the total number of SNe of type per unit mass which will ever occur. The SN rate , i.e. the number of SNe per unit time, is then given by

 ˙S(τi,τj,t)=∫t0Φ(τi,t−s)dMF(τj,s)dsds=Fjχ∗(τj,tH)νiμ2ij[˙fSN(τi,t)+˙χ∗(τj,t)τjτi(1μij−τjt)], (A4)

where ; we will be assuming and throughout. Finally, the cumulative number of SNe of type which have occurred by time is

 S(τi,τj,t)=∫t0dS(τi,τj,s)dsds=Fjχ∗(τj,tH)νiχ(τi,τj,t), (A5) χ(τi,τj,t)=μ2ij[fSN(τi,t)+1μijτjτiχ∗(τj,t)+τjτi(e−t/τj−1)]. (A6)

Thus it only remains to specify and for each type of SN.

For core collapse SNe, the delay time is short; we will adopt Gyr. We will employ a simple IMF with a slope of -0.5 for stellar masses in the range and Salpeter slope -1.35 for (BaldryG03). Assuming all stars in the range 8-50 end their lives as SNCC, then .

The picture for SNIa is more complicated. There is growing evidence for a prompt SNIa component with a short delay time after star formation (e.g. DellaVallePPCMT05; MannucciDVPCCMPT05; MannucciDVP06; Neill+06; Sullivan+06; AubourgJHSS08, and references therein). On the other hand, much observational evidence points towards much longer delay times of 2-4 Gyr (e.g. GallagherGBCJK05; MannucciDVPCCMPT05; Sullivan+06; ForsterSchawinski08z, and references therein). ScannapiecoBildsten05 proposed that two evolutionary channels exist, a prompt population that would be proportional to the recent star formation rate, and a tardy one with a much longer delay time that would be proportional to the total stellar mass. MannucciDVP06 found that the data are best matched by equal amounts of SNIa in two populations, one exploding of order yr after star formation and the other following an exponential distribution with a decay time of about 3 Gyr. We will follow this model for SNIa, setting a prompt delay time Gyr for one half, and a tardy Gyr for the other half.

Assuming SNIa result from 3-8 stars, then , where is the fraction of the stars in this mass range that go supernova. This fraction is not well constrained; Maoz08 examines a number of different estimates of this fraction and finds that they vary from 2% to 40%, but also that is consistent with all the estimates he considered. Thus we will set , or , such that there are equal numbers of prompt and tardy SNIa.

In a cosmology, this model gives predicted SN rates per unit mass at of 0.07 SNuM for both SNCC and SNIa. For comparison, MannucciMSBDGP08 find and SNuM ( statistical errors) for SNCC and SNIa, respectively. The predicted rates increase with redshift; the predicted SNIa rate is well within the (statistical) errors of the rate measured by SharonGMFG07 in the redshift range 0.06–0.19, and of the rate measured by Graham+08 in the range 0.2–1. From the supernova rate one can predict metallicity evolution. Using iron yields per event of 0.077 for SNCC and 0.749 for SNIa, the resulting history of iron production is also shown in Fig. 11. Even though there is little star formation at low redshift, iron production continues because of tardy SNIa. Note the predicted metallicity will be proportional to the amount of star formation; also note that we are assuming spatial homogeneity. Fig. 11 also shows the predicted cumulative ratio of SNIa to SNCC; there are roughly 2.5 SNCC for each SNIa, and SNCC produce 20% of the iron.

The rate of change in the actual stellar mass is the difference between the star formation rate and the rate at which mass is lost from stars. Codes exist for tracking the mass loss from a stellar population, e.g. the PEGASE.2 code of FiocRV99z. We will instead assume that all mass loss, whether from SNe or winds, can be modeled with Eqn. A4. Then the stellar mass follows

 ˙M∗(τj,t)=˙MF(τj,t)−∑iyi˙S(τi,τj,t), (A7)

where the mass yield is the mean initial mass of stars of type minus the mean remnant mass, if any. The delay times are adjusted to match the mass loss calculated by the PEGASE.2 code111Available at http://www2.iap.fr/users/fioc/PEGASE.html The choices of delay times and yields are summarized in Table 1. We assume SNCC leave behind neutron stars of mass 1.4 and SNIa leave no remnant. Stars larger than 50 are all assumed to collapse to 2 black holes, on the same time scale as SNCC. Stars in the range which do not go supernova instead remain as 0.55 white dwarfs. The delay time for mass loss from these stars is set to 0.05 Gyr, slower than the prompt SNIa component; this choice is required to match the mass loss on short ( Gyr) time scales. Stars in the range also become white dwarfs, with mass loss on a much longer time scale of 6 Gyr; the lower limit of is treated as a free parameter, set to match the mass in white dwarfs at long time scales. Given these parameters, the evolution of the recycled gas fraction and the mass in remnants resulting from a single starburst is shown in Fig. 12; also shown for the same IMF are the predictions of the PEGASE.2 code. Our relatively simple scheme follows the more detailed evolutionary code reasonably well.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters