where (z) = 3H_{0 m}(0) / 8 Ga³ is the proper background density of the universe.

The cosmological scale factor a(t) is obtained by integrating the Friedmann equation (Eq. 4). To complete the specification of the problem we need the ideal gas equation of state p = ( - 1)e, and the gas heating and cooling rates. When simulating the ICM, the simplest approximation is to assume and = 0; i.e., no heating or cooling of the gas other than by adiabatic processes and shock heating. Such simulations are referred to as adiabatic (despite entropy-creating shock waves), and are a reasonable first approximation to real clusters because except in the cores of clusters, the radiative cooling time is longer than a Hubble time, and gravitational heating is much larger than sources of astrophysical heating. However, as discussed in the paper by Cavaliere in this volume, there is strong evidence that the gas in cores of clusters has evolved non-adiabatically. This is revealed by the entropy profiles observed in clusters [29] which deviate substantially from adiabatic predictions. In the simulations presented below, we consider radiative cooling due to thermal bremsstrahlung, and mechanical heating due to galaxy feedback, details of which are described below.

4.3. Numerical methods overview

A great deal of literature exists on the gravitational clustering of CDM using N-body simulations. A variety of methods have been employed including the fast grid-based methods particle-mesh (PM), and particle-particle+particle-mesh (P³M) [30], spatially adaptive methods such as adaptive P³M [31], adaptive mesh refinement [24], tree codes [32, 33], and hybrid methods such as TreePM [34]. Because of the large dynamic range required, spatially adaptive methods are favored, with Tree and TreePM methods the most widely used today. Fig. 13 shows a high resolution N-body simulation of the substructure within a dark matter halo performed by Springel using the GADGET code [78].

Figure 13. Ultra-high resolution N-body simulations of the clustering of dark matter reveal a wealth of substructure. From [78].

When gas dynamics is included, only certain combinations of hydrodynamics algorithms and collisionless N-body algorithms are "natural". Dynamic range considerations have led to two principal approaches: P³MSPH and TreeSPH, which marries a P³M or tree code for the dark matter with the Lagrangian smoothed-particle-hydrodynamics (SPH) method [35, 21, 22], and adaptive mesh refinement (AMR), which marries PM with Eulerian finite-volume gas dynamics schemes on a spatially adaptive mesh [23, 26, 25, 36]. Pioneering hydrodynamic simulations using non-adaptive Eulerian grids [37, 38, 13] yielded some important insights about cluster formation and statistics, but generally have inadequate resolution to resolve their internal structure in large survey volumes. In the following we concentrate on our latest results using the AMR code Enzo [26]. The reader is also referred to the paper by Borgani et al. [39] which presents recent, high-resolution results from a large TreeSPH simulation.

Enzo is a grid-based hybrid code (hydro + N-body) which uses the block-structured AMR algorithm of Berger & Collela [40] to improve spatial resolution in regions of large gradients, such as in gravitationally collapsing objects. The method is attractive for cosmological applications because it: (1) is spatially- and time-adaptive, (2) uses accurate and well-tested grid-based methods for solving the hydrodynamics equations, and (3) can be well optimized and parallelized. The central idea behind AMR is to solve the evolution equations on a grid, adding finer meshes in regions that require enhanced resolution. Mesh refinement can be continued to an arbitrary level, based on criteria involving any combination of overdensity (dark matter and/or baryon), Jeans length, cooling time, etc., enabling us to tailor the adaptivity to the problem of interest. The code solves the following physics models: collisionless dark matter and star particles, using the particle-mesh N-body technique [41]; gravity, using FFTs on the root grid and multigrid relaxation on the subgrids; cosmic expansion; gas dynamics, using the piecewise parabolic method (PPM) [42]; multispecies nonequilibrium ionization and H₂ chemistry, using backward Euler time differencing [28]; radiative heating and cooling, using subcycled forward Euler time differencing [43]; and a parameterized star formation/ feedback recipe [44]. At the present time, magnetic fields and radiation transport are being installed. Enzo is publicly available at http://lca.ucsd.edu/projects/enzo.

4.4. Structure of nonradiative clusters: the Santa Barbara test cluster

In Frenk et al. [20] 12 groups compared the results of a variety of hydrodynamic cosmological algorithms on a standard test problem. The test problem, called the Santa Barbara cluster, was to simulate the formation of a Coma-like cluster in a standard CDM cosmology (_m = 1) assuming the gas is nonradiative. Groups were provided with uniform initial conditions and were asked to carry out a "best effort" computation, and analyze their results at z = 0.5 and z = 0 for a set of specified outputs. These outputs included global integrated quantities, radial profiles, and column-integrated images. The simulations varied substantially in their spatial and mass resolution owing to algorithmic and hardware limitations. Nonetheless, the comparisons brought out which predicted quantities were robust, and which were not yet converged. In Fig. 14 we show a few figures from Frenk et al. (1999) which highlight areas of agreement (top row) and disagreement (bottom row).

Figure 14. The Santa Barbara test cluster. Top row, left to right: profiles of dark matter density, gas density, and gas pressure. Bottom row, left to right: profiles of gas temperature, gas entropy, and X-ray emissivity. Different symbols correspond to different code results. From [20].

The top row shows profile of dark matter density, baryon density, and pressure for the different codes. All are in quite good agreement for the mechanical structure of the cluster. The dark matter profile is well described by an NFW profile which has a central cusp [45]. The baryon density profiles show more dispersion, but all codes agree that the profile flattens at small radius, as observed. All codes agree extremely well on the gas pressure profile, which is not surprising, since mechanical equilibrium is easy to achieve for all methods even with limited resolution. This bodes well for the interpretation of SZE observations of clusters, since the Compton y parameter is proportional to the projected pressure distribution. In section 5 we show results from a statistical ensemble of clusters which bear this out.

The bottom row shows the thermodynamic structure of the cluster, as well as the profile of X-ray emissivity. The temperature profiles show a lot of scatter within about one-third the virial radius (= 2.7 Mpc). Systematically, the SPH codes produce nearly isothermal cores, while the grid codes produce temperature profiles which continue to rise as r 0. The origin of this discrepancy has not been resolved, but improved SPH formulations come closer to reproducing the AMR results [51]. This discrepancy is reflected in the entropy profiles. Again, agreement is good in the outer two-thirds of the cluster, but the profiles show a lot of dispersion in the inner one third. Discounting the codes with inadequate resolution, one finds the SPH codes produce an entropy profile which continues to fall as r 0, while the grid codes show an entropy core, which is more consistent with observations [29]. The dispersion in the density and temperature profiles are amplified in the X-ray emissivity profile, since _{x b}² T^1/2. The different codes agree on the integrated X-ray luminosity of the cluster only to within a factor of 2. This is primarily because the density profile is quite sensitive to resolution in the core; any underestimate in the core density due to inadequate resolution is amplified by the density squared dependence of the emissivity. This suggests that quite high resolution is needed, as well as a good grasp on non-adiabatic processes operating in cluster cores, before simulations will be able to accurately predict X-ray luminosities.

4.5. Effect of additional physics

Within r = 0.15 r_vir, Vikhlinin et al. [50] found large variation in temperature profiles, but in all cases the gas is cooler than the cluster mean. This suggests that radiative cooling is important in cluster cores, and possibly other effects as well. It has been long known that ~ 60 percent of nearby, luminous X-ray clusters have central X-ray excesses, which has been interpreted as evidence for the presence of a cluster-wide cooling flows [64]. More recently, Ponman et al. [29] have used X-ray observations to deduce the entropy profiles in galaxy groups and clusters. They find an entropy floor in the cores of clusters indicative of extra, non-gravitational heating, which they suggest is feedback from galaxy formation. It is easy to imagine cooling and heating both may be important to the thermodynamic evolution of ICM gas.

To explore the effects of additional physics on the ICM, we recomputed the entire sample of clusters changing the assumed baryonic physics, keeping initial conditions the same. Three additional samples of about 100 clusters each were simulated: The "radiative cooling" sample assumes no additional heating, but gas is allowed to cool due to X-ray line and bremsstrahlung emission in a 0.3 solar metallicity plasma. The "star formation" sample uses the same cooling, but additionally cold gas is turned into collisionless star particles at a rate _SF = _{sf b} / max(_cool , _dyn), where _sf is the star formation efficiency factor ~ 0.1, and _cool and _dyn are the local cooling time and freefall time, respectively. This locks up cold baryons in a non-X-ray emitting component, which has been shown to have an important effect of the entropy profile of the remaining hot gas [56, 57]. Finally, we have the "star formation feedback" sample, which is similar to the previous sample, except that newly formed stars return a fraction of their rest mass energy as thermal and mechanical energy. The source of this energy is high velocity winds and supernova energy from massive stars. In Enzo, we implement this as thermal heating in every cell forming stars: _sf = _{SN SF} c². The feedback parameter depends on the assumed stellar IMF the explosion energy of individual supernovae. It is estimated to be in the range 10^-6 _SN 10^-5 [44]. We treat it as a free parameter.

Fig. 15 shows synthetic maps of X-ray surface brightness, temperature, and Compton y-parameter for a M = 2 × 10¹⁵ M cluster at z = 0 for the three cases indicated. The "star formation" case is omitted because the images are very similar to the "star formation feedback" case (see reference [46].) The adiabatic cluster shows that the X-ray emission is highly concentrated to the cluster core. The projected temperature distribution shows a lot of substructure, which is true for the adiabatic sample as a whole [47]. A complex virialization shock is toward the edge of the frame. The y-parameter is smooth, relatively symmetric, and centrally concentrated. The inclusion of radiative cooling has a strong effect on the temperature and X-ray maps, but relatively little effect on the SZE map. The significance of this is discussed in Section 5. In simulations with radiative cooling only, dense gas in merging subclusters cools to 10⁴ K and is brought into the cluster core intact [48]. These cold lumps are visible as dark spots in the temperature map. They appear as X-ray bright features. The inclusion of star formation and energy feedback erases these cold lumps, producing maps in all three quantities that resemble slightly smoothed versions of the adiabatic maps. However, an analysis of the radial temperature profiles (Fig. 15) reveal important differences in the cluster core. The temperature continues to rise toward smaller radii in the adiabatic case, while it plummets to ~ 10⁴ K for the radiative cooling case. While the temperature profile looks qualitatively similar to observations of so-called cooling flow clusters, our central temperature is too low and the X-ray brightness too high. The star formation feedback case converts the cool gas into stars, and yields a temperature profile which follows the UTP at r 0.15r_vir, but flattens out at smaller radii. This is consistent with the high resolution Chandra observations of Vikhlinin et al. [50].

Figure 15. Effect of physical processes on simulated galaxy cluster ICM observables. Left: Columns show X-ray surface brightness, projected temperature, and Compton y-parameter for a M = 2 × 1015 M cluster assuming different baryonic physics. Field of view is 5 h-1 Mpc. Right: Corresponding spherically averaged radial temperature profiles. From [46].