text:
# Test of a theoretical equation of state for elemental solids and liquids

## 1 Introduction

Over approximately the last sixty years, numerous models and techniques have been developed for creating equations of state (EOS) for a variety of materials that are valid up to very extreme pressures (tens of Mbar) and temperatures (several eV). In the EOS community at the national laboratories, for instance, we have often used models based on the Mie-Grüneisen EOS together with the Thomas-Fermi or Thomas-Fermi-Dirac model (or one of its modifications) to include the contributions from the electrons (see  for examples). The models usually contain enough independent parameters to adjust the EOS until it correctly reproduces the experimentally measured Hugoniot (and perhaps a few other data points), but it is generally an open question how accurate the EOS is away from the Hugoniot. In this paper we argue that for one class of materials, elemental solids and liquids, our understanding of the underlying condensed matter Hamiltonian for the nuclei and electrons has grown to the point that we can construct highly accurate EOS from essentially first principles, and we also propose a means for doing so. We also argue that, since the underlying physics is well understood, an EOS derived this way should have the right functional form, even if we are unsure of the values of some of its parameters; thus, if the resulting EOS is shown to be accurate in one thermodynamic region (say, along the Hugoniot), then we can be confident that it is roughly equally accurate elsewhere.

In this formalism, the EOS in the solid phase depends on a decomposition of the Hamiltonian due to Wallace (see Chapter 1 of ), extending the work of Born  to metals as well as insulators; the resulting free energy contains terms describing the harmonic motion of the nuclei about their lattice sites (phonons), thermal excitation of the electrons from their ground state, anharmonic corrections to the nuclear motion (represented as phonon-phonon interactions), and interactions between the electron excitations and the nuclear motion, represented as electron-phonon interactions. (Please note that this description is exact; all of the physics contained in the true Hamiltonian of the system is included here. Specific EOS models usually neglect the anharmonic and electron-phonon terms, arguing that anharmonicity is small and making reference to some form of the Born-Oppenheimer approximation; we will take a somewhat different route, commenting on approximations below.) A recently developed theory of the dynamics of monatomic liquids (see  for a review) uses the same Hamiltonian to derive a liquid free energy which is quite similar to the expression for a solid, with additional terms accounting for the fact that the liquid, as opposed to the solid, traverses many potential valleys and thus sees the boundaries between them. For both phases, the resulting free energies have been compared with experimental data in the low-pressure regime ($`P100`$ kbar), with the following results (Sections 17-19 and 23 of ):

(a) Molecular dynamics (MD) calculations of the anharmonic contribution to the entropy of several solids match experimental entropy data to the accuracy of the data themselves.

(b) Low-temperature ($`T20`$ K) calculations of the electron-phonon term for several solids lead to predictions that also match experimental entropy to the accuracy of the data.

(c) Theoretical arguments show that the electron-phonon contribution is entirely negligible except when the electronic contribution dominates the free energy, such as in metallic solids at low temperatures.

(d) For the 27 elemental solids for which accurate data are available from low $`T`$ (but not too low; see point (c)) to the melting temperature $`T_m`$, the free energy excluding the anharmonic and electron-phonon terms accounts for the experimental thermal energy and entropy to an accuracy of $`5\%`$ (in fact, an accuracy of $`2\%`$ for all but about five materials).

(e) For the 6 elements in the liquid phase for which accurate data are available at temperatures up to around $`3T_m`$, the effect of neglecting the anharmonic, boundary, and electron-phonon contributions to the energy and entropy is similarly small.

This tells us that at low pressures, we can neglect the anharmonic, boundary, and electron-phonon terms in both the solid and liquid free energy (which happen to be the hardest terms to calculate), and the resulting thermal energy and entropy are both simple in form and accurate at the $`5\%`$ level. It is for this reason, not an appeal to the Born-Oppenheimer or other approximations, that we know we can simplify our EOS and what the results of the simplification will be, at least at low pressures.

In this paper we do two things: (1) We describe in more detail this framework for constructing EOS and discuss the theoretical and experimental inputs needed to implement it, and (2) we construct a sample EOS, neglecting anharmonic, boundary, and electron-phonon terms, both to illustrate the method and to discover whether points (d) and (e) above continue to hold in the high-pressure regime. We use Aluminum as our sample because of the availability of extensive electronic structure calculations, up to a compression of three, and highly accurate shock Hugoniot data, which provide a test of our EOS through both phases to pressures of around 5 Mbar. In Subsection 2.1 we develop the general theory of the solid EOS, and in Subsection 2.2 we do the same for the liquid. In Subsection 3.1, we construct our sample EOS for Al, comparing it with other EOS work, and in Subsection 3.2 we compute the Hugoniot predicted by the EOS and compare it with experimental data. The results are encouraging. Finally, we review our work, discuss the advantages and disadvantages of this formalism (and how to address the disadvantages), and suggest directions for future development.

## 2 General theory

### 2.1 Solid phase

The condensed matter Hamiltonian, decomposed as described above, consists of terms describing the motion of the nuclei in a potential generated by the electrons in their ground state, plus additional terms that lead to the thermal excitation of the electrons and describe their interactions with the nuclear motion. With this Hamiltonian, the Helmholtz free energy per atom for a solid at temperature $`T`$ with volume $`V`$ per atom takes the form

$$F^\mathrm{s}(V,T)=\mathrm{\Phi }_0^\mathrm{s}(V)+F_{\mathrm{ph}}^\mathrm{s}(V,T)+F_{\mathrm{el}}^\mathrm{s}(V,T)+F_{\mathrm{anh}}^\mathrm{s}(V,T)+F_{\mathrm{ep}}^\mathrm{s}(V,T).$$
(1)
Here $`\mathrm{\Phi }_0^\mathrm{s}`$ is the static lattice potential (the electronic ground state energy when the nuclei are fixed at their lattice sites); it depends on the particular crystal structure. $`F_{\mathrm{ph}}^\mathrm{s}`$ is the contribution from the harmonic motion of the nuclei about their lattice sites, $`F_{\mathrm{el}}^\mathrm{s}`$ represents the thermal excitation of the electrons when the nuclei are fixed at their lattice sites, $`F_{\mathrm{anh}}^\mathrm{s}`$ accounts for the anharmonicity of the nuclear motion (which may be represented as phonon-phonon interactions), and $`F_{\mathrm{ep}}^\mathrm{s}`$ expresses the interactions between the electron excitations and the nuclear motion, represented as electron-phonon interactions. (We emphasize again that this free energy is exact; it includes all of the physics present in the Hamiltonian.) The discussion in the Introduction justifies our approximating the solid free energy as

$$F^\mathrm{s}(V,T)=\mathrm{\Phi }_0^\mathrm{s}(V)+F_{\mathrm{ph}}^\mathrm{s}(V,T)+F_{\mathrm{el}}^\mathrm{s}(V,T),$$
(2)
so let us now consider the forms of $`F_{\mathrm{ph}}^\mathrm{s}`$ and $`F_{\mathrm{el}}^\mathrm{s}`$ and the parameters on which they depend. The phonon term in the Hamiltonian describes harmonic motion, which leads uniquely to the free energy of lattice dynamics:

$$F_{\mathrm{ph}}^\mathrm{s}(V,T)=_0^{\mathrm{}}g^\mathrm{s}(\omega )\left[\frac{1}{2}\mathrm{}\omega +\mathrm{ln}(1e^{\beta \mathrm{}\omega })\right]𝑑\omega ,$$
(3)
where $`\beta =1/kT`$ and $`g^\mathrm{s}(\omega )`$ is the distribution of phonon frequencies in the Brillouin zone. (Note that $`g^\mathrm{s}(\omega )`$ is volume dependent.) Sometimes we require not the full Eq. (3) but only its high- and low-temperature limits, for which we need not the full $`g^\mathrm{s}(\omega )`$ but only three of its moments, expressed in terms of the characteristic temperatures $`\mathrm{\Theta }_0^\mathrm{s}`$, $`\mathrm{\Theta }_1^\mathrm{s}`$, and $`\mathrm{\Theta }_2^\mathrm{s}`$ defined by

$`\mathrm{ln}k\mathrm{\Theta }_0^\mathrm{s}`$ $`=`$ $`\mathrm{ln}\mathrm{}\omega _{\mathrm{BZ}}`$
$`k\mathrm{\Theta }_1^\mathrm{s}`$ $`=`$ $`{\displaystyle \frac{4}{3}}\mathrm{}\omega _{\mathrm{BZ}}`$
$`k\mathrm{\Theta }_2^\mathrm{s}`$ $`=`$ $`\left[{\displaystyle \frac{5}{3}}(\mathrm{}\omega )^2_{\mathrm{BZ}}\right]^{1/2},`$ (4)
where $`\mathrm{}_{\mathrm{BZ}}`$ indicates an average over all the frequencies in the Brillouin zone. Then the following limits hold:

$$F_{\mathrm{ph}}^\mathrm{s}(V,T)\frac{9}{8}k\mathrm{\Theta }_1^\mathrm{s}\mathrm{as}T0$$
(5)
and

$$F_{\mathrm{ph}}^\mathrm{s}(V,T)=3kT\left[\mathrm{ln}\left(\frac{T}{\mathrm{\Theta }_0^\mathrm{s}}\right)\frac{1}{40}\left(\frac{\mathrm{\Theta }_2^\mathrm{s}}{T}\right)^2+\mathrm{}\right]\mathrm{at}\mathrm{high}T.$$
(6)
The leading term in Eq. (6) describes purely classical nuclear motion, while the series of terms in powers of $`T^2`$ are quantum corrections. Keeping only the first quantum correction, the thermodynamic functions derived from Eq. (6) are accurate to $`1\%`$ at temperatures above $`\frac{1}{2}\mathrm{\Theta }_2^\mathrm{s}`$.

The electronic excitation free energy $`F_{\mathrm{el}}^\mathrm{s}`$ can be expressed generally as an integral function of the electronic density of states per atom, $`n^\mathrm{s}(ϵ)`$, and the Fermi distribution

$$f(ϵ)=\frac{1}{e^{\beta (ϵ\mu )}+1},$$
(7)
where $`\beta `$ is still $`1/kT`$ and $`\mu `$ is the chemical potential. If each atom contributes $`Z`$ electrons to the valence bands (notice that $`Z`$ is not necessarily the atomic number), with the lowest valence energy set to zero, then $`\mu `$ is a function of $`T`$ determined by the normalization condition

$$_0^{\mathrm{}}n^\mathrm{s}(ϵ)f(ϵ)𝑑ϵ=Z.$$
(8)
The electronic free energy is then

$`F_{\mathrm{el}}^\mathrm{s}(V,T)=`$ (9)
$`\mu Z{\displaystyle _0^{ϵ_F}}ϵn^\mathrm{s}(ϵ)𝑑ϵkT{\displaystyle _0^{\mathrm{}}}n^\mathrm{s}(ϵ)\mathrm{ln}[1+e^{\beta (ϵ\mu )}]𝑑ϵ,`$
where $`ϵ_F`$, the Fermi energy, is the value of $`\mu `$ when $`T=0`$. The second term on the right hand side of Eq. (9) is the subtraction of the electronic ground state energy, which ensures that $`F_{\mathrm{el}}^\mathrm{s}0`$ as $`T0`$. This property makes sense if $`F_{\mathrm{el}}^\mathrm{s}`$ represents purely thermal excitation of the electrons. (It also avoids double counting of the energy, as the electronic ground state energy is already represented as $`\mathrm{\Phi }_0^\mathrm{s}`$.)

We see from this discussion that to evaluate the terms in Eq. (2) for the solid free energy we require three unknown functions: $`\mathrm{\Phi }_0^\mathrm{s}`$, $`g^\mathrm{s}(\omega )`$ (or $`\mathrm{\Theta }_0^\mathrm{s}`$, $`\mathrm{\Theta }_1^\mathrm{s}`$, and $`\mathrm{\Theta }_2^\mathrm{s}`$ if we are concerned only with the high- and low-$`T`$ limits), and $`n^\mathrm{s}(ϵ)`$ (and the associated quantities $`Z`$ and $`ϵ_F`$). These can be determined in various ways: compressibility data and diamond anvil cell data can be used to construct $`\mathrm{\Phi }_0^\mathrm{s}(V)`$; neutron scattering experiments can determine $`g^\mathrm{s}(\omega )`$ or its various moments at $`P=1`$ bar; and for many elements all three of these functions can be computed reliably using electronic structure theory. (Or one could use results from multiple sources in combination, which is often an option with $`\mathrm{\Phi }_0^\mathrm{s}`$ and is basically a necessity with $`g^\mathrm{s}(\omega )`$.) One must keep in mind, however, that the accuracy of one’s answers will be limited by the accuracy and range of applicability of these functions, regardless of how they are determined.

### 2.2 Liquid phase and two-phase region

According to the theory of liquid dynamics reviewed in , the same Hamiltonian that gave us the solid free energy leads to a similar form for the free energy of a monatomic liquid. In this theory, the region of the many-body potential surface in which the system moves in the liquid phase is dominated by a large number of intersecting nearly-harmonic valleys, called “random” valleys because they correspond to particle configurations which retain no remnant crystal symmetry, and which are all macroscopically identical. In particular, the valleys all have the same distribution of normal mode frequencies, and they all have the same depth (which, as in the solid case, is the electronic ground state energy when the nuclei are fixed at the valley minimum). The resulting liquid free energy per atom is

$`F^l(V,T)`$ $`=`$ $`\mathrm{\Phi }_0^l(V)+F_{\mathrm{ph}}^l(V,T)+F_{\mathrm{el}}^l(V,T)+F_{\mathrm{ab}}^l(V,T)+`$ (10)
$`F_{\mathrm{ep}}^l(V,T)kT\mathrm{ln}w.`$
All of the terms correspond to their solid counterparts with the following exceptions:

(1) $`\mathrm{\Phi }_0^l`$, now called the static structure potential, is the depth of a typical valley in which the liquid system moves.

(2) The normal mode spectrum appearing in $`F_{\mathrm{ph}}^l`$ is that of a typical liquid potential valley, not the unique solid potential valley.

(3) The term $`F_{\mathrm{ab}}^l`$ includes corrections due both to anharmonicity and to the fact that the potential valleys have boundaries, which the liquid (as opposed to the solid) encounters as it transits from valley to valley.

(4) The extra term $`kT\mathrm{ln}w`$ corresponds to an increase in entropy of $`k\mathrm{ln}w`$ per atom; the value $`\mathrm{ln}w0.8`$ is estimated from entropy of melting data of the elements (again, see  for details). In liquid dynamics theory, this term is due to the hypothesis that the number of potential valleys among which the liquid moves is of order $`w^N`$, where $`N`$ is the number of atoms in the system.

We emphasize that the same Hamiltonian gives rise to both Eq. (1) and (10); the differences are that the potential is expanded about different equilibrium configurations in the two cases, and that the region of configuration space over which the liquid moves is obviously far larger than the space available to the solid (hence the $`kT\mathrm{ln}w`$ term).

Again making the approximations discussed in the Introduction, our form for the liquid free energy becomes

$$F^l(V,T)=\mathrm{\Phi }_0^l(V)+F_{\mathrm{ph}}^l(V,T)+F_{\mathrm{el}}^l(V,T)kT\mathrm{ln}w,$$
(11)
and the additional term $`kT\mathrm{ln}w`$ is fully determined by setting $`\mathrm{ln}w=0.8`$, as mentioned above. The form of the phonon term is dictated by a central hypothesis of liquid dynamics theory: The motion of the liquid consists of oscillations in the macroscopically similar valleys described above together with occasional transits between valleys; the transits are of such short duration that they do not contribute to the thermodynamics to lowest order. Thus we will take $`F_{\mathrm{ph}}^l`$ to have the same form as the solid phonon term, Eq. (3), with a possibly different phonon frequency distribution $`g^l(\omega )`$. The electronic excitation term for the solid was derived using only the assumption that the electrons are thermally distributed over the available states using Fermi statistics; all of the information about the configuration of the nuclei is contained in the density of states. Hence $`F_{\mathrm{el}}^l`$ also takes the same form as the corresponding solid term, Eq. (9), with a density of states $`n^l(ϵ)`$ appropriate for the liquid phase. (What this means is discussed briefly below.)

The liquid and solid EOS together determine the melting temperature as a function of pressure $`T_m(P)`$ by the requirement that the solid and liquid Gibbs free energies match along the melt curve, or

$$G^\mathrm{s}(P,T_m(P))=G^l(P,T_m(P)).$$
(12)
Once the melt curve is determined, one can use the solid and liquid EOS separately to compute $`V_m^\mathrm{s}(T)`$ and $`V_m^l(T)`$, the atomic volumes of the solid and liquid at melt as functions of temperature. (Of course, using $`T_m(P)`$ one can express $`V_m^\mathrm{s}`$ and $`V_m^l`$ as functions of pressure equally well.)

In the case $`V_m^\mathrm{s}(T)<V_m^l(T)`$ for all $`T`$, which we assume here, the computation of the full two-phase EOS ($`F`$, $`E`$, $`S`$, and $`P`$) for any $`V`$ and $`T`$ proceeds as follows. If $`VV_m^\mathrm{s}(T)`$ for the given $`T`$, the system is in the solid phase; $`F^\mathrm{s}`$ is computed as described in Subsection 2.1, and the other functions follow from

$`S^\mathrm{s}`$ $`=`$ $`\left({\displaystyle \frac{F^\mathrm{s}}{T}}\right)_V,`$
$`E^\mathrm{s}`$ $`=`$ $`F^\mathrm{s}+TS^\mathrm{s},`$
$`P^\mathrm{s}`$ $`=`$ $`\left({\displaystyle \frac{F^\mathrm{s}}{V}}\right)_T.`$ (13)
If $`VV_m^l(T)`$, the system is in the liquid phase; $`F^\mathrm{l}`$ is computed as described in this Subsection, and the other functions follow from expressions analogous to Eqs. (13). If $`V_m^\mathrm{s}(T)<V<V_m^l(T)`$, the system is in the two-phase region; defining $`\eta `$, the volume fraction of the system in the liquid phase, by

$$\eta =\frac{VV_m^\mathrm{s}(T)}{V_m^l(T)V_m^\mathrm{s}(T)},$$
(14)
the thermodynamic functions are

$`F(V,T)`$ $`=`$ $`\eta F^l(V_m^l(T),T)+(1\eta )F^\mathrm{s}(V_m^\mathrm{s}(T),T),`$
$`E(V,T)`$ $`=`$ $`\eta E^l(V_m^l(T),T)+(1\eta )E^\mathrm{s}(V_m^\mathrm{s}(T),T),`$
$`S(V,T)`$ $`=`$ $`\eta S^l(V_m^l(T),T)+(1\eta )S^\mathrm{s}(V_m^\mathrm{s}(T),T),`$
$`P(V,T)`$ $`=`$ $`{\displaystyle \frac{F^l(V_m^l(T),T)F^\mathrm{s}(V_m^\mathrm{s}(T),T)}{V_m^l(T)V_m^\mathrm{s}(T)}}.`$ (15)
Finally, we note that just as with the solid, to evaluate the terms in Eq. (11) we require three unknown functions: $`\mathrm{\Phi }_0^l`$, $`g^l(\omega )`$ (or $`\mathrm{\Theta }_0^l`$ and $`\mathrm{\Theta }_2^l`$, since we’re not likely to need the low-$`T`$ limit), and $`n^l(ϵ)`$. In contrast to the solid case, these functions are generally not available experimentally. (It is possible that one might be able to compute $`\mathrm{\Phi }_0^l`$ using liquid compressibility data, but we suspect that this will be very difficult.) However, for many materials these functions should be computable using electronic structure theory, proceeding much as one would in the solid case except that the nuclei would be arranged not in a crystal configuration but in a disordered structure characteristic of a “random” valley in the liquid potential surface . To our knowledge very few such calculations have been attempted; the only ones we are aware of are $`\mathrm{\Phi }_0^l`$ and $`g^l(\omega )`$ at a single volume for liquid sodium in  (the results are referred to in  and a graph of $`g^l(\omega )`$ using their results appears as Fig. 1 in ). Another function that is sometimes available is the melt curve $`T_m(P)`$, but this curve cannot be chosen independently of the others, since the solid and liquid EOS determine it jointly; but this can be an advantage, since if $`T_m(P)`$ is known from experiment, for example, it can be used to compute one of the other needed functions if it is not otherwise available. In fact, this is how we will determine $`\mathrm{\Phi }_0^l`$ in our example EOS, to which we now turn.

## 3 An example: Aluminum

To illustrate the application of the theory we’ve described, we will now construct an EOS for Aluminum, which has been the subject of extensive electronic structure calculations and for which a great deal of high-quality experimental data are available. We will then compare the Hugoniot predicted by our EOS with data up to pressures of approximately 5 Mbar; this will tell us whether the approximations we discussed in the Introduction (neglecting anharmonic, boundary, and electron-phonon effects), known to be accurate at low pressures, continue to be reasonable in the high-pressure domain.

### 3.1 Constructing the EOS

We recall from Subsection 2.1 that the solid EOS requires three functions: $`\mathrm{\Phi }_0^\mathrm{s}`$, $`g^\mathrm{s}(\omega )`$, and $`n^\mathrm{s}(ϵ)`$. Since we will be testing the EOS by comparison with Hugoniot data, we will always be in the high-$`T`$ region (except for one brief low-$`T`$ excursion; see below), so we use Eq. (6) for $`F_{\mathrm{ph}}^\mathrm{s}`$ instead of Eq. (3); this means that we require only $`\mathrm{\Theta }_0^\mathrm{s}`$, $`\mathrm{\Theta }_1^\mathrm{s}`$, and $`\mathrm{\Theta }_2^\mathrm{s}`$ in place of $`g^\mathrm{s}(\omega )`$.

To determine these functions, we began by consulting the results of density functional theory (DFT) calculations carried out in the local density approximation by Straub et al. . They worked with fcc and bcc Al at atomic volumes from $`37a_0^3`$ to $`160a_0^3`$, where $`a_0`$ is the Bohr radius, corresponding to densities from $`8.17`$ g/cm<sup>3</sup> to $`1.89`$ g/cm<sup>3</sup> (the density of Al at 293 K and 1 bar is $`2.700`$ g/cm<sup>3</sup>). Their calculations indicate a $`T=0`$ transition from fcc to bcc at $`51a_0^3`$, corresponding to $`\rho =5.93`$ g/cm<sup>3</sup>, but we will neglect this phase change and treat solid Al as an fcc crystal for two reasons: The DFT calculations themselves suggest that the effect of the phase change on the thermodynamic functions will be quite small; and we know from experiment that the solid-liquid transition on the Hugoniot takes place well before reaching the density of concern (see Subsection 3.2), so we are confident of our assumption of fcc along the Hugoniot until melting. However, this assumption may have an effect on the liquid EOS at high densities, which we will comment on below. (Other electronic structure work, discussed on pp. 89-90 of , suggests the possibility of an hcp phase between the fcc and bcc phases, but as  also mentions, no experimental signature of this phase has been found, so we will proceed under the assumption of a single solid phase.)

Straub et al. computed $`\mathrm{\Phi }_0^\mathrm{s}`$ for fcc by fitting their results to a Birch-Murnaghan form,

$$\mathrm{\Phi }_0^\mathrm{s}(V)=c_0+V_b\underset{n=2}{\overset{5}{}}\frac{c_n}{n!}\left\{\frac{1}{2}\left[\left(\frac{V}{V_b}\right)^{2/3}1\right]\right\}^n,$$
(16)
with coefficients

$`V_b=106.302a_0^3,c_0=287.7832\mathrm{mRy},`$
$`c_2=761.2029\mathrm{GPa},c_3=1319.036\mathrm{GPa},`$
$`c_4=\mathrm{13\hspace{0.17em}661.06}\mathrm{GPa},c_5=\mathrm{50\hspace{0.17em}315.53}\mathrm{GPa}.`$ (17)
The $`\mathrm{\Theta }_n^\mathrm{s}`$ were determined by computing the bulk modulus and four zone-boundary phonons at several volumes, and these results were used to calibrate a pseudopotential model at each volume. The pseudopotential was then used to calculate phonon frequencies throughout the Brillouin zone, allowing the determination of $`\mathrm{\Theta }_0^\mathrm{s}`$, $`\mathrm{\Theta }_1^\mathrm{s}`$, and $`\mathrm{\Theta }_2^\mathrm{s}`$. Their results are shown in Table 1 and Figure 1. (The full set of results was not reported in .)

To check these results, Straub et al. compared experimental phonon moments for Al at $`T=80`$ K and $`P=1`$ bar based on Born-von Karmen fits to neutron scattering data  with their predictions interpolated to the appropriate atomic volume of $`110.7a_0^3`$. The experimental points, also shown in Figure 1, are in very good agreement with their calculations; hence these results for the $`\mathrm{\Theta }_n^\mathrm{s}`$ are acceptable for use in our EOS without modification. To determine the $`\mathrm{\Theta }_n^\mathrm{s}`$ at any volume, we first constructed a functional fit to the $`\mathrm{\Theta }_0^\mathrm{s}`$ points, with the result

$$\mathrm{\Theta }_0^\mathrm{s}(V)=2852.69+\frac{\mathrm{17\hspace{0.17em}319.9}}{V}+2.33667V633.858\mathrm{ln}(V),$$
(18)
where $`\mathrm{\Theta }_0^\mathrm{s}`$ is in K and $`V`$ is in $`a_0^3`$. Then we noted that according to the DFT results both $`\mathrm{\Theta }_1^\mathrm{s}`$ and $`\mathrm{\Theta }_2^\mathrm{s}`$ approximately equal $`e^{1/3}\mathrm{\Theta }_0^\mathrm{s}`$, so we made the same approximation using Eq. (18) for $`\mathrm{\Theta }_0^\mathrm{s}`$ to determine $`\mathrm{\Theta }_1^\mathrm{s}`$ and $`\mathrm{\Theta }_2^\mathrm{s}`$ at any volume. These functions are also shown in Figure 1.

The DFT calculations also provided data on the electronic density of states $`n^\mathrm{s}(ϵ)`$. Graphs of $`n^\mathrm{s}(ϵ)`$ for fcc and bcc Al at atomic volume $`112.0a_0^3`$ (corresponding to $`P=0`$ and $`T=295`$ K) are shown in Figure 2, along with the free electron $`n^\mathrm{s}(ϵ)`$, for which

$$n^\mathrm{s}(ϵ)=\sqrt{\frac{ϵ}{ϵ_F}}\left(\frac{3Z}{2ϵ_F}\right)\mathrm{and}ϵ_F=\frac{\mathrm{}^2}{2m_e}\left(\frac{3\pi ^2Z}{V}\right)^{2/3},$$
(19)
at $`V=112.0a_0^3`$ and $`Z=3`$.

The Figure shows that the free electron model is a good approximation for either crystal structure, for electronic excitations to around $`\frac{1}{2}`$ Ry. The same match, at a slightly poorer level of approximation and for excitations to around $`1`$ Ry, is found at our smallest atomic volume of $`37a_0^3`$. For all volumes of our study and temperatures up to $`T_m`$, the total electronic excitation contribution to the energy, entropy, and pressure is at most $`5\%`$, so the error introduced by using the free electron $`n^\mathrm{s}(ϵ)`$ in our calculations is negligible. Making this approximation, the normalization condition from Eq. (8) becomes

$$F_{1/2}(\beta \mu )=\frac{2}{3}(\beta ϵ_F)^{3/2},$$
(20)
which determines $`\mu `$ as a function of $`V`$ and $`T`$, and Eq. (9) for the free energy becomes

$$F_{\mathrm{el}}^\mathrm{s}(V,T)=Z\left(\mu \frac{2}{3}kT\frac{F_{3/2}(\beta \mu )}{F_{1/2}(\beta \mu )}\frac{3}{5}ϵ_F\right)$$
(21)
where $`Z=3`$. The $`F_n(x)`$ are the standard Fermi integrals; their properties are discussed on pp. 332-334 of  and their values for various $`n`$ are tabulated in .

The solid EOS that results from assembling all of these functions is reliable over a large range of volumes and temperatures; however, it is not in perfect agreement with the highly accurate experimental data that are available at low pressures. Specifically, experiments on Al at $`T=0`$ and $`P=0`$ show that

$`V_0=110.6a_0^3,E_0=249\mathrm{mRy},`$
$`B_0=79.4\mathrm{GPa},{\displaystyle \frac{dB_0}{dP}}=4.7,`$ (22)
but the EOS yields $`V_0=107.3a_0^3`$, which is outside the experimental error. Therefore, we chose to make a small correction to our purely theoretical free energy to agree with experiment. These quantities are obviously determined by the $`T=0`$ form of the free energy, $`F_0^\mathrm{s}=\mathrm{\Phi }_0^\mathrm{s}+\frac{9}{8}k\mathrm{\Theta }_1^\mathrm{s}`$ (see Eq. (5)); since $`\mathrm{\Theta }_1^\mathrm{s}`$ is already in excellent agreement with available experiment, we chose to modify only $`\mathrm{\Phi }_0^\mathrm{s}`$. To proceed, we noted that the data determine $`P_0`$, the $`T=0`$ pressure, in the vicinity of $`V=V_0`$ by the relation

$`P_0(V)`$ $`=`$ $`P_0(V_0)+{\displaystyle \frac{dP_0}{dV}}|_{V_0}(VV_0)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^2P_0}{dV^2}}|_{V_0}(VV_0)^2+\mathrm{}`$ (23)
$`=`$ $`{\displaystyle \frac{B_0}{V_0}}(VV_0)+{\displaystyle \frac{B_0}{2V_0^2}}\left(1+{\displaystyle \frac{dB_0}{dP}}\right)(VV_0)^2+\mathrm{}`$
while at higher compressions we have no information to supplement the electronic structure result; so we decided to construct a $`\mathrm{\Phi }_0^\mathrm{s}`$ that correctly reproduces Eq. (23) near $`V_0`$ but smoothly interpolates to Eq. (16) at lower volumes. To do this, we computed $`P_0`$ at 10 volumes between $`110a_0^3`$ and $`111.25a_0^3`$ using Eq. (23), and we also computed $`P_0=F_0^\mathrm{s}/V`$ using the above form for $`F_0^\mathrm{s}`$, with $`\mathrm{\Phi }_0^\mathrm{s}`$ from Eq. (16), at 23 volumes between $`30a_0^3`$ and $`41a_0^3`$. We then performed a least-squares fit to these points using an expression similar to the Birch-Murnaghan form, but carried to a slightly higher order; after integrating, adjusting the constant of integration to correctly match $`E_0`$ from Eq. (22), and subtracting off $`\frac{9}{8}k\mathrm{\Theta }_1^\mathrm{s}`$, we had a new $`\mathrm{\Phi }_0^\mathrm{s}`$ given by

$`\mathrm{\Phi }_0^\mathrm{s}(V)`$ $`=`$ $`1.64615\times 10^6+{\displaystyle \frac{2.07608\times 10^7}{V^{2/3}}}{\displaystyle \frac{4.61515\times 10^8}{V^{4/3}}}`$ (24)
$`+{\displaystyle \frac{5.71249\times 10^9}{V^2}}{\displaystyle \frac{5.49998\times 10^{10}}{V^{8/3}}}+{\displaystyle \frac{3.71978\times 10^{11}}{V^{10/3}}}`$
$`{\displaystyle \frac{1.66284\times 10^{12}}{V^4}}+{\displaystyle \frac{4.41118\times 10^{12}}{V^{14/3}}}{\displaystyle \frac{5.25064\times 10^{12}}{V^{16/3}}}`$
$`+{\displaystyle \frac{2.26789\times 10^7}{V}}220.716V+\mathrm{23\hspace{0.17em}6788}\mathrm{ln}(V)`$
where $`\mathrm{\Phi }_0^\mathrm{s}`$ is in mRy and $`V`$ is in $`a_0^3`$. This $`\mathrm{\Phi }_0^\mathrm{s}`$, which reproduces the data in Eq. (22) and interpolates smoothly to the DFT curve at higher compressions, is what we use in our EOS instead of Eq. (16). The $`T=0`$ pressure-volume curves constructed using both the original and new $`\mathrm{\Phi }_0^\mathrm{s}`$ are shown in Figure 3.

Our choice of a Birch-Murnaghan-like form was dictated by the fact that the Straub et al. result provides most of our information about the shape of $`\mathrm{\Phi }_0^\mathrm{s}`$; so our goal was to preserve that form insofar as was possible, interpolating back to their result as quickly as we could without introducing enough curvature to compromise agreement with $`dB_0/dP`$. This correction to $`\mathrm{\Phi }_0^\mathrm{s}`$ constitutes a small change to the overall EOS; the effect of this change on the Hugoniot will be considered in the next Section. This modification completes the full solid free energy, so we can now consider the liquid.

For the liquid we need the same three functions that we needed for the solid, and we must also consider the melting curve $`T_m(P)`$. Having chosen to use Eq. (6) for $`F_{\mathrm{ph}}^\mathrm{s}`$, we certainly did the same for $`F_{\mathrm{ph}}^l`$, since the Hugoniot will obviously enter the liquid region only at rather high temperatures; thus we needed only $`\mathrm{\Theta }_0^l`$ and $`\mathrm{\Theta }_2^l`$. From experiment we know that Al is what is called in liquid dynamics a “normal melting element” (the entropy of melting at constant density is approximately $`0.8k`$ per atom), and we argue in  that $`\mathrm{\Theta }_0`$ in solid and liquid phases of such an element are approximately equal. (Experimental and computational work supporting this conjecture are also discussed in .) Thus we took the liquid to have the same $`\mathrm{\Theta }_0`$ as the solid. It is also true that in the liquid, $`T`$ is typically much larger than $`\mathrm{\Theta }_2^\mathrm{s}`$ (for example, in liquid Al at normal density $`T2\mathrm{\Theta }_2^\mathrm{s}`$), rendering the first quantum correction to $`F_{\mathrm{ph}}^l`$ negligible (roughly $`1\%`$ at normal density), so even if $`\mathrm{\Theta }_2^l`$ differs from $`\mathrm{\Theta }_2^\mathrm{s}`$ by $`25\%`$ or more, the impact on the phonon term will be very small; therefore we also used the same $`\mathrm{\Theta }_2`$ in the liquid as in the solid.

Since the free electron model approximates the DFT result for $`n^\mathrm{s}(ϵ)`$ so well for both fcc and bcc Al (Figure 2), which correspond to two valleys in the many-body potential surface with rather different structures, we also expect this model to be a good approximation for $`n^l(ϵ)`$, the density of states for the structure characteristic of a liquid. Since at all volumes and temperatures up to $`5T_m`$ (the relevance of this number will appear below), the electronic contribution to the thermodynamic functions does not exceed $`25\%`$, the error introduced by the free electron model is still acceptable.

We fixed the only remaining term in Eq. (11), $`\mathrm{\Phi }_0^l`$, by the requirement that the Gibbs free energies of the solid and liquid match along the Al melting curve. Boehler and Ross  suggested that

$$T_m(P)=933.45\mathrm{K}\left(\frac{P}{6.049\mathrm{GPa}}+1\right)^{0.531}$$
(25)
on the basis of their experimental work up to 80 GPa (0.8 Mbar), and experiments by McQueen et al.  and Hänström and Lazor  and theoretical work by Pélissier  suggest that this curve continues to be accurate up to 200 GPa. In the absence of evidence to the contrary, we took Eq. (25) to be valid to higher pressures as needed. (As we will see later, our EOS will assume Eq. (25) no higher than 400 GPa.) We computed $`\mathrm{\Phi }_0^l`$ as follows: We made a guess for $`\mathrm{\Phi }_0^l`$ not very different from $`\mathrm{\Phi }_0^\mathrm{s}`$, and then we used it and Eq. (25) to calculate the difference between the two Gibbs free energies,

$$\mathrm{\Delta }G(P)=G^\mathrm{s}(P,T_m(P))G^l(P,T_m(P)),$$
(26)
at several hundred values of $`P`$ over the entire pressure range considered in this study. We also calculated the liquid melt volume $`V_m^l(P)`$ at each $`P`$. If the rms average of Eq. (26) over all calculated points was not sufficiently small, we used the following easily verified fact: To first order, a small change $`\delta \mathrm{\Phi }_0^l`$ produces a small change $`\delta G^l(P,T_m(P))`$ given by $`\delta G^l(P,T_m(P))=\delta \mathrm{\Phi }_0(V_m^l(P))`$. Thus we performed the substitution

$$\mathrm{\Phi }_0^l(V)\mathrm{\Phi }_0^l(V)+\mathrm{\Delta }G(P_m^l(V)),$$
(27)
where $`\mathrm{\Delta }G`$ was computed by Eq. (26) and $`P_m^l(V)`$ is the inverse of $`V_m^l(P)`$, and calculated Eq. (26) again. We iterated until the rms average was sufficiently small (less than $`0.001`$ mRy in our case), giving us the needed $`\mathrm{\Phi }_0^l`$, which is shown in Figure 4 along with $`\mathrm{\Phi }_0^\mathrm{s}`$. We recorded $`\mathrm{\Phi }_0^l`$ as a list of points, and we did not create a functional fit for it; instead we used an interpolator to calculate it and its derivative when needed.

It is at this point that the existence of other solid phases in Al, discussed earlier, affects the EOS of the liquid. It is likely that the liquid region borders the fcc crystal only over part of its boundary, beyond which the liquid borders the bcc region or other solid phases. In such a case, at sufficiently high pressures we should use the free energy appropriate for that solid phase, not the fcc free energy, in Eq. (26). This suggests that $`\mathrm{\Phi }_0^l`$ may become inaccurate beyond densities in the neighborhood of 6 g/cm<sup>3</sup>, where the $`T=0`$ fcc-bcc phase transition occurs. We will take this fact into consideration when we discuss the limits of applicability of the EOS below.

Once we had the full solid and liquid EOS, we then solved Eq. (12) directly to compute $`T_m(P)`$, verifying that we had reproduced the Boehler-Ross curve; our result is shown in Figure 5, together with the data from  and some points from Pélissier’s theoretical curve. (According to , their data point at 125 GPa marks the onset of melting along the Hugoniot. We will comment on this in the next Subsection.)

Now that we have the full two-phase EOS, it is profitable to compare our work with another EOS for Al, due to Moriarty et al. . These authors also use a full lattice dynamics treatment of the crystal phonons, although they calculate their $`g^\mathrm{s}(\omega )`$ two separate ways, using both Moriarty’s Generalized Pseudopotential Theory (GPT) and a local pseudopotential model with parameters chosen to match solid-phase EOS data. We strongly prefer to rely on DFT results, as we believe DFT has reached such a stage of maturity that it more accurately captures the physics contained in the real Hamiltonian of the system, which as we have emphasized we believe to be understood in sufficient depth that it should underlie all of our work. Second, in their treatment of the liquid phase Moriarty et al. rely on fluid variational theory, described in detail in , to compute the least upper bound to the “real” liquid free energy (from a liquid Hamiltonian based on pseudopotentials) that can be obtained from the free energy of a reference system; Moriarty et al. investigate hard-sphere, soft-sphere, and one-component plasma reference systems before settling on the soft-sphere system as providing the best bound. Based on the investigations summarized in , we claim that we have the actual Hamiltonian of the liquid itself, not a Hamiltonian based on pseudopotential theory; furthermore, this Hamiltonian decomposes naturally into a dominant term that produces a free energy that can be used directly (instead of requiring approximation by the free energy of a reference system) and remaining terms whose contributions to the free energy are known to be small (see the Introduction). The same point we made above for the solid phase applies; we argue that it is a better strategy in developing EOS to try to understand the true Hamiltonian of the system, and then to use it when doing statistical mechanics. Almost inevitably, one must make approximations (which we certainly have done here), but we believe we are in a better position to understand and improve upon them if the physical foundation of the EOS is as solid as we can make it.

Finally, let us make some conservative estimates of the range of applicability of this EOS. Any limits will arise from two sources: the validity of the approximation that $`F_{\mathrm{ab}}^l`$ is negligible in the liquid (see the Introduction), and the limited ranges over which the functions $`\mathrm{\Phi }_0(V)`$, $`g(\omega )`$, $`n(ϵ)`$, and $`T_m(P)`$ are known. Let’s consider each in turn.

(1) We know from experiment that $`F_{\mathrm{ab}}^l`$ is negligible when $`T3T_m`$ (again, see the Introduction), and judging from trends in the data we suspect $`F_{\mathrm{ab}}^l`$ will still be small up to $`T5T_m`$, but clearly this term must become relevant as the nuclear motion becomes more gaslike. Thus we shall take care with any data at $`\rho `$ and $`T`$ such that $`T`$ approaches or exceeds $`5T_m(\rho )`$.

(2) At densities below approximately 6 g/cm<sup>3</sup>, we are confident that the solid is in the fcc phase, and the liquid free energy is based on this phase, so we trust the full EOS here. At higher densities, we must be more circumspect; the solid may have undergone a phase transition to bcc, and the liquid EOS at this density may be based on the wrong solid free energy. Further, as we have indicated earlier, Eq. (25) for the melt curve has received independent support only up to 200 GPa, so we must be cautious with the liquid EOS in regions beyond this point. We decided to be brave and accept the melt curve as valid up to 400 GPa; this corresponds to a liquid density of $`6.15`$ g/cm<sup>3</sup>, and since this is not far from the probable location of the solid fcc-bcc transition, we take it as the density limit of our EOS. (Even if we did not have this concern, we would be restricted to densities below $`8.17`$ g/cm<sup>3</sup>, where electronic structure results are available.) Also, the free electron approximation to $`n^\mathrm{s}(ϵ)`$ has been validated only for $`ϵϵ_F`$ up to 1/2 Ry, or 6.8 eV, at low compression and 1 Ry, or 13.6 eV, at high compression, but at higher temperatures the electronic energy and entropy are sensitive to the details of $`n^\mathrm{s}(ϵ)`$ to energies above these limits. We estimated the values of $`T`$ that begin to probe the unvalidated region of $`n^\mathrm{s}(ϵ)`$ (roughly $`3kT=ϵϵ_F`$), and we found that over our valid density range the $`T=5T_m`$ limit always took precedence. Hence this limit is not relevant for us, but we mention it for completeness, as it may become a concern if the EOS is extended to higher densities.

Figure 6 shows the limits $`\rho 6.15`$ g/cm<sup>3</sup> and $`T5T_m`$ of the EOS in $`TP`$ space, together with the melt curve and the Hugoniot (see the next Subsection), while Figure 7 shows the same three features in $`T\rho `$ space.

In this Figure, the melt curve becomes a two-phase region, which we will consider in more detail in the next Subsection.

We are confident that this EOS is valid within these limits, but we don’t know how far beyond them the inaccuracies begin to appear; thus we will not be shy about considering data not too far outside this range.

### 3.2 Comparison with Hugoniot data

If a shock wave travels at speed $`u_s`$ through a sample of material, accelerating its particles from rest to speed $`u_p`$ and changing its density, atomic volume, pressure, and internal energy per atom from $`\rho _0`$, $`V_0`$, $`P_0`$, and $`E_0`$ to $`\rho `$, $`V`$, $`P`$, and $`E`$, then (assuming thermal equilibrium before and after the shock) these quantities must satisfy the Rankine-Hugoniot relations,

$`\rho (u_su_p)`$ $`=`$ $`\rho _0u_s`$
$`PP_0`$ $`=`$ $`\rho _0u_su_p`$
$`EE_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}(P_0+P)(V_0V),`$ (28)
derived from considerations of mass, momentum, and energy conservation. (It is assumed that the wave is steady and strength effects are negligible.)

By solving these equations together with the EOS, which relates $`P`$, $`V`$, and $`E`$, we can determine the Hugoniot, the curve of all possible end states of the shocked material.

We used our EOS and Eqs. (28) to compute $`u_s`$ as a function of $`u_p`$ and $`P`$ as a function of $`\rho `$ along the Al Hugoniot; the results are shown in Figures 8 and 9 along with the intersection of the Hugoniot with the limit of validity of the EOS.

Hugoniot data from several sources  are also included.

The low-pressure region of the Hugoniot is highlighted in Figures 10 and 11, and the intermediate-pressure region, including the intersections with the phase boundaries, is shown in Figures 12 and 13.

Three important considerations in selecting which data to include are (1) the initial densities of the samples, (2) the quality of the experimental technique, and (3) whether the measurements were absolute or relative.

All of the available data were taken using Al alloys with densities that differ from the known pure metal value of $`2.70`$ g/cm<sup>3</sup> (predicted correctly by our EOS); some alloys are as close as $`2.71`$ g/cm<sup>3</sup> while others differ much more.

Since Hugoniots in general are quite sensitive to the initial density, we chose to compare only with the data for which $`\rho _0`$ clustered around $`2.71`$ g/cm<sup>3</sup>. (Thus we used only one data point from , which mainly concerns porous materials. All of the data from the other references were used.) We also avoided sources which gathered data using unusual shock wave geometries (such as ), and we also chose not to use the results of indirect or relative measurements, such as , preferring to rely on the absolute measurements that are available. Finally, we did not use the few data points available (primarily nuclear-driven) that lie very far beyond the limits of applicability of our EOS (but see below).

The theoretical Hugoniot compares well with both the $`u_s`$-$`u_p`$ and the $`P`$-$`\rho `$ data all the way up to the predicted limit of its validity, at approximately $`500`$ GPa ($`5`$ Mbar). More specifically, theory agrees with experiment at $`P\stackrel{<}{}40`$ GPa (Figures 10 and 11); at 40-125 GPa, theory falls below the experimental error bars by around 1% at most (Figures 10 through 13); and theory again lies within the experimental error bars through the liquid phase (Figures 8, 9, 12, and 13). (We recall that given percentage errors in $`u_s`$ and $`u_p`$ correspond to roughly the same percentage errors in the $`P`$-$`\rho `$ plane.) The presence of theoretical error only in the solid phase is likely due to strength effects, which are present in the solid but not in the liquid, and which are neglected in our Hugoniot calculations. Furthermore, as Figure 13 shows, we predict that the Hugoniot crosses the two-phase region between $`\rho =4.43`$ g/cm<sup>3</sup> and $`\rho =4.58`$ g/cm<sup>3</sup>, corresponding to a range in $`P`$ from 126 to 156 GPa; this agrees very well with , in which melting was found to occur between 125 and 150 GPa. (We note that  used Al 2024, an alloy whose density is sufficiently different from pure Al that we did not use Hugoniot data taken with that alloy in our Figures. We consider their melting results because, as we saw in the previous Subsection, their data are consistent with other experiments that did use pure Al.) The correction to $`\mathrm{\Phi }_0^\mathrm{s}`$ from the previous Subsection, shown in Figure 3, shifts the Hugoniot at pressures below 30 GPa, bringing it into excellent agreement with experiment, while at pressures above 60 GPa or so, the effect on the Hugoniot is insignificant.

We have also compared our results with data just beyond the EOS limits of validity; the points in  that match our initial density (one of which is a reanalysis of the single point in ), lying at about 10 Mbar, fall noticeably below our Hugoniot, and their consistency with the very-high-pressure points of Ragan  strongly suggests that they represent the true Hugoniot, which thus falls beneath our prediction at higher pressures. Possible errors in our EOS at such densities include, in what we estimate to be decreasing order of importance, (1) the shift from the fcc to the bcc crystal, with a corresponding change in $`\mathrm{\Phi }_0^l`$ as discussed in the last Subsection, (2) deviations in the melt curve $`T_m(P)`$ from the Boehler-Ross form at higher pressures (the densities of the points in  correspond to melt pressures around 620 GPa according to our EOS), (3) the fact that at such high $`T`$ the EOS is probing the high-energy region of $`n^\mathrm{s}(ϵ)`$, and (4) the neglect of the anharmonic and boundary contributions to the liquid EOS ($`T`$ is only slightly below $`5T_m`$ at these points according to our EOS).

## 4 Conclusions

Drawing upon theory developed in , we have described a framework for constructing EOS for elemental solids and liquids, and we have discussed experimental and theoretical results indicating that the framework remains highly accurate at low pressures when certain small effects (anharmonicity, boundaries, electron-phonon interactions) are neglected. After displaying the resulting formulas for the Helmholtz free energy, we considered the information one needs to evaluate them, and we discussed the combination of experiment and theoretical work that could be used to get this information. Finally, as an illustration we constructed an EOS for Al, established its range of validity based on the inputs to the EOS, and compared it with Hugoniot data to 5 Mbar; our EOS matched the data to the accuracy we expected based on the low-pressure results.

We consider the primary advantage of this method to lie in the fact that it incorporates into the decomposition of the Hamiltonian a great deal of accumulated knowledge of condensed matter physics both for the solid and liquid phases (for example, the fact that the electronic ground state energy is the most appropriate potential for the nuclear motion). If we have indeed captured the correct physics (and we expect no new physics to enter until the relativistic domain), then the EOS should have the right functional form, which means that if it is shown to agree with available data, then we have reason to believe that it will be equally accurate in regions where no data are available; and making predictions where we have no data is the point of having an EOS to begin with. Furthermore, the better the foundation we can build, the better our position for intelligently investigating and controlling our approximations.

This discussion bears on the second goal of this paper, which was to learn whether the approximation of neglecting anharmonic, boundary, and electron-phonon effects remains useful at higher pressures. We already knew, as discussed in the Introduction, that at low $`P`$ the anharmonic and electron-phonon terms are small, and we found this from direct calculations; we also knew that for several elements, over a range of $`T`$, at low $`P`$ the approximations in question yielded thermal energies and entropies that disagreed with experiments by $`5\%`$ at most. Our work here has shown that for one material at much higher $`P`$ the approximations yield results that match data along a single curve, the Hugoniot, to comparable accuracy. Based on our arguments above, that the EOS incorporates the correct physics and is thus of the correct functional form, we claim to have shown that this Al EOS is trustworthy throughout its range of validity, for all $`T`$ and $`P`$.

The main disadvantage of this method is that it relies on many inputs ($`\mathrm{\Phi }_0`$, $`g(\omega )`$, and $`n(ϵ)`$ for each phase) which may be available only over limited ranges, and each of these limits also restricts the range of validity of the EOS. Our Al example amply illustrates this problem; with a compression range from a little under one to just over two, and a temperature range that reaches only to slightly under 4 eV, this EOS is inadequate for many applications at the national laboratories. We argue, however, that this problem does not indicate a deficiency in the approach; it only underscores the need for many more DFT calculations of these quantities for more materials with ever greater accuracy over ever larger ranges.

In the meantime, though, we would like to be able to say something about elemental solid and liquid EOS at higher compressions. We do know that as density increases, the electrons come to dominate the free energy, and it is also known that TFD correctly describes the electrons in the limit of infinite density. This suggests the following possibility: Construct an EOS using the present techniques to compressions as high as the available experimental or DFT results allow, and then interpolate between these results and the predictions of TFD for higher compressions. This raises an important question: At what pressures does TFD begin to become accurate? Conventional wisdom, usually traced back to Feynman et al. , has held that TFD becomes reliable starting at $`P10`$ Mbar, but other work  suggests that TFD (or TF in their case, but TF and TFD converge at high pressures) deviates noticeably from electronic structure results until 100 Mbar at least. This suggests that the pressure threshold at which TFD is trustworthy has not yet been adequately established; it would be of great interest to settle this question more definitively.

Acknowledgment

This work was supported by the U. S. DOE through contract W-7405-ENG-36.


format: text

text:
# Equivalent forms of the Bessis-Moussa-Villani conjecture

(January 17, 2003)

## Abstract

The BMV conjecture for traces, which states that $`\mathrm{Tr}\mathrm{exp}(𝖠\lambda 𝖡)`$ is the Laplace transform of a positive measure, is shown to be equivalent to two other statements: (i) The polynomial $`\lambda \mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ has only non-negative coefficients for all $`𝖠,𝖡0`$, $`p`$ and (ii) $`\lambda \mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ is the Laplace transform of a positive measure for $`𝖠,𝖡0`$, $`p>0`$.

$\*$$\*$footnotetext: Work partially supported by U.S. National Science Foundation grant PHY 01 39984.$\dagger$$\dagger$footnotetext: Erwin Schrödinger Fellow, supported by the Austrian Science Fund.
© 2003 by the authors. This paper may be reproduced, in its entirety, for non-commercial purposes.

Dedicated with best wishes to Giovanni Jona-Lasinio on his 70th birthday

An intriguing conjecture , which is more than a quarter century old, was formulated by Bessis, Moussa and Villani (BMV) in an attempt to simplify the calculation of partition functions of quantum mechanical systems. It refers to a positivity property of traces of matrices which, if true, would permit the calculation of explicit error bounds in a sequence of approximations known as the Padé approximants.

The BMV conjecture is easy to state. Let $`𝖠`$ and $`𝖡`$ be any two Hermitian $`n\times n`$ matrices with $`𝖡0`$, and let $`\lambda `$. Then the function of $`\lambda `$ given by

$$\lambda \mathrm{Tr}\mathrm{exp}(𝖠\lambda 𝖡)$$
(1)
(with $`\mathrm{Tr}=`$ trace) is the Laplace transform of a positive measure on $`[0,\mathrm{})`$.

Oddly, the positivity can be easily seen to be true for quantum mechanics defined by the Schrödinger equation for bosons without magnetic fields because the partition function in that case can be represented by Wiener integrals. This representation fails for fermions, and hence the importance of the conjecture for condensed matter physics.

The BMV conjecture appears naturally in other areas of matrix analysis as well and there is a significant literature devoted to it. For $`2\times 2`$ matrices there is an easy proof of its correctness, but for $`3\times 3`$ matrices neither a proof nor a counterexample is known to this day. Some recent work that shows the validity of the conjecture in some ‘average’ or ‘typical’ sense is in . Ref.  contains a review.

The purpose of the present paper is to show that the BMV conjectured positivity is equivalent to two other positivity statements that appear to be somewhat different. One of them is, in fact, easier than the original conjecture because it involves only a statement about certain polynomials. Namely for every pair of positive matrices $`𝖠`$ and $`𝖡`$ and every positive integer $`p`$, the polynomial $`\lambda \mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ has all its coefficients positive. It is relatively easy to see that this statement implies the BMV conjecture but the striking fact, which is the main point of our paper, is that the BMV conjecture implies the polynomial statement. The current consensus is that the BMV conjecture is true, in which case our equivalence leads to a non-obvious and interesting theorem about traces. It also presents an alternative route to proving or disproving the BMV conjecture. We are indebted to D. Petz for urging us to write this note about the equivalences.

###### THEOREM 1.

For fixed $`n`$ let $`𝖠`$ and $`𝖡`$ denote arbitrary hermitian $`n\times n`$ matrices over $``$, and let $`\lambda `$. The following statements are equivalent:

* For all $`𝖠`$ and $`𝖡`$ positive, and all $`p`$, the polynomial $`\lambda \mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ has only non-negative coefficients.
* For all $`𝖠`$ hermitian and $`𝖡`$ positive, $`\lambda \mathrm{Tr}\mathrm{exp}(𝖠\lambda 𝖡)`$ is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$.
* For all $`𝖠`$ positive definite and $`𝖡`$ positive, and all $`p0`$, $`\lambda \mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$.

###### Remark 1.

$`𝖠`$ positive means that $`𝖠`$ is Hermitian and $`(v,𝖠v)0`$ for all vectors $`v`$. Positive definite means positive and invertible.

###### Remark 2.

By Bernstein’s theorem , a function of $`\lambda `$ defined for $`\lambda 0`$ is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$ if and only if it is completely monotone, i.e., its $`r`$-th derivative has the sign $`(1)^r`$ for all values of $`\lambda 0`$.

###### Remark 3.

In the case of $`2\times 2`$ matrices, it is easy to see that statement (i) of Theorem 1 (and therefore all three statements) are true: there exists a basis in which both $`𝖠`$ and $`𝖡`$ have only non-negative entries, and therefore all coefficients of $`\mathrm{Tr}(𝖠+\lambda 𝖡)^p`$, which is independent of the choice of the basis, are necessarily non-negative. Some other special cases are treated in .

###### Remark 4.

Note the logical structure in the proof that statement (iii) for $`p`$ implies item (i), which then implies (iii) for all $`p>0`$.

###### Remark 5.

The proof will show (see (13)) that item (iii) is really a very simple corollary of item (ii). Indeed, the proof extends to the statement that $`\lambda \mathrm{Tr}f(𝖠+\lambda 𝖡)`$ is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$ for $`𝖠0,𝖡0`$ whenever $`f`$ is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$.

###### Remark 6.

An important remark about item (i) is addressed in . When $`\mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ for $`p`$ is multiplied out, the coefficient of $`\lambda ^r`$ is a sum of terms.  shows that some of these terms can be negative. The first interesting case occurs for $`p=6`$ and $`r=3`$. The term $`\mathrm{Tr}𝖠^2𝖡^2\mathrm{𝖠𝖡}`$ can be negative even though $`𝖠`$ and $`𝖡`$ are positive. Nevertheless, the BMV conjecture implies that the sum of all terms of order $`\lambda ^3`$ is non-negative. We cannot prove this, even for $`n=3,r=3,p=6`$.

The proofs of (i)$``$(ii) and (ii)$``$(iii) will turn out to be very easy. The difficult direction is (iii)$``$(i) and the following lemma is needed for this case.

###### LEMMA 1.

Let $`𝖺`$ and $`𝖻`$ be hermitian $`n\times n`$ matrices over $``$, with $`𝖺`$ positive definite. Define $`𝖠=𝖺^1`$ and $`𝖡=𝖺^{1/2}\mathrm{𝖻𝖺}^{1/2}`$, and let $`\lambda `$. For all $`p,r`$

$$\frac{d^r}{d\lambda ^r}\mathrm{Tr}\frac{1}{\left(𝖺+\lambda 𝖻\right)^p}|_{\lambda =0}=\frac{p}{p+r}(1)^r\frac{d^r}{d\lambda ^r}\mathrm{Tr}\left(𝖠+\lambda 𝖡\right)^{p+r}|_{\lambda =0}.$$
(2)
###### Proof.

By induction it is easy to show that

$$\frac{d^r}{d\lambda ^r}\left(𝖠+\lambda 𝖡\right)^{p+r}=r!\underset{\begin{array}{c}0i_1,\mathrm{},i_{r+1}p\\ _ji_j=p\end{array}}{}(𝖠+\lambda 𝖡)^{i_1}𝖡\mathrm{}𝖡(𝖠+\lambda 𝖡)^{i_{r+1}}.$$
(3)
By taking the trace at $`\lambda =0`$ we obtain

$$I_1\frac{d^r}{d\lambda ^r}\mathrm{Tr}\left(𝖠+\lambda 𝖡\right)^{p+r}|_{\lambda =0}=r!\underset{\begin{array}{c}0i_1,\mathrm{},i_{r+1}p\\ _ji_j=p\end{array}}{}\mathrm{Tr}𝖠^{i_1}𝖡\mathrm{}\mathrm{𝖡𝖠}^{i_{r+1}}.$$
(4)
Moreover, by similar arguments,

$$\frac{d^r}{d\lambda ^r}\frac{1}{\left(𝖺+\lambda 𝖻\right)^p}=(1)^rr!\underset{\begin{array}{c}1i_1,\mathrm{},i_{r+1}p\\ _ji_j=p+r\end{array}}{}\frac{1}{(𝖺+\lambda 𝖻)^{i_1}}𝖻\mathrm{}𝖻\frac{1}{(𝖺+\lambda 𝖻)^{i_{r+1}}}.$$
(5)
By taking the trace at $`\lambda =0`$ and using cyclicity, we get

$$I_2\frac{d^r}{d\lambda ^r}\mathrm{Tr}\frac{1}{\left(𝖺+\lambda 𝖻\right)^p}|_{\lambda =0}=(1)^rr!\underset{\begin{array}{c}0i_1,\mathrm{},i_{r+1}p1\\ _ji_j=p1\end{array}}{}\mathrm{Tr}𝖠𝖠^{i_1}𝖡\mathrm{}\mathrm{𝖡𝖠}^{i_{r+1}}.$$
(6)
We have to show that

$$I_2=\frac{p}{p+r}(1)^rI_1.$$
(7)
To see this we rewrite $`I_1`$ in the following way. Define $`p+r`$ matrices $`𝖬_j`$ by

$$𝖬_j=\{\begin{array}{cc}𝖡\hfill & \mathrm{for}1jr\hfill \\ 𝖠\hfill & \mathrm{for}r+1jr+p.\hfill \end{array}$$
(8)
Let $`𝒮_n`$ denote the permutation group. Then

$$I_1=\frac{1}{p!}\underset{\pi 𝒮_{p+r}}{}\mathrm{Tr}\underset{j=1}{\overset{p+r}{}}𝖬_{\pi (j)}.$$
(9)
Because of the cyclicity of the trace we can always arrange the product such that $`𝖬_{p+r}`$ has the first position in the trace. Since there are $`p+r`$ possible locations for $`𝖬_{p+r}`$ to appear in the product above, and all products are equally weighted, we get

$$I_1=\frac{p+r}{p!}\underset{\pi 𝒮_{p+r1}}{}\mathrm{Tr}𝖠\underset{j=1}{\overset{p+r1}{}}𝖬_{\pi (j)}.$$
(10)
On the other hand,

$$I_2=(1)^r\frac{1}{(p1)!}\underset{\pi 𝒮_{p+r1}}{}\mathrm{Tr}𝖠\underset{j=1}{\overset{p+r1}{}}𝖬_{\pi (j)},$$
(11)
so we arrive at the desired equality. ∎

###### Proof of Theorem 1.

(i)$``$(ii): By continuity of the trace,

$$\mathrm{Tr}\mathrm{exp}(𝖠\lambda 𝖡)=e^𝖠\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{k!}\mathrm{Tr}(𝖠+𝖠𝕀\lambda 𝖡)^k,$$
(12)
where $`𝕀`$ denotes the identity matrix. It follows from Bernstein’s theorem and (i) that the right side is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$.

(ii)$``$(iii): This follows from taking the trace of the matrix equation

$$\frac{1}{(𝖠+\lambda 𝖡)^p}=\frac{1}{\mathrm{\Gamma }(p)}_0^{\mathrm{}}\mathrm{exp}\left[t(𝖠+\lambda 𝖡)\right]t^{p1}𝑑t.$$
(13)
(iii)$``$(i): It suffices to assume (iii) only for $`p`$. For invertible $`𝖠`$ we observe that the $`r`$-th derivative of $`\mathrm{Tr}(𝖺+\lambda 𝖻)^p`$ at $`\lambda =0`$ is related to the coefficient of $`\lambda ^r`$ in $`\mathrm{Tr}(𝖠+\lambda 𝖡)^p`$ as given by (2) with $`𝖠,𝖺,𝖡,𝖻`$ related as in Lemma 1. The left side of (2) has the sign $`(1)^r`$ because it is the derivative of a function that is the Laplace transform of a positive measure supported in $`[0,\mathrm{})`$. Thus the right side has the correct sign as stated in item (i). The case of non-invertible $`𝖠`$ follows from continuity. ∎


text:
# Universal Behavior of Charged Particle Production in Heavy Ion Collisions at RHIC Energies

## 1 Introduction

The PHOBOS experiment has measured $`dN_{ch}/d\eta `$ and the average multiplicity of charged particles $`N_{ch}`$ produced in heavy ion collisions for center of mass energies in the nucleon-nucleon center of mass system, $`\sqrt{s_{NN}}`$, of 19.6, 130 and 200 GeV. The data is also binned as a function of event centrality (impact parameter) characterized by the number of participating nucleons, $`N_{part}`$, allowing comparisons to elementary systems, like $`pp/\overline{p}p`$ and $`e^+e^{}\mathrm{hadrons}`$.

The PHOBOS multiplicity detector consists of several arrays of silicon detectors which cover nearly the full solid angle for collision events. The event centrality is characterized by the multiplicity of charged particles measured by two sets of 16 paddle counters covering $`3<|\eta |<4.5`$. The methods used for measuring the multiplicity of charged particles as well as for extracting $`N_{part}`$ has been described in more detail in Ref. .

## 2 Limiting Behavior in Pseudorapidity Distributions

Fig. 1 shows $`dN_{ch}/d\eta ^{}/N_{part}/2`$ ($`\eta ^{}=\eta y_{\mathrm{beam}}`$) measured at three different RHIC energies for peripheral ($`N_{part}100`$) and central events ($`N_{part}355`$), in the left and right panels, respectively. These show a clear “limiting behavior” in the fragmentation region. That is, the distributions are independent of beam energy in a substantial range in $`\eta ^{}`$. As the beam energy increases, $`dN/d\eta ^{}`$ follows the universal trend until it reaches 85-90% of it’s maximum value at midrapidity, at which point it stops following the trend. Similar behavior has been observed in elementary collisions as well, both in $`\overline{p}p`$ collisions  and in $`e^+e^{}`$ collisions over a large range of energies .

The limiting curve constrains the energy dependence of the charged particle multiplicity. It also varies with centrality in such a way that the increases seen at low $`\eta ^{}`$ (which is midrapidity in $`\eta `$) as $`N_{part}`$ increases, are accompanied by decreases near $`\eta ^{}0`$ (forward rapidities), as seen in Fig. 1. It is not clear why this behavior occurs, e.g. whether it is from energy conservation or a true long-range correlation.

## 3 Comparison with Elementary Systems

The angular distributions of charged particles for different strongly-interacting systems are shown in Fig. 2, where central Au+Au (divided by $`N_{part}/2`$), $`\overline{p}p`$  and $`e^+e^{}`$  data are compared, all at $`\sqrt{s}=200`$ GeV. In this figure, the Au+Au and $`\overline{p}p`$ data are shown as $`dN_{ch}/d\eta `$, while the $`e^+e^{}`$ data (with cuts applied to reject events with large initial-state photon radiation) is shown as $`dN/dy_T`$, the rapidity distribution along the event thrust axis, calculated assuming the pion mass. The shapes of Au+Au and $`e^+e^{}`$ are similar (within 10%) in shape and magnitude, especially within $`|\eta |<4`$. In the lower panel, we also observe that the shapes of the Au+Au and $`\overline{p}p`$ distributions are also very similar over a large range in $`\eta `$, but the integrals differ by $`40\%`$.

In Fig. 3, we compare $`N_{ch}/N_{part}/2`$ in heavy ion collisions  to $`e^+e^{}`$ and $`pp/\overline{p}p`$ data over a large range in $`\sqrt{s}`$ . It is observed that $`N_{ch}/N_{part}/2`$ lies below $`pp`$ at low energies, passes through the $`pp`$ data around $`\sqrt{s}10`$ GeV, and then gradually joins with the $`e^+e^{}`$ trend above CERN SPS energies. These comparisons can be seen more clearly by dividing all of the data by a fit to the $`e^+e^{}`$ data .

The $`pp/\overline{p}p`$ data follows the same trend as $`e^+e^{}`$, but it can be shown that it matches very well if the “effective energy” $`\sqrt{s}_{\mathrm{eff}}=\sqrt{s}/2`$ is used, which accounts for the leading particle effect seen in $`pp`$ collisions . Ref.  finds that bulk particle production in $`pp`$ and $`e^+e^{}`$ data does not depend in detail on the collision system but rather the energy available for particle production. In this scenario, the Au+Au data suggests a substantially reduced leading particle effect in central collisions of heavy nuclei at high energy.

The alleviation of the leading particle effect might not be so surprising in nuclear collisions. Each participating nucleon is typically struck $`\overline{\nu }=\frac{N_{coll}}{N_{part}/2}>3`$ times on average as it passes through the oncoming gold nucleus for $`N_{part}>65`$. One could speculate that the multiple collisions transfer much more of the initial longitudinal energy into particle production. This naturally leads to the scaling of total particle production in heavy ion collisions with $`N_{part}`$, as seen in Fig. 4, reminiscent of the “wounded nucleon model”  but with the scaling factor determined by $`e^+e^{}`$ rather than $`pp`$.

In conclusion, PHOBOS has observed two kinds of universal behavior. The first is an energy-independent, but centrality-dependent, universal limiting distribution of charged particle production away from midrapidity. The second is that the total charged particle multiplicity per participant pair in heavy ion collisions above CERN SPS energies scales with sqrt(s) in a similar way as $`e^+e^{}`$ collisions. These features may offer a new perspective on particle production in heavy ion collisions.

Special thanks to the ICHEP organizers and Heavy Ion convenors for the invitation to speak at the conference. This work was partially supported by U.S. DOE grants DE-AC02-98CH10886, DE-FG02-93ER40802, DE-FC02-94ER40818, DE-FG02-94ER40865, DE-FG02-99ER41099, and W-31-109-ENG-38 as well as NSF grants 9603486, 9722606 and 0072204. The Polish groups were partially supported by KBN grant 2 PO3B 10323. The NCU group was partially supported by NSC of Taiwan under contract NSC 89-2112-M-008-024.


text:
# The Effect of Radiative Cooling on the Sunyaev-Zel’dovich Cluster Counts and Angular Power Spectrum: Analytic Treatment

## 1 Introduction

With the rapid development of observing techniques (Birkinshaw, Hughes & Arnaud 1991; Jones et al. 1993; Birkinshaw & Hughes 1994; Carlstrom, Joy & Grego 1996; Myers et al. 1997; Hughes & Birkinshaw 1998; etc.), the Sunyaev-Zel’dovich (SZ) effect has become one of the most powerful tools for the detections of high redshift clusters (Joy et al. 2001) and cosmic microwave background (CMB) anisotropy on small scales (Mason et al. 2002). Indeed, the redshift-independence is the major advantage of non-targeted SZ surveys over traditional optical and X-ray observations. This arises from the fact that the SZ effect depends uniquely on the intrinsic properties of the warm-hot gas associated with cosmic structures, whilst the photons interacting with the gas come from CMB at very high redshift $`z1000`$. Because robust constraints on cosmological models are provided by the most massive and distant clusters, growing interest over the past years has been focused upon how well the fundamental cosmological parameters can be constrained with non-targeted SZ cluster surveys (Molnar, Birkinshaw & Mushotzky 2002) and SZ power spectrum (Bond et al. 2002). The sensitivity of the expected SZ cluster counts and SZ power spectrum to the underlying cosmological model is quite impressive.

On arcminute scales, the strength of thermal SZ signals is directly proportional to the total thermal energy of the hot gas confined in clusters. Unlike the dark matter component of clusters whose dynamical behavior is governed purely by gravity, the intracluster gas is easily disturbed by many complicated physical processes in addition to gravity and thermal pressure, which include (non)gravitational heating, radiative cooling, star formation and energy feedback, magnetic fields, etc. The reliability of cosmological applications of SZ surveys on arcminute scales is closely connected with the question of how well one can handle these non-gravitational mechanisms. At present, a sophisticated treatment of the problem must rely on hydrodynamical simulations coupled with some semi-analytic approximations (e.g. da Silva et al. 2000, 2001; Seljak, Burwell & Pen 2001; Zhang, Pen & Wang 2002; White et al. 2002). Yet, the effective probe of cosmological models using non-targeted SZ surveys can be more easily achieved by semi-analytic approaches because of the requirement of the continuity in cosmological parameter space. It is thus desirable to understand the essential physics that dominates the dynamical evolution of intracluster gas.

While intracluster gas is mainly driven by gravitational shocks and adiabatic compression, two completely different mechanisms have been suggested thus far that may significantly affect the distribution of the intracluster gas: nongravitational preheating and radiative cooling. This is primarily motivated by the steepening of the X-ray luminosity-temperature relation and the entropy excess of groups and clusters reported by a number of investigators over the past years (e.g. David et al. 1993; Wu, Xue & Fang 1999; Ponman, Cannon & Navarro 1999). However, it has been realized that the prevailing preheating scenario is facing the so-called energy crisis – an unreasonably high efficiency of energy injection into the intracluster medium from supernovae must be required (Wu, Fabian & Nulsen 2000). Moreover, a uniform preheating of the cosmic baryons to a temperature of $`10^6`$ K would make the Ly$`\alpha `$ forest to disappear at high redshifts. In contrast, radiative cooling is a natural process of the hot intracluster gas, in which no exotic physics is needed. Several recent studies have shown that radiative cooling alone can allow one to successfully reproduce not only the observed distributions of global X-ray luminosity and entropy but also the internal structures of hot gas in groups and clusters (Bryan 2000; Pearce et al. 2000; Muanwong et al. 2001, 2002; Voit & Bryan 2001; Wu & Xue 2002a,b; Borgani et al. 2002; Voit et al. 2002). If confirmed, this would have profound implications for our understanding of the evolution of hot gas and the formation of stars in the most massive systems in the universe.

So far, only can one employ hydrodynamical simulations to assess the effect of radiative cooling on the SZ counts and power spectrum (da Silva et al. 2000; White et al. 2002), while in the preheating scenario several (semi)analytic models have been applied for the prediction of SZ cluster surveys (Cavaliere & Menci 2001; Holder & Carlstrom 2001; Benson, Reichardt & Kamionkowski 2002). Now, a fully analytic treatment of radiative cooling effect on the prediction of SZ cluster counts and power spectrum may become possible if the observed central entropy floor of groups and clusters can be attributed to radiative cooling of hot intracluster/intragroup gas (Voit & Brayn 2001). In such a scenario, the minimum entropy of gas distribution in each cluster can be uniquely determined by the conservation of energy, which is equivalent to specifying the equation of state for the intracluster gas. Consequently, we will be able to derive the gas distribution of clusters at various redshifts as a result of cooling, in combination with the hydrostatic equilibrium hypothesis, provided that the underlying dark matter profile can be approximated by some kinds of analytic form, e.g., the universal density profile (Navarro, Frenk & White 1997; NFW). Together with the Press-Schechter (1994; PS) formalism for the abundance of dark halos at different cosmic epoch, we can eventually compare the predicted SZ counts and power spectra with and without radiative cooling, and demonstrate the uncertainty in the determination of the cosmological parameters arising from radiative cooling.

## 2 Gas distribution with and without cooling

### 2.1 Dark matter distribution

Intracluster gas with and without cooling is always assumed to be in hydrostatic equilibrium with the underlying gravitational potential dominated by dark matter component $`\rho _{\mathrm{DM}}`$:

$$\frac{1}{\mu m_\mathrm{p}n_\mathrm{e}}\frac{d(n_\mathrm{e}k_\mathrm{B}T)}{dr}=\frac{GM_{\mathrm{DM}}(r)}{r^2},$$
(1)
where $`n_\mathrm{e}`$ and $`T`$ are the electron density and temperature, respectively, and $`\mu =0.585`$ is the mean molecular weight. We use the universal density profile suggested by numerical simulations (NFW) for $`\rho _{\mathrm{DM}}`$

$$\rho _{\mathrm{DM}}(r)=\frac{\delta _c\rho _{\mathrm{crit}}}{(r/r_s)(1+r/r_s)^2},$$
(2)
where $`\delta _c`$ and $`r_s`$ are the characteristic density and length of the halo, respectively, and $`\rho _{\mathrm{crit}}`$ is the critical density of the universe at cosmic time $`t`$. We follow the prescription of Eke, Navarro & Steinmntz (2001) to fix the two free parameters, $`\delta _c`$ and $`r_s`$, in the NFW profile. To do this, the concentration parameter $`c=r_{\mathrm{vir}}/r_s`$ of a dark halo identified at redshift $`z`$ is related to the collapsing redshift $`z_{\mathrm{coll}}`$ through

$$c^3=\frac{\mathrm{\Delta }_\mathrm{c}(z_{\mathrm{coll}})}{\mathrm{\Delta }_\mathrm{c}(z)}\frac{\mathrm{\Omega }_\mathrm{M}(z)}{\mathrm{\Omega }_\mathrm{M}(z_{\mathrm{coll}})}\left(\frac{1+z_{\mathrm{coll}}}{1+z}\right)^3,$$
(3)
where $`\mathrm{\Delta }_\mathrm{c}`$ is the overdensity of dark halo within virial radius $`r_{\mathrm{vir}}`$ with respect to $`\rho _{\mathrm{crit}}`$, for which we take $`\mathrm{\Delta }_\mathrm{c}=18\pi ^2+82[\mathrm{\Omega }_\mathrm{M}(z)1]39[\mathrm{\Omega }_\mathrm{M}(z)1]^2`$ for a flat universe, and $`\mathrm{\Omega }_\mathrm{M}(z)`$ is the cosmic density parameter. The collapsing redshift $`z_{\mathrm{coll}}`$ is determined by

$$D(z_{\mathrm{coll}})\sigma _{\mathrm{eff}}(M_s)=\frac{1}{C_\sigma },$$
(4)
where $`D`$ is the normalized linear growth factor, $`C_\sigma =25`$ for $`\mathrm{\Lambda }`$CDM model, and $`\sigma _{\mathrm{eff}}(M_s)`$ is the so-called modulated rms linear density at mass scale $`M_s`$

$$\sigma _{\mathrm{eff}}(M_\mathrm{s})=\sigma (M_\mathrm{s})\left[\frac{d\mathrm{ln}\sigma (M_\mathrm{s})}{d\mathrm{ln}M_\mathrm{s}}\right],$$
(5)
and $`M_s`$ corresponds to the mass contained within $`r=2.17r_s`$ where the circular velocity reaches the maximum. Finally, we specify a temperature $`T_{\mathrm{vir}}`$ to the hot gas in cluster of mass $`M`$ at redshift $`z`$ in terms of cosmic virial theorem (Bryan & Norman 1998)

$$k_\mathrm{B}T_{\mathrm{vir}}=1.39f_\mathrm{T}(h^2\mathrm{\Delta }_\mathrm{c}E^2)^{1/3}\mathrm{keV}\left(\frac{M}{10^{15}M_{}}\right)^{2/3},$$
(6)
in which we will choose the normalization factor to be $`f_T=0.92`$, $`E^2=\mathrm{\Omega }_\mathrm{M}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }+(1\mathrm{\Omega }_\mathrm{M}\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2`$, and $`h`$ is the Hubble constant in units of 100 km s-1 Mpc-1.

### 2.2 Gas distribution without cooling

In the absence of cooling, we adopt two models for the gas properties:

Model I – Gas is assumed to follow the same distribution as the dark matter in clusters:

$$n_\mathrm{e}^0=\frac{f_\mathrm{b}}{\mu _\mathrm{e}m_\mathrm{p}}\rho _{\mathrm{DM}}$$
(7)
in which $`f_\mathrm{b}`$ is the universal baryon fraction, and $`\mu _e=2/(1+X)`$ is the mean electron weight with $`X=0.768`$ being the hydrogen mass fraction in the primordial abundance of hydrogen and helium. We solve the equation of hydrostatic equilibrium to get the temperature profile under the boundary restriction $`T^0(r\mathrm{})0`$:

$$k_\mathrm{B}T^0(r)=k_\mathrm{B}T^{}\frac{r}{r_s}\left(1+\frac{r}{r_s}\right)^2_{r/r_s}^{\mathrm{}}\frac{(1+x)\mathrm{ln}(1+x)x}{x^3(1+x)^3}𝑑x,$$
(8)
where $`k_\mathrm{B}T^{}=4\pi G\mu m_p\delta _c\rho _{\mathrm{crit}}r_s^2`$.

Model II – Gas is assumed to be isothermal and $`T^0(r)=T_{\mathrm{vir}}`$. In this case, the electron number density in terms of the equation of hydrostatic equilibrium reads (Makino, Sasaki & Suto 1998)

$$n_\mathrm{e}^0(r)=n_{\mathrm{e0}}e^\alpha \left(1+\frac{r}{r_\mathrm{s}}\right)^{\alpha /(r/r_\mathrm{s})},$$
(9)
where $`\alpha =4\pi G\mu m_\mathrm{p}\delta _c\rho _{\mathrm{crit}}r_\mathrm{s}^2/k_\mathrm{B}T`$. The normalization parameter $`n_{\mathrm{e0}}`$ can be fixed through

$$_0^{r_{\mathrm{vir}}}4\pi \mu _\mathrm{e}m_\mathrm{p}r^2n_\mathrm{e}^0(r)𝑑r=M_{\mathrm{vir}}f_\mathrm{b}.$$
(10)
### 2.3 Gas distribution with cooling

The cooling time scale for a steady, highly subsonic flow is determined by the conservation of energy. Setting the energy loss rate due to bremsstrahlung emission to equal the change in the specific energy of gas yields

$$t_\mathrm{c}=2.869\times 10^{10}\mathrm{yr}\left(\frac{1.2}{g(T^0)}\right)\left(\frac{k_\mathrm{B}T^0}{\mathrm{keV}}\right)^{1/2}\left(\frac{n_\mathrm{e}^0}{10^3\mathrm{cm}^3}\right)^1,$$
(11)
where $`n_\mathrm{e}^0`$ and $`T^0`$ are the electron number density and temperature without cooling, respectively, and $`g`$ is the total Guant factor. If the cooling time $`t_c`$ is chosen to be the age of clusters, or approximately the age of the universe, the above equation would allow us to set up a link between $`n_\mathrm{e}^0`$ and $`T^0`$. Note that this may lead to an overestimate of the cooling effect, and meanwhile $`t_c`$ becomes to be cosmological model dependent. Perhaps, a more reasonable approach is to define the age of a cluster as the cosmic time between its collapsing redshift $`z_{\mathrm{coll}}`$ and the redshift $`z`$ when it is identified (cf. Section 2.1). Once the gas density is specified, we will be able to work out the cooling radius and the amount of gas that cools out of the hot phase by the time $`t_c`$ (e.g. Wu & Xue 2002b). The fate of the cooled materials within cooling radius may be associated with star formation. Now, our task is to work out the new distribution of the intracluster gas with cooling under different assumptions of the equation of state.

Model III – Before and after cooling, the equation of state for gas is always assumed to be isothermal and $`T(r)=T_{\mathrm{vir}}`$. Consequently, the functional form of equation (9) also applies to the cooling case. However, the normalization after cooling is made through

$$_0^{r_{\mathrm{vir}}}4\pi \mu _\mathrm{e}m_\mathrm{p}r^2n_\mathrm{e}(r)𝑑r=M_{\mathrm{vir}}f_\mathrm{b}M_{\mathrm{cool}}$$
(12)
where the total cooled material $`M_{\mathrm{cool}}`$ is given by

$`M_{\mathrm{cool}}`$ $`=`$ $`{\displaystyle _0^{r_{\mathrm{cool}}}}4\pi \mu _em_\mathrm{p}n_\mathrm{e}^0(r)r^2𝑑r,`$ (13)
$`=`$ $`4\pi f_\mathrm{b}\delta _\mathrm{c}\rho _{\mathrm{crit}}r_\mathrm{s}^3\left[\mathrm{ln}\left(1+{\displaystyle \frac{r_{\mathrm{cool}}}{r_\mathrm{s}}}\right){\displaystyle \frac{r_{\mathrm{cool}}}{r_\mathrm{s}+r_{\mathrm{cool}}}}\right],`$ (14)
and $`r_{\mathrm{cool}}`$ is the cooling radius which can be obtained by combining equations (9) and (11).

Model IV – Gas traces dark matter (Model I) before cooling, and the equation of state for the new gas distribution after cooling is given by the entropy before cooling plus a constant entropy floor $`S_c`$ (cf. Holder & Carlstrom 2001; Voit & Bryan 2001)

$$S=\frac{k_\mathrm{B}T}{n_\mathrm{e}^{2/3}}=S_\mathrm{c}+\frac{k_\mathrm{B}T^0}{(n_\mathrm{e}^0)^{2/3}},$$
(15)
where the entropy floor is determined by the cooling time $`t_c`$

$`S_c`$ $`=`$ $`{\displaystyle \frac{k_\mathrm{B}T_{\mathrm{vir}}}{(n_\mathrm{e}^0)^{2/3}}}`$ (16)
$`=`$ $`100\mathrm{keV}\mathrm{cm}^2\left({\displaystyle \frac{t_\mathrm{c}}{2.869\times 10^{10}\mathrm{yr}}}\right)^{2/3}\left({\displaystyle \frac{g}{1.2}}\right)^{2/3}\left({\displaystyle \frac{k_\mathrm{B}T_{\mathrm{vir}}}{\mathrm{keV}}}\right)^{2/3}.`$
This “entropy floor” increases with cosmic time, resulting in a deposition of cooled materials in the central regions of clusters (Voit & Bryan 2001). Solving the equation of hydrostatic equilibrium yields

$$k_\mathrm{B}T(r)=\frac{2}{5}G\mu m_\mathrm{p}S^{3/5}(r)_r^{\mathrm{}}\frac{M_{\mathrm{DM}}(r)}{r^2}S^{3/5}(r)𝑑r,$$
(17)
and $`n_\mathrm{e}=(k_\mathrm{B}T/S)^{3/2}`$.

Model V – The same as Model IV except a varying metallicity $`Z=0.3Z_{}(t/t_0)`$ is assumed instead of $`Z=0.3Z_{}`$ for the rest four models, where $`t_0`$ denotes the present cosmic epoch. The parameters of the five models are summarized in Table 1.

## 3 Expectation for SZ cluster counts

The total SZ flux observed at frequency $`\nu `$ by a cluster of mass $`M`$ at redshift $`z`$ is (e.g. Barbosa et al. 1996)

$$S_\nu (x,M,z)=\frac{g_\nu (x)}{D_a^2(z)}\left(\frac{\sigma _\mathrm{T}}{m_\mathrm{e}c^2}\right)k_\mathrm{B}T(r)n_\mathrm{e}(r)\mathrm{\hspace{0.33em}4}\pi r^2𝑑r,$$
(18)
where

$`g_\nu (x)`$ $`=`$ $`2{\displaystyle \frac{(k_\mathrm{B}T_{\mathrm{CMB}})^3}{(h_\mathrm{p}c)^2}}{\displaystyle \frac{x^4e^x}{(e^x1)^2}}f_\nu (x);`$ (19)
$`f_\nu (x)`$ $`=`$ $`x\mathrm{coth}{\displaystyle \frac{x}{2}}4,`$ (20)
$`x=h_\mathrm{p}\nu /k_\mathrm{B}T_{\mathrm{CMB}}`$, $`T_{\mathrm{CMB}}=2.728`$ K is the temperature of CMB (Fixsen et al. 1996), and $`D_a`$ is the angular diameter distance to the cluster. In the isothermal case, the integral in the right hand of equation (18) can be replaced by the total mass of the cluster if the universal baryon fraction is introduced:

$$S_\nu =\frac{g_\nu (x)}{D_a^2(z)}\left(\frac{k_\mathrm{B}T}{m_\mathrm{e}c^2}\right)\left(\frac{f_b\sigma _\mathrm{T}}{\mu _\mathrm{e}m_\mathrm{p}}\right)M.$$
(21)
This is the standard method adopted in the literature for the theoretical prediction of SZ cluster counts (Model II). While the assumption of isothermality is more or less reasonable in terms of current X-ray observations, the dependence of the gas fraction on temperature claimed by many observations throws doubt upon the direct utilization of equation (21). Nevertheless, if the gas fraction $`f_b`$ is allowed to vary according to mass or temperature, or the total gas mass is replaced by the hot gas component alone, $`M_{\mathrm{gas}}(T)`$, we may modify the above equation to be

$$S_\nu =\frac{g_\nu (x)}{D_a^2(z)}\left(\frac{k_\mathrm{B}T}{m_\mathrm{e}c^2}\right)\left(\frac{\sigma _\mathrm{T}}{\mu _\mathrm{e}m_\mathrm{p}}\right)M_{\mathrm{gas}}(T).$$
(22)
As shown in the above section, $`M_{\mathrm{gas}}(T)`$ can be analytically determined for any cluster within the framework of radiative cooling. This provides a simple approach to estimating the effect of radiative cooling on SZ cluster counts without knowing the density distribution of the hot, isothermal gas inside clusters, which will be applied to Model III. For rest three models, I, IV and V, we will adopt the exact formula of $`S_\nu `$, equation (18), in the theoretical prediction of SZ cluster counts.

The expected number of SZ selected clusters with flux greater than $`S_\nu `$ and in redshift interval ($`z`$, $`z+dz`$) is

$$\frac{dN(>S_\nu )}{dzd\mathrm{\Omega }}=\frac{dV}{dzd\mathrm{\Omega }}_{M_{\mathrm{min}}(z,S_\nu )}^{\mathrm{}}\frac{dn}{dM}𝑑M,$$
(23)
where the mass threshold $`M_{\mathrm{min}}(z,S_\nu )`$ is given by the definition of equation (18) or (21) or (22), depending on what approximation we would use. We adopt the PS mass function to describe the distribution of clusters

$$dn=\sqrt{\frac{2}{\pi }}\frac{\overline{\rho }}{M}\frac{\delta _\mathrm{c}(z)}{\sigma ^2}\frac{d\sigma }{dM}\mathrm{exp}\left(\frac{\delta _c^2(z)}{2\sigma ^2}\right)dM,$$
(24)
where $`\overline{\rho }`$ is the mean cosmic density, $`\delta _c`$ is the linear over-density of perturbations that collapsed and virialized at redshift $`z`$, and $`\sigma `$ is the linear theory variance of the mass density fluctuation in sphere of mass $`M`$:

$$\sigma ^2(M)=\frac{1}{2\pi ^2}\underset{0}{\overset{\mathrm{}}{}}k^2P(k)|W(kR)|^2𝑑k,$$
(25)
and $`W(kR)=3(\mathrm{sin}xx\mathrm{cos}x)/x^3`$ is the Fourier representation of the window function. The power spectrum, $`P(k)k^nT^2(k)`$, is normalized by the rms fluctuation on an $`8`$ $`h^1`$ Mpc scale, $`\sigma _8`$, and we take the transfer function $`T(k)`$ from an adiabatic CDM model given by Bardeen et al. (1986) for the Harrison-Zel’dovich case n=1. Note that our results may be moderately changed if a modified mass function of dark halos is adopted (e.g. Jenkins et al. 2001; Sheth & Tormen 2001).

We work with a flat cosmological model ($`\mathrm{\Lambda }`$CDM) of $`\mathrm{\Omega }_bh^2=0.019`$, $`\mathrm{\Omega }_\mathrm{M}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. We adopt $`\sigma _8=0.90`$ and a Hubble constant of $`h=0.65`$. We choose the SZ flux limit to be $`S_\nu =15`$ mJy at frequency $`\nu =30`$ GHz. We perform numerical calculations for the five models of gas distribution with and without cooling, and demonstrate the resulting differential SZ cluster counts in Figure 1. While there are some differences between the redshift distribution of clusters predicted by the gas-traces-mass assumption and the one by the isothermal assumption, the effect of radiative cooling leads to a moderate decrease in the expected number of SZ clusters although the peak locations with and without cooling remain roughly the same. Alternatively, our prediction alters only slightly if the cosmic evolution of metallicity according to $`0.3Z_{}(t/t_0)`$ is included.

Recall that in the standard treatment, the expectation for SZ cluster counts is made by presumably taking the intracluster gas to be isothermal and without cooling. If this prediction is directly applied to the determination of cosmological parameters in future SZ surveys, large uncertainty may be introduced because of the ignorance of radiative cooling correction. Now we can evaluate the uncertainty in the determination of cosmological parameters by comparing the theoretical predictions with and without radiative cooling. To do this, we still work with a flat cosmological model of $`\mathrm{\Omega }_\mathrm{M}+\mathrm{\Omega }_\mathrm{\Lambda }=1`$, but allow $`\mathrm{\Omega }_\mathrm{M}`$ to vary according to $`\sigma _8`$ until the predicted redshift distribution of SZ cluster counts by the standard model (II) matches those by the cooling models III–V.

Figure 2 shows the resulting $`\sigma _8`$ versus $`\mathrm{\Omega }_\mathrm{M}`$, in which the contours indicate $`68\%`$ joint confidence intervals on the two parameters, corresponding to $`\mathrm{\Delta }\chi ^2=2.30`$. Uncertainties in $`\mathrm{\Omega }_\mathrm{M}`$ and $`\sigma _8`$ due to radiative cooling correction can be easily demonstrated by noticing that the true input values in the cooling models are $`\mathrm{\Omega }_\mathrm{M}=0.30`$ and $`\sigma _8=0.90`$, while our best-fit results are $`(\mathrm{\Omega }_\mathrm{M},\sigma _8)=(0.24_{0.11}^{+0.13},0.89_{0.10}^{+0.08})`$, $`(0.23_{0.11}^{+0.11},0.88_{0.10}^{+0.08})`$ and $`(0.23_{0.11}^{+0.10},0.89_{0.10}^{+0.07})`$ for the three cooling models III, IV, and V, respectively. Namely, without the correction of radiative cooling we might underestimate both $`\mathrm{\Omega }_\mathrm{M}`$ and $`\sigma _8`$ parameters for $`\mathrm{\Lambda }`$CDM cosmological model, although within the uncertainties of our ’observations’, the best-fit parameters are still consistent with the input values.

## 4 SZ Power spectrum

The angular power spectrum of temperature fluctuation on the CMB sky due to the SZ effect of clusters can be separated into the Poisson term $`C_l^{(P)}`$ and clustering term $`C_l^{(C)}`$ (Cole & Kaiser 1988; Komatsu & Kitayama 1999):

$`C_l^{(P)}`$ $`=`$ $`f_\nu ^2(x){\displaystyle _0^{z_{\mathrm{dec}}}}𝑑z{\displaystyle \frac{dV}{dz}}`$ (27)
$`{\displaystyle _{M_{\mathrm{min}}}^{\mathrm{}}}𝑑M{\displaystyle \frac{dn(M,z)}{dM}}|y_l(M,z)|^2,`$
and

$`C_l^{(C)}`$ $`=`$ $`f_\nu ^2(x){\displaystyle _0^{z_{\mathrm{dec}}}}𝑑z{\displaystyle \frac{dV}{dz}}P(k=l/D_0,z)`$ (28)
$`\left[{\displaystyle _{M_{\mathrm{min}}}^{\mathrm{}}}𝑑M{\displaystyle \frac{dn(M,z)}{dM}}b(M,z)y_l(M,z)\right]^2,`$
where $`z_{\mathrm{dec}}1000`$ is the CMB photon decoupling redshift, $`D_0`$ is the comoving distance to cluster of mass $`M`$ at $`z`$, $`y_l(M,z)`$ is the Fourier transform of the Compton $`y`$-parameter in the thermal SZ effect: $`\mathrm{\Delta }T/T_{\mathrm{CMB}}=f_\nu (x)y(\theta )`$, $`f_\nu (x)`$ reflects the spectral dependence, and $`b(M,z)`$ is the so-called bias parameter, for which we use the analytic approximation of Mo & White (1996). In our numerical computation, the minimum cluster mass is taken to be $`M_{\mathrm{min}}=1\times 10^{13}M_{}`$, and our final results are unaffected by this choice if we focus on large-scale ($`>`$arcminute) fluctuations on the CMB sky.

We illustrate in Figure 3 the SZ power spectra with and without cooling for the five models in Table 1, together with the primary CMB signal produced by CMBFAST. It appears that the power spectra predicted by the gas-traces-mass and isothermal assumptions become indistinguishable, and the total SZ power spectrum is dominated by the Poisson distribution of groups and clusters. Inclusion of radiative cooling leads to a significant decrease of SZ power spectrum, especially at small angular scales or high-$`\mathrm{}`$. This is simply because cooling removes more efficiently the hot gas at the central regions of clusters and in low-mass groups, and therefore suppresses the SZ signal at small scales. This result is in good agreement with the one found by hydrodynamical simulations (e.g. da Silva et al. 2001).

## 5 Discussion and conclusions

As a natural process, radiative cooling plays an important role in the formation of galaxies. In the central regions of clusters and groups, the typical cooling time scale of the hot gas is usually less than the age of the universe. It follows that a considerable amount of the intracluster/intragroup gas must have cooled out of the hot phase since the formation of clusters and groups unless the cooling efficiency can be suppressed by other non-gravitational heating processes such as energy feedback from supernovae and AGN activity. Since cooling removes the hot gas in groups and clusters, it meanwhile reduces the SZ signal, leading to a decease of the expected SZ cluster counts and power spectrum. This scenario has been confirmed by numerical simulations (e.g. da Silva et al. 2001) and our calculations based on a simple analytical model. If the excess entropy detected in the central cores of groups and clusters can be attributed to radiative cooling (Voit & Bryan 2001), some conclusions drawn from a phenomenological preheating model based on the observed entropy distribution can be equally applied to the cooling scenario. Indeed, it has been shown that the preheating model also results in a decrease of SZ cluster counts and power spectrum (e.g. Komatsu & Kitayama 1999; Holder & Carlstrom 2001; da Silva et al. 2001). Recall that current X-ray observations alone are still unable to distinguish between preheating and cooling models (e.g. Voit et al. 2002; Borgani et al. 2002).

It is well known that radiative cooling is a runaway process. Consequently, the expected cooled materials within the framework of radiative cooling from both hydrodynamical simulations and analytical approaches exceed the observed stellar mass fraction in the local universe (Balogh et al. 2001; Wu & Xue 2002b and references therein). Inclusion of star formation and energy feedback, which should in principle resolve the so-called cosmic cooling crisis, may significantly reduce the effect of radiative cooling on the X-ray properties and SZ signal of clusters. Indeed, recent numerical simulations by White et al. (2002) have shown that as a consequence of the combined effect of energy injection by star formation and radiative cooling, the SZ power spectrum remains roughly unchanged.

Our simple analytic approach to estimating the effect of radiative cooling on the non-targeted SZ cluster survey and SZ power spectrum yields a result consistent with those found by hydrodynamical simulations (e.g. da Silva et al. 2001). This analytic method allows an effective evaluation of the significance of various factors in the gas cooling process. For example, we have shown that within the framework of radiative cooling, both SZ cluster counts and power spectrum are insensitive to temperature variation. In other words, the isothermality hypothesis is a good approximation in the theoretical study of cluster SZ effect. Also, one may choose to use the gas-traces-mass assumption, which provides a reasonable description of the gas distribution outside the X-ray core. Recall that unlike the X-ray emission which depends on $`n_\mathrm{e}^2`$, the SZ signal is proportional to $`n_\mathrm{e}`$. So, the gas outside the X-ray core makes an important contribution to SZ signal too. Finally, in order to eliminate the cooling crisis in the present model and include another natural process, i.e., heating by supernovae, in the evolution of hot gas, it is worth exploring how energy injection into the intracluster/intrgroup gas from star formation is incorporated into cooling scenario in a simple analytic way.

This work was supported by the National Science Foundation of China, and the Ministry of Science and Technology of China, under Grant No. NKBRSF G19990754.