Results from a large set of hydrodynamical smoothed particle hydrodynamics (SPH) simulations of galaxy clusters in a flat ΛCDM cosmology are used to investigate the metal enrichment and heating of the intracluster medium (ICM). The physical modelling of the gas includes radiative cooling, star formation, energy feedback and metal enrichment which follow from the explosions of supernovae of type II and Ia. The metallicity dependence of the cooling function is also taken into account. The gas is metal-enriched from star particles according to the SPH prescriptions. The simulations have been performed to study the dependence of final metal abundances and heating of the ICM on the numerical resolution and the model parameters.
For a fiducial set of model prescriptions the results indicate radial iron profiles in broad agreement with observations; global iron abundances are also consistent with data. It is found that the iron distribution in the intracluster medium is critically dependent on the shape of the metal deposition profile. At large radii the radial iron abundance profiles in the simulations are steeper than those in the data, suggesting a dynamical evolution of simulated clusters different from those observed. For low-temperature clusters simulations yield iron abundances below the allowed observational range, unless a minimum diffusion length of metals in the ICM is introduced. The simulated emission-weighted radial temperature profiles are in good agreement with data for cooling flow clusters, but at very small distances from the cluster centres (∼2 per cent of the virial radii) the temperatures are a factor of ∼2 higher than the measured spectral values.
The luminosity—temperature relation is in excellent agreement with the data; cool clusters (TX ∼ 1 keV) have a core excess entropy of ∼200 keV cm2 and their X-ray properties are unaffected by the amount of feedback energy that has heated the ICM. The findings support the model proposed recently by Bryan, where the cluster X-ray properties are determined by radiative cooling. The fraction of hot gas fg at the virial radius increases with TX, and the distribution obtained from the simulated cluster sample is consistent with the observational ranges.
Galaxy clusters are the largest virialized structures known in the Universe and are considered useful probes to constrain current cosmological theories of structure formation. X-ray observations of galaxy clusters show that most of the baryonic cluster mass is in the form of a hot ionized intracluster medium (ICM) at temperatures ≃107–108 K, with the bulk of the emission in the X-ray band from bremsstrahlung processes (Sarazin 1986). The dependence of the X-ray emission on the square of the gas density allows us to construct cluster samples without the biases that may arise in the optical band. The final physical state of the ICM is primarily determined by the gravitational processes, which have driven the dynamical evolution of the gas and dark matter mass components during the cluster collapse. Under the assumption of hydrostatic equilibrium, the ICM gas distribution can be modelled to connect the gas temperature TX to the cluster virial mass or to the bolometric X-ray luminosity LX. The X-ray temperature and luminosity functions are then predicted for a given cosmological model according to the standard theoretical Press—Schechter (Press & Schechter 1974) mass function. The observed evolutionary history of these functions (Edge et al. 1990; Henry & Arnoud 1991; Henry 1997) can then be used to put severe constraints on the allowed cosmological models (Henry & Arnoud 1991; White, Efstathiou & Frenk 1993; Eke, Cole & Frenk 1996; Bahcall & Fan 1998; Kitayama, Sasaki & Suto 1998).
If the ICM has been shock-heated solely by gravitational processes, the cluster scaling relations are predicted to obey a self-similar behaviour. For instance, the LX−TX relation should scale as LX∝TX2 (Kaiser 1986). For clusters with TX≳ 2 keV there is a wide range of observational evidence (David et al. 1993; Allen & Fabian 1998; Markevitch 1998) that the observed bolometric X-ray luminosity scales with temperature with a slope steeper than expected (LX∝TX3). This implies that low-temperature clusters have central densities lower than those predicted by the self-similar scaling relations (Bower 1997; Ponman, Cannon & Navarro 1999; Lloyd-Davies, Ponman & Cannon 2000). This break of self-similarity is usually taken as strong evidence that non-gravitational heating of the ICM has played an important role in the ICM evolution (Evrard & Henry 1991; Kaiser 1991; White 1991), at least for cool clusters.
The energy source that has been most widely considered as a heating mechanism for the ICM is supernova (SN) driven winds, which inject energy into the ICM through SN explosions of type Ia and II (White 1991; Loewenstein & Mushotzky 1996). The energy input required to heat the gas at the entropy level necessary to reproduce the observed departure from self-similarity is estimated to lie in the range ≃0.5 –3 keV per particle (Balogh, Babul & Patton 1999; Wu, Fabian & Nulsen 2000; Tozzi & Norman 2001). Several authors have argued that it is unlikely that SNe can provide the required energy to heat the gas, and have suggested active galactic nuclei (Valageas & Silk 1999; Wu et al. 2000) as the extra energy source for ICM heating. Bryan (2000) has proposed the alternative view that radiative cooling and the subsequent galaxy formation (Pearce et al. 2000) can explain the observed LX−TX relation because of the removal of low-entropy gas at the cluster cores.
The estimates of the energy injected from SNe into the ICM are based on the measured abundances of metals in the ICM. The elements that have enriched the ICM have been synthesized in the SN explosions of the stellar population of the cluster. Two enrichment mechanisms have been proposed: SN-driven galactic winds from early-type galaxies (Mushotzky & Lowenstein 1997), and ram pressure stripping of the enriched gas from galaxies (Gunn & Gott 1972). Another possibility is that of a significant contribution to the metal enrichment and heating of the ICM, associated with an early-epoch generation of massive Population III stars (Loewenstein 2001).
A large set of observations confirms that the ICM of galaxy clusters is rich in metals (Arnaud et al. 1992; Mushotzky & Lowenstein 1997; Fukazawa et al. 1998; Dupke & White 2000a; Finoguenov, David & Ponman 2000; Fukazawa et al. 2000; Matsumoto et al. 2000). Analyses of X-ray spectra show that the abundance of heavy elements in the ICM is nearly ≃1/3 solar. These measurements provide strong support for the SN scenario as a heating source for the ICM.
Abundance gradients have also been measured (Ezawa et al. 1997; Dupke & White 2000b; White 2000; Finoguenov et al. 2000; De Grandi & Molendi 2001; Irwin & Bregman 2001). For instance, Ezawa et al. (1997) found a decline in the iron abundance of AWM7 from ≃0.5 solar in the centre to ≃0.2 solar at a distance of ≃500 kpc. Measurements of the relative abundance of the heavy elements can be used to constrain the enrichment mechanism of the ICM and the energy input from SNe. Analysis of the spatial distribution of metallicity gradients is also important to discriminate between the proposed enrichment scenarios (Dupke & White 2000b). From the observed metallicities one can estimate the SN energy released given a shape for the initial mass function (IMF) and a nucleosynthesis yield model for the SN explosions. Finoguenov, Arnaud & David (2001) have obtained Si and Fe abundances using the X-ray data of a selected sample of X-ray clusters. From the measured abundances, they find a significant contribution to the energy per particle associated with SN explosions at temperatures ≃3 keV.
These kinds of estimates of the SN energy input suffer from theoretical uncertainties in the yield models and in the form of the IMF (Gibson, Loewenstein & Mushotzky 1997). Moreover, these estimates can vary widely with the assumed spatial distribution of the ICM metals and the transfer efficiency of the kinetic energy released in an SN explosion to the ambient gas (Kravtsov & Yepes 2000). Hydrodynamical simulations of cluster evolution have the advantage over analytical methods that they take into account the dynamical evolution of the gas. From the results of numerical simulations of galaxy formation, Kravtsov & Yepes (2000) conclude that it is unlikely that SNe can provide the required energy input, even assuming the existence of radial metallicity gradients in the ICM.
In order to implement self-consistently the metal enrichment of the ICM in hydrodynamical simulations of galaxy clusters, one first must consider the effects of non-gravitational processes on the cluster gas distribution. The effects of radiative cooling in simulated clusters have been investigated by a number of authors (Anninos & Norman 1996; Lewis et al. 2000; Pearce et al. 2000; Yoshikawa, Jing & Suto 2000; Valdarnini 2002). One of the main conclusions of these simulations is that the modelling of radiative processes for the gas cannot be decoupled from a prescription for turning cold, dense gas into stars. This is done in order to avoid unphysically high densities (cooling catastrophe). The luminosities LX of the simulated clusters are found to be physically plausible, provided that a suitable prescription for the treatment of the cold gas has been taken into account.
In a previous paper (Valdarnini 2002), results from a set of hydrodynamical smoothed particle hydrodynamics (SPH) simulations of galaxy clusters were used to investigate how final X-ray properties of the simulated clusters depend upon the numerical resolution of the simulation and the chosen star formation (SF) prescription. For a chosen SF model, final X-ray luminosities were found to be numerically stable and consistent with data, with uncertainties of a factor of ≃2. In these simulations the metal enrichment of the ICM was included with a minimal number of prescriptions and it was shown that, for the simulated clusters, the final iron abundances of the ICM are in broad agreement with measured values.
Chemical evolution in hydrodynamical SPH simulations has already been considered in a variety of contexts (Steinmetz & Müller 1994; Raiteri, Villata & Navarro 1996; Carraro, Lia & Chiosi 1998; Buonomo et al. 2000; Mosconi et al. 2001; Churches, Nelson & Edmunds 2001; Aguirre et al. 2001; Lia, Portinari & Carraro 2002). Metzler & Evrard (1994) have investigated the metal enrichment of the ICM in particle—particle—mesh smoothed particle hydrodynamics (P3MSPH) simulations of galaxy clusters in a standard cold dark matter (CDM) scenario. The authors have not included radiative cooling in the simulations and have adopted a phenomenological prescription to model the chemical enrichment.
In a more recent paper, Kravtsov & Yepes (2000) have estimated SN heating of the ICM using fixed-grid Eulerian hydrodynamical simulations. However, the numerical resolution of the simulations was not adequate to study the evolution of simulated clusters. The number of SNe occurring in a given cluster was estimated statistically from many small-box galaxy formation simulations. The implementation of a self-consistent metal enrichment model for the ICM in hydrodynamical simulations is therefore important for investigations of the ICM metal evolution. This is particularly relevant in connection with the observed metallicity gradients and for assessing the reliability of the calculated amount of ICM heating from SNe, inferred from measured metallicity abundances.
The main aim of this paper is to analyse, in hydrodynamical SPH simulations of galaxy clusters, the dependence of the final iron abundance on a number of model parameters that control the ICM metal enrichment. This is done in order to obtain for the simulated clusters a final ICM distribution that can consistently fit a set of observational constraints, such as the observed iron abundances and at the same time the luminosity—temperature relation. Implications for the ICM heating from SNe are also discussed. This paper constitutes a generalization of a previous work by Valdarnini (2002, hereafter V02), where the investigation was mainly concerned with the numerical stability of different SF models. Section 2 presents the hydrodynamical simulations of galaxy clusters that have been performed. Section 3 describes the SF algorithm and the modelling of ICM metal enrichment implemented in the simulations. The results are presented in Section 4 and the conclusions in Section 5.
Here a short description of the simulations is given. Further details can be found in V02.
The cosmological model considered is a flat CDM model, with a vacuum energy density ΩΛ = 0.7, matter density parameter Ωm = 0.3 and Hubble constant h = 0.7 in units of 100 km s−1 Mpc−1. The primeval spectral index of the power spectrum n is set to 1 and Ωb = 0.015 h−2 is the value of the baryonic density. The power spectrum of the density fluctuations has been normalized in order to match at the present epoch the measured cluster number density (Eke et al. 1996; Girardi et al. 1998). Initial conditions for the cluster simulations are constructed as follows. A collisionless cosmological N-body simulation is first run in an L = 200 h−1 Mpc comoving box using a P3M code with 843 particles, starting from an initial redshift zin. At z = 0 clusters of galaxies are located using a friends-of-friends algorithm, so to detect densities ≃200 Ω−0.6m times the background density within a radius r200. The corresponding mass M200 contained within this radius is defined as M200 = (4 π/3) ΩmρcΔcr3200, where Δc = 187 Ωc−0.55 for a flat cosmology and ρc is the critical density. The 40 most massive clusters within the simulation box are identified according to this procedure and the most massive and least massive cluster (labels 00 and 39, respectively) of this sample are selected for the hydrodynamical simulations. This procedure has been already followed in V02 and the initial conditions of the cosmological simulation are the same. In addition to these two clusters, two other clusters are considered. The original sample is enlarged to include 80 more clusters of decreasing virial mass. This amounts to a total of 120 clusters. The additional clusters extracted from the new sample are the cluster highest in mass after L39 (L40), and the least massive of the sample (L119). Table 1 lists the properties of the four selected clusters.
For each of the four test clusters, hydrodynamical TREESPH simulations are performed in physical coordinates. The initial conditions of the hydrodynamical simulations are determined as follows. The cluster particles at z = 0 within r200 of the test cluster are located in the original simulation box back at zin. A cube of size Lc≃ 15–25 h−1 Mpc ∝M1/3200 enclosing these particles is placed at the cluster centre. A lattice of grid points is set in the cube. At each grid point are associated a dark matter particle and a gas particle, of corresponding mass and coordinates. The particle positions are then perturbed, using the same initial conditions of the cosmological simulations. Finally, the particles for which the perturbed positions lie inside a sphere of radius Lc/2 from the cube centre are kept for the hydrodynamical simulations.
To model the effects of the external gravitational fields, the inner sphere is surrounded out to a radius Lc by dark matter particles with a mass eight times larger than the sum of the masses of a gas particle and a dark matter particle of the inner sphere. For each particle, gravitational softening parameters are set according to the scaling εi∝m1/3i. The numerical parameters of the hydrodynamical simulations are given in Table 2.
The simulations with index L00 and L39 have a number of gas particles Ng≃ 22 600. For this mass resolution, the corresponding runs in V02 have been found to give fairly stable final X-ray luminosities for the simulated clusters. Additional runs have been considered, with an increased resolution with respect to the standard runs. The high-resolution (H) runs have Ng≃ 70 000 and the very high-resolution (VH) runs have Ng≃ 210 000. The other numerical parameters have their values scaled accordingly.
The gravitational forces of the hydrodynamical simulations are computed using a hierarchical tree method with a tolerance parameter θ = 1 and taking into account quadrupole corrections. The hydrodynamical variables of the gas are followed in time according to the SPH Lagrangian method (Hernquist & Katz 1989, and references cited therein). In SPH, local fluid variables are estimated from the particle distribution by smoothing over a number of neighbours. A common choice is the B2 -spline smoothing kernel Ws(r, h) (Monaghan & Lattanzio 1985), which has compact support and is zero for interparticle distances |r|⩾ 2h. The smoothing kernel is normalized according to ∫W(r, h) d r = 1. The smoothing length h fixes the spatial resolution of the simulation. The resolution is greatly increased in high-density regions when individual smoothing lengths hi are allowed, so that the number of neighbours of a gas particle is nearly constant (≃32). A lower limit to the smoothing lengths hi is set by hi⩾hmin≡εg/4, where εg is the gravitational softening parameter of the gas particles. The time integration is done allowing each particle its own time-step. The accuracy of the time integration is controlled by a number of constraints that the individual time-steps must satisfy (see, e.g. Valdarnini, Ghizzardi & Bonometto 1999). The minimum allowed time-step for the gas particles is 6.9 × 105 yr.
The thermal energy equation for the gas particles includes a term that models the radiative processes of an optically thin plasma in ionization equilibrium. The total cooling function Λc depends on the gas temperature and metallicity. The cooling function takes into account contributions from recombination and collisional excitation, bremsstrahlung and inverse Compton cooling. Heating from an ionizing ultraviolet background has not been considered. The cooling rate of the gas in the simulations is then dependent on the gas metallicity. This is an important difference with respect to the previous simulations (V02) and is essential in order to analyse low-temperature (TX≃ 2 keV) clusters consistently. The dependence on the metallicity indeed has larger effects, as it increases the cooling rate. The small scales resolved by the simulations will cool faster than larger ones. The inclusion of cooling with its metallicity dependence has then a strong impact also on the formation of larger clusters because of the hierarchical growth of structure.
Tables of the cooling rates as a function of the temperature and gas metallicities have been constructed from Sutherland & Dopita (1993) and stored in a file. During the simulations a cubic spline interpolation is then used to calculate from the tabulated values the cooling function Λc(T, Z) of a gas particle of given temperature T and metallicity Z. Here Z is the mass fraction of metals of the gas particle. Conversion from the metallicity Z to the corresponding value of the iron-to hydrogen ratio is done as in Sutherland & Dopita. A good fit to this relation is given by 2 of Tantalo, Chiosi & Bressan (1998). The X-ray luminosities are computed from the gas emissivities within the cluster virial radius according to the standard SPH estimator [see 2 of V02]. The X-ray emissivity ε(r) associated with a gas particle is calculated with a Raymond—Smith code (Raymond & Smith 1977) as a function of the gas temperature and metallicity.
3 Star Formation and Icm Enrichment
Cold gas in high-density regions will be thermally unstable and subject to SF. In SPH simulations, SF processes have been implemented using a variety of algorithms. Here conversion of cold gas particles into stars is performed according to Navarro & White (1993). In V02 it has been found that for this SF method final profiles of the simulated clusters are robust against the numerical resolution of the simulation. According to Navarro & White, any gas particle in a convergent flow and for which the gas density exceed a threshold,Navarro & White (1993) defines the gas density threshold ρc,g now reads 2 must be solved numerically. A plot of ρc,g as a function of Z shows that ρc,g(Z) is nearly flat for Z≲ 10−3 and has a fast decay above Z ≃ 5 × 10−3.
Once a star particle i is created at the time ts(i) it will release energy into the surrounding gas through SN explosions. SNe of type II (SNe II) originate from the explosions of stars of mass mu ⩾ m ⩾ 8 M⊙ at the end of their lifetime τ(m); here τ(m) is defined as in Navarro & White (1993).
Each SN explosion produces εSNII≡εSN ≃ 1051 erg, which is added to the thermal energy of the gas, and leaves a ≃1.4 M⊙ remnant. The number of SN II explosions associated with the star particle i in the time interval [t − Δt, t] is determined as1 The normalization of the IMF is set to ; with this normalization Npop = ms(i) is the number of stellar populations of the star particle i. Several forms of the IMF have been considered. A Miller—Scalo (Miller & Scalo 1979) IMF has been chosen for a consistent comparison with the simulations of V02. For this IMF mu = 100. A standard IMF is the one of Salpeter (1955), where ϕ(m) ∝m−(1+x), with x = 1.35. Finally, a less steep IMF is given by Arimoto & Yoshii (1987) for elliptical galaxies, for which x≃ 1. The latter two IMFs have mu = 40. The indices of the simulations with different IMFs are presented in Table 3.
The energy produced in the time interval [t−Δt, t] by SN explosions of type II associated with the star particle i is εSNΔNII. This feedback energy is returned entirely to the nearest neighbour gas particles of the star particle i. The velocity field of the neighbouring gas particles is left unperturbed by the SN explosion, since the typical SPH simulation resolutions are much larger than the size of the shell expansion (Carraro et al. 1998). The energy is smoothed among the internal energies of the gas particles according to the SPH smoothing prescription. The internal energy increment Δuj of the gas particle j is thenSection 4.2, where the dependence of the simulated cluster profiles on a number of parameters is analysed.
The number of SNe of type Ia (SNe Ia) associated with the star particle i in the time interval [t − Δt, t] has been determined according to Greggio & Renzini (1983). The scheme implemented follows Lia et al. (2002). The number of SN Ia events at epoch t isPortinari, Chiosi & Bressan 1998). Inverting the order of integrations, one obtains Woosley & Weaver 1986). This energy is smoothed as in 2 over the nearest-neighbour gas particles of the star particle i.
SN explosions also inject enriched material into the ICM, thus increasing its metallicity with time. The mass of the kth heavy element produced in an SN explosion is defined as the stellar yield yS,k(m), with S = II or Ia. The total yield is the sum over the masses of the heavy elements: yS(m) = ∑ kyS,k(m). For SNe Ia the yield yIa,k is a constant which is independent of the progenitor mass. The adopted SN Ia yields are those of Iwamoto et al. (1999, W7 model), with yIa = 1.4 and yIa,Fe = 0.7. The yields yII,k of type II SNe depend on the progenitor mass m, and it is useful to define an average yieldGibson et al. (1997). The theoretical model chosen here is model B of Woosley & Weaver (1995). For this model the explosion energy of massive stars is boosted and the produced metals are not reabsorbed in the SN explosion. The yields of this model are therefore the most favourable from the point of view of the amount of iron synthesized from SNe II. Final metal abundances for a different theoretical SN explosion model can be obtained from the SN II component of the gas metallicity profile by a rescaling of the yields. This is a valid approximation as long as there are small differences between the yields yII which determine the gas metallicity Z, and therefore the local cooling rate and the SF threshold (2).
Woosley & Weaver (1995) yields are available down to ml = 11. The yields between ml = 8 and m = 11 have been obtained with a linear extrapolation; the associated errors are, however, small, since SN nucleosynthesis is negligible below m = 11. From table 1 of Gibson et al. (1997), 〈yII,Fe〉≃ 0.14 for the Woosley & Weaver B model with a Salpeter IMF, ml = 10 and mu = 50. I obtain the same value for an Arimoto—Yoshi IMF (Arimoto & Yoshii 1987) with ml = 8 and mu = 40. The mass of heavy elements that has been injected into the ICM at age t as a result of the SN II explosions associated with the star particle i is determined as (Tinsley 1980)2) for returning the SN feedback energy to the gas, i.e.
Since the choice of the smoothing kernel WZ is somewhat arbitrary, a standard choice is the SPH smoothing kernel (Mosconi et al. 2001; Lia et al. 2002). However, if the ejected material is deposited at the end of the shock expansion, a less steep deposition profile is a better description of the deposition mechanism. For a uniform distribution, WZ ∝ constant. Aguirre et al. (2001) have investigated the metal enrichment of the diffuse intergalactic medium in cosmological simulations. They have considered a power-law shape W(r, h) ∝rα for the deposition kernel. They assumed a default value α = 3; the uniform case considered here corresponds to their α = 2. The simulations have been performed using for WZ the SPH smoothing kernel as the default kernel. Table 4 reports the indices of the runs for which WZ = constant. Each simulation has a label obtained by merging together the labels of Tables 2 and 3, or Tables 2 and 4, which determine the simulation parameters.
A measure of the regularity of the numerical cluster sample can be assessed from a comparison of the mean gas fractions of the four test clusters fg = Mg/MT within the radius r = 0.5 h−1 Mpc with the corresponding values of the Matsumoto et al. (2000) sample (Fig. 1). This sample has been chosen for a comparison of the global cluster iron abundances with those predicted by several runs (see below). The sample values of fg shown in Fig. 1 have a large scatter around the cosmological value of the model fb = Ωb/Ωm = 0.1, while those of the simulated clusters are systematically smaller. This indicates that the simulated clusters are regular objects (Evrard, Metzler & Navarro 1996) in a quiet dynamical state, without having undergone recent mergers, a conclusion supported also by previous substructure analyses of X-ray maps (Valdarnini et al. 1999) for a ΛCDM cluster sample with the same cosmological initial conditions. Therefore the profiles of simulations with different model parameters can be consistently compared without biases that can follow from the presence of substructure.
4.1 Simulations with different IMF and cooling
Results from simulations with different IMF and cooling parameters are discussed first. The relevant parameters for these runs are presented in Table 3. The simulation results are shown here only for the cluster ΛCDM39; for the other three test clusters there are not qualitatively relevant differences in the final profiles. The simulations have been performed keeping the numerical resolution, which is given by the index L39 of Table 2, fixed. The first three runs of Table 3 (L39NZIa1, L39NZIa2 and L39NZIa3) correspond to the following choices of the IMF: Miller—Scalo (MS), Salpeter (S) and Arimoto—Yoshi (AY), respectively. All the other parameters have been left unchanged. These simulations do not consider the possible dependence of the cooling function on the gas metallicity and the contribution of SNe of type Ia to the gas metal enrichment. For this choice of parameters, the run L39NZIa1 is just case cl39-10 of V02, in order to compare final profiles with previous findings. L39Ia4 has the metallicity dependence of the cooling function switched on. The other two indices of Table 3 (S and A) are for runs with SNe Ia as a heating source of the ICM, as well as of metal enrichment.
Fig. 2 shows the final profiles of temperature, gas and iron abundances in solar units for the different runs. The SF rate (SFR) as a function of age is also plotted. While a comparison of the iron abundances obtained from simulations with available data will be discussed later, several conclusions about the choice of the IMF can already be drawn from the iron profiles seen in Fig. 2. For a cluster like ΛCDM39, the iron abundance expected at a distance ≃0.4 Mpc (≃ 0.2 r200) from the centre is in the range ≃0.3–0.35 (De Grandi & Molendi 2001). Hereafter iron abundances are given in solar units (Fe/H = 4.68 × 10−5 by number). The iron profiles of Fig. 2 show that a MS IMF is completely ruled out as a possible IMF and cannot produce the amount of iron required by observations. The same conclusion holds also for an S IMF, for which Fe/H ≲ 0.1 at r≃ 0.4 Mpc. The only IMF for which the simulation gives a significant amount of iron is Arimoto—Yoshi (AY). The iron profile of the S run L39NZIa2 is always a factor of 2 below that of L39NZIa3 (AY). At r≃ 0.4 Mpc the latter profile gives Fe/H ≃ 0.2–0.15. These values are still below the ones required by observations; nevertheless, it appears that the AY IMF gives the best results from the point of view of the amount of iron required to fit the data. These conclusions are in agreement with those of Loewenstein & Mushotzky (1996).
These iron profiles originate from the metal ejecta of SNe II in the ICM, with the bulk of the SN explosions which have enriched the ICM already at z≃ 0.7. Therefore the gradients in the final iron profiles indicate that the cluster was already dynamically relaxed at this epoch, without major mergers of substructure which could have remixed the ICM and erased the original gradients. The most important differences in the profiles arise when the cooling function also depends on the gas metallicity (L39NIa4). For this run, the SFR is much higher than in the previous cases. For the range of temperatures considered here, this is a direct consequence of the higher cooling rate. As a consequence, the density threshold criterion for SF is lower and the SF activity is much higher than in the no-metal cases. The top right panel of Fig. 2 shows that final gas density profiles are not strongly affected by the metallicity dependence of the cooling rate.
One of the most important consequences of considering metallicity effects for the cooling rate is seen in the temperature profiles. From the top left panel of Fig. 2 there is clear evidence that the temperature profile of L39NIa4 has a different shape from the other profiles already at r≃ 200 kpc from the cluster centre. Between this distance and the cluster centre, the profile shows values of the temperature higher than in the runs with no-metal cooling. The metallicity dependence of the cooling rate influences the final temperature profiles at the cluster centre through two main effects. The first is the term Λc (T, Z), which enters in the gas thermal equation. Even for gas temperatures above ≃2 keV the increase in the cooling rate with respect to the no-metal case can be significant for high metallicities (Fe/H ≳ 1), which can be present at the cluster centre if there are strong metallicity gradients. The second effect is indirect: because of the increased cooling rate, SF has been higher in the past than in the no-metal case, and as a consequence there is a larger amount of cold gas that has been converted into stars and a deeper potential well at the cluster centre. As a result, there is a higher inflow at the cluster centre of the surrounding high-entropy gas than in the runs that do not consider the metallicity dependence of the cooling function. This in turn implies higher temperatures at the cluster centre (see Section 4.2). An important result that therefore follows from these simulations is that the temperature profiles in the cluster inner regions cannot be considered as approximately flat. This is in disagreement with what has been found in V02, where the simulations included radiative cooling and SF, but did not take into account the metal contribution to the cooling. These conclusions are rather general and are valid not only for ΛCDM39, which has a mass-weighted temperature Tm ≃ 3 keV, but also for a cluster like ΛCDM00, for which Tm ≃ 6 keV (see Fig. 5, later).
The iron abundance profile of the run L39NIa4 is lower by a factor of ∼2 with respect that of L39NZIa3, where the cooling function does not depend on the gas metallicity. This is also confirmed by contrasting the metal density profiles in the top panel of Fig. 2. If the cooling rate increases when the metallicity dependence is taken into account, then the final amounts of stars and metals at the cluster centre are expected to be higher than in the runs with no-metal cooling. This is in contrast with what is found. It is not clear how this anticorrelation between the final metal ICM abundance and the metal dependence of the cooling rate originates. A possible explanation lies in the strong SF activity that occurs earlier for the run under consideration. As a consequence, most of the ICM is metal-enriched by stars at earlier epochs. An iron profile that is in better agreement with data is recovered when SNe of type Ia (L39A) are also considered as sources of metal enrichment for the ICM. For this run, the iron abundance at ≃300 kpc is ≃0.2, about two times that of the parent run with the same parameters but without SNe Ia (L39NIa4). This is obtained for A = 0.07 (2). According to Portinari et al. (1998), this choice of A corresponds to a number of SN Ia ≃20 per cent of the total number of SNe in the Galaxy. A numerical simulation with A = 0.04 has shown that the amount of iron that is injected from SNe Ia in the ICM is roughly proportional to A. The above results indicate that for an AY IMF as much as ≃50 per cent of the iron in the ICM comes from SNe Ia, provided that the current rate is normalized to match that estimated from spiral galaxies.
The stability of the final iron profiles against the numerical resolution of the simulation can be estimated from Fig. 3, where for a chosen IMF (S) the iron profiles of the two clusters ΛCDM00 and ΛCDM40 are displayed for the standard resolution and the very high-resolution runs (VH). There is a good agreement between the profiles with different resolution, with small differences localized within ≲ 0.1 Mpc.
4.2 Simulations with different metal ejection parameters
Values of the global iron abundances AFe = MFe/MH are listed in solar units in Table 5 for the various simulations. These values have been calculated at a fiducial radius r = 0.5 h−1 Mpc and are in the range ≃0.15 –0.2 for the simulations with index A. From the Matsumoto et al. (2000) sample of nearby clusters (z < 0.1) the estimated abundances at 0.5 h−1 Mpc are all above ≃0.2 and in the range ≃0.2 –0.4 (Fig. 5, later). Therefore it is important to investigate the effects on the final iron abundances of adopting a prescription of metal ejection different from that of the A runs, keeping fixed the other parameters of the simulations.
The largest amounts of iron ejected in the ICM are obtained for the AY IMF, and this is the IMF that is chosen in the remainder of the paper in order to account for the measured iron abundances. Hence, for this IMF, different model parameters of metal ejection have been tested. Table 4 presents the indices of the runs and the associated relevant parameters.
The metal enrichment of the gas is modelled according to 2, with the mass that is ejected by a star particle i in the interval Δt distributed over the gas neighbours; the mass fraction is weighted according to the smoothing kernel WZ(r, h). For the A runs the kernel WZ is the standard SPH B2 spline. The radial profile of the deposition kernel is one of the parameters that governs the distribution of the ejected metals among the gas particles surrounding the star particle. According to Aguirre et al. (2001), the form of the distribution function WZ(r, h) can be generically assumed with a radial power-law ∝rα behaviour. A radial profile shallower than that of the B2 spline clearly implies that more metals are deposited near the limiting radius. Therefore gas particles that are not part of SF activity can be metal-enriched and diffuse metals into the ICM. The metal enrichment of the ICM is expected to be higher in this case because metals are more effectively mixed. This has also been pointed out by Mosconi et al. (2001). As an alternative to the B2 spline, the shape of the deposition profile chosen here is WZ = constant, which corresponds to the case of a uniform metal distribution.
The simulations with index A3 in Table 4 have WZ = const and are the mirror simulations of the runs with index A, the only different parameter being the choice of the deposition kernel. Final iron profiles are compared against those of models A in Fig. 4. The four panels of Fig. 4 refer to the four test clusters. In each panel final iron profiles are plotted for the runs with different metal ejection parameters and different numerical resolutions. The plots of Fig. 4 show that, for models A3, the profiles are higher than that of models A by a factor that is in the range ≃1–2. The differences are largest for ΛCDM00 and ΛCDM39, and modest for ΛCDM40. The global iron abundances AFe are now in the range ∼0.23–0.35 (Table 5), in better agreement with data than the values of models A. These results demonstrate that the choice of the deposition profile is a key parameter in determining iron abundances in the ICM. The differences in the final profiles between the runs A3 and A give a measure of the scatter in the final abundances that is inherent to the choice of different kernels.
It must be stressed that other choices are clearly possible: in particular a steeper IMF may even require a deposition kernel with a positive radial derivative. A quantitative comparison with several data of the profiles obtained from model A3 for the four test clusters is performed later. The results show that for models A3 the simulations are in good agreement with the measured values and the parameters of the A3 runs are then taken here as those of a ‘fiducial model’ against which to compare the results from the other simulations.
The stability of the simulated iron profiles against the numerical resolution is tested in Fig. 4 by plotting (filled symbols) profiles of high-resolution simulations (H in Table 2) against the corresponding ones with standard resolution. As can be seen, there are no large differences in the final profiles between the high- and standard-resolution runs. The only important exception is for the coldest cluster ΛCDM119 (Tm≃ 1.5 keV). The iron profile of L119HA3 is much lower than that of the standard-resolution run L119A3. This strongly suggests that there are numerical resolution problems that affect the profile of ΛCDM119, when the numerical resolution is that of the standard runs. Increasing the numerical resolution implies that smoothing lengths are smaller. As a consequence the metals are distributed among gas neighbours in smaller volumes; this implies a smaller mixing of metals in the ICM (see also Mosconi et al. 2001). Clearly, this effect should not be seen (cf. Fig. 3), and if present implies that low-resolution runs are undersampling the gas distribution. In order to predict metallicity profiles reliably, it is then safe to assume that simulations of low-temperature (T ≲ 3 keV) clusters require at least the mass resolution of high-resolution runs. To explain the low values of AFe for the run L119HA3, an alternative hypothesis is that the metal enrichment of the ICM by SNe is characterized by a minimum diffusion length. Simulations have been performed without a lower limit hmins for the smoothing lengths his of the diffusion kernel WZ of ejected metals. For cool clusters high-resolution runs may therefore imply values of his in the central high-density regions below hsmin. This in turn would imply lower values for the metallicity profiles because of the reduced mixing of the ejected metals in the ICM. This effect is very similar to the behaviour expected when the numerical resolution is insufficient and high-resolution runs yield lower profiles than those of low-resolution runs. The profiles of Fig. 3 show that this effect is relevant only for the less massive, cool clusters (≲ 3 keV). For the numerical parameters of the high-resolution runs (H) this implies that hmins must be of the order of a few kpc. In order to investigate the dependence of the final metallicity profiles on the value of hmins, two simulations have been performed for the cluster ΛCDM119. Models A5 and A6 have the same parameters as model A3, but with hmins = 12 and 6 kpc, respectively. Iron profiles for these two runs are plotted in the bottom right panel of Fig. 4 for the cluster ΛCDM119. The results demonstrate that requiring a minimum diffusion length is highly effective to obtain for cool clusters a large-scale mixing of iron in the ICM. The iron abundances AFe are now in the range ∼ 0.4–0.5, in better agreement with the observational values for cool clusters.
Another parameter that is important to determine the amount of metal enrichment of the ICM is the upper limit hM (see Table 4) of the smoothing lengths hsi. These are defined in 2 as the half-radius of a sphere enclosing ≃32 gas neighbours of the generic star particle i. The standard value assumed for the simulations is hM = 30 kpc.2 The sensitivity of the final profiles to the assumed values can be estimated from the plots of Fig. 4. For the runs A2, hM = 10 kpc. Contrasting the profiles with those of the A3 runs shows clearly that models A2 yield final abundances well below those of models A3 and comparable to the abundances of models A. For the runs A2, the abundances AFe are of the order of ≃ 0.1–0.13, a factor of ≃2 smaller than the values of runs A3. The value of hM = 10 kpc is clearly ruled out by the measured iron abundances. The choice hM = 30 kpc gives much better results for the metal distribution in the ICM. In fact, the values of his are rarely fixed by this upper limit. A radial binning of the squared distribution (his)2 shows that on average the rms values of the smoothing lengths hsi are below 10 kpc in high-density regions (r≲ 0.1 Mpc), and grow up to ≃20 kpc in low-density regions (r ≃ 500 kpc) for which ρ/ρc ≃ 10. This shows that the choice of the value of hM has a relevant impact on the final metallicity profiles at radial distances r ≳ 0.1 Mpc, in low-density regions where it is most likely that gas particles that are not in an SF region can get metal-enriched and diffuse metals. These values of the smoothing lengths correspond to maximum radii ≃40 kpc and are much higher than the upper limit of ≃10 kpc estimated by Ezawa et al. (1997) for the diffusion of ions in a Hubble time in a low-density plasma (ρ ≃ 10 ρc) at a temperature of ≃4 keV. On the other hand, the Ezawa et al. (1997) limit corresponds here to hM = 5 kpc; as has been found for models A2, this choice of hM would imply for the simulated clusters final iron abundances in the ICM much lower than the measured values. The plots of Figs 3 and 4 show also that for clusters with temperatures above ≳ keV the profiles of model A3, with hM = 30 kpc, do not depend on the numerical resolution of the simulations. A possible way of reconciling these discrepancies lies in the fact that analytical estimates of the maximum diffusion length of metals in the ICM do not take into account the mixing of metals that can occur because of the dynamical interactions between cluster galaxies and the ICM (Dupke & White 2000b). Simulation profiles show that to model the mixing of metals the best results are obtained for the parameters of model A3.
4.3 Comparison with data
For the four test clusters, observational variables from simulations with different model parameters are compared in Fig. 5 against a number of data. For model A3 the projected emission-weighted temperature profiles are plotted as a function of radius in panel (a). Data points are the mean error-weighted temperature profiles of 11 cooling flow clusters from De Grandi & Molendi (2002). The profiles have been calculated according to equation (A3) of De Grandi & Molendi. The cluster centre is defined as the maximum of the gas density, and a peak emission criterion for defining the centre does not modify the calculated profiles in a significant way. Each smoothed profile has been rescaled to match the last data point. There is a remarkable agreement of the simulated profiles with data, the only important exception being for the two innermost bins. The temperature profiles show a radial increase toward the cluster central region, followed by a strong drop at the cluster centre. This feature is common to all the runs and is not shared by the data, for which the innermost bin has a lower temperature than the nearby ones. This behaviour of the temperature profile is robust and is not sensitive to an increase of the numerical resolution. On the other hand, this increase of the cluster temperature at the centre is a consequence of the entropy conservation during galaxy formation and the subsequent removal of low-entropy gas (Wu & Xue 2002a). A possible explanation for this discrepancy lies in the fact that what are being measured are spectral temperatures. Mathiesen & Evrard (2001) argue that spectral fit temperatures are weighted by the number of photons, so line cooling from small clumps can bias the spectrum toward lower temperatures. This issue can be settled by performing a spectral fit analysis as in Mathiesen & Evrard (2001), whose simulations did not include radiative cooling.
The projected metallicity profiles are displayed as a function of distance in panel (b). Data points represents the mean profile from the nine cooling flow clusters of De Grandi & Molendi (2001). For model A3, the iron profile of L00A3 is the one in better agreement with data. The profile of L39A3 also has a shape similar to that of L00A3, but with a much higher abundance at r∼ 0.02 r200 (Fe/H ≃ 0.7 –0.8). The profile of L40A3 agrees fairly well with data only for the first three radial bins. The overall shape is similar to that of the other two runs, but with a lower amplitude. The systematic difference between the profiles of runs L39A3 and L40A3 is mostly due to the different dynamical histories of the two clusters. Therefore there are uncertainties in the final profiles which can be as high as ∼50 per cent and are related to the cluster dynamical evolution. At outer radii all the profiles show a radial decay steeper than that of data points, which in fact can be considered to have an almost constant profile. It appears very difficult to modify the model parameters of the simulations in such a way that the simulated profiles match the data points at outer radii, without also increasing central abundances.
The flatness of the observed profiles suggests that, at early epochs, a series of merger events has erased the existing metallicity gradients. Accordingly, the metal abundances from SNe II do not show significant spatial gradients, while the iron abundance gradient can be attributed to SNe Ia (Dupke & White 2000b). In this scenario the metal excess distribution is an indicator of the optical light distribution of an early-type galaxy, as has been at least partially confirmed by De Grandi & Molendi (2001). This is not the behaviour of the simulated cluster sample. Fig. 3 shows for two clusters the iron radial distribution originating from both SNe II and Ia. The iron distributions are very similar; because the time-scale of metal ejection is much higher for SNe Ia than for SNe II, this is indicative that the dynamical evolution of the two clusters has been very smooth. As already stressed, this cluster sample has been chosen for the regularity of its members. In order to assess in a significant way the effects of dynamical evolution on the shape of the metallicity profiles, it is necessary to perform a statistical analysis over a large (say, ≳40) cluster sample. This task is left to a future paper, where a number of issues will be investigated with a statistically robust cluster sample. For the sake of clarity, in panel (b) are also shown the profiles of the runs L40A and L40A2, to be compared with that of L40A3. These profiles are clearly below the measured values and confirm the parameters of model A3 as the ones yielding the best agreement with data.
The iron profile of L119A3 is inconsistent with the data points; the iron abundances are smaller at all radii. The profile is very similar to those of L40A2 and L40A. This is a failure of model A3 that is hard to reconcile within the framework of the adopted prescriptions, unless the diffusion of metals in the ICM is characterized by a minimum diffusion length. The most important difference of ΛCDM119 is that this is the coldest of the four test clusters, with a mass-weighted temperature of ∼1 keV. For this cluster a proper comparison of the simulated iron profiles with the mean metallicity profile of the sample is not possible. The averaged profile is that of nine cooling flow clusters, with minimum temperatures ≳ 4 keV. Without measured profiles for cool clusters it is therefore difficult to put observational constraints on different model parameters from the simulated profiles.
Global values, AFe = MFe/MH, of the iron abundances for the simulated clusters are compared in panel (c) for models A and A3 against the estimated values for the nearby cluster sample of Matsumoto et al. (2000). These values have been taken from column 4 of table 1 of Matsumoto et al. (2000) and are estimated at a radius of 0.5 h−1 Mpc. For model A3 there is a fair agreement with data, while model A yields values of AFe that are marginally consistent with the allowed uncertainties (90 per cent confidence level). The run L119HA3 has a value AFe ≃ 0.2, about a factor ∼2 smaller than that expected by extrapolating the sample average below the minimum (≃2 keV) cluster sample temperature. This is connected with the results shown in the previous panel. Without a minimum smoothing length there is a clear tendency in the simulations to produce a lesser amount of iron than that inferred from observations for cool clusters.
The open diamond of Fig. 5(c) is the value of AFe (∼0.47, see Table 5) for the run L119HA5. For this simulation, a minimum value hmins = 12 kpc has been assumed for the smoothing length of metals in the ICM. Such a high value of AFe is clearly more indicated to account for the iron abundance of low-temperature clusters. The observational evidence of an iron abundance decreasing with TX is, however, statistically weak (Mushotzky & Lowenstein 1997; Finoguenov et al. 2001). Furthermore, an increase of AFe for cool clusters can be an artefact arising from the presence of a dominant galaxy in the cluster central region. For example, Fukazawa et al. (2000) have removed the contribution of the central region for those clusters with a cD galaxy, and for their cluster sample found that there is not a significant correlation AFe−kTX. This is also in agreement with what was found by Finoguenov et al. (2001). A statistical comparison between the AFe distribution generated by the simulations and the one from real data has been performed for the models discussed above. The statistical tests applied to the AFe data and the AFe distribution of the simulated cluster sample are the Student t-test for the means, the F-test for the variances and the Kolmogorov—Smirnov (KS) test for the distributions. The corresponding probabilities give the significance level that, according to the analysed quantities, the two sets originated from the same process. The results are presented in Table 6. Model A is clearly ruled out; model A3 is marginally consistent. This is because the value AFe of L119HA3 is significant to lower the mean of the sample. A much better agreement is obtained if the global iron abundance is given by the run L119HA5 for the cluster ΛCDM119.
Finally, the final values of the bolometric X-ray luminosity are shown in panel (d) as a function of the cluster temperature. Mass-weighted temperatures have been used as unbiased estimators of the spectral temperatures (Mathiesen & Evrard 2001). The luminosity LX of the simulated clusters is calculated as described in Section 2. Data points are those of fig. 11 of Tozzi & Norman (2001). For the sake of clarity, not all the points in the figure are plotted in the panel. For a consistent comparison with data, a central region of size 50 h−1 kpc has been excised (Markevitch 1998) in order to remove the contribution to LX of the cooling flow central region. For model A3, the LX s of the simulations are in excellent agreement with data over the entire range of temperatures. An additional run has been performed for the cluster with the lowest temperature (ΛCDM119). The parameters of this run (A4) are the same as of model A3, expect for an SN feedback energy of 1050 erg being used for both SNe II and Ia. This value is 10 per cent of that of model A3 and has been considered in order to investigate the effects on final X-ray properties of the amount of heating returned to the ICM by SNe. As can be seen, for model A4 the final LX of the run L119HA4 is very similar to that of L119HA3. The result demonstrates that the final X-ray luminosities of the simulations are not sensitive to the amount of SN feedback energy that has heated the ICM. In fact, a simulation run with a zero SN energy returned to the ICM (εSN = 0) yields very similar results. These results are particularly relevant in connection with the recent proposal (Bryan 2000) that the X-ray properties of the ICM are driven by the efficiency of galaxy formation, rather than by heating caused by non-gravitational processes. Final profiles for models A3 and A4 are investigated in more detail in Figs 6 and 7.
4. SN heating of the ICM and cooled gas fractions
In Fig. 6, the final profiles of runs L119A3 and L119A4 are compared against those of L119HA3 and L119HA4. The main differences are between the profiles of the standard-resolution run and the high-resolution run. The profiles of the runs for the models A3 and A4 are very similar, with the only important difference between the two models being the SFR at early redshifts. The temperature profiles are almost identical, so energy feedback from SNe is not relevant to determine final gas properties. This conclusion is also supported by the entropy profiles of Fig. 7.
In Fig. 7, final distributions of the two runs L119HA3 and L119HA4 are compared. In panel (a) are plotted the Si abundances synthesized in the explosions of SNe II and SNe Ia. There are not large differences between the corresponding profiles. The Si abundances of SNe Ia are almost identical. For SNe II, the Si abundances profile of L119HA4 is higher than that of L119HA3 in the cluster central region (≲0.1 Mpc). These differences follow because run L119HA4 has at early redshifts a higher SFR than that of L119HA3 (see also c). This is a consequence of the reduced thermal pressure and higher gas densities in the former simulation. This effect is not relevant for the SN Ia explosions because the involved time-scales are much higher than those of the early generations of SNe II.
The cumulative injected energy per baryon within the radius r is defined as the ratio ESN(<r)μmp/Mg(<r), where ESN(<r) is the total SN energy injected within that radius, and Mg(<r) is the corresponding gas mass. This ratio is shown in panel (d) for the two runs. The ratio has a radially decreasing behaviour, for the run L119HA3 drops from ∼1 keV particle−1 at ≃0.1 Mpc to ≃0.5 keV particle−1 at r≃ 1 Mpc ≃r200. There is a small drop in the cluster inner region, presumably arising from a biased estimate because of the small number of simulation particles for r≲ 40 kpc. For the run L119HA4 the energy per particle ratio has a shape similar to that of the corresponding ratio for the run L119HA3; at r∼ 1 Mpc it takes the value of ∼0.04 keV particle−1, which is about ∼10 per cent of the value of the parent run with εSN = 1051 erg. The binned distribution (panel c) is much more noisy than the cumulative one, but it qualitatively confirms the expected behaviour.
In order to assess the amount of heating from the SNe, this ratio is often inferred from the measured abundance of metals (Finoguenov et al. 2001). Si abundances are particularly useful because of the similar yields for different SN types. The ratio estimated from the Si abundance ASi is εSNASiμmp/ySi with ySi∼ 0.12 being the average SN yield of Si. The estimated ratio (thin lines) is in good agreement at all radii with the values obtained from the simulation. For εSN = 1051 erg the average energy per particle of the cluster ΛCDM119 is then ∼0.5 keV particle−1 at the virial radius (r200∼ 1 Mpc). This cluster has a virial temperature ∼1.5 keV, from a sample of 18 relaxed clusters with temperature below ∼4 keV. Finoguenov et al. (2001 see fig. 9) find an average SN injected energy of ∼0.5 keV particle−1 for a cluster with a temperature ∼1 keV. There is thus a good agreement between the measured amount of thermal energy per particle associated with SN feedback and that predicted by the simulations with the model parameters of Fig. 7.
The gas entropy is defined as S(r) =kBT(r)/ne(r)2/3, where ne is the electron number density. For the two runs L119HA3 and L119HA4 the entropy profiles in units of keV cm2 are plotted in panel (b) of Fig. 7. An important result that follows from the radial dependence of the two profiles is that they are almost identical. A simulation with a zero SN energy yields a very similar profile. This means that final ICM properties are largely unaffected by the amount of energy injected by SNe into the ICM. This follows because most of the energy is injected in the cluster central region where the gas density is higher and cooling is very efficient at radiating away the energy of the reheated gas. Energy feedback from SNe can modify the ICM state for small clusters or groups with TX≲ 1 keV. Another important result from the entropy profiles of Fig. 7 is the entropy level of the gas at r∼ 0.1RV≃ 0.1r200≃ 0.1 Mpc. At this radial distance the gas entropy of the simulated cluster is ∼200 keV cm2, a value consistent with a set of observations (Ponman et al. 1999; Lloyd-Davies et al. 2000; Wu & Xue 2002a) for a system with TX∼ 1.5 keV. Together with the good agreement between the LX−TX relation obtained from simulations and that from observations, these findings provide strong support for the radiative cooling model proposed by Bryan (2000).
According to this scenario (Bryan 2000; Voit & Bryan 2001; Wu & Xue 2002a; Muanwong et al. 2002; Davé, Katz & Weinberg 2002), the efficiency of galaxy formation is higher in groups and cool clusters than in hot clusters. The low-entropy cooled gas is removed because of galaxy formation, and is replaced by the inflow of the surrounding high-entropy gas. The higher efficiency of galaxy formation for cool systems explains the central entropy excess over the self-similar predictions (Ponman et al. 1999). One of the most important observational consequences of the cooling model is that the fraction of hot gas fg = Mg/MT increases going from cool clusters to hot clusters. Conversely, the fraction fstar of cooled gas turned into stars should decrease with the mass of the system.
Observational evidence for a dependence of fg with TX is controversial (David et al. 1990; David, Jones & Forman 1995; Mohr, Mathiesen & Evrard 1999; Arnaud & Evrard 1999; Roussel, Sadat & Blanchard 2000; Balogh et al. 2001). The main difficulty is the extrapolation to virial radii of the X-ray data for cool clusters, which can bias the estimates of the cooled gas fraction. According to Roussel et al. (2000), there is not strong observational support for an increase of fgas with TX. This is in disagreement with the findings of David et al. (1995), who show that the gas fraction fgas increases with the X-ray temperature. According to Mohr et al. (1999) there is a weak dependence of fgas on TX in their X-ray ROSAT sample of 45 galaxy clusters. A constant fgas is inconsistent with the data at the 95 per cent confidence level.
Arnaud & Evrard (1999) have analysed the temperature dependence of fgas for a sample of 24 clusters with accurately measured temperatures and low cooling flows. They have determined fgas for two radii enclosing gas overdensities δ = 500 and 200. Total masses have been estimated according to two different methods. The β-model (BM) assumes an isothermal gas density profile in hydrostatic equilibrium with the shape determined according to the X-ray surface brightness. The virial theorem (VT) estimates the cluster total mass according to virial equilibrium at a fixed density contrast. The relation is calibrated from a set of numerical experiments (Evrard et al. 1996) and is not sensitive to the assumed cosmological model. According to Arnaud & Evrard (1999), it is difficult to draw conclusions on the dependence of fgas on the cluster temperature. The results depend on the chosen model. For the BM model, the difference in the value of fgas between subsamples of cool clusters and hot (>4 keV) clusters is not statistically significant. For the VT model this difference is at the 3σ level at a radius corresponding to a gas overdensity δ = 500. Arnaud & Evrard (1999) indicate a conservative upper limit of ∼30 per cent in the 1σ fractional error of fgas.
For the cluster sample of Arnaud & Evrard (1999), the gas fraction fgas is shown as a function of the total cluster mass in the top panels of Fig. 8. The sample distributions are derived from the VT model at a density contrast of δ = 500 (top left panel) and δ = 200 (top right panel). The data points are those in the top panels of fig. 3 of Arnaud & Evrard (1999), with an assumed ≃30 per cent uncertainty (top right panel). For the model A3 the corresponding values of fgas at δ = 500 and 200 from the numerical cluster sample are also shown as open squares. There is a clear tendency for the two simulated cluster distributions to follow the fgas distribution of the data, but with a reduced amplitude. The result of a statistical comparison are reported in Table 6. The fgas distribution of model A3 at δ = 500 is inconsistent with data with a high significance level (95 per cent confidence level), and the situation at δ = 200 is even worse. For this overdensity, an extrapolation of X-ray data up to the required radius has been performed (Arnaud & Evrard 1999) in order to evaluate the corresponding data points. Therefore possible biases can undermine the estimate of fgas at δ = 200. For fgas(δ = 500) a possible source of disagreement between the numerical distribution and that of data points lies in the assumed value of the cosmological baryonic density Ωb in the simulations. Here it has been assumed that Ωbh2 = 0.015, but recent measurements (Burles & Tytler 1998) favour Ωbh2 = 0.019. In correspondence to this value of Ωb in the initial conditions of the simulations, the values of fgas(δ = 500) are displayed in Fig. 8(a) for a subsample of three simulated clusters (ΛCDM00, ΛCDM119 and a new cluster with virial mass ∼2.5 × 1014, h−1 M⊙ ). The distribution is now in a better agreement with the data, and the confidence levels for rejecting the null hypothesis are now below 95 per cent. Therefore it is fair to say that the ICM gas fractions obtained here are consistent with available observational estimates.
The predictions of the radiative cooling model for ICM evolution have been recently investigated in a number of papers, either through analytical methods (Bryan 2000; Voit & Bryan 2001; Wu & Xue 2002a; Voit et al. 2002; Wu & Xue 2002b) or with numerical simulations (Muanwong et al. 2001, 2002; Davé et al. 2002). Results from the employed methods clearly show a decreasing fstar with cluster temperature. The dependence of fstar(δ = 200) on TX is shown in Fig. 8(c) for the four test clusters. Open squares are for model A3 with standard resolution and filled squares denote the high-resolution runs. For comparative purposes the fstar predicted by the analytical model of Bryan (2000) is also shown. Because of the different values of Ωb and Ωm in the paper by Bryan (2000), fstar has been approximately rescaled by the factor fb (here) /fb (Bryan) = 0.1/0.16 = 0.625. There is a fair agreement of the model with the simulations at high temperatures, but for TX∼ 1 keV the fstar of the simulations is below the model predictions by a factor of ∼2. The disagreement at low temperatures is ameliorated for the runs with Ωbh2 = 0.019; in such a case fstar(Bryan) must be rescaled by the factor ∼0.81. The predicted value of fstar is now ∼0.075 at TX≃ 1 keV, which is ≃50 per cent higher than the value ∼0.05 of fstar at the same temperature for the runs with Ωbh2 = 0.019.
The comparison of fstar with data is a controversial issue (Balogh et al. 2001), also because of the lack of firm bounds on the allowed uncertainties (Roussel et al. 2000); however, it is worth stressing that Roussel et al. (2000) found a weak dependence of the gas fraction on cluster temperature when the estimates of virial masses are calibrated from numerical simulations rather than inferred from the BM. The dependence is not statistically significant because of the large size of the sample error bars. A comparison of fstar predicted by an analytical model has been performed by Wu & Xue (2002b) against the Roussel et al. (2000) data. The theoretical predictions are about twice as high as the sample values. From Fig. 8(c) the values of fstar(δ = 200) for the three runs with Ωbh2 = 0.019 can be compared with the analytical predictions of the Wu & Xue (2002b) model (cf. fig. 2 of their paper). The values of Ωb and h are slightly different, but there is a broad agreement. At TX∼ 1 keV fstar( here) ≃ 0.05 is lower than the theoretically predicted value ≃0.07. As noted by Wu & Xue (2002b), theoretical models of cooling yield a dependence of fstar on TX steeper than that found in hydrodynamical simulations and observations. At low temperatures the fstar of the simulations is in better agreement with the Roussel et al. (2000) data, but the observed uncertainties do not allow firm conclusions. The simulation results for fstar are also in qualitative agreement with those of the radiative model of Muanwong et al. (2002). In units of Ωb/Ωm the values of fstar(fg) range from ∼0.38(0.46) from M200 = 0.8 × 1014h−1 M⊙ to ∼0.34 (0.56) for a cluster with a virial mass of M200 ∼ 2.4 × 1014h−1 M⊙. These values are not very different from those that can be inferred from the distributions plotted in fig. 3(b) of Muanwong et al. (2002) for their radiative model. The values of the cosmological parameters of the simulations of Muanwong et al. (2002) are close to those assumed here (Ωm is ∼20 per cent higher), and the numerical resolution is broadly similar.
Finally, the effects of the numerical resolution on the predicted final amount of cooled gas are shown in panel (d) of Fig. 8. As a function of the gas temperature, the fractional contribution to the total X-ray luminosity is plotted from the right for runs of three clusters (ΛCDM00, ΛCDM40 and ΛCDM119) with different resolutions. For the first two clusters the ratio of LX(<T)/LT for very high-resolution runs is compared against that of standard resolution. The ratio LX(<T)/LT of the run L119A3 is instead compared against that of L119HA3. From the left, the gas masses above a given gas temperature are plotted for the corresponding runs. The distribution of the gas mass versus the gas temperature shows that resolution effects are most important in the low-temperature region of the gas temperature distribution, whereas X-ray luminosities are largely determined by the high-temperature part of the distribution. This demonstrates that the fraction of cold gas depends on the numerical resolution, as can be seen from Fig. 8(c), but is not a source of worry as far as X-ray properties are concerned.
For the cluster ΛCDM119, the quantity Mgas(>T) is weakly dependent on the numerical resolution of the simulations. From Fig. 8(c), one can see that fstar at low temperatures also seems to converge as the resolution is increased. This issue is strictly related to the global value of the baryonic fraction of cooled gas Ωcold/Ωb = fb (global). Observational upper limits indicate per cent (Balogh et al. 2001). According to Balogh et al. (2001), the global fraction fb (global) is expected to increase as the resolution limit of the simulation is increased, unless a feedback model is incorporated in the simulation. The results of the high-resolution runs suggest that at low temperatures convergence is being achieved for fstar and that the feedback model implemented here effectively regulates the amount of cooled gas. This is strongly indicated by the large value of fstar (∼0.05) when the SN feedback energy is reduced (open diamond of Fig. 8c, which corresponds to the run L119HA4). For clusters with temperatures above ≳ 2 keV, the values of fstar appear to converge for simulations with number of gas particles Ng ≳ 200 000. This is shown in Fig. 8(c), where for the cluster ΛCDM00 the value of fstar for the very high-resolution run L00VHA3 (filled triangle) is ∼0.03 and is close to that of fstar for L00HA3. Note that the mass of the gas particle for the run L00VHA3 is only ∼10 per cent higher than that of L119HA3 (≃ 7 × 108 M⊙). These issues can be clarified in more detail with a set of cosmological simulations with different resolutions that incorporate the SN feedback scheme implemented here. Such analysis is beyond the scope of this paper and will be considered in a future work.
5 Summary and Conclusions
In this paper, results from a large set of hydrodynamical smoothed particle hydrodynamics simulations of galaxy clusters have been used to investigate the dependence of iron abundances and heating of the intracluster medium on a number of model parameters. The simulations have been performed with different numerical resolutions and a numerical cluster sample covering nearly a decade in cluster mass. The modelling of the gas physical processes in the simulations incorporates radiative cooling, star formation and energy feedback. The gas is metal-enriched from supernova ejecta of type II and Ia, and the cooling rate depends on the local gas metallicity. This allows us to follow self-consistently the ICM evolution also for cool clusters, whose departures of cluster scaling relations from self-similarity are most important.
The metal enrichment of the ICM is governed by a number of model parameters. Theoretical uncertainties on the shape of the initial mass function and nucleosynthesis stellar yields lead to final iron abundances which can differ by a factor ∼2. These two parameters have then been kept fixed for a large fraction of the performed simulations. This is done in order to restrict the full range of the parameter space. They have been chosen so that the corresponding final amounts of iron in the ICM are those that most closely agree with that indicated by observations. The ejected metals are distributed among gas neighbours according to the SPH formalism. Final iron abundances have been found to be sensitive to the shape of the smoothing kernel of the ejected metals WZ and to the constraints on the corresponding smoothing lengths hs. For a certain choice of the above parameters, a fair agreement has been found between the observed radial metallicity profiles and those obtained from the simulated clusters. Global iron abundances are also consistent with measured values. There are, however, a few issues that are still unresolved.
For cool clusters the ejection parameters are not well determined because of the lack of measured iron profiles; also constraints from global iron abundances suffer from selection biases in the cluster samples. The simulated profiles have a radial decay steeper than observed. This difference could be mainly due to the lack in the simulated cluster sample of mergers that have efficiently remixed the metal content of the ICM. The small size of the numerical cluster sample does not allow us to reach firm conclusions. In order to investigate the correlation between the shape of the final metallicity profiles and the cluster dynamical evolution, it is then necessary to analyse simulation results from a statistically robust (say, ≳40) cluster sample.
A discrepancy with observations is given by the radial behaviour of the projected emission-weighted temperature profiles, which have a steep rise at r≲ 0.1r200. This feature is not shared by observations; for the innermost bin the simulated temperatures are higher by a factor of ∼2 than the measured values. This difference cannot be easily accommodated by the gas physical modelling of the simulations. Because of radiative cooling, an increase of the gas temperature toward the cluster centre is expected as a consequence of entropy conservation (Wu & Xue 2002a). A certain amount of cooled gas is present in the core of the clusters and is responsible of the strong decline of the temperature profiles at distances very close to the cluster centres (r≲ 0.02 r200). This feature could help to reconcile the differences with the shapes of the observed profiles. As found by Mathiesen & Evrard (2001), the measured spectral fit temperatures can be biased toward lower values because of line cooling. An X-ray spectral fit analysis is then necessary to resolve the discrepancies between the simulated and measured temperatures profiles in the cluster core regions.
The ICM X-ray properties of cool clusters are largely unaffected by the amount of feedback energy injected by SNe, with the luminosity—temperature relation being in excellent agreement with data. The final gas entropy distribution has been found to be almost independent of the heating efficiency; for TX∼ 1 keV the core entropy excess is in agreement with estimated values. These findings support the radiative cooling model of Bryan (2000), where the ICM X-ray scaling relations are driven by the efficiency of galaxy formation. The model predicts that at the cluster virial radius the fraction of hot gas fg is positively correlated with the cluster temperature TX. A number of observational uncertainties prevent tight constraints, but the fg−TX distribution of the simulated cluster sample is not inconsistent with available estimates from X-ray data.
To explain the steepness of the luminosity—temperature relation, an alternative to the cooling model is the pre-heating scenario, where the gas has been heated by an energy injection that occurred at early epochs. Simulation results (Bialek, Evrard & Mohr 2001; Muanwong et al. 2002; Borgani et al. 2002) show that a number of observed X-ray properties can be reproduced in the pre-heating model. Future X-ray observations of cool clusters and groups with Chandra will help to discriminate between the two models. It is worth noting that simulations of galaxy clusters in the pre-heating scenario must also include radiative cooling and SN feedback. Therefore the agreement with observations of X-ray properties reproduced here by the simulated clusters in the cooling model suggests that additional heat sources, if required, are most likely to affect the ICM properties of systems with temperatures below ≲1 keV.
The author is grateful to M. Arnaud, S. De Grandi and P. Tozzi for providing their data files. G. Carraro and S. Recchi are also gratefully acknowledged for clarifying comments on chemical evolution. The author also thanks an anonymous referee for useful comments and suggestions which improved the presentation of the paper.