text:
# Cooling and heating the ICM in hydrodynamical simulations

## 1 Introduction

The simplest picture to describe the thermal properties of the intra–cluster medium (ICM) is based on the assumption that gas heating occurs only by the action of gravitational processes, such as adiabatic compression from gravitational collapse, and by hydrodynamical shocks from supersonic accretion (Kaiser 1986). Since gravity does not have characteristic scales, this model predicts that galaxy systems of different mass look like scaled versions of each other. Under the assumptions of thermal bremsstrahlung emissivity and hydrostatic equilibrium, this model provides precise predictions for $`X`$–ray scaling properties of galaxy systems: (a) $`L_XT^2(1+z)^{3/2}`$ for the shape and evolution of the relation between $`X`$–ray luminosity and gas temperature; (b) $`ST(1+z)^2`$ for the entropy-temperature relation, where $`S=T/n_e^{2/3}`$ is the gas entropy and $`n_e`$ is electron number density; (c) $`MT^{3/2}`$ for the relation between total cluster virial mass and temperature, with normalization determined by the parameter $`\beta =\mu m_p\sigma _v^2/k_BT`$ ($`\mu =0.59`$ mean molecular weight for primordial composition; $`m_p`$: proton mass; and $`\sigma _v`$: line-of-sight velocity dispersion). Numerical simulations that only include gravitational heating showed that $`\beta 1`$–1.3 (e.g. Navarro, Frenk & White 1995; Evrard, Metzler & Navarro 1996; Bryan & Norman 1998; Eke et al. 1998b; Borgani, Governato, Wadsley et al. 2002, BGW hereafter).

A number of observational facts demonstrate that this picture is too simplistic, thus calling for the consideration of extra physics in the description of the ICM. The $`L_X`$$`T`$ relation is found to be steeper than predicted, with $`L_XT^3`$ at $`T_X>2`$ keV (e.g., White, Jones & Forman 1997; Markevitch 1998; Arnaud & Evrard 1999; Ettori, De Grandi & Molendi 2002), possibly approaching the self–similar scaling only for the hottest systems with $`T>\text{ }8`$ keV (Allen & Fabian 1998). Evidences also emerged for this relation to further steepen for colder groups, $`T<\text{ }1`$ keV (e.g., Ponman et al. 1996; Helsdon & Ponman 2000; Mulchaey 2000). Furthermore, no evidence for a strong positive evolution of the $`L_X`$$`T`$ relation has been found to date out to $`z1`$ (e.g., Mushotzky & Scharf 1997; Reichart et al. 1999; Fairley et al. 2000; Borgani et al. 2001a; Holden et al. 2002; Novicki, Sornig & Henry 2002; cf. also Vikhlinin et al. 2002). As for the $`S`$$`T`$ relation, Ponman, Cannon & Navarro (1999) found from ROSAT and ASCA data an excess of entropy within the central regions of $`T<\text{ }2`$ keV systems (see also Lloyd–Davis et al. 2000, Finoguenov et al. 2002a), possibly approaching the value $`S100`$ keV cm2 for the coldest groups. Finally, a series of evidences, based on ASCA (e.g., Horner, Mushotzky & Scharf 1999; Nevalainen, Markevitch & Forman 2000; Finoguenov, Reiprich & Böhringer 2001b), Beppo–SAX (Ettori et al. 2002) and Chandra (Allen, Schmidt & Fabian 2001) data, shows that the observed $`M`$$`T`$ relation has a $`40`$ per cent lower normalization than predicted by simulations that only include gravitational heating.

In the attempt of interpreting these data, theoreticians are currently following two alternative routes, based either on introducing non–gravitational heating of the ICM or on alluding to the effects of radiative cooling.

An episode of non–gravitational heating, occurring before or during the gravitational collapse, has the effect of increasing the entropy of the gas, preventing it from reaching high densities in the central cluster regions and suppressing its $`X`$–ray emissivity (e.g., Evrard & Henry 1991, Kaiser 1991; Bower 1997). For a fixed amount of specific heating, the effect is larger for poorer systems, i.e. when the extra energy per gas particle is comparable to the halo virial temperature. This produces both an excess entropy and a steeper $`L_X`$$`T`$ relation (e.g., Cavaliere, Menci & Tozzi 1998; Balogh, Babul & Patton 1999; Tozzi & Norman 2001). Arguments based on semi–analytical work (e.g., Tozzi & Norman 2001) and numerical simulations (Bialek, Evrard & Mohr 2001; Brighenti & Mathews 2001; Borgani et al. 2001b, 2002) suggest that a specific heating energy of $`E_h1`$ keV/part or, equivalently, a pre–collapse entropy floor of $`S100`$ keV cm2, can account for the observed $`X`$–ray properties of galaxy systems (cf. also Babul et al. 2002, Finoguenov et al. 2002a for arguments suggesting a stronger pre–heating). Yet, the origin for this energy has still to be determined. Energy release from supernovae feedback has been advocated as a possibility (e.g., Bower et al. 2001; Menci & Cavaliere 2001). Using the abundance of heavy elements of the ICM as a diagnostic for the past history of the star formation within clusters (e.g., Renzini 1997; Kravtsov & Yepes 2000; Pipino et al. 2002; Valdarnini 2002), a number of studies concluded that SN may fall short in providing the required extra–energy budget (cf. also Finoguenov, Arnaud & David 2001a). The other obvious candidate is represented by energy from AGN (e.g., Valageas & Silk 1999; Wu, Fabian & Nulsen 2000; Mc Namara et al. 2000; Nath & Roychowdhury 2002; Cavaliere, Lapi & Menci 2002). In this case, the large amount of energy that is available requires some degree of tuning of the mechanisms responsible for its conversion into thermal energy of the gas. While a suitable amount of non–gravitational heating can account for the observed $`L_X`$$`T`$ relation and entropy excess, the $`M`$$`T`$ relation is only marginally affected by extra heating (e.g. BGW), thus leaving the discrepancy between observed and predicted relation unresolved.

As for cooling, its effect is to selectively remove those low–entropy particles from the diffuse $`X`$–ray emitting phase which have cooling times shorter than the Hubble time (e.g., Voit & Bryan 2002; Wu & Xue 2002). Conversion of cooled gas into collisionless stars decreases the central gas density and, at the same time, the resulting lack of pressure support causes higher–entropy shocked gas to flow in from the outskirts of the cluster or group. As a result, the $`X`$–ray luminosity is suppressed, while the entropy increases, much like in a pre–heating scenario (Pearce et al. 2001; Muanwong et al. 2002; Davé, Katz & Weinberg 2002). However, by its nature, cooling is known to be a runaway process: cooling causes gas to be accumulated into dense structures, and the efficiency of cooling increases with gas density. As a result, most simulations consistently predict a significant fraction of gas to be converted into cold “stars”, $`f_{\mathrm{cold}}>\text{ }30`$ per cent (e.g., Suginohara & Ostriker 1998; Lewis et al. 2000; Yoshida et al. 2002; BGW), while observations indicate a considerably lower value of $`f_{\mathrm{cold}}<\text{ }10`$ per cent (e.g., Balogh et al. 2001; Wu & Xue 2002).

This suggests that in real clusters some source of extra heating is increasing the entropy of the gas, preventing overcooling. Voit et al (2002) have developed a semi–analytical approach to derive $`X`$–ray observable properties of the ICM in the presence of both cooling and extra heating. Based on this approach, these authors found that cooling and a modest amount of extra heating are able to account for basically all the $`X`$–ray ICM observables. Oh & Benson (2002) pointed out that pre–heating is needed to increase the cooling time and prevent overcooling, by suppressing the gas supply to galaxies (see also Finoguenov et al. 2002b). It is however clear that, as for any analytical approach, suitable assumptions and approximations are needed to choose criteria for removing cooled gas from the hot diffuse phase, and to follow the complex dynamics of cooling/heating of gas during the process of cluster formation.

Muanwong et al. (2002) and Kay, Thomas & Theuns (2002) used hydrodynamical simulations within a cosmological box to study the interplay of gas cooling and a few prescriptions for non–gravitational heating. As a general result, they found that increasing the heating can suppress the amount of cooled gas. While the choice of simulating a whole cosmological box has the advantage of providing a large statistics of groups and clusters, it also severely limits the available mass and force resolution. On the other hand, by the very nature of cooling, increasing the mass resolution allows to follow the formation of smaller halos at progressively larger redshift, where cooling and, potentially, star formation are particularly efficient. As a consequence, unless very high mass resolution is achieved, cooling in simulations can be significantly underestimated (e.g., Balogh et al. 2001).

In this paper, we follow the alternative approach of simulating at very high resolution a limited number of group– and cluster–sized halos selected from a cosmological box, and we widen the explored range of possible patterns for non–gravitational heating (see also BGW). While this limits our ability to precisely calibrate shape and scatter of $`X`$–ray scaling relations, we are able to increase the resolution in the most interesting regions of the gas distribution. Indeed, the simulations presented in this paper are among the highest resolution attempts realized so far to follow the structure of gas cooling within groups and clusters in the presence of a variety of schemes for extra gas heating. Furthermore, we also investigate how the cooling efficiency depends both on numerical resolution and on details of the SPH implementation.

The structure of this paper is as follows. After providing a short description of the code, we present in Section 2 the procedure to simulate individual halos at high resolution and discuss the main characteristics of the four selected halos. In Section 3, we discuss the results on the cold fraction. Here we will concentrate on showing how this fraction depends on numerical resolution, integration scheme and removal of cold dense particles from the SPH computation (star formation). Finally, we present the adopted schemes for non–gravitational gas heating and discuss their impact on the resulting cold fraction and pattern of star formation. In Section 4, we present the predictions on $`X`$–ray properties of clusters and groups from our simulations, namely the entropy–temperature, the luminosity–temperature and the mass–temperature relations. Finally, we discuss our main results and draw conclusions in Section 5.

## 2 The simulations

### 2.1 The code

Our simulations are realized with GADGET111http://www.mpa-garching.mpg.de/gadget, a parallel tree N–body/SPH code (Springel, Yoshida & White 2001), with fully adaptive time–step integration. Gas cooling in the SPH part of the code is implemented following Katz, Weinberg & Hernquist (1996, KWH hereafter). Specifically, the abundances of ionic species are computed by assuming collisional equilibrium for a gas of primordial composition (mass–fraction $`X=0.76`$ of hydrogen and $`1X=0.24`$ of helium). Since we not follow metal production from star–formation, we do not include the effect of metals on the cooling function. We include the effect of a time–dependent uniform UV background (e.g., Haardt & Madau 1999), although its effect is only very small for the massive objects we focus on in this study. We set the number of neighbors for SPH computations to 32, allowing the SPH smoothing length to drop at most to the value of the gravitational softening length of the gas particles.

### 2.2 The simulated structures

We simulate four halos at high resolution, which are extracted from a low–resolution DM only simulation within a box of $`70h^1\mathrm{Mpc}`$ on a side, for a cosmological model with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, Hubble constant $`H_0=70`$ km s-1 Mpc-1 and normalization $`\sigma _8=0.8`$, consistent with recent determinations of the number density of nearby clusters (Pierpaoli et al. 2002, and references therein). As for the baryon content, we assume $`\mathrm{\Omega }_{\mathrm{bar}}=0.019h^2`$ (e.g., Burles & Tytler 1998). This choice of $`\mathrm{\Omega }_{\mathrm{bar}}`$ corresponds to $`f_{\mathrm{bar}}0.13`$ for the cosmic baryon fraction, which, for the assumed cosmology, is consistent with the value measured from cluster observations (e.g., Ettori 2002, and references therein).

The most massive halo we selected corresponds to a Virgo–like cluster, with virial mass of about $`4\times 10^{14}M_{}`$ (as usual, we call “virial” the mass within the radius encompassing the virial overdensity computed for the simulated cosmology; e.g. Eke et al. 1998a). This turns out to be the most massive system extracted from the simulation box. In the following, we will refer to this system as the “Virgo” cluster. The other three halos, which have been extracted from a single Lagrangian region, correspond to groups in the mass range (2–6)$`\times 10^{13}M_{}`$. In the following, we will refer to these three structures as “Group-1”, “Group-2” and “Group-3”. We provide in Table 1 the main characteristics of the simulated structures.

We follow the technique originally presented by Katz & White (1993) to increase the mass resolution and to add short wavelength modes within Lagrangian regions that contain the structures of interest. In these high–resolution regions, particles are split into a dark matter and a gaseous part, with mass ratio reflecting the value of the cosmic baryon fraction. Force and mass resolution are then gradually degraded in the outer regions, so as to limit the computational cost, while providing a correct representation of the large–scale tidal field. The size of the regions selected at $`z=0`$, to be resimulated at high resolution, typically corresponds to 10-20 Mpc in Lagrangian space, and is always chosen to be large enough that no low-resolution heavy particles contaminate the virial region of the simulated halos.

In order to assess numerical effects, structures have been simulated at different mass and force resolutions. We fix three different mass resolutions, which correspond to $`m_{\mathrm{gas}}2.5\times 10^9M_{}`$, $`3.2\times 10^8M_{}`$ and $`3.9\times 10^7M_{}`$ for the mass of the gas particles, respectively. In the following, the group runs with the smallest (intermediate) value of $`m_{\mathrm{gas}}`$, and the Virgo runs with the intermediate (largest) $`m_{\mathrm{gas}}`$ will be indicated as high–resolution (low–resolution) runs and labeled with HR (LR). We do not discuss Virgo runs with the smallest $`m_{\mathrm{gas}}`$ and Group runs with the largest $`m_{\mathrm{gas}}`$ among this list of three mass resolutions. With these choices for the mass resolution, the HR runs resolve the virial regions of the simulated structures with a number of gas particles ranging from about 70,000 to about 185,000 (see Table 1). The redshift $`z_i`$ at which initial conditions are generated is chosen such that the r.m.s. fluctuation in the density field of the high–resolution region is $`\sigma =0.1`$ (on the scale of the smallest resolved masses). With this requirement, we have $`z_i65`$. As for the choice of the softening scale for the computation of the gravitational force, we assume it to have a constant value in comoving units down to $`z=2`$, and a constant value in physical units at later epochs. The corresponding values of the Plummer–equivalent softening scale at $`z=0`$ are $`ϵ_{\mathrm{Pl}}=10`$, 5 and 2.5 kpc for the three different choices of $`m_{\mathrm{gas}}`$. This choice has been dictated by the requirement of resolving halos down to scales of about one percent of their virial radii, so as to correctly follow the gas clumpiness and, therefore, to have convergent estimates of the X–ray luminosity (e.g., Borgani et al. 2002).

We show in Figure 1 the gas density and entropy maps for the HR runs of the cluster and of the group regions at $`z=0`$ for the runs including cooling and star formation (see below). In the entropy map of the Virgo cluster (lower left panel), we note a tail of low entropy gas pointing toward the center. This feature is generated by a merging sub-group, whose gas has been tidally stripped during the first passage through the cluster virial region. The persistence of low entropy for this gas indicates that it has been only recently stripped and has still to thermalize within the cluster environment. As apparent from the gas–density map (upper left panel), this merging sub-halo gives rise only to a minor disturbance of the gas density, thus marginally disturbing the relaxed dynamical status of the cluster. As for the simulation of the region containing the three groups, we note that they are placed along a filamentary structure. Although they are still relatively isolated and separated from each other by a few virial radii, their motion shows that they are approaching each other and will merge to form a cluster–sized structure in a few Gyrs. In general, these maps witness that a rich variety of structures, emerging when high resolution is achieved, are naturally expected to characterize the ICM, much like shown by high resolution Chandra observations.

## 3 Computing the collapsed gas fraction in cluster simulations

### 3.1 Introducing radiative cooling

An important aspect when dealing with simulations that include cooling concerns the detailed scheme of SPH implementation. Most standard implementations integrate the specific thermal energy as an independent variable, differing however in the detailed method used for symmetrizing the pairwise hydrodynamic forces between gas particle pairs, where either a simple arithmetic or a geometric mean form the most common choices (e.g., Weinberg, Hernquist & Katz 1997; Davè et al. 1999; White, Hernquist & Springel 2001). While these SPH implementations conserve energy and momentum, Springel & Hernquist (2002, SH02 hereafter) have shown that several of the commonly used SPH implementations are characterized by a spurious loss of specific entropy in strongly cooling regions, an effect which can be particularly severe at low resolution, and which is stronger when the geometrical scheme for hydrodynamical force symmetrization is adopted. This problem is essentially due to spurious coupling between cool dense particles, which should have virtually left the collisional phase, and neighboring hot gas particles, which still belong to the diffuse phase. In order to avoid the resulting spurious overcooling, different techniques have been suggested by several authors (e.g., Pearce et al. 2001; Marri et al. 2002).

SH02 proposed a new SPH implementation based on integrating the specific entropy as an independent thermodynamic variable, an approach which explicitly conserves entropy in non–shock regions. Using a variational principle to derive the SPH equations of motion, they also showed that this new formulation removes any ambiguity in the choice of symmetrization and conserves energy, even when adaptive smoothing lengths are used.

Since one of the main purposes of this paper is to investigate the properties of gas cooling in galaxy clusters, we will study below by how much differences in the SPH implementation can change the resulting fraction of cold gas. Adopting the naming convention of SH02, we refer to a standard SPH implementation with geometric symmetrization as “geometrical”, and to one with arithmetic symmetrization as “arithmetic”, while the new formulation of SH02 will be referred to as “entropy–conserving”.

### 3.2 Introducing star formation

Star formation is introduced as an algorithm to remove dense cold gas particles from the SPH computation, treating them as collisionless “stars”. We follow the recipe originally proposed by KWH. According to this recipe, a gas particle is eligible to form stars if the following conditions are met: (i) locally convergent flow, $`𝐯<0`$; (ii) Jeans unstable, i.e. locally determined sound crossing time longer than dynamical crossing time; (iii) gas overdensity exceeding a critical overdensity value, $`\delta _g>55`$; (iv) local number density of hydrogen atoms $`n_H>0.1`$ cm-3.

Once a particle is eligible to form stars, its star formation rate (SFR) is given by $`\mathrm{d}\mathrm{ln}\rho _g/\mathrm{d}t=c_{}/t_g`$, where $`t_g`$ is the minimum between the local gas–dynamical time–scale, $`t_{\mathrm{dyn}}=(4\pi G\rho _g)^{1/2}`$, and the local cooling time–scale. We assumed $`c_{}=0.1`$ for the parameter regulating the rate of conversion of cold gas into stars, and verified with a low–resolution simulation of the Virgo cluster that basic results are left essentially unchanged by taking instead $`c_{}=0.01`$. A gas particle eligible for star formation is assumed to be gradually converted into a star particle, according to the above SFR. Instead of creating a new star particle for every star–formation (SF) instance, each gas particle undergoing SF behaves in a “schizophrenic” way, with its stellar part feeling only gravity (see, e.g., Mihos & Hernquist 1994). Once the SPH mass fraction decreases to 10 per cent, it is dissolved into SPH neighbors, thus leaving a purely stellar particle.

We also follow the recipe by KWH to compute the energy feedback from the SN associated with the star–formation produced in the simulations. Assuming a Miller–Scalo (1979) initial mass function (IMF), we compute the number of stars with mass $`>8M_{}`$, which we identify with instantaneously exploding SN. After assuming that each SN releases $`10^{51}`$ ergs, the resulting amount of energy per formed stellar mass turns out to be $`7\times 10^{48}`$ ergs $`M_{}^1`$. While the approximation of instantaneous explosion can be justified for type-II SN, due to the short life–time of their progenitor stars, it is not valid for type-Ia SN, which have stellar progenitors of smaller masses and much longer life–times (e.g. Lia, Portinari & Carraro 2001; Pipino et al. 2002; Valdarnini 2002). The resulting feedback energy is assigned as thermal energy to the star–forming gas particles. This scheme for SN feedback is known to thermalize a negligible amount of energy in the diffuse medium, since it acts mostly on cold dense particles which rapidly radiate away the feedback energy as a consequence of their short cooling time. In the following, we show results based on including this scheme for SN feedback while bearing in mind that it causes only negligible differences compared to simulations that lack any stellar feedback. In Section 3.3 we shall discuss a different SN feedback scheme, based on the predictions of semi–analytical modelling of galaxy formation.

The effect of including star formation on the fraction of collapsed gas in the Virgo cluster run is shown in Figure 2. Besides the population of collisionless star particles, we also define as belonging to the cold phase all the SPH particles which have overdensity $`\delta _{\mathrm{gas}}>500`$ and temperature $`T<3\times 10^4`$ K (see also Croft et al. 2001, Borgani et al. 2002). At $`z=1`$, the star formation simulation produce $`f_{\mathrm{cold}}25`$–35 per cent of collapsed gas, with a weak dependence on the integration scheme. However, in the cooling-only simulations, where collapsed gas is not converted to stars, the “geometric” scheme leads to substantially larger values, indicating a numerical overcooling problem in this method. This effect is absent in the entropy–conserving scheme, which proves effective in suppressing spurious cooling in the absence of an explicit SF scheme. Note that $`f_{\mathrm{cold}}`$ is seen to slightly decrease at later epochs, which is as a consequence of a reduction of the rate at which cooling and star formation proceeds with respect to diffuse gas accretion.

In Figure 3 we show the trend of $`f_{\mathrm{cold}}`$ against the emission–weighted temperature, which is defined here as

$$T_{ew}=\frac{\underset{i}{}\rho _iT_i^{3/2}}{_i\rho _iT_i^{1/2}},$$
(1)
where $`\rho _i`$ is the gas density carried by the $`i`$-th SPH particle. This definition is valid only for bremsstrahlung emissivity, although final values of $`T_{ew}`$ are left essentially unchanged if we account for the contribution from metal lines.

For the groups, we plot results for both the LR and the HR runs. The trend toward a higher $`f_{\mathrm{cold}}`$ in colder systems is a consequence of the shorter cooling times of the associated DM halos, which makes cooling to proceed faster in lower–mass systems. Quite apparently, increasing the resolution causes a $`20`$ per cent or a $`50`$ per cent increase of $`f_{\mathrm{cold}}`$ at the group and at the cluster scales, respectively. The effect of numerical resolution is also shown in the left panel of Figure 4. In this figure we plot the density of star–formation rate (SFR) within the virial region of the Virgo cluster and of the Group-1. Once the same mass resolution is used for the simulation of these two structures, the resulting SFRs are quite similar. Increasing the resolution of the group simulation allows to resolve smaller halos forming at higher $`z`$, where additional star formation can take place. As a result, the peak in the SFR moves from $`z2`$ to $`z3`$ and then declines more gently, while recovering the same shape at lower redshift. Note that the integrated star formation rate is dominated by the contribution from these low redshifts, where most of the physical time is being spent. The resolution achieved in the HR runs is sufficient to resolve “galaxy” halos well below $`L_{}`$. Therefore, we are confident that we are obtaining nearly converged estimates of the collapsed gas fraction, at least when the highest mass resolution is used. At the same time, our results should be considered as a warning on the interpretation of simulations that lack the resolution to adequately follow gas cooling.

In summary, our simulations demonstrate that cold fractions as large as $`f_{\mathrm{cold}}=25`$–35 per cent should be expected when radiative cooling and star formation are considered. These values are larger than the observed $`10\%`$ fraction of cold gas in clusters (e.g., Balogh et al. 2001). This calls for the need to introduce a suitable scheme of non–gravitational energy injection, allowing a regulation of the runaway cooling process.

### 3.3 Introducing extra heating

The SN feedback recipe that we discussed so far, where thermal energy is deposited into cold gas, does not produce any sizable effect. In order to overcome this problem, many different schemes have been proposed. All these schemes attempt to prevent feedback energy from being quickly radiated away, for example by suitably parameterizing “sub–grid” physics, such as the multi-phase structure of the interstellar medium or galactic winds (e.g., Kay et al. 2000; Springel & Hernquist 2002; Marri & White 2002).

Here we present different phenomenological approaches for non–gravitational heating. Rather than predicting feedback from the star formation actually produced in the simulations, these schemes are designed to shed light on how much extra energy is required and how it should be distributed in redshift and as a function of the local gas density, to prevent overcooling and, at the same time, to reproduce $`X`$–ray observables of galaxy clusters and groups. A summary of the characteristics of the heating recipes that we explore here is provided in Table 2. In the rest of the paper we will present results based only on the high–resolution (HR) runs.

#### 3.3.1 Impulsive heating

In our first class of heating schemes, we assume that all the energy is dumped into the diffuse baryons in an impulsive way, with a single heating episode occurring at some redshift $`z_h`$.

Entropy floors $`S_{\mathrm{fl}}=50`$ keV cm2 at redshift $`z_h=9`$ (S50-9 runs) and at $`z_h=3`$ (S50-3 runs), and $`S_{\mathrm{fl}}=25`$ keV cm2 at $`z_h=9`$. In this scheme, the entropy associated with each gas particle, $`s=T/n_e^{2/3}`$ ($`T`$: temperature in keV; $`n_e`$ electron number density in cm-3), is either increased to the value $`S_{\mathrm{fl}}`$ if smaller than that, or otherwise left unchanged (see also Navarro et al. 1995; Bialek et al. 2001; BGW). The choice of $`z_h=9`$ corresponds to a heating epoch well before a substantial amount of gas in simulations cools and forms stars and, therefore, heavily suppresses star formation. The existence of a pristine SN generation (from the so–called Pop III stars) has been invoked to account for the IGM metal enrichment (e.g., Madau, Ferrara & Rees 2001). However, were this heating able to rise the entropy to the above levels, it would prevent the later formation of the Ly-$`\alpha `$ forest, which is known to have about one order of magnitude lower entropy. Furthermore, the amount of heating energy would also correspond to a too high production of heavy elements. For these reasons, we consider this choice for $`S_{\mathrm{fl}}`$ to be motivated by the phenomenology of X–ray ICM properties alone, rather than by expectations from star formation processes at high redshift.

A fixed amount of heating energy per particle, $`E_h=0.75`$ keV/particle at $`z_h=3`$. This amount of energy is roughly the same as the average specific energy dumped by the S50-3 scheme within the halo virial radius (see Table 2). Therefore, it allows to check for differences induced in the final results by distributing the same energy budget in a different way as a function of gas density. The $`z_h=3`$ heating epoch is close to that at which the star-formation rate within a proto–cluster region peaks (e.g., Menci & Cavaliere 2000; Bower et al. 2001; BGW). An energy budget $`E_h0.6`$–0.8 keV/part has been also suggested by Finoguenov et al. (2001a) to be consistent with the Si abundance detected in groups and clusters.

#### 3.3.2 SAM–predicted SN feedback

This heating scheme is based on computing the star–formation rate (SFR) within clusters using a semi–analytic model (SAM) of galaxy formation (e.g., Kauffmann, White & Guiderdoni 1993, Somerville & Primack 1999, Cole et al. 2000, and references therein). Here we employ a variation of the scheme described by Menci & Cavaliere (2000, see also Bower et al. 2001), and we refer to their paper for a detailed description of the method, while we refer to BGW for further details on the its implementation in cluster simulations.

The hierarchical merging of DM halos is followed by means of the extended Press–Schechter formalism (e.g. Lacey & Cole 1993), while model parameters describing the gas physics, such as cooling, star formation and stellar feedback, are chosen so as to reproduce observed properties of the local galaxy population, such as the Tully–Fisher relation, or optical luminosity functions and disk–sizes (e.g. Poli et al. 2001). The model prediction we are interested in here is the integrated star formation history, $`\dot{m}_{}(z,M_0)`$ ,of all the condensations which are incorporated into a structure of total mass $`M_0`$ by the present time. For a halo of size similar to our Virgo–like cluster, the SFR peaks at $`z4`$ in this semi-analytic model, while it is $`z2.5`$–3 for the group–sized halos (see BGW, for a plot of the $`M_0`$–dependence of the cluster SFR).

The rate of total energy feedback released by type-II SN is then computed as

$$\frac{\mathrm{d}E_{SN}}{\mathrm{d}t}=\mathrm{\hspace{0.17em}10}^{51}\mathrm{ergs}\eta _0\dot{m}_{}(z,M_0),$$
(2)
where $`\eta _0`$ is the number of SN per solar mass of formed stars. This value depends on the assumed initial mass function (IMF), and is obtained by integrating over the IMF for stellar masses $`>8M_{}`$. In the following, we will use the values $`\eta _0=3.2\times 10^3M_{}^1`$, which follows from a Scalo IMF (Scalo 1986), $`\eta _0=7\times 10^3M_{}^1`$, from the Salpeter IMF (Salpeter 1955) and $`\eta _0=1.5\times 10^2M_{}^1`$, as an extreme case.

Since our simulations include cooling, radiative losses of SN energy do not need to be assumed a priori, rather they are self–consistently computed by the code. However, we need to specify the gas overdensity, $`\delta _g`$, at which the SN heating energy is assigned to the gas. In the following we take $`\delta _g=50`$ or 500, and assume that $`E_{\mathrm{SN}}`$ is shared in equal parts among all the gas particles at overdensity larger than $`\delta _g`$. The choice $`\delta _g=50`$ corresponds to assuming that the virial region of the whole halo is heated and, therefore, that physical processes like galactic winds, for example, are rather efficiently transferring energy to the IGM. Increasing $`\delta _g`$ implies two competitive effects: on one hand, it decreases the number of heated gas particles, therefore it increases the amount of extra energy assigned to each of them; on the other hand, the energy is assigned to denser particles, which have shorter cooling time and, therefore, larger radiative losses.

### 3.4 The effect of extra heating on the cold fraction

As we have already discussed, introducing cooling causes a too large fraction of gas to be converted into stars. This is a well known feature of hydrodynamical simulations, which has been widely discussed in the literature (e.g., Suginohara & Ostriker 1998, BGW). Even worse, the runaway nature of the cooling process causes its efficiency to be highly sensitive to numerical resolution (see Fig.3, see also Balogh et al. 2001). Therefore, one should be very cautious in the interpretation of results from simulations that do not resolve halos with luminosity well below $`L_{}`$.

In Figure 5, we show the effect of the different heating schemes on the resulting cold fraction within the virial radius of our simulated structures. As expected, we find a decrease of $`f_{\mathrm{cold}}`$ when non–gravitational heating is included. However, the efficiency of this suppression of star–formation does not exclusively depend on the amount of dumped energy. For instance, imposing an entropy floor of 50 keV cm2 at $`z_h=9`$ (left panel of Fig.5) is far more efficient than at $`z_h=3`$. The reason for this is illustrated by the different patterns of SFR history, that we show in Figure 4. The impulsive heating at $`z_h=3`$ causes a suppression of the SFR at later epochs, but a fair amount of stars are already in place at $`z_h`$ (central panel of Fig. 4). Quite interestingly, the results for the two runs with heating at $`z_h=3`$ produce quite similar SFR. This indicates that, once the heating epoch is fixed, the degree of SFR suppression depends only on the amount of heating energy, while being largely independent on its distribution as a function of the local gas density.

However, heating at $`z=9`$ with a comparable amount of energy does not allow gas to reach high densities within DM halos and to cool before $`z1`$. Once cooling takes place, it converts less than $`10`$ per cent of the gas into stars, within a short episode of star formation. The resulting SFR peaks at very low redshift, $`z0.3`$, which is highly discrepant with observational determinations of the SFR history in clusters (e.g., Kodama & Bower 2001). Of course, this is not the only feature which rules out the picture of a strong heating occurring at such a high redshift. For instance, at $`z=0`$, stars are all concentrated in one single object located at the center, a cD–like galaxy, while no other galaxy–sized DM halos contains significant amounts of collapsed gas.

As for the SN heating (right panel of Fig. 4), suppressing the star fraction below the 10 per cent level requires a high, probably unrealistic value for $`\eta _0`$. Also, taking $`\eta _0=1.5\times 10^2M_{}^1`$ generates an implausible SFR history, resembling that found for the runs based on setting the entropy floor at $`z_h=9`$: the large amount of extra heating at high redshift prevents the occurrence of star formation down to $`z1.5`$. Taking $`\eta _0`$ in the range 3–7$`\times 10^3`$ predicts more realistic SFRs, but it is not able to suppress $`f_{\mathrm{cold}}`$ below the $`20`$ per cent level.

A general conclusion of our analysis is that heating schemes producing plausible SFR histories are not efficient in suppressing the fraction of cold gas below the 20 and 25 per cent values at the cluster and group scales, respectively. Vice-versa, a more efficient suppression is obtained by preventing gas to cool at high redshift, at the expense of delaying star formation to unreasonably low redshifts.

## 4 $`X`$–ray properties of simulated clusters

### 4.1 The entropy of the ICM

Measurements of the excess entropy in central regions of poor clusters and groups are considered to provide direct evidence for the lack of self–similarity of the ICM properties (e.g., Ponman et al. 1999; Finoguenov et al. 2002a). In a separate paper (Finoguenov et al. 2002b), a self–consistent comparison is realized between the entropy properties of the simulations with impulsive heating, that we present here, and the observational data for groups and clusters by Finoguenov et al. (2002a). The main result of this comparison is that, although cooling and star formation tend to somewhat increase entropy in central cluster regions, they still fall short in producing the entropy excess which is observed at the group scale. While preheating at $`z_h=3`$ is shown to increase the entropy to the observed values, runs with $`z_h=9`$ are characterized by a low entropy level in central regions of clusters and groups.

Instead of attempting any further comparison with observations, we want to discuss here the dynamical reasons for such a behavior. To this end, we show in Figure 6 the effect of cooling and non–gravitational heating on the entropy profiles for our whole set of “Virgo” simulations. As expected, when cooling and SF are included, low entropy gas is selectively removed in central cluster regions, thus inducing a flattening of the profile. This is explicitly shown in Figure 7: while the run including only gravitational heating has a population of high–density low entropy gas particles, such particles are removed from the diffuse phase once cooling and star formation are introduced. This result is consistent with the expectation from analytical arguments based on the comparison between cooling time–scale and typical cluster age (e.g., Voit et al. 2002, Wu & Xue 2002). The inclusion of extra heating has a non–trivial effect on the efficiency of cooling in removing particles from the lower left side of the $`S`$$`\delta _g`$ phase diagram. For instance, imposing the same entropy floor at $`z_h=9`$ and at $`z_h=3`$ has quite different effects on the entropy pattern (see lower panels of Fig. 7). Heating at $`z_h=9`$ has the effect of increasing the cooling time for most of the gas particles, so as to allow star formation to take place only quite recently (see Fig. 4). The increased time-scale for cooling causes this process to proceed in a more gradual way. For this reason, the entropy of gas particles undergoing cooling decreases slowly, thus making their removal from the hot diffuse phase less efficient.

### 4.2 The luminosity–temperature and luminosity–mass relations

The slope of the $`L_X`$$`T`$ relation also provides important observational evidence for the lack of self–similar behaviour of the ICM. Since the first measurements of ICM temperatures for sizable sets of clusters, it has been recognised that $`L_XT^\alpha `$ with $`\alpha 3`$, although with a considerable scatter (e.g., White et al. 1998, and references therein). Better quality observations established that a significant contribution to this scatter is associated with the different strength of cooling flows detected in different clusters. Either excluding clusters with pronounced signatures of cooling flows or correcting for their effect (e.g., Markevitch 1998; Allen & Fabian 1998; Arnaud & Evrard 1999; Ettori et al. 2002) results in a much tighter $`L_X`$$`T`$ relation, albeit still with a rather steep slope. At the same time, hints have also been found for a further steepening of this relation at $`T<\text{ }1`$ keV (e.g., Helsdon & Ponman 2000, and references therein), possibly indicating that the mechanism responsible for the $`L_X`$$`T`$ scaling should act in a different way for clusters and groups.

Simulations that allow for non–gravitational heating (e.g., Bialek et al. 2001; BGW) and radiative cooling (e.g. Pearce et al. 2001; Davè et al. 2002; Muanwong et al. 2002) have been shown to be able to account for the observed $`L_X`$$`T`$ relation. However, a sometimes overlooked issue in determining the $`X`$–ray luminosity of clusters in simulations concerns the contribution of metal lines to the emissivity. While this contribution is negligible above 2 keV, it becomes relevant at the scale of groups. For instance, neglecting the contribution from line emissivity for an ICM enriched to a metallicity of $`Z=0.3Z_{}`$ leads to an underestimate of the X–ray luminosity by almost 50 per cent at 1 keV, and by more than a factor 2 at 0.5 keV (e.g. BGW).

A correct procedure would require simulations that include a treatment of metal enrichment and a self–consistent estimate of the contribution of line cooling to the X–ray emissivity. However, only preliminary attempts have been realized so far to include the treatment of ICM metal enrichment from SN ejecta (Lia et al. 2002; Valdarnini 2002). Pearce et al. (2001) include the contribution of metals to the cooling function adopted in their simulation by assuming $`Z=0.3Z_{}`$ at the present epoch, linearly decreasing with time towards the past. Davé et al. (2002) did not include the metal contribution in their cooling function, but estimated $`X`$–ray luminosities by assuming a phenomenological relation between metal abundance and temperature of the galaxy system. However, while the ICM metallicity at the scale of rich clusters is quite well established from observations, the situation is less clear for poor clusters and groups (e.g., Davis, Mulchaey & Mushotzky 1999; Renzini 2000, and references therein).

The cooling function used in our simulations assumes zero metallicity, but we compute the $`X`$–ray luminosity by adding to the bremsstrahlung emissivity the contribution from lines for a $`Z=0.3Z_{}`$ plasma. This represents a reasonable approximation as long as gas spends most of the time at low metallicity, being enriched to high metallicity only recently. Owing to the uncertainties connected to these assumptions, the reliability of $`L_X`$ values at $`T<\text{ }1`$ keV is unclear, however. Precise predictions will require a fully self–consistent treatment of metal enrichment of the ICM from star formation activity.

In Figure 8, we show the profiles of emissivity (energy released per unit time and unit volume) for the different Virgo runs. As expected, including only cooling and star formation has the effect of flattening the profiles in the central cluster regions as a consequence of gas removal from the hot phase. When extra heating is included, the profiles change according to the amount of gas left at relatively low entropy in the central cluster regions. For instance, the fairly large population of low entropy particles in the run with $`S_{\mathrm{fl}}=50`$ keV cm2 at $`z_h=9`$ (see Fig.7) is responsible for the spike in the $`X`$–ray emissivity. In the same way, the efficient removal of low–entropy gas for the run where an entropy floor was imposed at $`z_h=3`$ explains the flattening of the luminosity profile in the central cluster region. These results confirm the existence of a non–trivial interplay between the effects of cooling and extra heating. In some cases, one reaches the apparently paradoxical conclusion that combining heating and cooling increases the $`X`$–ray luminosity, although their separate effects are that of suppressing $`L_X`$.

Figure 9 shows the comparison between the simulated and the observed $`L_X`$$`T`$ relation for clusters and groups. As expected, cooling causes a sizeable suppression of the $`X`$–ray luminosity. At the same time, the emission–weighted temperature is increased as a consequence of the steepening of the temperature profiles in the central halo regions (see below). While the mere introduction of cooling and star–formation brings the “Virgo” cluster into agreement with observations, the simulated groups are somewhat overluminous with respect to data. The inclusion of pre–heating at $`z_h=3`$ has a smaller effect on the Virgo cluster, consistent with the result from the luminosity profile, while it further suppresses $`L_X`$ at the scale of groups. A similar result is also found for the runs with SAM–predicted SN feedback.

Quite interestingly, the $`L_X`$ value for Group-2 in the runs with no extra heating appears to be systematically in excess with respect to that inferred from the $`L_X`$$`T`$ scaling of the other three simulated structures. This deviation is due to the occurrence of a recent merger shock in the Group-2 run, which produced a sudden increase in the $`X`$–ray emission. When extra heating is included, its effect is that of decreasing the strength of the shock, thus also reducing the jump in luminosity.

A similar constraint is provided by the relation between $`X`$–ray luminosity and mass. Reiprich & Böhringer (2002) have estimated this relation by applying the equation of hydrostatic equilibrium to a fairly large ensemble of clusters and groups, under the assumption of isothermal gas. They used ICM temperatures based on ASCA data, in combination with ROSAT-PSPC data for the surface brightness profile. Ettori, De Grandi & Molendi (2002) used the better quality data from Beppo–SAX observations to resolve the temperature profiles for a smaller ensemble of clusters. Although the analysis by Ettori et al. explicitly includes temperature gradients when solving the equation of hydrostatic equilibrium, it is restricted to clusters with $`T>\text{ }3`$ keV, thus hotter than those simulated here. For this reason, we here compare our simulation results to the data by Reiprich & Böhringer (see Figure 10). In this analysis, the cluster masses, $`M_{500}`$, are computed within the radius encompassing an average density $`\overline{\rho }=500\rho _{\mathrm{crit}}`$, while observed luminosities are provided in the 0.1–2.4 keV ROSAT energy band. We use the MEKAL spectral model to correct bolometric luminosities from simulations by assuming $`Z=0.3Z_{}`$ for the global ICM metallicity. Consistent with the results from the analysis of the $`L_X`$$`T`$ relation, we find that the runs with heating at $`z_h=3`$ and that with SN feedback, based on a Salpeter IMF, are able to follow the steep slope of the observed $`L_{0.12.4}`$$`M_{500}`$ relation.

In principle, the $`L_X`$$`T`$ and the $`L_X`$$`M`$ relations do not provide independent information, since masses are anyway estimated using temperature data. Still, both relations are obtained by using largely different observational data sets and analysis procedures. Therefore, the fact that the same simulations are able to account for both scalings lends support to the robustness of our results and indicates that our conclusions are not affected by observational biases or systematics.

Owing to the uncertainties mentioned above in modelling the luminosities of groups, it appears prudent not to make strong claims about how much extra heating is needed to reproduce the observations. Overall, we note that all the runs that produce a delayed star formation, such as those with $`z_h=9`$ and the one with SN–feedback and large $`\eta _0=1.5\times 10^2M_{}^1`$ (see Fig.4), are quite inefficient in suppressing $`L_X`$. On the contrary, runs with pre–heating at $`z_h=3`$ or with SN–feedback combined with more reasonable values for $`\eta _0`$ succeed to account for the steep slopes of the $`L_X`$$`T`$ and $`L_X`$$`M`$ relations.

Having warned about the reliability of the emissivity modeling for gas at $`T<1`$ keV, a word of caution should also be spent on the reliability of the interpretation of current observational data. Estimating temperature and luminosity for small groups from pre–Chandra and pre–XMM data is not a trivial task, mostly due to the difficulty of separating the contribution of the diffuse intra–group medium from that of member galaxies, and of detecting $`X`$–ray emission out to a significant fraction of the virial radius (see, e.g., Mulchaey 2000, for a review on the $`X`$–ray properties of groups). The situation is likely to improve as newer and better quality data will be accumulated, although we will probably have to wait for a few more years before a critical amount of Chandra and Newton–XMM observations of groups will be available.

### 4.3 The mass–temperature relation

Under the assumptions of spherical symmetry and an isothermal gas distribution, the condition of hydrostatic equilibrium predicts a precise relationship between the virial mass of a cluster and its temperature: $`k_BT=1.38\beta ^1M_{15}^{2/3}[\mathrm{\Omega }_m\mathrm{\Delta }_{\mathrm{vir}}(z)]^{1/3}(1+z)`$ keV for a gas of primordial composition, with $`M_{15}`$ being the virial mass in units of $`10^{15}h^1M_{}`$ and $`\mathrm{\Delta }_{\mathrm{vir}}`$ being the ratio between the virial density and the average cosmic matter density at redshift $`z`$. Under the above assumptions, the $`\beta `$ parameter gives the ratio between the specific kinetic energy of dark matter particles and the thermal energy of the gas. Simulations including only gravitational heating demonstrated that this relation is reproduced quite well, with $`\beta 1`$–1.2 (e.g., Evrard et al. 1996; Bryan & Norman 1998; Frenk et al. 1999; BGW). For these reasons, the $`M`$$`T`$ relation has been considered for several years as a fairly robust prediction of hydrostatic equilibrium: gas temperature, unlike $`X`$–ray emissivity, is primarily determined by the action of gravity and, as such, depends on global cluster properties, and only weakly on local structure of the ICM.

However, data based on ASCA and ROSAT observations show an $`M`$$`T`$ relation which is about 40 per cent lower than predicted (Horner et al. 1999; Nevalainen et al. 1999; Finoguenov et al. 2001b), a result which has been confirmed by Beppo–SAX (Ettori, De Grandi & Molendi 2002) and Chandra (Allen et al. 2001) data for relatively hot systems ($`T>\text{ }4`$ keV).

Non–gravitational heating could be naively expected to solve this discrepancy by increasing the ICM temperature at fixed cluster mass. However, BGW have shown that for a broad class of pre–heating models similar to those discussed here the $`M`$$`T`$ relation is left almost unchanged by the injection of extra–energy (cf. also Lin et al. 2002). In fact, as long as gas has time after being heated to settle back into hydrostatic equilibrium within the gravitational potential well, its temperature is mainly determined by the amount of collapsed dark matter, which is unaffected by the heating process.

An alternative explanation for the observed low amplitude of the $`M`$$`T`$ relation, based on the effect of radiative cooling, has been shown by Thomas et al. (2002) to be much more promising. In this case, gas left in the diffuse phase flows towards the central cluster region, where it is compressed, thus increasing its temperature. As a result, the overall mass–weighted temperature remains almost unchanged, but the emission weighted temperature significantly increases. Our results, as shown in Figure 11, actually confirm this picture and generalise it to a large range of schemes for extra heating: while the value of $`M_{500}`$ is left unchanged by the cooling/heating processes, $`T_{\mathrm{ew}}`$ increases as a consequence of the temperature increase in the central cluster regions.

In order to better understand the effect of cooling on the central temperature structure of the ICM, we plot in Figure 12 the gas pressure, $`P=\rho _{\mathrm{gas}}k_BT/(\mu m_p)`$, as a function of gas density for the simulations of the “Virgo” cluster. We introduce here the effective polytropic index $`\gamma =\mathrm{d}\mathrm{log}P/\mathrm{d}\mathrm{log}\rho _{\mathrm{gas}}`$ to describe the run of pressure as a function of gas density. In the external cluster region the gas is characterised by $`\gamma >\text{ }1`$, thus consistent with the slowly outward-declining temperature profiles, almost independent of the presence of cooling and extra heating. However, cooling leads to a loss of pressure in central high–density regions. As cooling partially removes low–entropy gas from the diffuse phase, gas of higher entropy flows in from more external regions. As long as this gas has sufficiently long cooling time, its entropy is conserved and the gas is adiabatically compressed during the inflow. In this regime, the effective polytropic index increases towards $`\gamma =5/3`$, thus indicating an adiabatic behaviour of the ICM. This result is essentially independent of whether gas is preheated. The only effect of imposing an entropy floor at $`z_h=9`$ is that of making the cooling process more gradual. This allows a larger amount of gas to remain in the diffuse phase, so as to reach higher density and higher pressure in central regions (see also Fig. 8).

### 4.4 The temperature profiles

The way in which cooling acts in reconciling the observed and the simulated $`M`$$`T`$ relations implies that temperature profiles should steepen in central cluster regions. From an observational viewpoint, the possibility of realizing spatially resolved spectroscopy has recently opened the possibility to determine temperature profiles for fairly large samples of clusters. Interestingly, observations based on the ASCA (e.g., Markevitch et al. 1998) and Beppo–SAX (De Grandi & Molendi 2002) satellites show declining temperature profiles in the outer regions, at cluster-centric distances $`>\text{ }0.2`$–0.3$`R_{\mathrm{vir}}`$ (cf. also Irwin & Bregman 2000). This behaviour is generally reproduced by simulations that do not include cooling (e.g., BGW). Furthermore, both Beppo–SAX (De Grandi & Molendi 2002), Chandra (e.g., Ettori et al. 2002; Allen et al. 2001; Johnstone et al. 2002) and XMM (e.g., Tamura et al. 2001) data for fairly hot systems, $`T_X>\text{ }4`$ keV, show temperature profiles declining towards the very central regions of clusters, thus indicating the presence of cooling cores. This behaviour is grossly at variance with respect to that found for the “Virgo” runs, as reported in Figure 13: the only case where a somewhat declining profile is produced is the one with gravitational heating, while cooling always gives rise to steeply increasing profiles with no evidence for any decline, independent of the presence of extra heating.

A more comprehensive comparison with the observations would require simulations to be realized for a set of clusters with higher temperature. On the other hand, our simulated Virgo cluster has been chosen as a fairly relaxed system. Therefore, as long as observations suggest profiles to be universal for such systems (Allen et al. 2001), such a discrepancy should be taken quite seriously. A steepening of the temperature profiles caused by cooling has been already noticed by Lewis et al. (2000), Muanwong et al. (2002) and Valdarnini (2002). The temperature profiles in Fig. 13 generalise this result also in the presence of a variety of extra–heating mechanisms.

We also note that the steep temperature profiles predicted by simulations are also at variance with respect to those predicted by the semi–analytical model for ICM heating/cooling by Voit et al. (2002). A detailed comparison between the predictions of semi–analytical models and simulations is beyond the scope of this paper. However, a full understanding of the physical processes taking place in the ICM will only be obtained if the reasons for such differences can be understood and eventually sorted out.

If the discrepancy between observed and simulated temperature profiles will be confirmed, it may indicate that we are missing some basic physical mechanism which affects the thermal properties of the gas in the high density cooling regions. For instance, thermal conduction has been advocated by some authors as a mechanisms that, in combination with central heating, may regulate gas cooling (e.g. Voigt et al. 2002) while providing acceptable temperature profiles for a suitable choice of the conductivity parameter (e.g., Zakamska & Narayan 2002; Ruszkowski & Begelman 2002). In this scenario, one expects the outer layers to heat gas in the innermost regions, so as to increase its cooling time, allowing it to stay in the diffuse phase at a relatively low temperature. However, the detection of sharp features in the temperature map of several clusters, as observed by the Chandra satellite, led some authors to suggest that thermal conduction is suppressed in the ICM (e.g., Ettori & Fabian 2000). Magnetic fields are naturally expected to produce such a suppression (e.g., Sarazin 1988). Still, it is not clear whether this mechanism can act in an ubiquitous way inside clusters or whether the turbulence associated with the presence of magnetic fields is actually able to maintain a relatively efficient thermal conduction (e.g., Narayan & Medvedev 2001).

## 5 Discussion and conclusions

We presented results from high resolution Tree+SPH simulations of a moderately poor “Virgo”–like cluster and of three group–sized halos, including the effects of radiative cooling and non–gravitational gas heating. The numerical accuracy reached in these simulations was aimed at following in detail the pattern of gas cooling and its effect on the $`X`$–ray properties of groups and clusters of galaxies. The main results that we obtained can be summarised as follows.

Including cooling and star formation causes a fraction $`f_{}0.25`$ of baryons to be converted into a collisionless “stellar” phase in the Virgo cluster and $`f_{}0.35`$–0.40 in the simulated groups. Given the sensitivity of cooling on numerical resolution, it is likely that the result for the “Virgo” run should still be interpreted as a lower limit on $`f_{}`$.

The cold fraction is reduced by including non–gravitational heating. The degree to which overcooling is suppressed depends not only on the amount of feedback energy, but also on the redshift and on the gas overdensity at which it is released into the diffuse medium. For instance, heating at $`z_h=9`$ is very efficient in decreasing $`f_{}`$ below the 10 per cent level, at the expense of delaying the bulk of star formation to $`z<\text{ }1`$. A more realistic star formation history, peaking at $`z3`$, consistently requires that most of the non–gravitational heating takes place at a similar redshift, with at least 20 per cent of the baryons still being converted into stars.

Heating at $`z_h=3`$ with $`E_h0.75`$ keV/part is shown to produce scalings of $`X`$–ray luminosity, mass and entropy vs. temperature which agree in general with observational data. This result holds independent of whether an equal amount of energy is assigned to all gas particles or whether an entropy floor is created. A similar agreement is also found for the SAM–predicted SN feedback, once realistic models for the IMF are used. Both, heating at $`z_h=9`$, or using an IMF which produces a large number of SN, are not efficient in suppressing the $`X`$–ray luminosity, which is a consequence of the fairly large amount of gas that, while avoiding cooling, is concentrated in central cluster regions.

Including cooling and star formation increases the ICM temperature in the central regions. While this helps in reconciling simulations with the observed $`M`$$`T`$ relation, it steepens the temperature profiles, which show no evidence for any decline at small cluster-centric distances. This result, which holds independent of the scheme for non–gravitational heating, is discrepant with recent observations.

Over the last year or so, different groups have presented simulations aimed at studying the effect of radiative cooling and feedback on the $`X`$–ray properties of the ICM. Most of these studies are based on simulations which follow the gas hydrodynamics within the full volume of a cosmological box (e.g., Muanwong et al. 2002; Davé et al. 2002; Kay et al. 2002). One common result of these simulations, which agrees with what we find in our analysis, is that the effect of cooling is able to alleviate or even solve the discrepancy between simulated and observed $`X`$–ray scaling properties of clusters and groups, but the fraction of baryons converted into stars is too large. To remedy this problem, Muanwong et al. (2002) pre–heated the gas by adding 1.5 keV thermal energy to all the gas particles at $`z_h=4`$. As a result, they found that the cold fraction in groups and clusters is decreased from 15 per cent to 0.4 per cent, which is somewhat smaller than the values found in our simulations. Kay et al. (2002) implemented a feedback mechanism in their simulations, which accounts for the rate of both type Ia and type II SN. By assuming an energetics twice as large as that provided by standard supernova computations, they were able to reproduce the observed $`X`$–ray scaling properties, while obtaining only 3 per cent of the gas to be converted into stars. The main limitation of this type of simulations is that one is restricted to relatively poor numerical resolution in order to limit the computational cost. For instance, the simulations by Muanwong et al. (2002) have a mass resolution which is about one order of magnitude worse than that of our “Virgo” runs and almost two orders of magnitude worse than that of our group runs. A better mass resolution within a smaller box was used by Kay et al. (2002), for which the mass of gas particles are a factor 2.8 and 22 smaller than for our Virgo and Group runs, respectively.

The results that we presented in Section 3 demonstrate that the cooling efficiency is quite sensitive to mass resolution. For this reason, one has to be careful in drawing conclusions about overcooling and how it is suppressed by extra heating, in the presence of limited numerical resolution. In fact, our simulations demonstrate that the two main problems caused by the introduction of radiative cooling, namely the overproduction of stars and the steeply increasing temperature profiles in central cluster regions, may not be easily solved by the introduction of non–gravitational heating.

Does this imply that none of our heating schemes is a realistic approximation to what happens in real clusters? The energy release in all these schemes misses, although to different degrees, to faithfully follow the simulated rate of star production. A realistic scheme for SN feedback should dump thermal energy with a rate that accurately follows the star formation rate, properly accounting for the typical life–times of different stellar populations. Furthermore, our schemes for energy release demonstrate that for feedback to have a sizeable effect on the ICM thermodynamics, it has to act in a non–local way, so as to assign most of the energy on gas particles which have a sufficiently long cooling time. Such non–local feedback mechanisms may arise from AGN activity, cosmic rays or galactic winds, for example.

While further work is clearly needed to study such feedback mechanisms self-consistently in simulations, a better understanding is also required as to whether optical/$`X`$–ray data really implies a stellar fraction as small as $`<\text{ }10`$ per cent within clusters and groups. Balogh et al. (2001) used the 2MASS results on the K–band luminosity function by Cole et al. (2001) to estimate the cosmic fraction of baryons converted into stars. After assuming a Kennicutt IMF (Kennicutt 1983), they find $`f_{}0.05`$ for our choice of $`\mathrm{\Omega }_m`$ and $`h`$, and argued that no much evidence exists for $`f_{}`$ to increase inside clusters, or to depend on the cluster mass (cf. also Bryan 2000). However, this estimate of the cosmic value of $`f_{}`$ increases by about a factor 2 if a Salpeter IMF (Salpeter 1955) were used instead. Furthermore, it is worth reminding that the estimate inside clusters relies to some degree of extrapolation. For instance, Balogh et a. (2001) obtained the stellar mass in clusters from the B–band luminosity data by Roussel et al. (2000), using $`M/L_B=4.5`$, and correcting for undetected galaxies by extrapolating the luminosity function to the faint end slope. It is clear that a more robust determination of $`f_{}`$ in clusters should rather rely on K– or H–band luminosity, which is more directly related to stellar mass (e.g., Gavazzi et al. 1996), rather than to B–band luminosity whose conversion to stellar mass is quite sensitive to galaxy morphology. More recently, Huang et al. (2002) used the Hawaii-AAO K-band redshift survey to estimate the K-band luminosity function in the local Universe. They found that the K–band luminosity density is twice as large as that from 2MASS, thus implying a twice as large $`f_{}`$ value. In light of this discussion, a $`f_{}`$ value somewhat larger than 10 per cent, possibly as large as 20 per cent, may still be viable at present, which would tend to alleviate the problem of ICM overcooling.

As for the temperature profile, our results indicate that the discrepancy between observations and simulations is unlikely to be solved by the inclusion of feedback mechanisms that are similar to the ones explored here. If this is the case, it would demonstrate that our simulations are missing some basic physical mechanisms. For instance, as we discussed, thermal conduction has been proposed to be an important effect in clusters. Another piece of physics which is currently missing from most simulation work is the effect of magnetic fields (e.g., Dolag, Bartelmann & Lesch 2002). Their introduction might give rise to non–trivial structures in the gas distribution if they can locally suppress thermal conduction, or it they provide a non–thermal contribution to the gas pressure.

There is little doubt that including such more complex physics will represent a significant, non-trivial challenge for cluster simulations of the next generation. Most of the processes involved require both, a rather sophisticated numerical method, and a treatment of sub-grid physics. Still, the inclusion of more physics in numerical codes is mandatory if the reliability and the predictive power of cluster simulations want to keep pace with the increasing quality of observational data.

## Acknowledgements

Simulations were run at the CINECA Supercomputing Center, with CPU time provided by a grant of the National Institute for Astrophysics (INAF), at the Computing Center of the Astronomical Observatory of Catania and at the Computing Center of the University of Trieste. We wish to thank Hans Böhringer, Alexis Finoguenov, Fabio Governato, Paolo Tozzi and Xian–Ping Wu for enlightening discussions.


format: text

text:
# A formula for the number of tilings of an octagon by rhombi

## Introduction

The enumeration of tilings of a centro-symmetric polygon by rhombi is a notoriously difficult problem that concerns discrete mathematics and theoretical computer science, as well as theoretical physics, in relation with quasicrystallography. In the latter community, these tilings are usually called “random tilings with octagonal symmetry”. We address the following issue: given a centro-symmetric octagon $`O_{a,b,c,d}`$, of integral sides lengths $`a,b,c`$ and $`d`$ (read clockwise; see figure 2, left), in how many ways is it possible to fill it entirely, without any gap or overlap, with the following six species of tiles: two differently oriented squares, and four differently oriented 45$`^{\text{}}`$ rhombi, the six of them with unitary side lengths? So far, this question has been solved in very particular instances only. We denote by $`𝒯_{a,b,c,d}`$ the set of all the tilings of $`O_{a,b,c,d}`$, and by $`T_{a,b,c,d}=|𝒯_{a,b,c,d}|`$ the cardinality of $`𝒯_{a,b,c,d}`$. For example, figure 1 displays the eight tilings of the set $`𝒯_{1,1,1,1}`$.

Small systems have been studied in references  up to sizes of some hundred tiles (see table 1). However, the technique employed cannot reasonably provide tiling enumerations for bigger octagons.

On the other hand, Elnitsky gave in ref.  two formulas when two sides of the octagon are set to 1:

$$T_{a,1,c,1}=\underset{r+s=a}{}\underset{t+u=c}{}\left(\genfrac{}{}{0pt}{}{r+t}{r}\right)\left(\genfrac{}{}{0pt}{}{s+t}{s}\right)\left(\genfrac{}{}{0pt}{}{r+u}{r}\right)\left(\genfrac{}{}{0pt}{}{s+u}{s}\right),$$
(1)
and

$$T_{a,b,1,1}=\frac{2(a+b+1)!(a+b+2)!}{a!b!(a+2)!(b+2)!}.$$
(2)
The first formula has been later partially simplified :

$$T_{a,1,c,1}=\frac{(a+c+1)!}{a!c!(2a+1)(2c+1)}\left[\frac{2(a+c+1)!}{a!c!}+\underset{k=0}{\overset{a}{}}\frac{1}{2k1}\left(\genfrac{}{}{0pt}{}{a}{k}\right)\left(\genfrac{}{}{0pt}{}{c}{k}\right)\right]$$
(3)
where the last sum can be written in terms of a hypergeometric function

$$\underset{k=0}{\overset{a}{}}\frac{1}{2k1}\left(\genfrac{}{}{0pt}{}{a}{k}\right)\left(\genfrac{}{}{0pt}{}{c}{k}\right)=_3F_2[1/2,a,c;1/2,1;1].$$
(4)
We propose a generalization of the first formula (1) to any side lengths, where $`T_{a,b,c,d}`$ is written as a sum of products of determinants. Even if the complexity of our formula increases with the system size, it is the first explicit algebraic expression to count tilings of an octagon (see eq. (12)), which can be in principle calculated for any system size.

As it is discussed below into detail, Elnitsky’s proof uses a “square grid representation” of tilings (see figure 2), which is closely related to the “de Bruijn dualization”, a wide-spread technique in quasicrystal science. This dualization has proved powerful to handle rhombus tilings in several circumstances. The present paper also uses this technique, thus generalizing Elnitsky’s proof.

We now state our main result. Given the side lengths $`a`$, $`b`$, $`c`$ and $`d`$, we denote by $`X`$ (resp. $`Y`$) the set of families of integers $`(x_{k,l})`$ (resp. $`(y_{k,l})`$), $`k=1,\mathrm{},b`$, $`l=1,\mathrm{},d`$, satisfying the relations:

$`0x_{k,l}a,`$ (5)
$`x_{k,l}x_{k^{},l^{}}`$ if $`kk^{}\text{ and }ll^{},`$ (6)
$`0y_{k,l}c,`$ (7)
$`y_{k,l}y_{k^{},l^{}}`$ if $`kk^{}\text{ and }ll^{}.`$ (8)
Note that conditions (6) and (8) are not exactly similar. In the following, these integers will be the coordinates of the white disks in the grid representation (see figure 2, right, and section 1). In addition, we set by convention for $`k=1,\mathrm{},b`$ and $`l=1,\mathrm{},d`$

$$\begin{array}{cc}x_{k,0}=0,\hfill & y_{k,0}=0,\hfill \\ x_{0,l}=0,\hfill & y_{0,l}=c,\hfill \\ x_{b+1,l}=a,\hfill & y_{b+1,l}=0,\hfill \\ x_{k,d+1}=a,\hfill & y_{k,d+1}=c.\hfill \end{array}$$
(9)
The reasons for this convention will be explicited below. For any two such sequences $`x=(x_{k,l})`$ and $`y=(y_{k,l})`$, we define the matrices $`M^{(u)}(x,y)`$ and $`P^{(v)}(x,y)`$ as follows: $`M^{(u)}(x,y)`$ is a $`b\times b`$ matrix of coefficients

$$m_{ij}=\left(\genfrac{}{}{0pt}{}{x_{j,u}x_{i,u1}+y_{j,u}y_{i,u1}}{x_{j,u}x_{i,u1}+ji}\right)$$
(10)
for $`1i,jb`$ and $`P^{(v)}(x,y)`$ is a $`d\times d`$ matrix of coefficients

$$p_{ij}=\left(\genfrac{}{}{0pt}{}{x_{v,j}x_{v1,i}+y_{v1,i}y_{v,j}}{x_{v,j}x_{v1,i}+ji}\right).$$
(11)
for $`1i,jd`$. Note that, by convention, we set $`\left(\genfrac{}{}{0pt}{}{A}{B}\right)=0`$ whenever $`A<0`$ or $`B<0`$ or $`B>A`$.

###### Theorem 1

With the previous notations, the number of rhombus tilings of a centro-symmetric octagon of sides $`a,b,c,d`$ reads:

$$T_{a,b,c,d}=\underset{(x,y)X\times Y}{}\underset{u=1}{\overset{d+1}{}}detM^{(u)}(x,y)\underset{v=1}{\overset{b+1}{}}detP^{(v)}(x,y).$$
(12)
It is demonstrated below that the determinants come from the enumeration, by the Gessel-Viennot method (presented below), of tilings of independent sub-domains of the octagon delimited by some points of coordinates $`(x_{k,l},y_{k,l})`$ in the square grid representation.

When $`b=d=1`$, the previous expression is reduced to Elsnitsky’s relation (1). Note that by contrast, relation (2) is not a spacial case of this formula. Beyond this simple case, the number of terms in the formula grows with the octagon size. For example, for $`(a,b,c,d)=(2,2,2,1)`$, the formula contains $`6\times 6=36`$ terms to count the 480 tilings. For $`(a,b,c,d)=(2,2,2,2)`$, there are $`20\times 20=400`$ terms and 5383 tilings. More generally, the number of terms grows exponentially with the number of tiles, but it nevertheless grows exponentially more slowly than the number of tilings. As a consequence, this formula is exponentially more compact than the crude enumeration of tilings. This point is discussed in the conclusion.

## 1 Octagonal tilings and the square grid representation

In this section, we show that octagonal tilings are conveniently represented by families of directed paths running on a rectangular patch of square grid. This representation was used by Elnitsky  and it is reminiscent of the prior “de Bruijn dualization”  and derived representations . We first expose briefly the de Bruijn dualization process. Figures 2 and 3 will help the reader. To begin with, we notice that tile edges can have four possible orientations. We define a family of de Bruijn lines for each orientation: de Bruijn lines are made up of adjacent rhombi sharing an edge with this orientation. It is always possible to extend these lines through the whole tiling up to a boundary tile. Two examples of lines are presented in both figures 2 and 3, belonging to two different families. A rhombic tile is situated at the intersection of two lines of different families (see the figures). By construction, lines of a same family never intersect; there are respectively $`a`$, $`b`$, $`c`$ and $`d`$ lines in each family.

Now we show how to translate the de Bruijn’s representation of a tiling into its square grid representation. In figure 3, we show this correspondence in the simplest case $`b=d=1`$ . The idea is to shrink the de Bruijn lines of two families among four, so that they become paths on a square grid, as displayed in figures 2 or 3. Because all tiles of a de Bruijn line have an edge of a given orientation, these paths are directed. The $`b`$ paths of the first family (denoted by $`SW`$) go from the south-west corner to the north-east one (dark gray); they can follow eastward and northward edges only; the $`d`$ paths of the second family (denoted by $`NW`$) go from the north-west corner to the south-east one (light gray); they can follow eastward and southward edges only.

In the simplest case $`b=d=1`$, to avoid ambiguity due to path tangency and to make this correspondence bijective , we keep track of the intersection of the de Bruijn lines thanks to a distinguished vertex, represented by a white disk in the right figure. It marks the position of the unique tilted square (medium gray). When $`b>1`$ or $`d>1`$ as in figure 2, there are $`bd`$ intersections and therefore $`bd`$ tilted squares. Each of them must be located by a distinguished vertex on the square grid. Paths do not cross in a same family even though they can be locally adjacent (see figure 2). We denote the paths of $`SW`$ (resp. $`NW`$) by $`SW_1,\mathrm{},SW_b`$ (resp. $`NW_1,\mathrm{},NW_d`$) from left to right. As a consequence, distinguished vertices are indexed by two integers $`k`$ and $`l`$, and are denoted by $`DV_{k,l}`$ (see figure 4).

## 2 The Gessel-Viennot method

The Gessel-Viennot method  is a combinatorial technique for the counting of configurations of directed non-intersecting paths on oriented graphs. This technique has already proved very useful for the enumeration of rhombus tilings (see  for examples, as well as section 3 of the present paper). It has been extensively described in the literature  and we shall only briefly explain it in the present paper, focusing on the underlying ideas and not on technical details. The method is rather general and can be applied to any acyclic oriented graph $`𝒢`$, in which are selected two families of vertices, $`d_i`$ (“departure” vertices) and $`a_j`$ (“arrival” vertices), $`i,j=1,\mathrm{},n`$. We consider directed paths, running on $`𝒢`$, starting from one vertex $`d_i`$ and arriving at one vertex $`a_j`$. By “directed”, we naturally mean that the paths must follow the edge orientations. In addition, this graph is supposed to satisfy the property of compatibility: if two directed paths on $`𝒢`$ are going respectively from $`d_{i_1}`$ to $`a_{j_1}`$ and from $`d_{i_2}`$ to $`a_{j_2}`$, if these paths do not cross, and if $`i_1<i_2`$ then $`j_1<j_2`$. This property is very specific to two-dimensional graphs.

We are interested in the number $`D_n`$ of configurations of $`n`$ non-intersecting directed paths on $`𝒢`$, where the $`i`$-th path goes from $`d_i`$ to $`a_i`$: two paths are said to be non-intersecting if they share no vertex; $`n`$ paths are said to be non-intersecting if any two paths are non-intersecting. If we denote by $`\lambda _{ij}`$ the number of paths going from $`d_i`$ to $`a_j`$, then the Gessel-Viennot theorem states that

$$D_n=det(\lambda _{ij})_{1i,jn}.$$
(13)
The idea of the proof is that in this determinant, all configurations of $`n`$ paths, whether intersecting or not, the $`i`$-th path going from $`d_i`$ to $`a_{\sigma (i)}`$, for any permutation $`\sigma `$, are counted, with a $`+`$ or $``$ sign. Because of these signs, all configurations with one or more intersections cancel two by two. Only the non-intersecting configurations remain. They are exactly the configurations under interest thanks to the property of compatibility. The interested reader is referred to Stembridge  for more detailed explanations.

## 3 Proof of theorem 1

We are now ready to prove theorem 1. First of all, we need to endow the square grid with integer coordinates in order to locate the positions of distinguished vertices. They are defined according to the usual conventions, so that the south-west and north-east corners have respective coordinates $`(0,0)`$ and $`(a,c)`$. The coordinates of $`DV_{k,l}`$ are denoted by $`(x_{k,l},y_{k,l})`$. By extension, we also define the vertices $`DV_{k,0}=(0,0)`$, $`DV_{k,d+1}=(a,c)`$, $`DV_{0,l}=(0,c)`$ and $`DV_{b+1,l}=(a,0)`$. They are the ends of paths of families $`SW`$ and $`NW`$. If their coordinates are also denoted by $`(x_{k,l},y_{k,l})`$, these last definitions are compatible with the conventions (9) given in introduction. All these coordinates naturally obey relations (5) and (7).

Furthermore, because of the directed character of de Bruijn lines, distinguished vertices $`DV_{k,l}`$ are constrained by some conditions when they belong to the same paths, and they must obey relations (6) and (8) as well. These four conditions define the sets $`X`$ and $`Y`$, as it was stated in the introductory part.

Let $`x=(x_{k,l})X`$ and $`y=(y_{k,l})Y`$ be an admissible set of coordinates of the distinguished vertices. We denote by $`𝒯_{x,y}`$ the subset of tilings of $`𝒯_{a,b,c,d}`$ in the square grid representations of which the distinguished vertices have these coordinates. The subsets $`𝒯_{x,y}`$ are two-by-two disjoint so that $`T_{a,b,c,d}=_{(x,y)X\times Y}|𝒯_{x,y}|`$. Our purpose is now to calculate each $`|𝒯_{x,y}|`$. This calculation is feasible because for a given $`(x,y)`$, the subset $`𝒯_{x,y}`$ can be factorized into simple sets (see eq. (15)). Each of them can in turn be counted by the Gessel-Viennot method, which leads to relation (12).

First we need to introduce two definitions. Given a configuration of vertices $`DV_{k,l}=(x_{k,l},y_{k,l})`$, we fix $`l`$ and we consider in isolation the vertices $`DV_{k,l1}`$ as well as $`DV_{k,l}`$, $`k=1,\mathrm{},b`$ (see figures 4 and 5). Then we define the set $`sw(l)`$ of all the configurations of $`b`$ directed non-crossing paths, the $`k`$-th path going from $`DV_{k,l1}`$ to $`DV_{k,l}`$, with $`k=1,\mathrm{},b`$. These paths are directed from south-west to north-east. They have no constraint except that they are directed and non-crossing (these paths can have tangencies). In a similar way, we define the sets $`nw(k)`$ for any $`k`$: they are the sets of all configurations of $`d`$ directed non-crossing paths, going from north-west to south-east. The $`l`$-th path goes from $`DV_{k1,l}`$ to $`DV_{k,l}`$.

Now in order to prove theorem 1, we start from the following observation, illustrated in figure 4: in $`𝒯_{x,y}`$ the distinguished vertices $`DV_{k,l}`$ are held fixed and one can consider independent patches of the families $`SW`$ or $`NW`$, as follows. Without loss of generality, we focus on $`SW`$. We cut each path $`SW_k`$ into $`d+1`$ sections, denoted by $`SW_k(l)`$, where $`l=1,\mathrm{},d+1`$. The section $`SW_k(l)`$ goes from $`DV_{k,l1}`$ to $`DV_{k,l}`$. Then all the $`b`$ sections $`SW_k(l)`$, $`k=1,\mathrm{},b`$ form a local path configuration denoted by $`P_{SW}(l)`$. It belongs to $`sw(l)`$. In a similar way, the corresponding grid patches defined with respect to the family $`NW`$ are denoted by $`P_{NW}(k)`$ and belong to $`nw(k)`$. Therefore, when $`x`$ and $`y`$ are fixed, we have the natural inclusion

$$𝒯_{x,y}\underset{u=1}{\overset{d+1}{}}sw(u)\underset{v=1}{\overset{b+1}{}}nw(v),$$
(14)
where the products are direct. We prove below that

###### Lemma 1

The previous inclusion is an equality:

$$𝒯_{x,y}=\underset{u=1}{\overset{d+1}{}}sw(u)\underset{v=1}{\overset{b+1}{}}nw(v).$$
(15)
It follows from eq. (15) that

$$𝒯_{a,c,b,d}=\underset{(x,y)X\times Y}{\overset{\text{disjoint}}{}}\underset{u=1}{\overset{d+1}{}}sw(u)\underset{v=1}{\overset{b+1}{}}nw(v),$$
(16)
and that

$$T_{a,c,b,d}=\underset{(x,y)X\times Y}{}\underset{u=1}{\overset{d+1}{}}|sw(u)|\underset{v=1}{\overset{b+1}{}}|nw(v)|.$$
(17)
The remainder of the proof consists in calculating the cardinalities $`|sw(u)|`$ and $`|nw(v)|`$ by the Gessel-Viennot method. Indeed, it is also demonstrated below that

###### Lemma 2

When $`x`$ and $`y`$ are fixed,

$$|sw(u)|=detM^{(u)}(x,y)\text{;}|nw(v)|=detP^{(v)}(x,y).$$
(18)
### Proof of lemma 1:

We need to prove the reverse inclusion

$$𝒯_{x,y}\underset{u=1}{\overset{d+1}{}}sw(u)\underset{v=1}{\overset{b+1}{}}nw(v).$$
(19)
Configurations from the sets $`sw(u)`$ provide sections of paths from $`DV_{k,l1}`$ to $`DV_{k,l}`$. When concatenated, these sections provide complete directed non-crossing paths from $`(0,0)`$ to $`(a,c)`$, which form a family $`SW`$. In a similar way, sections from the $`nw(v)`$ provide directed non-crossing paths from $`(0,c)`$ to $`(a,0)`$, forming a family $`NW`$. We only need to check that any two paths from $`SW`$ and $`NW`$ only cross at the distinguished vertices $`DV_{k,l}`$. This point is ensured by the directed character of path sections (see figure 5). This last observation is crucial and all the demonstration relies on it: it ensures the reverse inclusion and therefore the direct character of the product (15), from which our enumerating formula ensues.

### Proof of lemma 2:

So far we have used the terminologies “non-intersecting paths” in section 2 and “non-crossing paths” in section 3. Now it is time to precise what aspects these two terms cover. We have seen that non-crossing paths can have tangencies, that is to say they can share vertices or edges of the grid, but they cannot step over one another. In particular, non-crossing paths of families $`SW`$ (or $`NW`$) share their ends, but can be indexed from west to east without ambiguity.

On the contrary, non-intersecting paths cannot share any vertex or edge. Therefore, if we want to use the Gessel-Viennot method, we need to transform configurations of non-crossing paths on the square grid into configurations of non-intersecting paths. The trick consists in shifting non-crossing paths, as it is illustrated in figure 6. The trick is standard and was already used by Elnitsky  for example.

We use a unitary shift vector $`𝒖=(1,1)`$ and we shift the $`k`$-th path section $`SW_k(l)`$ by a vector $`(k1)𝒖`$ (see the figure). The new paths still belong to the square grid. The $`k`$-th new path section goes from the new vertex $`DV_{k,l1}^{}`$ to the new vertex $`DV_{k,l}^{}`$, of coordinates

$$x_{k,l}^{}=x_{k,l}+(k1)\text{and}y_{k,l}^{}=y_{k,l}(k1).$$
(20)
The so-obtained path configuration is non-intersecting by construction. This correspondence between configurations of non-crossing paths and configurations of non-intersecting paths is bijective.

Now we use the Gessel-Viennot technique by setting the departure vertices $`d_i=DV_{i,l1}^{}`$ and the arrival ones $`a_j=DV_{j,l}^{}`$. The number of directed paths on the square grid going from a vertex $`D=(x_D,y_D)`$ to a vertex $`A=(x_A,y_A)`$ is simply given by the binomial coefficient

$$\lambda =\left(\genfrac{}{}{0pt}{}{x_Ax_D+y_Ay_D}{x_Ax_D}\right).$$
(21)
Then one computes the coefficients $`\lambda _{ij}`$ involved in the Gessel-Viennot method: $`\lambda _{ij}=m_{ij}`$. Thus one obtains the matrices $`M^{(u)}(x,y)`$, the determinants of which count the elements of $`sw(u)`$. In the same way, to count the elements of the sets $`nw(v)`$, one must shift the sections of paths $`NW_l`$ by $`(l1)𝒗`$ where $`𝒗=(1,1)`$. One gets the matrices $`P^{(v)}(x,y)`$ and takes their determinant, which completes the proof of lemma 2.

## 4 Conclusion

We have demonstrated how Elnitsky’s technique can be generalized to octagons of any size, leading to an explicit enumerative formula (theorem 1).

We also notice that the conditions (5) and (6) (resp. (7) and (8)) that define the set $`X`$ (resp. $`Y`$) are identical to the conditions defining plane partitions of height $`a`$ (resp. $`c`$) on a $`b\times d`$ grid . This point is remarkable because such plane partitions are known to be equivalent to rhombus tilings filling a centro-symmetric hexagon of sides lengths $`b`$, $`d`$ and $`c`$ (resp. $`b`$, $`d`$ and $`a`$. We have derived a partial combinatorial interpretation of our formula (12) in terms of these tilings of hexagons. It is related to a natural decomposition of the configuration sets of tilings of octagons, as described in ref. . But it goes beyond the scope of the present paper and will be described elsewhere .

If $`T_{a,b,c}^{\text{hex}}`$ denotes the number of tilings of the centro-symmetric hexagon of sides $`a`$, $`b`$ and $`c`$, the previous remark leads to the lower bound

$$T_{a,b,c,d}T_{b,d,c}^{\text{hex}}T_{b,d,a}^{\text{hex}},$$
(22)
the number of terms the formula. By construction, the sets $`nw(k)`$ and $`sw(l)`$ are not empty and all terms are positive. In statistical physics and more specifically in quasicrystal science, people are interested in thermodynamic quantities such as the configurational entropy (per tile): $`S=\mathrm{ln}(T_{a,b,c,d})/N_T`$ where $`N_T`$ is the number of tiles. With our polygonal boundary conditions, this quantity has a finite limit when $`N_T`$ goes to infinity provided the relative ratios of the side lengths also have a finite limit . In the so-called “diagonal” case where all side lengths are equal, taking into account the number of tiles, the previous relation becomes $`SS^{\text{hex}}=(3/2)\mathrm{ln}32\mathrm{ln}20.262`$ . The actual value of $`S`$ is numerically known to be close to 0.36 . The previous lower bound is manifestly loose and its improvement requires a better knowledge of the asymptotic behavior of the determinants in (12) at the large size limit.

But the main advantage of our formula precisely lies on the fact that the previous bound is weak: the formula realizes an exponential reduction of the number of terms as compared to a crude enumeration of tilings. Indeed, as it was just discussed in the previous paragraph, the number of terms grows exponentially like $`\mathrm{exp}(0.26N_T)`$ wheres the number of tilings grows like $`\mathrm{exp}(0.36N_T)`$. Even if in practice we cannot compute numerically the number of tilings of octagons bigger than in table 1, the progress is already significant. Moreover there exists some hope to simplify our formula, at least partially, as in eq. (3).

In addition, the formula brings a new insight into the structure of tiling sets: it emphasizes a natural decomposition of the sets into smaller disjoint subsets, the cardinality of which is simply given by evaluation of determinants.

## Acknowledgments

We are indebted to one of the referees for a careful reading of our manuscript and sound suggestions of improvements.