text:
# The Impact of Galaxy Formation on the X-ray Evolution of Clusters

## 1 Introduction

One of the fundamental puzzles of the X-ray universe concerns the relation between the X-ray luminosity and gas temperature of clusters of galaxies. A simple scaling analysis (Kaiser 1986) suggests that the temperature and luminosity should be related by $`TL^{1/2}`$. Temperatures have now been measured for the diffuse X-ray emission for an extensive range of groups and clusters (David et al. 1993; Ponman et al. 1996; Allen & Fabian 1998; Markevitch 1998; Mulchaey & Zabludoff 1998; Arnaud & Evrard 1999; Helsdon & Ponman 2000). In contrast to the theoretical prediction, the observations show a much shallower trend, approximately $`TL^{1/3}`$.

A closely related problem is the evolution of the cluster X-ray luminosity function. Kaiser’s (1986) analysis of the evolution of the X-ray properties of clusters suggested that dense, X-ray luminous associations of galaxies should be more numerous in the intermediate and high redshift universe. This possibility was soon ruled out by the initial results of the EMSS cluster survey (Gioia et al. 1990; Henry et al. 1992), which quickly established that clusters in the distant universe have a comparable space density to those of the local universe. This has been confirmed in more recent ROSAT surveys (eg., Jones et al. 2000).

Initially, one might have hoped that including radiative cooling of the gas (omitted from Kaiser’s analysis) might resolve this discrepancy. Unfortunately, it is extremely difficult to obtain numerically convergent results from simulations once gas cooling is included. The difficulty is inherent to the problem. Because the universe is dense at early times, cooling is initially very efficient. This leads to an unrealistically large fraction of the halo baryon content cooling at high redshifts to form very small galaxies. As White & Rees (1978), White & Frenk (1991), Cole (1991), Suginohara & Ostriker (1998) and Pearce et al. (2000) amongst others have shown, some form of heating is required to overcome this catastrophe. In Appendix A, we examine how the cooled gas fraction depends on cluster mass, and find that the cooled fraction depends too weakly on cluster temperature to explain the discrepancy. Therefore radiative cooling cannot by itself solve the problem with the temperature-luminosity relation (see Bryan et al. 2000 and Balogh et al. 2001 for an extended discussion).

One approach to this problem that has given encouraging results is to assume that the gas is “preheated” before collapsing into the cluster (Evrard & Henry 1991; Kaiser 1991; Navarro, Frenk & White 1995). This creates an entropy floor in the gas, ensuring that it remains at low densities in low mass systems, and results in a much improved match to the T-L relation (Balogh, Babul & Patton 1999; Valageas & Silk 1999; Tozzi & Norman 2000). This model also provides an encouraging match to the surface brightness profiles of low mass groups (Ponman et al. 1999). The problem with this model is in explaining the origin of this diffuse heating and its apparent uniformity.

Prolonged heat input from galaxy formation has been suggested as a solution by Wu et al. (1998, 1999a) and Cavaliere et al. (2000). They calculated the response of the gas profile to the energy input from supernovae by using a simple energetic approach. They start by assuming a 1-parameter form for the gas profile or equation of state and then calculate how the total energy of the halo gas depends on this parameter. By then solving for the dependence of the parameter on the energy excess relative to an initial profile corresponding to the case of no heat input, the contributions to the energy balance from gravity, radiative cooling and supernova heating can be taken into account in calculating the new gas distribution. This approach successfully accounts for the shallow present-day temperature-luminosity (T-L) relation if galaxy formation has a roughly uniform efficiency in haloes of different masses. Since the binding energy per particle increases with halo mass, while the additional heating remains roughly constant, high mass clusters are almost unaffected while the gas in low mass groups becomes unbound.

In this paper, we develop a new model for the evolution of the masses and density profiles of the hot gas in galaxy, group and cluster haloes, using the semi-analytic galaxy formation scheme of Cole et al. (2000) to predict the evolution of the supernova heating rate in dark haloes of different masses, and calculating the response of the gas profile to this heating using a method related to those of Wu et al. and Cavaliere et al. The semi-analytic scheme is an elaboration of that described by Baugh et al. (1998), and is based on similar principles to the models described by Kauffmann, White & Guiderdoni (1993) and reviewed by Somerville & Primack (1999). We apply our model to study the evolution of the X-ray luminosity function and the temperature-luminosity relation.

The structure of this paper is as follows. Our method for relating the gas distribution in the halo to the supernova energy input is presented in §2. The predicted X-ray properties are detailed in §3. In section §3.1, we show that supernova heating is able to produce the observed slope and normalization of the present-day T-L relation only if the efficiency with which the supernova explosion energy is transferred to the diffuse intracluster medium (ICM) is very high, or if the stellar initial mass function (IMF) is tilted to produce an overabundance of high-mass stars relative to the IMF in the solar neighbourhood. Alternatively, heat input from AGN activity may be required. In §3.2, we apply this model to the X-ray luminosity function of galaxy clusters. We compare the evolution predicted by the model within a flat, $`\mathrm{\Omega }_0=0.3`$, CDM cosmology with the available observations of intermediate redshift clusters. In §3.3, we consider the X-ray properties of the universe at very high redshifts, and in §3.4, we compare the expectations based on our galaxy formation model with those from two extreme models for the redshift evolution of the heat input. A further discussion of the problems and a restatement of our conclusions are given in §4 and §5.

## 2 The Model

Wu et al. (1998, 1999a) have suggested a simple approach that allows non-gravitational heating to be incorporated into the calculation of the properties of cluster gas. Starting from a default distribution, gas is redistributed to larger and larger radii until the total energy increase matches the energy input from galaxy formation. The effect of heat input may affect the distribution of gas within clusters in a variety of ways. Our approach differs from that of Wu et al. both in the way we determine the default gas distribution and in the way we modify the gas distribution in response to the excess energy input. Firstly, while Wu et al. adopt a complex prescription for the default gas distribution based on the cluster’s gravitational binding energy, our default profile is based explicitly on the observed properties of high temperature rich clusters. We are able to do this because the ranges of excess energy that we consider have little impact on the gas distribution in these systems. Secondly, Wu et al. explore a variety of models in which heating occurs either by uniformly varying the gas temperature, or by varying the polytropic index of the gas. In contrast, our approach is empirical and motivated by the observations of Arnaud & Evrard (1999) and Lloyd-Davies, Ponman & Cannon (2000) who find that the gas profiles of clusters become systematically shallower at lower temperatures. We therefore assume that the overriding effect of heating is to reduce the slope of the radial density profile of the gas. Our empirical approach does not require us to choose explicitly between the isothermal and polytropic regimes. Instead, for our given density profile, we solve for hydrostatic equilibrium in order to determine the gas temperature.

We assume that the dark matter density of the halo follows a Navarro, Frenk & White (1997, NFW) profile:

$$\rho _{\mathrm{DM}}(r)=\frac{\rho _sa^3}{r(r+a)^2}$$
(1)
The density normalization $`\rho _s`$ and dependence of scale radius $`a`$ on halo mass are calculated as described in Cole et al. (2000). We parameterize the gas distribution using a conventional $`\beta `$ model (Cavaliere & Fusco-Femiano 1976):

$$\rho _{\mathrm{gas}}(r)=\frac{\rho _c}{(1+(r/r_c)^2)^{3\beta /2}}$$
(2)
Both of these distributions are assumed to apply within the virial radius $`R_{\mathrm{vir}}`$, where $`R_{\mathrm{vir}}`$ is calculated from a spherical collapse model as described by Cole et al. The first step is to fix the parameters of the default radial gas profile that applies in the absence of any energy input from supernovae. We initially distribute the gas with a core radius $`r_c`$ that is a fixed fraction (7%) of the virial radius $`R_{\mathrm{vir}}`$, and set $`\beta =0.7`$ in order to match observations of the most massive clusters (eg. Lloyd-Davies, Ponman & Cannon 2000). The temperature of the gas at the virial radius is set to $`0.6T_{\mathrm{vir}}`$, as suggested by the numerical simulations of Eke, Navarro & Frenk (1998) and Frenk et al. (1999). Here the virial temperature is defined as

$$T_{\mathrm{vir}}=\frac{1}{2}\frac{\mu m_H}{k}\frac{GM}{R_{\mathrm{vir}}}$$
(3)
where $`M`$ is the total mass within $`R_{\mathrm{vir}}`$, $`\mu `$ is the mean molecular mass (we take $`\mu =0.59`$ for fully ionized gas) and $`m_H`$ is the mass of a hydrogen atom. The temperature of the gas at smaller radii is then found by solving for hydrostatic equilibrium in the gravitational potential of the dark matter. This technique accurately reproduces the luminosity weighted temperature of the cluster simulated by Frenk et al. (1999). We adjust the normalization $`\rho _c`$ of the default gas profile so that the baryonic mass fraction (ie. gas plus galaxies) enclosed within the virial radius is equal to the cosmic baryon fraction. Treating the total mass in this way takes into account the effect of cooled gas and stars that are locked into galaxies, thus reducing the hot gas fraction of the cluster. X-ray luminosities are calculated from the gas within the cluster virial radius, since material at larger radii is unlikely to be in hydrostatic equilibrium. In practice, this cut-off has little influence on the X-ray luminosity, since that is dominated by the densest material in the cluster core.

Having established the default profile for a given cluster, we reduce the slope $`\beta `$ of the gas profile while keeping the core radius $`r_c`$ constant, until the total energy (thermal plus gravitational) of the gas is increased by the required amount. As the profile changes, we keep the pressure and density (and thus temperature) at the cluster virial radius fixed at the value found for the default profile. This results in the mass of gas within the virial radius being less than for the default profile. The excess gas is assumed to be expelled. It is displaced to the virial radius and included in the energy balance calculation, but not in the calculation of the X-ray luminosity. The temperature of this material is assumed to be the same as that of the gas at $`R_{\mathrm{vir}}`$. This corresponds to the lowest plausible temperature for the expelled gas to be both in pressure equilibrium with its surroundings and buoyant with respect to the remaining cluster material.

We have chosen the virial radius as the point at which to normalize our density profiles because this approximately delineates the region of the cluster that is in dynamical equilibrium and separates it from the outer parts of the cluster that are dominated by bulk inflow. Outside the virial radius, the gas is unlikely to be in hydrostatic equilibrium. Close to the virial radius, the infalling gas is shocked so that its bulk motion is converted to internal energy. In the one dimensional simulations of Knight & Ponman (1997), where the infalling material has uniform initial entropy, the shock radius occurs at 1-$`1.5R_{\mathrm{vir}}`$, in line with the boundary radius we assume here. In three dimensional simulations, the shock radius is more poorly defined because the infalling material already has a range of initial entropies and this tends to smooth out the shock, but the same general picture applies. The gas pressure at the virial radius is thus regulated by the dynamical pressure of the infalling gas.

There is a limit to the overall energy increase that can be accommodated by flattening the gas profile. Rather than letting $`\beta `$ become arbitrarily low, we impose a minimum value $`\beta _{\mathrm{min}}=0.2`$. If the required slope falls below this value, we set $`\beta =\beta _{\mathrm{min}}`$, and instead allow the temperature of the gas (both inside and outside $`R_{\mathrm{vir}}`$) to rise to accommodate the excess energy, giving up the condition $`T(R_{\mathrm{vir}})=0.5T_{\mathrm{vir}}`$, but maintaining the condition that $`T=T(R_{\mathrm{vir}})`$ for the expelled gas. Since the pressure at $`R_{\mathrm{vir}}`$ is kept fixed, the gas density must then fall and a greater fraction of the gas mass is expelled. Our results are not sensitive to the exact choice of $`\beta _{\mathrm{min}}`$, since the total energy of the cluster depends only very weakly on $`\beta `$ for $`\beta <0.4`$. The lowest values in observed systems are $`\beta 0.35`$ (Lloyd-Davies, Ponman & Cannon 2000).

Figure 1 shows the relation between energy input (ie. the excess energy) per baryon, $`E_\mathrm{x}`$, and X-ray luminosity $`L_\mathrm{X}`$ for clusters with virial temperatures of 1 and 5 keV (Panel (a)), and the fraction of the original gas mass that remains within the cluster virial radius (Panel (b)). Note that the decline in X-ray luminosity is much more rapid than the decline in the gas mass within $`R_{\mathrm{vir}}`$. Experimenting with different schemes for modelling the effects of heating, such as keeping the mass within $`R_{\mathrm{vir}}`$ constant, shows that the fixed pressure assumption is the most effective at reducing the X-ray luminosity for a given energy input. For comparison, Panel (c) shows the dependence of the slope parameter $`\beta `$ on the injected energy.

The effects on the temperature-luminosity relation of radiative cooling are discussed in Appendix A. We argue there that radiative cooling will cause some flattening of the T-L relation compared to the case of no cooling and no energy input, but not enough by itself to match the observed relation. This is because, although the fraction of gas inside the cooling radius increases with decreasing halo mass, this dependence is too weak in the relevant mass range to flatten the relation to $`TL^{1/3}`$.

X-ray luminosities and luminosity-weighted temperatures for individual dark matter haloes are calculated using Peacock’s (1996) analytic fit to the Raymond-Smith cooling function, assuming a metal abundance of $`1/3`$ solar. The results change by only a few percent if we use the tabulated cooling function directly. This includes both bremsstrahlung and recombination processes and is adequate for the range of halo masses considered here. Representative halo formation and merging histories are generated using a Monte-Carlo method based on the extended Press-Schechter model as described by Cole et al. (2000). This ensures that our model includes the correct halo mass distribution and assigns collapse redshifts to individual haloes. We use the properties of the halo at its collapse time for determining the X-ray properties.

We adopt the cosmological parameters $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }_0=0.7`$, $`h=0.7`$, $`\sigma _8=0.8`$ and $`\mathrm{\Gamma }=0.19`$, where $`\mathrm{\Lambda }_0`$ is the cosmological constant measured in units of $`3H_0^2/c^2`$, $`\sigma _8`$ is the linear theory mass variance in spheres of radius $`8h^1`$Mpc at the present, and $`\mathrm{\Gamma }`$ is the shape parameter of the initial spectrum of density fluctuations defined by Efstathiou, Bond & White (1992). With these parameters, our model X-ray temperature function matches the data of Eke et al. (1996, 1998). Note our value of $`\sigma _8`$ differs slightly from that inferred by Eke et al. because our luminosity-weighted gas temperatures are $``$15% higher than the cluster virial temperatures that they use. This temperature offset is consistent with the results of hydrodynamical simulations of clusters (eg. Frenk et al. 2000), and depends on the profile adopted for the gas distribution in the central regions of the clusters (which dominate the X-ray luminosity). In order to match the observed temperature function, we have lowered $`\sigma _8`$ from 0.93 to 0.80. We retain the $`\mathrm{\Gamma }=0.19`$ power spectrum shape preferred by Eke et al.

We normalize the model to fit the observed temperatures and luminosities of the most luminous X-ray clusters by adjusting the cosmic baryon fraction. These clusters are almost unaffected by energy injection. Their predicted luminosities scale as $`L_X^{\mathrm{model}}\mathrm{\Omega }_b^2h(1+z_f)^{3/2}T^2`$ for given values of $`\mathrm{\Omega }_0`$ and $`\mathrm{\Lambda }_0`$ (see Appendix A), while the observationally inferred values scale as $`L_X^{\mathrm{obs}}h^2`$. Matching models to observations, we find that $`\mathrm{\Omega }_\mathrm{b}=0.025h^{3/2}`$ gives a good fit to the observed X-ray luminosities of clusters with virial temperatures greater than $`7\text{keV}`$. Once the model is normalized in this way, the gas fractions within $`1.5h^1`$ Mpc are consistent with more detailed observational estimates.

We calculate the amount of energy injected into the halo gas using the detailed semi-analytic model of galaxy formation of Cole et al. (2000, GALFORM). This model follows the formation of dark matter halos through hierarchical clustering and merging, using merger histories generated by a Monte Carlo scheme, and then calculates the cooling and collapse of gas within halos to form galaxies, and the formation of stars from the cool gas. The model includes the effects of feedback from supernova explosions in expelling cold gas from galaxy disks, and traces the mergers of galaxies within common dark matter halos. From the star formation history of each galaxy, we can then calculate the rate of supernovae as a function of time, once we assume an IMF. The model therefore predicts for each dark halo both the total mass of baryons in galaxies, and the total number of supernovae that have occurred in galaxies in that halo and in its progenitor halos at earlier times. We calculate the excess energy of the hot gas in a particular halo by summing the energies of all the supernovae that have occurred in the progenitors of the halo up to its formation time, after allowing for some efficiency for this energy to be transferred to the halo gas. Both the excess energy and the fraction of baryons in galaxies have a systematic variation with halo mass and a random scatter in haloes of a given mass due to different cooling and star formation histories.

We ran the GALFORM model with the same parameter values (for star formation, feedback, etc) as in Cole et al. (2000). These values were chosen to match observed properties of present-day galaxies, in particular luminosity functions, colours and stellar mass-to-light ratios. The GALFORM model does not include the effects of the modification of the halo gas profile due to energy injection when it calculates the rate of gas cooling, so our modelling is not fully self-consistent. A fully self-consistent treatment would require us to reconsider the form of the star formation law and investigate afresh what combination of parameters gives the best fit of predicted to observed galaxy properties in the present-day universe, once we include the effects of energy injection on gas cooling. This self-consistent treatment is postponed to a future paper. In the present paper, we have a more limited aim, which is to investigate the consequences for the X-ray properties of the ICM of including energy injection based on a specific ab initio model of galaxy formation which has already been shown to reproduce a wide range of observational data on galaxy properties. In a self-consistent calculation, the effect of injection will be to reduce the amount of gas which cools onto galaxies, mimicking the effect of reducing $`\mathrm{\Omega }_b`$. For this reason, we allow the value of $`\mathrm{\Omega }_b`$ used in GALFORM to be smaller than the value used in calculating the ICM properties. Specifically, GALFORM was run with $`\mathrm{\Omega }_b=0.02`$, as in Cole et al. (2000), to calculate galaxy masses and supernova rates. The total baryonic mass was calculated assuming $`\mathrm{\Omega }_\mathrm{b}=0.025h^{3/2}`$ ($`\mathrm{\Omega }_b=0.043`$ for $`h=0.7`$), and the ICM mass calculated as the difference of the total baryon mass and the mass in galaxies.

We treat the energy injected into the ICM per unit mass of stars formed as a free parameter, which we will adjust in order to fit the present-day form of the T-L relation. We adopt the parameterisation that an energy $`ϵ_{\text{sn}}10^{49}\text{erg}`$ goes into heating the ICM per M of stars formed. There are two sources of uncertainty in trying to predict the value of $`ϵ_{\text{sn}}`$ from first principles: (a) The number of Type II supernovae per $`M_{}`$ of stars formed, $`\eta _{sn}`$, depends on the IMF and on the minimum stellar mass $`m_{sn}`$ for core-collapse. For a Salpeter IMF with an upper mass limit of 125M, lower mass limit of 0.1 M and $`m_{sn}=8M_{}`$, $`\eta _{sn}=0.007`$. The Cole et al. (2000) semi-analytic model uses a somewhat different IMF, that of Kennicutt (1983), which is a better fit to that observed in the solar neighbourhood at $`m<M_{}`$, but this predicts the same number of supernovae per unit mass as in the Salpeter case, once the fraction of brown dwarfs (with $`m<0.1M_{}`$) is normalized to match the observed galaxy luminosity function as in Cole et al. A higher SN rate applies if the IMF is skewed towards high mass stars, or if the lower mass limit for the progenitors of supernovae is reduced (e.g. Chiosi et al. 1992). Lower supernova rates are suggested by recent analyses of the metal abundance of the intracluster medium (Renzini, 1997, Kravtsov & Yepes, 2000). (b) Each Type II supernova explosion releases an energy of about $`10^{51}`$ ergs (eg. Woosley & Weaver 1986), but some fraction of this energy is lost by radiative cooling as the remnant is expanding into the interstellar medium of the host galaxy, and so is not available to heat the ICM. For instance, Thornton et al. (1998) find that 70-90% of the energy is lost in this way. Thus, for a Salpeter IMF, we predict $`ϵ_{\text{sn}}=0.7`$ if none of the supernova energy is lost by radiation, but in practice a much smaller value seems likely.

We can convert the heating efficiency $`ϵ_{\text{sn}}`$ into an excess energy per baryon $`E_\mathrm{X}`$ once the fraction of baryons converted into stars $`f_{\mathrm{gal}}`$ is known. For clusters with $`1<T<10\mathrm{k}\mathrm{e}\mathrm{V}`$, our models give $`f_{\mathrm{gal}}0.13`$ for $`h=0.7`$, with only a weak dependence on $`T`$ in this range. This value scales with the Hubble constant roughly as $`f_{\mathrm{gal}}0.16h^{1/2}`$, if the parameters in the semi-analytic model are adjusted to match observed galaxy luminosities and mass-to-light ratios at each $`h`$, resulting in the stellar mass scaling as $`h^1`$, and if the total baryon fraction is scaled as $`\mathrm{\Omega }_bh^{3/2}`$ to match the X-ray luminosities. The excess energy per baryon in the ICM is then

$$E_\mathrm{X}0.50ϵ_{\text{sn}}\left(\frac{f_{\mathrm{gal}}}{0.16h^{1/2}}\right)h^{1/2}\text{keV per particle}$$
(4)
The value of $`ϵ_{\text{sn}}`$ that is required to make the model clusters fit the observed T-L relation then scales with $`h`$ approximately as $`h^{1/2}`$, since a certain ($`h`$-independent) value of $`E_\mathrm{X}`$ is required at each $`T`$ to shift the clusters onto the observed relation from the relation that applies in the absence of heating or cooling.

As explained above, we expect the heating efficiency $`ϵ_{\text{sn}}`$ due to supernovae to be significantly less than unity, if the IMF has the standard form and radiative losses are significant. However, additional energy may be available from active galactic nuclei. AGNs may emit mechanical energy in the form of winds or jets that can heat the gas in the surrounding dark halo. Although the details of the fuelling of AGN activity are unclear (see Nulsen & Fabian 2000, for a recent discussion), the requirements for this fuelling are similar to those for star formation, and the two processes may be closely linked. We will assume that the AGN activity parallels the star formation activity in the galaxies. If all galaxies harbour black holes with masses close to those suggested by Magorrian et al. (1998), we can estimate the available energy as follows. If we assume that the total energy released by the formation of each black hole of mass $`M_{\mathrm{BH}}`$ is approximately $`0.1M_{\mathrm{BH}}c^2`$ (eg., Rees 1984). Magorrian’s relation suggests $`M_{\mathrm{BH}}0.06M_{\mathrm{stars}}`$ where $`M_{\mathrm{stars}}`$ is the mass in stars (strictly, the bulge mass). Combining these relations shows that the available energy is $`10^{52}`$ ergs per $`M_{}`$ of stars, or $`ϵ_{\text{sn}}=1000`$, compared to $`ϵ_{\text{sn}}1`$ from supernovae. Thus, the energy contribution to the ICM from AGN could, in principle, exceed that from galaxies by several orders of magnitude (depending on the fraction of the AGN luminosity emitted as kinetic energy). For this reason we will allow for the possibility that $`ϵ_{\text{sn}}>1`$.

## 3 Results

### 3.1 The Temperature-Luminosity Relation

As expected, if the ICM is assumed to be heated only by gravitational collapse, with no energy injection from galaxies, then the model clusters fail to match the observed slope of the T-L relation. Data from David et al. (1993) show a slope close to $`TL^{1/3}`$, a result that has been confirmed by the analysis of more recent ASCA observations (Arnaud & Evrard 1999). Although, the brightest clusters may follow a steeper slope than this when the luminosities are corrected for contributions from cooling flows (Markevitch 1998; Allen & Fabian 1998), the $`L^{1/3}`$ slope extends down to groups of much lower luminosity (Ponman et al. 1996; Mulchaey & Zabludoff 1998; Helsdon & Ponman 2000). The X-ray properties of our model clusters are compared with observational data in Figure 2, in which the dotted line shows the predicted median relation for the case when there is no energy injection. We prefer to plot this relation with temperature on the vertical axis as the observational uncertainties are far greater for X-ray temperatures than luminosities.

In order to match the observed form of the T-L relation, it is necessary to introduce very substantial heating of the ICM. In the upper panel of Fig 2, we show the T-L relation at $`z=0`$ in a model with a heating efficiency $`ϵ_{\text{sn}}=1.3h^{1/2}`$. (Note that we calculate all our models for $`h=0.7`$, but then assume the scaling of $`ϵ_{\text{sn}}`$ with $`h`$ that was derived in the previous section.) This value of $`ϵ_{\text{sn}}`$ is already larger than can be accounted for by supernova feedback alone, if the IMF has the conventional solar neighbourhood form, even if there are no radiative energy losses. This suggests that a significant contribution from AGN is probably also required. If the heating produced by galaxies is smaller than this, the model predictions at the bright end fall too steeply with decreasing luminosity. Even with an efficiency of $`ϵ_{\text{sn}}=1.3h^{1/2}`$, the predicted T-L relation seems somewhat too steep for the most luminous clusters ($`L_\mathrm{X}>10^{44}h^2\text{ergs}\text{s}^1`$). These clusters are not much affected by this amount of heating and tend to follow the self-similar slope. To bring the most luminous clusters into line with the observed T-L slope requires that the energy injection efficiency be increased to $`ϵ_{\text{sn}}=2.0h^{1/2}`$. However, the model then fails to reproduce the presence of X-ray luminous clusters with temperatures below 1 keV (Fig. 2, lower panel). Thus, it seems that the heating efficiency needs to be slightly greater in the progenitor haloes of the most massive clusters, in order to produce the best match to the T-L relation. This might be the case if galaxy formation and/or AGN activity were even more strongly biased to high density regions than in the Cole et al. model.

The model results show considerable scatter around the median relation, which arises from two sources. Firstly, haloes collapse over a range of redshifts, leading to some variation in core density. Secondly, the efficiency of galaxy formation varies from halo to halo, resulting in considerable scatter in the level of heating. The scatter in the model is in reasonably good agreement with the observational data, although it fails to encompass a small number of clusters with high temperature and low X-ray luminosity. The transient effects of cluster mergers are not included in the present model.

The free parameters of the model have now been fixed to match the present-day T-L relation, and so the evolution to higher redshift provides a test of the model. As discussed in the previous section, the evolution of the T-L relation is determined by a competition between the increasing density of collapsed structures, the temperature distribution of the clusters and the relative importance of the excess energy. The thick dashed line in Fig. 2 shows the predicted median T-L relation at $`z=0.5`$. There is little evolution in this relation, consistent with presently available data on distant clusters (Mushotzky & Scharf 1997; Fairley et al. 2000). There is a tendency in the model for clusters of a given temperature to become more X-ray luminous at high redshift, but the trend is too weak to be rejected on the basis of currently available data. Fairley et al. (2000) have analysed the evolution of the T-L relation in a large sample of clusters from $`z=0.2`$ to 0.8. They fit their results to a parameterised form, $`LT^{3.15}(1+z)^\eta `$, and derive $`\eta =0.60\pm 0.38`$ for an open $`\mathrm{\Omega }_0=0.3`$ universe. This corresponds to $`\eta =0.75\pm 0.48`$ in our flat cosmology, since the luminosities inferred from the data are then greater. At $`T=5\text{keV}`$ our model produces a factor of 1.86 increase in the median cluster luminosity over the redshift interval 0.0 to 0.5, corresponding to $`\eta =1.54`$. Thus, the evolution predicted by our model is compatible (at $`1.6\sigma `$) with that observed by Fairley et al.

### 3.2 The X-ray Luminosity Function

The heating model provides a good description of the present-day T-L relation, and can account for its observed lack of evolution. We now consider the X-ray luminosity function (XLF). Since the galaxy formation model generates a statistical sample of haloes, the X-ray luminosity function can be readily calculated by combining the different halos with appropriate weights. In Fig. 3 we show the predicted luminosity functions at $`z=0`$ (solid line) and $`z=0.5`$ (dotted line). We show the luminosity function that is derived without heating (ie., $`ϵ_{\text{sn}}=0`$) in the top panel. The middle and lower panels correspond to the values of the efficiency parameter, $`ϵ_{\text{sn}}=1.3h^{1/2}`$ and $`ϵ_{\text{sn}}=2h^{1/2}`$ respectively, chosen to match the observed temperature-luminosity relation. These predictions are compared to the observed local luminosity function derived by Ebeling et al (1997) from the ROSAT all-sky “Bright Cluster” survey (BCS).

Without heating, the model is an extremely poor fit to the observed data: this is expected since we have chosen the CDM power spectrum to match the observed cluster temperature function. In this case the luminosities of low temperature haloes are too high, and this is reflected in the luminosity function which is too steep. Since the available XLF data are restricted to relatively bright clusters, we expect to obtain the best fit with $`ϵ_{\text{sn}}=2.0h^{1/2}`$ rather than with $`ϵ_{\text{sn}}=1.3h^{1/2}`$. This is indeed the case, although even for $`ϵ_{\text{sn}}=2.0h^{1/2}`$ the model luminosity function is still somewhat too steep. The discrepancy can be traced back to the slight bend in the T-L relation seen in Fig. 2, at the temperature at which the effect of the injected energy becomes significant. The fit could be fine-tuned by making the energy input increase more strongly with halo mass (eg. if galaxy formation were more efficient in proto-cluster regions), or by adjusting the cosmological parameters. For example, adopting $`\sigma _8=0.73`$ and $`\mathrm{\Gamma }=0.07`$ reduces the number of small mass haloes while boosting the abundance of the highest mass objects. This gives a significantly improved match to the luminosity function, but such a small value of $`\mathrm{\Gamma }`$ may not be compatible with measurements of large-scale galaxy clustering (Peacock & Dodds, 1994, Hoyle et al., 1999, Eisenstein & Zaldarriaga, 2000).

Below the limits probed by the BCS data, the model predicts a significant flattening of the luminosity function. This is an unavoidable consequence of energy injection: in low mass haloes, most of the gas is ejected, resulting in very low luminosities and a ‘stretching’ of the luminosity function in this region. The space density of low-luminosity ($`L_\mathrm{X}<10^{42}h^2\text{ergs}\text{s}^1`$) systems is therefore a strong test of this model. The absence of luminous haloes around spiral galaxies reported by Benson et al. (2000) supports this aspect of the model.

The evolution of the luminosity function is another important test of the model. The dotted line in Fig 3 shows the XLF at $`z=0.5`$. This evolves very little relative to the present-day function. The results from the weak evolution of the temperature function in this cosmological model (Eke et al. 1998) combined with the weak trend of increasing luminosities with higher redshift at fixed temperature seen in Fig 2. The model predictions compare very favourably with the available measurements from deep ROSAT surveys (eg. Scharf et al. 1997; Rosati et al. 1998; Vikhlinin et al. 1998; Nichol et al. 1999; Jones et al. 2000) which show no significant evolution of the luminosity function over the redshift range $`0<z<0.8`$. The evolution seen at the bright end is, however, sensitive to the CDM power spectrum adopted. For instance, the $`\sigma _8=0.73`$, $`\mathrm{\Gamma }=0.07`$ model discussed above would imply that the most massive clusters ($`L_\mathrm{X}>5\times 10^{44}h^2\text{ergs}\text{s}^1`$) should have significantly lower space density at $`z=0.5`$ than at the present day. It is unclear whether this is supported by current X-ray data (see Jones et al. 2000 for a discussion).

### 3.3 X-ray Emission in the High-Redshift Universe

We can use the model to predict the evolution of the X-ray emission from haloes out to high redshifts ($`z>2`$). The Cole et al. (2000) model of galaxy formation and evolution matches reasonably well observations of the evolution of the cosmic star formation rate over these long look-back times. We can thus predict the evolution of the supernova heating rate out to very high redshift, as is required in order to model the evolution of the XLF at high redshifts. We focus on the $`ϵ_{\text{sn}}=2.0h^{1/2}`$ model in what follows, because this produces the best fit to the present-day XLF.

The model predictions for the T-L relation and the XLF at $`z=0`$ and $`z=2`$ are shown in Fig. 4. At a given temperature, high-$`z`$ clusters are substantially more luminous than their present-day counterparts. However, because of hierarchical clustering, high temperature systems are increasingly rare at high redshift. At $`z=2`$, this decline in abundance offsets the modest increase in X-ray luminosity at given $`T`$. As a result, even at $`z=2`$ the luminosity function is predicted to be close to that observed at the present-day. This is consistent with observational limits on the contribution of clusters to the X-ray background (e.g. Burg et al. 1993 and Wu et al. 1999b).

### 3.4 The Epoch of Galaxy Formation

We have argued that the slope of the temperature-luminosity relation reflects the energy input from galaxies and AGNs. Now we examine whether the evolution of clusters can be used to constrain the epoch at which this heating occurs. We contrast the galform model (with $`ϵ_{\text{sn}}=2.0h^{1/2}`$) with two simple ad hoc heating models. In the first, the heating occurs at a constant rate over cosmic time (model A); in the second, the heating occurs only at high redshift so that the excess energy remains constant below $`z=2.0`$ (model B). Model B is intended to mimic the effect of “pre-heating” of the intergalactic medium as in the model proposed by (eg.) Balogh et al. (1999). The total energy injection has been adjusted to match the present-day XLF of the galform model. Models A and B give XLFs at $`z=0`$ which are almost identical to that from GALFORM, when normalized by this procedure.

We contrast these two simple models with our fiducial model based on hierarchical galaxy formation in Figure 5. The upper panel shows the median T-L relations predicted by each of the models at $`z=2`$. Similar but less pronounced differences exist at $`z=1`$ and at $`z=0.5`$. The models diverge at low luminosities since the relative effect of the injected energy is greatest for small clusters. It is not surprising that the differences between the models at the bright end are small, as the heat input is fairly unimportant for these clusters. As expected, the two simple heating models bracket the galform model, although the latter seems closer to model B in which the heating occurs at high redshift. This reflects the fact that galaxy formation in proto-cluster regions is accelerated relative to that in an average region of the universe, simply due to the higher density there.

The lower panel in Fig. 5 shows the differences between the luminosity functions at $`z=2`$ for the three models. As expected from the upper panel, the two simple models again bracket the behaviour of galform. The luminosity function for the constant heating model shows stronger positive evolution (ie. a higher number density at higher redshift) than the model in which the heating has already occurred before this epoch. These differences offer an interesting possibility for determining the epoch of galaxy formation: if it becomes possible to distinguish between different models for the redshift dependence of the heating using X-ray observations, this will then provide a strong constraint on the epoch at which most of the stars in the universe formed.

## 4 Discussion

As we have shown, a model in which the intracluster gas is heated as galaxy formation proceeds provides a good explanation for the slope of the T-L relation and for the evolution of the X-ray luminosity function. The problem with associating this heating with supernovae is the large amount of energy that is required, between 1.3 and $`2.0\times 10^{49}h^{1/2}\text{erg}`$ per solar mass of stars formed. This corresponds to an energy of 0.6 – 1.0 keV per particle in the intracluster medium. This is comparable to the energy injection requirement (1 – 2 keV per particle) derived by Wu et al. (1999a), showing that the conclusions about the energetics do not depend greatly on the details of the heating model. Even with optimistic assumptions about the supernova rate, this amount of heating would require that the energy of the supernova explosions be transferred to the intergalactic plasma with an efficiency close to unity. This seems unrealistic.

An alternative source for heating the ICM is AGNs and quasars. If AGN activity is closely linked to the fuelling of star formation, then such activity will effectively enhance the value of $`ϵ_{\text{sn}}`$. In this case, the effect of AGN heating can easily be incorporated into our model. Our conclusions would be unchanged apart from the interpretation of $`ϵ_{\text{sn}}`$. If, on the other hand, the energy input comes predominantly from the most powerful AGN early in the history of the universe, it would be more appropriate to treat the energy injection as a uniform preheating of the gas prior to gravitational collapse of the dark matter halos. Assuming the energy sources were sufficiently uniformly distributed, the effects of the heating might be better modeled by assuming that the gas entropy is raised to some uniform value before collapse (eg. Evrard & Henry 1991; Navarro, Frenk & White 1995; Bower 1997; Balogh, Babul & Patton 1999; Valageas & Silk 1999). A possible problem of this scheme is the high temperature it implies for the diffuse IGM. For example, Balogh, Babul & Patton (1999) require a temperature of $`1.8\times 10^6`$ K for a preheating epoch of $`z=3`$ in our $`\mathrm{\Lambda }`$CDM cosmology. This is in stark contrast to the IGM temperature estimated from the Doppler widths of Ly-$`\alpha `$ forest lines. For example, Theuns et al (2000) estimate $`T_{\mathrm{IGM}}15,000`$ K at this redshift. Thus, unless the clouds giving rise to the Ly-$`\alpha `$ forest or the precursor gas of the ICM are atypical, a model in which the heating occurs within already virialised haloes seems preferable.

Dark haloes build up by mergers of smaller progenitor halos. Although we include energy input from the complete history of star formation in each halo, we assume that the way this heating is distributed between the earlier haloes is unimportant. In particular, we ignore any dependence of the energy released during the gravitational collapse of a halo on the distribution of the gas in the progenitor halos. If the gas has already been heated by supernovae or AGN in the progenitor halos, then the dynamics of the collapse of the gas and its shock heating will be modified. It is quite possible that this effect could give rise to an “amplification” of an initial energy excess, thus easing the requirements on the heating efficiency. Another limitation of our approach is that we have only included the effects of cooling in an approximate way, by removing from the ICM the gas which should have cooled to low temperature, based on a calculation which does not explicitly include the effects of the energy injection. However, the cooling rates will be modified by the effects of the energy injection, and conversely some of the excess energy may be lost by radiative cooling. It is clearly vital to understand all of these processes better, which can best be done through well-targeted numerical simulations (eg. Pearce et al. 2000).

Finally, we must recall that galaxy formation and X-ray evolution have not been treated in a fully self-consistent fashion in this paper. We have taken the successful galform model of galaxy formation, with the same parameters as in Cole et al. (2000), and used it to predict the energy injection and thus the evolution of cluster and group X-ray properties. In practice, we should use the methods developed here to calculate the gas density profiles of all haloes at each epoch, as modified by the energy injection, compute gas cooling rates using these modified profiles, and then calculate the energy injection from the resultant star formation histories. This represents a large computational overhead on the standard galform model, but is clearly an important next step to take.

## 5 Conclusions

In this paper, we have addressed the problem of why the observed properties of X-ray clusters do not conform to simple scaling relations. In particular, we have considered why the observed correlation between X-ray temperature and X-ray luminosity is significantly shallower than the adiabatic scaling solution, while the X-ray luminosity function evolves less rapidly than predicted in popular cold dark matter cosmologies. First, we argued that the effects of gas cooling in clusters (which break the scaling relations) do not resolve the problem. We then considered the heating of the intracluster gas by the energy released during galaxy formation, by combining the semi-analytic model of Cole et al. (2000) with a simple model for the radial profile of the intracluster gas. Our main conclusions, applicable in the $`\mathrm{\Lambda }`$CDM cosmology, may be summarized as follows:

* Heat input into the intracluster gas by processes associated with the formation of cluster galaxies, such as supernovae and/or AGN winds, will flatten the slope of the temperature-luminosity relation. The combined model gives a reasonable match to the observations if energy is injected at a level of 1.3–2.0 $`\times 10^{49}h^{1/2}`$ ergs per solar mass of stars formed (or, equivalently, 0.6–1 keV per particle in the intracluster medium). Values within this range produce broadly acceptable models, but lower values result in a better match to groups with $`T1\text{keV}`$, while higher values produce a better match to the most massive clusters.
* The interplay between the ongoing energy injection from galaxies and the growth of clusters by hierarchical clustering causes the $`TL`$ relation to evolve little to moderate redshifts. This is consistent with recent data based on ASCA observations.
* The present-day X-ray luminosity function in the model approximately matches observations, but the model over-produces low luminosity groups and under-produces very luminous clusters. Fine tuning the cosmological parameters or other details of the model may remove these discrepancies.
* Similar factors to those that regulate the evolution of the $`TL`$ relation result in only weak evolution of the luminosity function to $`z=0.5`$. This too is consistent with current data.
* The properties of clusters at high redshift provide a test of the model, since all of the free parameters are fixed to achieve agreement with present-day data. In particular, the model predicts little evolution in the X-ray luminosity function even out to $`z=2`$. The predicted near constancy of the luminosity function is testable with the current generation of X-ray satellites.
* The main difficulty of our model is that it requires an amount of energy injected per unit mass of stars formed which is comparable to the total energy available from supernovae. This would require the supernova explosion energy to couple to the intracluster gas with very high efficiency, with minimal losses by radiative cooling during the expansion of the supernova remnants through the ISM of the host galaxies. However, additional energy sources associated with galaxy formation may also contribute, such as the mechanical energy liberated by AGN winds. Alternatively (or additionally), an initial heat input to the intracluster medium might be amplified during the build-up of the cluster by mergers. Detailed numerical simulations are required to test the effectiveness of this process.

Our work demonstrates that the shape and evolution of the X-ray luminosity function and T-L relation are potentially powerful probes of the mode and efficiency of galaxy formation. Future observations with Newton and Chandra should be able to test these ideas.

## Acknowledgements

Thanks to Ed Lloyd-Davies and Trevor Ponman for extensive discussions and to Simon White for valuable suggestions. This project has made extensive use of Starlink computing facilities and was supported by the PPARC rolling grant on ‘Extra-galactic astronomy and cosmology’ at Durham, and by the EC TMR Network on ‘Galaxy Formation and Evolution’. CSF acknowledges the support of a Leverhulme Research Fellowship. CGL acknowledges support at SISSA from COFIN funds from MURST and funds from ASI.

## Appendix A The Effect of Cooling on the T-L Relation

In the absence of radiative cooling and energy input from galaxy formation, and assuming that all clusters have density profiles which are simply rescaled versions of each other, then the bolometric X-ray luminosities of clusters should vary as

$$L_\mathrm{X}f_g^2\rho _{\mathrm{vir}}M\mathrm{\Lambda }(T)$$
(5)
where $`f_g`$ is the fraction of the cluster mass in the form of hot gas, $`\rho _{\mathrm{vir}}M/R_{\mathrm{vir}}^3`$ is the mean total density within the virial radius, and the cooling rate per unit volume is proportional to $`\rho ^2\mathrm{\Lambda }(T)`$, with $`\mathrm{\Lambda }(T)T^{1/2}`$ for bremsstrahlung radiation. Assuming that the density depends on the collapse redshift $`z_\mathrm{f}`$ of the cluster as $`\rho _{\mathrm{vir}}\rho _0(1+z_\mathrm{f})^3`$ ($`\rho _0`$ being the present mean density of the universe), and using equation (3) for the temperature, we obtain the scaling law

$$L_\mathrm{X}f_g^2\rho _0^{1/2}(1+z_\mathrm{f})^{3/2}T^2$$
(6)
(eg. Kaiser 1986, 1991; Evrard & Henry 1991; Bower 1997; Kay & Bower 1999). We have explicitly included the dependence on the collapse redshift of the cluster, $`z_\mathrm{f}`$, to make it clear that the scaling depends on this rather than on the redshift at which the cluster is observed.

As we have discussed, the T-L relation implied by eqn. (6) is too steep compared with the observed luminosities and temperatures of clusters. Eqn. (6) suggests that the relation might be made shallower if lower temperature clusters had systematically lower collapse redshifts. In hierarchical models, however, smaller mass clusters are expected to collapse at higher redshifts — the opposite to the required trend.

In this appendix, we will use simple scaling arguments to argue that the effects of radiative cooling by itself are not sufficient to bring the predicted T-L relation into line with the observed one. The effects of radiative cooling on the density and temperature profiles, and thus X-ray luminosities, of spherical clusters have been the subject of various analytical (e.g. Bertschinger 1989) and numerical (e.g. Lufkin et al. 2000) investigations. These studies show that outside the cooling radius $`r_{\mathrm{cool}}`$, defined such that the local radiative cooling time is equal to the age of the system, the density and temperature are almost unchanged from their initial values. Inside $`r_{\mathrm{cool}}`$, the gas flows in and then drops out of the flow completely due to radiative cooling. In the case that the initial density profile is steep, this results in a reduction in the gas density within $`r_{\mathrm{cool}}`$, and thus also a reduction in $`L_\mathrm{X}`$, On the other hand, if the initial density profile is very shallow, then the gas density and $`L_\mathrm{X}`$ may be boosted.

Consider first the simple case that the gas density profile is that of a singular isothermal sphere, $`\rho (r)\rho _0(1+z_\mathrm{f})^3(r/R_{\mathrm{vir}})^2`$, with $`T=T_{\mathrm{vir}}`$. We define $`r_{\mathrm{cool}}`$ as the radius where the local cooling time equals the age of the universe $`t_\mathrm{H}`$ (we choose the age of the universe rather than that of the halo in order to derive the maximum effect of cooling). Thus:

$$t_{\mathrm{cool}}(r_{\mathrm{cool}})\rho (r_{\mathrm{cool}})^1T(r_{\mathrm{cool}})/\mathrm{\Lambda }(T)t_\mathrm{H}.$$
(7)
(Note that we suppress the dependence on the gas fraction $`f_g`$ in this and the following equations.) The fraction of the gas which is able to cool is then

$$f_{\mathrm{cool}}=\frac{r_{\mathrm{cool}}}{R_{\mathrm{vir}}}\left(\rho _0(1+z_\mathrm{f})^3T^{1/2}t_\mathrm{H}\right)^{1/2}.$$
(8)
The self-similar cooling flow solutions of Bertschinger (1989) show that the density profile flattens to $`\rho r^{3/2}`$ within $`r_{\mathrm{cool}}`$, so the X-ray luminosity scales as

$$L_\mathrm{X}r_{\mathrm{cool}}^3\rho (r_{\mathrm{cool}})^2\mathrm{\Lambda }(T)\frac{T^2\rho _0^{1/2}(1+z_\mathrm{f})^{3/2}}{f_{\mathrm{cool}}}.$$
(9)
The greater the fraction of gas that is able to cool, the more the luminosity is reduced below that of eqn. (6). Substituting equation (8) for $`f_{\mathrm{cool}}`$, we then find

$$L_\mathrm{X}T^{9/4}t_\mathrm{H}^{1/2}.$$
(10)
This only slightly improves the match to the data compared to the case of no cooling. For high temperature clusters, the relation becomes slightly shallower than before, but the effect is not sufficient. Moreover, in cooler clusters, the emissivity is enhanced by recombination radiation and the slope of the T-L relation becomes steeper again.

The above equations apply in the case that the initial gas density profile is singular. If, as is more realistic, the gas profile initially has a core of radius $`r_c`$ within which the density is constant, then the behaviour is modified. The cooling time will be constant within the core, so the density profile will remain almost unchanged until the age of the system is equal to this cooling time, and $`L_\mathrm{X}`$ will approximately follow equation (6). As a cooling flow starts up within the core, the X-ray luminosity may be enhanced, but at most by a factor of a few. The cooling radius $`r_{\mathrm{cool}}`$ will then grow beyond the core radius $`r_c`$, and $`L_\mathrm{X}`$ will converge towards the behaviour in equation (10). Thus, $`L_\mathrm{X}`$ will scale approximately as (6) for $`t_\mathrm{H}<t_{\mathrm{cool}}(r_c)`$ and as (10) for $`t_\mathrm{H}>t_{\mathrm{cool}}(r_c)`$, and our previous conclusions about the T-L relation not being reproduced remain unchanged.

An estimate of the effects of cooling based on our $`\beta `$-model approach is shown by the dot-dashed line in Figure 2. Starting from the default halo gas profile, we have calculated the gas mass within the cooling radius for each of the simulated haloes. This gas is removed, and the $`\beta `$-profile is then adjusted so that the remaining gas mass is distributed within the virial radius using the same boundary conditions as discussed in §2, keeping $`r_c`$ constant but reducing $`\beta `$. (If instead $`r_c`$ is varied and $`\beta `$ kept constant, the resulting T-L relation is very similar.) The T-L relation that results is close to that predicted by the scaling arguments discussed above, and fails to match the observed data.


format: text

text:
# References

Breit-Wigner to Gaussian transition in strength functions

V.K.B. Kotaa and R. Sahua,b

aPhysical Research Laboratory, Ahmedabad 380 009, India

bPhysics Department, Berhampur University, Berhampur 760 007, India

ABSTRACT: Employing hamiltonians defined by two-body embedded Gaussian orthogonal ensemble of random matrices(EGOE(2)) plus a mean-field producing one-body part, strength functions (for states defined by the one-body part) are constructed for various values of the strength of the chaos generating two-body part. Numerical calculations for six and seven fermion systems clearly demonstrate Breit-Wigner to Gaussian transition, in the chaotic domain, in strength functions as found earlier in nuclear shell model and Lipkin-Meshkov-Glick model calculations.

PACS number(s): 21.10.Pc, 24.60.Lz, 21.60.Cs

Strength functions (also called local density of states (LDOS) in literature) of simple modes are basic ingredients of many particle systems such as atomic nuclei . In the last few years, with developments in quantum chaos , nature of strength functions for finite isolated interacting many particle systems are being investigated using a variety of models: (i) Zelevinsky et al using nuclear shell model ; (ii) Wang et al using the three-orbital Lipkin-Meshkov-Glick model ; (iii) Borgonovi et al using a symmetrized coupled two-rotor model ; (iv) Benet et al using a chaotic model of two coupled quartic oscillators . Results from the shell model and Lipkin-Meshkov-Glick model studies clearly showed that strength functions exhibit, for interacting particle systems, Breit-Wigner to Gaussian transition. Moreover it is seen that the transition takes place much after the onset of chaos in energy levels (i.e. much after level fluctuations start following GOE). In order to establish that this transition is generic to interacting particle systems, we carried out one plus two-body embedded Gaussian orthogonal ensemble (EGOE(1+2)) calculations and the results are reported in this short communication (see also ).

Given a compound state $`\varphi _k`$, the probability of its decay into stationary states $`\psi _E`$ (generated by $`H`$) is given by $`\left|\varphi _k\psi _E\right|^2`$. Then the strength function $`F_k(E)`$ is,

$$|\varphi _k=\underset{E}{}C_k^E|\psi _EF_k(E)=\underset{E^{}}{}\left|C_k^E^{}\right|^2\delta (EE^{})=\delta (HE)^k$$
(1)
The standard form , normally employed in many applications in nuclear physics, for strength functions is the Breit-Wigner (BW) form characterized by a spreading width $`\mathrm{\Gamma }_k`$,

$$F_{k:BW}(E)=\frac{1}{2\pi }\frac{\mathrm{\Gamma }_k}{(E\overline{E_k})^2+\mathrm{\Gamma }_k^2/4}$$
(2)
where $`\overline{E_k}=kHk=_{\mathrm{}}^{\mathrm{}}F_k(E)E𝑑E`$. Starting with $`H=h(1)+\lambda V(2)`$ (a one plus two-body interaction as in the shell model) with $`h(1)`$ defining the $`\varphi _k`$ states, it is easily seen that the assumptions that give $`F_{k:BW}(E)`$ will break down when the mixing is strong. Thus it is expected that the form of $`F_k(E)`$ will be different from BW for large $`\lambda `$. In order to find the form of strength functions in large $`\lambda `$ limit, it is useful to think of $`\varphi _k`$ as a compound state generated by the action of a transition operator $`𝒪`$ on a state $`\psi _{E_k}`$ (for example ground state),

$$\begin{array}{c}|\varphi _{E_k}=\frac{𝒪|\psi _{E_k}}{\left[\psi _{E_k}𝒪^{}𝒪\psi _{E_k}\right]^{1/2}};\varphi _{E_k}\varphi _{E_k}=1\\ |\varphi _{E_k}=\underset{E}{}C_{E_k}^E|\psi _EF_{k;𝒪}(E)=\frac{𝒪^{}\delta (HE)𝒪\delta (HE_k)}{𝒪^{}𝒪\delta (HE_k)}\end{array}$$

The numerator in the expression for $`F_k(E)`$ in (3) is nothing but the bivariate strength density generated by the operator $`𝒪`$  and then it is easy to recognize that $`F_k(E)`$ is the conditional density of this bivariate strength density. It is known  that EGOE($`k`$), the embedded GOE of $`k`$-body interactions (see Ref.  for the definition of EGOE(k)), in general gives bivariate Gaussian strength densities. Therefore, for EGOE($`2`$), which models a generic two-body hamiltonian, the strength functions take Gaussian form (conditional density of a bivariate Gaussian is a Gaussian). However, for generic interacting particle systems, it is more appropriate to consider one plus two-body EGOE(1+2) hamiltonians $`\{H\}=h(1)+\lambda \{V(2)\}`$ where $`h(1)`$ is the one-body mean-field producing part and $`\{V(2)\}`$ is the chaos generating two-body part; note that $`\{V(2)\}`$ is GOE in two-particle space with unit matrix elements variance and $`\lambda `$ is the interaction strength. Starting with EGOE(1+2) it is to be expected that for sufficiently large values of $`\lambda `$, EGOE(2) description should be valid and therefore (3) gives the shape of the strength function to be Gaussian for large $`\lambda `$; this result is indeed seen in numerical calculations. Now the important questions are: (i) how $`F_k(E)`$ changes as $`\lambda `$ is varied and (ii) from which value of $`\lambda `$ strength functions take Gaussian form.

For further understanding of the nature of $`F_k(E)`$, we performed EGOE(1+2) calculations in 924 dimensional space generated by six fermions ($`m=6`$) in twelve single particle states ($`N=12`$) and similarly in the 3432 dimensional $`N=14`$ and $`m=7`$ space. The single particle energies employed are $`ϵ_i=i+(1/i)`$, $`i=1,2,\mathrm{},N`$ as in . For various values of $`\lambda `$ in $`h(1)+\lambda \{V(2)\}`$, strength functions are constructed by choosing the $`\varphi _k`$ states to be the mean-field $`h(1)`$ states defined by the distribution of $`m`$ particles in the $`N`$ single particle states; their energies $`_k`$ are $`_k=k|h(1)+\lambda V(2)|k`$. In the calculations $`E`$ and $`_k`$ are zero centered for each member and scaled by the spectrum ($`E`$’s) width $`\sigma `$; $`\widehat{E}_k=(_kϵ)/\sigma `$ and $`\widehat{E}=(Eϵ)/\sigma `$. For each member $`|C_k^E|^2`$ are summed over the basis states $`|k`$ in the energy window $`\widehat{E}_k\pm \mathrm{\Delta }`$ and then ensemble averaged $`F_{\widehat{E}_k}(\widehat{E})`$ vs $`\widehat{E}`$ is constructed as a histogram; the value of $`\mathrm{\Delta }`$ is chosen to be 0.025 for $`\lambda 0.1`$ and beyond this $`\mathrm{\Delta }=0.1`$. Results for $`\widehat{E}_k=0`$ are shown for $`\lambda =`$ 0.05, 0.1, 0.15, 0.2, 0.3 and 0.5 in Fig. 1; in the plots $`F__k(\widehat{E})𝑑\widehat{E}=1`$. From the figures it is clearly seen that there is BW to Gaussian transition in $`F_k(E)`$. A measure for this transition is,

$$R(\lambda )=\frac{{\displaystyle \underset{i}{}}\left[F_k^{(\lambda )}(E_i)F_{k:BW}(E_i)\right]^2}{{\displaystyle \underset{i}{}}\left[F_{(k:BW)}(E_i)F_{k:ED}(E_i)\right]^2}$$
(3)
In (4) $`E_i`$ are defined by the center of each bin in the histogram representing $`F_k^{(\lambda )}(E)`$ and $`ED`$ corresponds to Gaussian with Edgeworth corrections; the ED incorporates  skewness ($`\gamma _1`$) and excess ($`\gamma _2`$) corrections. Strength function $`F_k^{(\lambda )}(E)`$ is BW for $`R=0`$ and Gaussian (or ED) for $`R=1`$. The value of the interpolating parameter $`\lambda =\lambda _{F_k}`$ for onset of the transition from BW to Gaussian is taken to be $`R(\lambda _{F_k})=0.7`$. For $`\lambda =`$ 0.05, 0.1, 0.15, 0.2, 0.3, 0.5 the values of $`R`$ are 0.13, 0.06, 0.21, 0.7, 0.9, 1.03 for the six particle example. For the seven particle example the corresponding numbers are 0.15, 0.02, 0.45, 0.75, 0.93, 1.03 respectively. Thus, with the measure $`R`$ in (4), $`\lambda _{F_k}0.2`$ for both the 6 and 7 particle cases.

In Fig. 2 shown are EGOE(1+2) results for various values of $`\lambda `$ for the nearest neighbour spacing dustribution and $`\overline{\mathrm{\Delta }_3}(L)`$, $`0L40`$ statistic. The $`\lambda _{F_k}`$ deduced from the results in Fig. 1 should be compared with $`\lambda =\lambda _c0.08`$ and 0.06, derived via Fig. 2, for six and seven particle examples repectively for onset of chaos in level fluctuations. The $`\lambda _c`$ numbers are consistant with $`(Nm^2)^1`$ variation as expected from the theory given in . It should be noted that in the $`(2s,1d)^{m=12,J^\pi T=0^+0}`$ shell model example considered in , $`\lambda _c0.3`$ and $`\lambda _{F_k}0.6`$. Thus the BW form for $`F_k(E)`$, which begins some what before $`\lambda `$ approches $`\lambda _c`$ (for $`\lambda <<\lambda _c`$ the strength functions are delta functions with perturbative corrections), extends much into the chaotic domain (defined by $`\lambda >\lambda _c`$) and the transition to Gaussian shape takes place in the second layer defined by $`\lambda _{F_k}`$ ($`\lambda _{F_k}>>\lambda _c`$). As the $`\lambda _{F_k}`$ in our two examples did not show much variation with $`(N,m)`$, it appears that $`\lambda _{F_k}`$ will have weak dependence on $`m`$ unlike $`\lambda _c`$. Perhaps for sufficiently large $`m`$, the transition to Gaussian form takes place in the thermalization regime discussed by Flambaum and Izrailev  (see also ). It should be remarked that EGOE($`2`$) gives Gaussian form for state densities for sufficiently large $`m`$ and its extension to partial state densities (which are sums of strength functions) is indeed the basis of statistical nuclear spectrocopy .

In conclusion, the results of the present EGOE(1+2) study (Figs. 1,2) and the earlier shell model  and Lipkin-Meshkov-Glick model  analysis establish firmly the Breit-Wigner to Gaussian transition for strength functions, in the chaotic domain, in interacting particle systems such as atomic nuclei.

—————

Acknowledgements This work has been partially supported by DST(India).

FIGURE CAPTIONS

Fig. 1 Strength functions for EGOE(1+2) for various values of the interaction strength $`\lambda `$: (i) for a system of 6 fermions in 12 single particle states with 25 members; (ii) for a system of 7 fermions in 14 single particle states (due to computational constraints, here only one member is considered just as in ). In the figure, the histograms are EGOE(1+2) results and continuous curves are BW fit. For the 6 fermions case, the dotted curves are Gaussian for $`\lambda 0.15`$ and Edgeworth corrected Gaussians (ED) for $`\lambda >0.15`$. Similarly for the 7 fermions case, the dotted curves are Gaussian for $`\lambda 0.1`$ and Edgeworth corrected Gaussians (ED) for $`\lambda >0.1`$. See text for further details.

Fig. 2 Nearest neighbour spacing distribution $`P(S)`$ vs $`S`$ and Dyson-Mehta $`\overline{\mathrm{\Delta }_3}(L)`$ statistic for $`0L40`$ for various values of the interaction strength $`\lambda `$ in EGOE(1+2). For $`P(S)`$, histograms are EGOE(1+2) results, dashed curves are Poisson and continuous curves are Wigner distribution. For $`\overline{\mathrm{\Delta }_3}(L)`$, filled circles are EGOE(1+2) results, dashed curves are Poisson and continuous curves are for GOE. Results are shown for the six and seven particle examples considered in Fig. 1.


text:
# Edwards measures for powders and glasses

\[

## Abstract

Can one construct a thermodynamics for compact, slowly moving powders and grains? A few years ago, Edwards proposed a possible step in this direction, raising the fascinating perspective that such systems have a statistical mechanics of their own, different from that of Maxwell, Boltzmann and Gibbs, allowing us to have some information while still ignoring dynamic details.

Recent developments in the theory of glasses have come to confirm these ideas within mean-field. In order to go beyond, we explicitly generate Edwards’ measure in a model that, although schematic, is three dimensional. Comparison of the results thus obtained with the irreversible compaction data shows very good agreement. The present framework immediately suggests new experimental checks.

\]

The classical way to go from the microscopic dynamics to statistical mechanics proceeds in two steps: one first identifies a distribution that is left invariant by the dynamics (e.g. the microcanonical ensemble), and then assumes that this distribution will be reached by the system, under suitable conditions of ’ergodicity’. For granular systems this approach seems doomed from the outset: because energy is lost through internal friction, and gained by a non-thermal source such as tapping or shearing, the dynamical equations do not leave the microcanonical or any other known ensemble invariant. Moreover, the compaction dynamics is extremely slow and does not approach any stationary state on experimental time scales. This raises the question of characterizing the typical configurations or the region of phase space visited dynamically.

The proposal of Edwards and collaborators  is to use an alternative distribution for very gently vibrated or sheared granular systems, with the static situation as a limiting case. It may be summarized as follows: given a certain situation attained dynamically, physical observables are obtained by averaging over the usual equilibrium distribution at the corresponding volume, energy, etc. but restricting the sum to the ‘blocked’ configurations defined as those in which every grain is unable to move. This definition leads immediately to an entropy (in the glass literature a ‘complexity’) $`S_{edw}`$, given by the logarithm of the number of blocked configurations of given volume, energy, etc., and its corresponding density $`s_{edw}S_{edw}/N`$. Associated with this entropy are the state variables such as ‘compactivity’ $`X_{edw}^1=\frac{}{V}S_{edw}(V)`$ and ‘temperature’ $`T_{edw}^1=\frac{}{E}S_{edw}(E)`$.

That configurations with low mobility should be relevant in a jammed situation is rather obvious, the strong assumption here is that, apart from the usual statistical weights, all blocked configurations are treated as equivalent — any extra weight of dynamical origin that might distinguish them is disregarded. The purpose of this letter is to argue that this ‘flatness’ assumption characterizing Edwards’ distributions is neither capricious (it leads to correct predictions for the compaction dynamics of a given class of systems), nor obvious (it does not apply to other classes of systems). To do this we devise a method to count the blocked configurations and compute averages over them.

Let us briefly summarize the state of the art. A first clue comes from exploiting the analogy between the settling of grains and powders, as when we gently tap a jar with flour to make space for more, and the aging of glassy systems : in both cases, the system remains out of equilibrium on all accessible time-scales, and displays very slow relaxations.

In the late eighties, Kirkpatrick et al.  recognised that a class of mean-field models contains, although in a rather schematic way, the essentials of glassy phenomena. When the aging dynamics of these systems was solved analytically, a feature that emerged was the existence of a temperature $`T_{dyn}`$ for all the slow modes (corresponding to structural rearrangements) . For our purposes here, $`T_{dyn}`$ can be defined by comparing the random diffusion and the mobility between two widely separated times $`t`$ and $`t_w`$ of any particle or tracer in the aging glass. Surprisingly, one finds in all cases an Einstein relation $`(r(t)r(t_w))^2=T_{dyn}\frac{\delta r(t)r(t_w)}{\delta f}`$, where $`r`$ is the position of the particle and $`f`$ is a constant perturbing field. While in an equilibrium system the fluctuation-dissipation theorem guarantees that the role of $`T_{dyn}`$ is played by the thermodynamic temperature, the appearance of such a quantity out of equilibrium is by no means obvious. $`T_{dyn}`$ is different from the external temperature, but it can be shown to have all other properties defining a true temperature .

As it turned out, despite its very different origin, this temperature matches exactly Edwards’ ideas: $`T_{edw}`$ and $`T_{dyn}`$ happen to coincide for mean-field glass models aging in contact with an almost zero temperature bath . In fact, given the energy $`E(t)`$ at long times, the value of any other macroscopic observable is also given by the flat average over all blocked configurations of energy $`E(t)`$. Within the same approximation, one can also treat systems that like granulars present a non-linear friction and different kinds of energy input, and the conclusions remain the same .

A first partial conclusion is then that Edwards’ scenario is at the very least correct within mean-field schemes and for very weak vibration or forcing. The problem that remains is to what extent it carries through to more realistic models.

In this direction, there have been recently studies  of Lennard-Jones glass formers from the perspective of the so-called ‘inherent structures’ (a partition of the phase-space in terms of the blocked configurations ). In this context a ‘flat weight’ assumption — similar in spirit but not quite equivalent to Edwards’ — also comes into question and is tested in various ways. Though there are caveats, the results are encouraging.

The path we follow is instead to construct the Edwards measure explicitly in the case of representative (non mean-field) systems, together with the corresponding entropy and expectation values of observables. We thus obtain results that are distinctly different from the equilibrium ones, and we can compare both sets with those of the irreversible compaction dynamics.

The first model we consider is the so-called Kob-Andersen (KA) model  that, though very schematic, reproduces rather well several aspects of glasses  and of granular compaction ; most important, this model is non mean-field. A particle can move to a neighbouring empty site, on a three dimensional lattice, only if it has less than four neighbours in the initial and in the final position. (In these ‘hard particle’ models the temperature is irrelevant, and we set it to one.) The dynamic rule guarantees that the equilibrium distribution is trivially simple since all the configurations of a given density are equally probable. However, at densities close to $`\rho _\mathrm{g}`$ ($`.88`$), the particle diffusion becomes extremely slow due to the kinetic constraints. In order to mimic a compaction (or aging) process without gravity, we simulate a ‘piston’ by freely creating and destroying particles only on the topmost layer with a chemical potential $`\mu `$ .

i) The dynamic measurements are taken as follows: starting from low density, we perform a slow compression by raising the chemical potential up to a high value $`\mu =3`$. Since the equilibrium density at $`\mu =3`$ is much larger than the jamming density $`\rho _\mathrm{g}`$, the system falls out of equilibrium and very slow compaction ensues. We record the density $`\rho (t)`$ and the spatial structure function $`g_{dyn}(r,t)`$ defined as the probability that two sites at distance $`r`$ are occupied. We also compute the dynamic temperature $`T_{dyn}`$ by comparing induced and spontaneous displacements. This is the set of observables we use for testing the different measures, which are obtained independently.

ii) In the equilibrium measure all configurations (whether they are blocked or not) have equal weight. It is easy to obtain the exact equilibrium entropy density per particle $`s_{equil}(\rho )=\rho \mathrm{ln}\rho (1\rho )\mathrm{ln}(1\rho )`$. Since $`T=1`$ as mentioned above, $`\frac{ds_{equil}(\rho )}{d\rho }=\mu `$. The equilibrium structure factor is easily seen to be a constant $`g_{equil}(r)=\rho ^2`$: indeed, one main advantage of this model is that it is particularly easy to compare small deviations from $`g_{equil}(r)`$, a notoriously difficult task in glassy systems.

iii) Finally, we obtain Edwards’ measure results as follows: we introduce an ‘auxiliary model’ in which particles have energy equal to one if the dynamic rule of the original model would allow them to move, and to zero otherwise. Performing simulated annealing of the auxiliary model at fixed number of particles is an efficient way to sample over the configurations with vanishing fraction of moving particles the Edwards ensembles structure function $`g_{edw}(r)`$, and to obtain (Figure 1) $`S_{edw}(\rho )`$ as the logarithm of the number of such configurations by thermodynamic integration of the energy of the auxiliary model with respect to its temperature. We then compute $`T_{edw}^1=\frac{1}{\mu }\frac{ds_{edw}(\rho )}{d\rho }`$.

We are now in a position to compare the long-time results of the out of equilibrium dynamics i) with those obtained with measures ii) and iii). Figure 2 shows a plot of the mobility $`\chi (t,t_w)=\frac{1}{3N}_{a=1}^3_{k=1}^N\frac{\delta (r_k^a(t)r_k^a(t_w))}{\delta f}`$ vs. the mean square displacement $`B(t,t_w)=\frac{1}{3N}_{a=1}^3_{k=1}^N(r_k^a(t)r_k^a(t_w))^2`$, testing the existence of a dynamical temperature $`T_{dyn}`$, in the compaction data ($`N`$ is the number of particles and $`a`$ runs over the spatial dimensions). The agreement between $`T_{dyn}`$ and the Edwards temperature $`T_{edw}`$, obtained from the blocked configurations as in Figure 1, is clearly excellent.

In Figure 3 we plot i) the long-time dynamic $`g_{dyn}(r,t)`$, ii) the equilibrium $`g_{equil}(r)=\rho ^2`$, and iii) the Edwards $`g_{edw}(r)`$ structure factors. The agreement between i) and iii) is good.

From the results shown so far, a picture emerges where the Edwards measure is able to correctly reproduce the sampling of the phase space generated by the out of equilibrium dynamics of this non mean-field model. We have found, however, that at short times or for excessively fast compressions, the quality of the agreement becomes worse, possibly due to heterogeneities. We refer to a longer, more technical paper for a discussion of these issues, as well as a study of other models (in particular the so-called Tetris model  for which one recovers the same conclusions as for the KA model) and more technical details on our numerical methods .

As already mentioned, Edwards’ construction can be inappropriate for certain models, even though they may have a logarithmically slow dynamics. As a representative example of this we consider the low temperature domain growth dynamics of a 3D Ising model in a weak random magnetic field, a model relevant to many physical problems . At large times the domain walls are pinned by the field, and the dynamics proceeds by thermal activation. The mean energy decreases slowly towards the ground state energy. In a large system, the long-time configurations are made of domains of ‘up’ and ‘down’ spins having similar volumes, the global magnetization being zero. This is quite different from the equilibrium configurations at the same energy, which are instead magnetized (since the energy is near the ground state energy).

The question in the present context is therefore whether a long-time configuration of (low) energy $`E_0`$ is well reproduced by the typical ‘blocked’ configuration of the same energy. By simulating the corresponding ‘auxiliary’ model, (with auxiliary energy equal to the number of spins not aligned with their local field, i.e. to the number of ’mobile’ spins), we have checked that this is not the case: the blocked configurations consituting Edwards’ distribution at energy $`E_0`$ are also magnetised. Therefore, neither Gibbs’ nor Edwards’ distributions describe the typical configurations obtained dynamically.

When is then the flatness assumption characterizing Edwards’ argument justified? A natural criterion, suggested by glass theory , consists in studying how a system explores its phase space, i.e. its ‘chaoticity’ properties. After aging for a time $`t_w`$, two copies (clones) are made of the system, and allowed to evolve subsequently with different realisations of the randomness in the updating procedure. We have checked that in the KA model the two clones always diverge (the slower the larger $`t_w`$, see figure 4), while for the 3D Random Field Ising model they do not. It is thus tempting to conjecture that this form of chaoticity is a necessary condition to have flat statistical weights for the blocked configurations. Note that for this criterion to make sense, it should always be applied at non-zero (though weak) tapping or shearing. The condition of chaoticity is however not sufficient: Bouchaud’s ‘trap model’  is chaotic but its fluctuation-dissipation properties are not directly related to the density of states.

To summarize, our study suggests that the proposal made by Edwards does indeed make sense and opens a door towards a statistical (thermodynamic) description of compact granular matter under very weak driving. In order to generalise these ideas to stronger forcing, lower chemical potential or higher temperatures (as required to analyse the experiments in ), one has to learn how to go from the concept of ‘blocked configuration’ to that of ‘metastable state’, and this requires other tools . The inherent structure construction could provide a practical shortcut.

The present setting of the problem immediately suggests experiments to check these ideas, e.g. by studying diffusion and mobility of tracer particles within driven granular media.

Finally, let us note that even in the simplest cases, the correspondence between Edwards’ distribution and long-time dynamics is at best checked but does not follow from any principle. The situation is thus as if one would have checked that the microcanonical distribution gives good results for gases, without knowing Liouville’s theorem that proves that such a distribution is indeed left invariant by the equations of motion. Such more refined arguments would be very welcome.

ACKNOWLEDGEMENTS: MS is supported by a Marie Curie fellowship of the European Commission (contract ERBFMBICT983561). This work has also been partially supported from the European Network-Fractals under contract No. FMRXCT980183.