text:
# Stability of the iterative solutions of integral equations as one phase freezing criterion

## I Introduction

When studying the structure and thermodynamics of classical fluids one is often faced with the task of solving the nonlinear integral equation which stems out of the combination of the Ornstein-Zernike equation and an approximate relation between pair potential and correlation functions (the closure) Hansen and McDonald (1986). Integral equations can be generally written in the form

$`\gamma (r)=A\gamma (r),`$ (1)
where $`\gamma (r)S`$ may be the total correlation function $`h(r)`$, the direct correlation function $`c(r)`$, or a combination of the two, $`S`$ is a set of a metric space of functions, and $`A:SS`$ is a non linear operator mapping $`S`$ into itself.

Numerical analysis of integral equations suggests the use of the following combination

$`\gamma (r)=h(r)c(r),`$ (2)
since $`\gamma `$ is a much smoother function than $`h`$ or $`c`$, especially in the core region.

It has been pointed out by Malescio et. al. G. Malescio, P. V. Giaquinta, and Y. Rosenfeld (1998, 2000); G. Malescio and P. V. Giaquinta (2000) that, amongst the different numerical schemes that one may choose to solve (1), the simple iterative scheme of Picard plays a special role. Picard scheme consists in generating successive approximations to the solution through the relationship

$`\gamma _{n+1}=A\gamma _n,`$ (3)
starting from some initial value $`\gamma _0`$. If the sequence of successive approximations $`\{\gamma _n\}`$ converges toward a value $`\gamma ^{}`$, then $`\gamma ^{}`$ is a fixed point for the operator $`A`$, i.e. it is a solution of Eq. (1), $`\gamma ^{}=A\gamma ^{}`$. Banach’ s fixed point theorem (see chapter 1 in Zeidler (1986) especially theorem 1.A) states that, given an operator $`A:SS`$, where $`S`$ is a closed nonempty set in a complete metric space, the simple iteration (3) may converge toward the only fixed point in $`S`$ ($`A`$ is $`k`$ contractive) or it may not converge ($`A`$ is non expansive). So the simple iterative method can be used to signal a fundamental change in the properties of the underlying operator.

The operator $`A`$ will in general depend on the thermodynamic state of the fluid. In order to determine the properties of the operator at a given state we can proceed as follows. First, we find the fixed point $`\gamma ^{}`$ using a numerical scheme (more refined then Picard’ s) capable of converging in the high density region. Next, we perturb the fixed point with an arbitrary initial perturbation $`\delta _0(r)`$ so that

$`A(\gamma ^{}+\delta _0)A\gamma ^{}+{\displaystyle \frac{A}{\gamma }}|_\gamma ^{}\delta _0=\gamma ^{}+M\delta _0,`$ (4)
where we have introduced the Floquet matrix $`M`$. Now $`\delta _1=M\delta _0`$ may be considered as the new perturbation. We then generate the succession $`\{\delta _n\}`$ where

$`\delta _n=M\delta _{n1}.`$ (5)
If the succession converges to zero then the operator $`A`$ is $`k`$ contractive, if it diverges the operator is non expansive. Malescio et. al. call $`\{\delta _n\}`$ fictitious dynamics and associate to the resulting fate of the initial perturbation the nature of the structural equilibrium of the fluid. If the succession converges to zero they say that the fluid is structurally stable and structurally unstable otherwise. We will call $`\rho _{inst}`$ the density where the transition between a structurally stable and unstable fluid occurs.

Following Malescio et. al. it is possible to define a measure for the structural stability of the system as follows. We define

$`S_i={\displaystyle \frac{M\delta _i(r)}{\delta _i(r)}},`$ (6)
where $`f(r)=\sqrt{_{i=1}^Nf^2(r_i)}`$ is the norm of a function $`f`$ defined over a mesh of $`N`$ points. We assume that the norm of the perturbation depends exponentially on the number of iterations

$`\delta _n=\delta _02^{\lambda n},`$ (7)
where $`\lambda `$ is the Lyapunov exponent related to the fictitious dynamics. Then one can write the average exponential stretching of initially nearby points as

$`\lambda =\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{n}}\mathrm{log}_2\left({\displaystyle \underset{i=0}{\overset{n1}{}}}S_i\right).`$ (8)
Malescio et. al. have calculated the dependence of $`\lambda `$ on the density for various simple three dimensional liquids (and various closures): hard spheres G. Malescio, P. V. Giaquinta, and Y. Rosenfeld (1998), Yukawa, inverse power and Lennard-Jones potentials G. Malescio, P. V. Giaquinta, and Y. Rosenfeld (2000). For all these systems they found that $`\lambda `$ increases with the density and the density at which $`\lambda `$ becomes positive, $`\rho _{inst}`$, falls close to the freezing density $`\rho _f`$ of the fluid system. This occurrence lead them to propose this kind of analysis as a one-phase criterion to predict the freezing transition of a dense fluid and to estimate $`\rho _f`$. However, we think that there are some practical and conceptual difficulties with such one-phase criterion.

First of all, it does not depend only on the closure adopted but also on the kind of algorithm used to solve the integral equation. Indeed, different algorithms give different $`\rho _{inst}`$ and Malescio et. al. choose to use as instability threshold for their criterion the one obtained using Picard algorithm, thus giving to it a special status. However, it is hard to understand why the particular algorithm adopted in the solution of the integral equation should be directly related to a phase boundary.

Moreover, one would expect that the estimate of $`\rho _{inst}`$ would improve in connection with improved closures. This is not the case, at least in the one component hard sphere fluid.

Even a more serious doubt about the validity of the proposed criterion comes from its behavior in one dimensional systems. In this paper we present the same Lyapunov exponent analysis on a system of hard rods in one dimension treated using either the Percus-Yevick (PY) or the hypernetted chain (HNC) approximations. What we find is that the Lyapunov exponent as a function of density has the same behavior as that for the three dimensional system (hard spheres): it becomes positive beyond a certain $`\rho _{inst}`$. Since it is known Hove (1950) that a one dimensional fluid of hard rods does not have a phase transition, our result sheds some doubts on the validity of the proposed criterion.

## II Technical details

As numerical scheme to calculate the fixed point we used Zerah’ s algorithm Zerah (1985) for the three dimensional systems and a modified iterative method for the hard rods in one dimension. In the modified iterative method input and output are mixed at each iteration

$`\gamma _{n+1}=A_{mix}\gamma _n=\alpha A\gamma _n+(1\alpha )\gamma _n,`$ (9)
where $`\alpha `$ is a real parameter $`0<\alpha <1`$. Note that while for a non expansive operator $`A`$ the Picard iterative method (3) needs not converge, one can prove convergence results on an Hilbert space for the modified iterative method with fixed $`\alpha `$ (see proposition 10.16 in Zeidler (1986)). In all the computations we used a uniform grid of $`N=1024`$ points with a spacing $`\delta r=0.025`$. Generally, we observed a marginal increase of $`\rho _{inst}`$ by lowering $`N`$.

A method to find a Lyapunov exponent, equivalent but more accurate than the one of Malescio et. al. (8), goes through the diagonalization of the Floquet matrix. Note that in general this matrix is non symmetric, thus yielding complex eigenvalues. A Lyapunov exponent can then be defined as R. Z. Sagdeev, D. A. Usikov, and G. M. Zaslavsky (1988)

$`\lambda ^{}=\mathrm{log}\left[\underset{i}{\mathrm{max}}\left(\sqrt{er_i^2+ei_i^2}\right)\right],`$ (10)
where $`er_i`$ and $`ei_i`$ are respectively the real and imaginary part of the $`i`$-th eigenvalue. In our numerical computations we always used recipe (10) to calculate the Lyapunov exponents since it is explicitly independent from the choice of an initial perturbation.

We constructed the Floquet matrix in the following way Gillan (1979). In a Picard iteration we start from $`\gamma (r)`$ we calculate $`c(r)`$ from the closure approximation, we calculate its Fourier transform $`\stackrel{~}{c}(k)`$, we calculate $`\stackrel{~}{\gamma }(k)`$ from the OZ equation, and finally we anti transform $`\stackrel{~}{\gamma }`$ to get $`\gamma ^{}(r)`$. For example for a three dimensional system a PY iteration in discrete form can be written as follows

$`c_i`$ $`=`$ $`(1+\gamma _i)\left(e^{\beta \varphi _i}1\right),`$ (11)
$`\stackrel{~}{c}_j`$ $`=`$ $`{\displaystyle \frac{4\pi \delta r}{k_j}}{\displaystyle \underset{i=1}{\overset{N1}{}}}r_i\mathrm{sin}(k_jr_i)c_i,`$ (12)
$`\stackrel{~}{\gamma }_j`$ $`=`$ $`\rho \stackrel{~}{c}_j^2/(1\rho \stackrel{~}{c}_j)`$ (13)
$`\gamma _i^{}`$ $`=`$ $`{\displaystyle \frac{\delta k}{2\pi ^2r_i}}{\displaystyle \underset{j=1}{\overset{N1}{}}}k_j\mathrm{sin}(k_jr_i)\stackrel{~}{\gamma }_j,`$ (14)
where $`r_i=i\delta r`$ are the $`N`$ mesh points in $`r`$ space, $`k_j=j\delta k`$ are the $`N`$ mesh points in $`k`$ space, with $`\delta k=\pi /(N\delta r)`$, $`c_i=c(r_i)`$, $`\gamma _i=\gamma (r_i)`$, $`\stackrel{~}{c}_j=\stackrel{~}{c}(k_j)`$, $`\stackrel{~}{\gamma }_j=\stackrel{~}{\gamma }(k_j)`$, and $`\varphi _i=\varphi (r_i)`$ is the interparticle potential calculated on the grid points. The Floquet matrix will then be

$`M_{ij}`$ $`=`$ $`{\displaystyle \frac{\gamma _i^{}}{\gamma _j}}={\displaystyle \underset{m=1}{\overset{N1}{}}}{\displaystyle \frac{\gamma _i^{}}{\stackrel{~}{\gamma }_m}}{\displaystyle \frac{\stackrel{~}{\gamma }_m}{\stackrel{~}{c}_m}}{\displaystyle \frac{\stackrel{~}{c}_m}{c_j}}{\displaystyle \frac{c_j}{\gamma _j}}`$ (15)
$`=`$ $`{\displaystyle \frac{\delta r\delta k}{\pi }}\left({\displaystyle \frac{r_j}{r_i}}\right)\left(e^{\beta \varphi _j}1\right)(D_{ij}D_{i+j}),`$
where

$`D_l={\displaystyle \underset{m=1}{\overset{N1}{}}}\mathrm{cos}(k_mr_l)\left[{\displaystyle \frac{2\rho \stackrel{~}{c}_m}{1\rho \stackrel{~}{c}_m}}+\left({\displaystyle \frac{\rho \stackrel{~}{c}_m}{1\rho \stackrel{~}{c}_m}}\right)^2\right].`$ (16)
The HNC case can be obtained replacing in (15) $`[\mathrm{exp}(\beta \varphi _j)1]`$ with $`[\mathrm{exp}(\beta \varphi _j+\gamma _j)1]`$.

To derive the expression for the Floquet matrix valid for the one dimensional system and consistent with a trapezoidal discretization of the integrals, we need to replace (12) and (14) with

$`\stackrel{~}{c}_j`$ $`=`$ $`2\delta r\left({\displaystyle \underset{i=1}{\overset{N1}{}}}\mathrm{cos}(k_jr_i)c_i+{\displaystyle \frac{1}{2}}c_0\right),`$ (17)
$`\gamma _i^{}`$ $`=`$ $`{\displaystyle \frac{\delta k}{\pi }}\left({\displaystyle \underset{i=1}{\overset{N1}{}}}\mathrm{cos}(k_jr_i)\stackrel{~}{\gamma }_j+{\displaystyle \frac{1}{2}}\stackrel{~}{\gamma }_0\right).`$ (18)
## III Numerical results

We checked our procedure for a three dimensional hard spheres fluid and a Lennard-Jones fluid at a reduced temperature $`T^{}=2.74`$. Our results, obtained using recipe (10), were in good agreement with those of Malescio et. al. G. Malescio, P. V. Giaquinta, and Y. Rosenfeld (1998, 2000) which used recipe (8) instead (another difference between our analysis and theirs is that we used for $`\gamma `$ the indirect correlation function (2) while they were using the total correlation function $`h`$). For the Lennard-Jones fluid our results where indistinguishable from those of Malescio et. al. G. Malescio, P. V. Giaquinta, and Y. Rosenfeld (2000). We found a reduced instability density $`\rho _{inst}^{}`$ around 1.09 in the PY approximation and around 1.06 in the HNC approximation. For the three dimensional hard sphere fluid we found slightly larger ($`4\%`$) values for $`\rho _{inst}`$. We found a $`\eta _{inst}=\rho _{inst}\pi d^3/6`$ of about 0.445 in the PY approximation and around 0.461 in the HNC approximation. We also checked the value corresponding to the Martynov-Sarkisov (MS) Martynov and Sarkisov (1983) closure and we found $`\eta _{inst}=0.543`$.

We feel that the differences are within what we can expect on the basis of small numerical differences in different implementations. We think that it is more worth of notice that closures providing better structural and thermodynamic properties, like PY or MS do not provide a better value of $`\eta _{inst}`$.

We looked at the structure of the Floquet matrix too but from direct inspection we can conclude that it is not diagonally dominated.

Then, we have calculated the Lyapunov exponent (10) as a function of the density for a fluid of hard rods in one dimension using both PY and HNC closures. The results of the calculation are shown in Fig. 1 and Fig. 2. The curves show the same qualitative behavior as the ones for the three dimensional fluid. From Fig. 1 we can see how the slope of the curves starts high at low densities and decreases rapidly with $`\rho `$. At high densities the Lyapunov exponent becomes zero at $`\rho _{inst}`$. As expected, to find the fixed point for $`\rho \rho _{inst}`$ it is necessary to choose $`\alpha <1`$ in the modified iterative scheme (9). Before reaching the instability threshold the curves show a rapid change in their slope at $`\rho _c<\rho _{inst}`$. Figure 2 shows a magnification of the region around $`\rho _c`$ from which we are lead to conclude that, within the numerical accuracy of the calculations, the slope of the curves $`d\lambda ^{}/d\rho `$ undergoes a jump at $`\rho _c`$.

## IV Conclusions

The fictitious dynamics associated to the iterative solution of an integral equation can signal the transition of the map of the integral equation from $`k`$ contractive to non expansive. If the Lyapunov exponent is negative the map is $`k`$ contractive, if it is positive the map is non expansive.

Since it is possible to modify in an arbitrary way the fictitious dynamics keeping the same fixed point, it is difficult to understand a deep direct connection between the stability properties of the map and a one-phase criterion for a thermodynamic transition.

Admittedly the correlations shown by Malescio et al. are striking. We calculated the Lyapunov exponent as a function of the density for various fluids (hard spheres in one and three dimensions and three dimensional Lennard-Jones fluid) both in the HNC and PY approximations. For the three dimensional fluids the instability density falls close to the freezing density $`\rho _f`$. For example, the Lennard-Jones fluid studied with HNC should undergo a freezing transition at $`\rho ^{}1.06`$ or at $`\rho ^{}1.09`$, if studied with PY , rather close to the freezing density $`\rho _f^{}1.113`$. For hard spheres $`\rho _{inst}^{}`$ is about $`10\%`$ smaller than $`\rho _f^{}0.948`$. The Hansen-Verlet “rule” states that a simple fluid freezes when the maximum of the structure factor is about 2.85 Hansen and Verlet (1969). According to this rule the three dimensional hard spheres fluid studied with HNC should undergo a freezing transition at $`\rho 1.01`$ while when studied with PY the transition should be at $`\rho 0.936`$. The corresponding estimates obtained through $`\rho _{inst}^{}`$, $`0.879`$ (HNC) and $`0.850`$ (PY) are poorer and, more important, are not consistent with the well known better performance of PY in the case of hard spheres.

In one dimension, a fluid of hard spheres (hard rods), cannot undergo a phase transition Hove (1950). From Fig. 1 we see that the system still becomes structurally unstable. This can be explained by observing that the structural stability as defined by Malescio et. al. is a property of the map $`A`$ and in particular of the algorithm used to get solution of the integral equation under study. As such, it is not directly related to the thermodynamic properties even at the approximate level of the theory (there is no direct relation between the contractiveness properties of $`A`$ and the thermodynamics). It looks more reasonable that the increase of the correlations would be the common origin of the numerical instability of the Picard iteration and, whenever it is possible, of thermodynamic phase transitions.

###### Acknowledgements.

G.P. would like to acknowledge preliminary exploratory work on this subject carried on in collaboration with Matteo Mosangini and Waheed Adeniyi Adeagbo.

LIST OF FIGURES

* For a fluid of hard rods in one dimension, we show the Lyapunov exponent as a function of the reduced density ($`\rho ^{}=\rho \sigma `$ where $`\sigma `$ is the rods width) as calculated using the PY and the HNC closures.
* We show a magnification of Fig. 1 in a neighborhood of the instability threshold.


format: text

text:
# Hadronization in heavy ion collisions: recombination or fragmentation?

## Abstract

We show that hadron production in relativistic heavy ion collisions at transverse momenta larger than 2 GeV/$`c`$ can be explained by the competition of two different hadronization mechanisms. Above 5 GeV/$`c`$ hadron production can be described by fragmentation of partons that are created perturbatively. Below 5 GeV/$`c`$ recombination of partons from the dense and hot fireball dominates. This can explain some of the surprising features of RHIC data like the constant baryon-to-meson ratio of about one and the small nuclear suppression for baryons between 2 to 4 GeV/$`c`$.

: JPG

The Relativistic Heavy Ion collider (RHIC) has provided exciting data about hadron production at transverse momenta of a few GeV/$`c`$ in central Au+Au collisions. The production of pions at high $`P_T`$ was found to be suppressed compared to the scaled yield from $`p+p`$ collisions . This jet quenching effect can be understood by final state interaction of fast partons with the dense and hot medium produced in central heavy ion collisions. Fast partons lose energy via induced bremsstrahlung before they can fragment into high $`P_T`$ hadrons . The suppression effect is dramatic and can be as large as a factor of 6 above $`P_T=5`$ GeV/$`c`$.

On the other hand the suppression of protons and antiprotons seems to be much less . Experimental data from PHENIX show a proton/pion ratio of about 1 between 1.5 GeV/$`c`$ and 4 GeV/$`c`$ . This is surprising since the production of protons and antiprotons is usually suppressed compared to the production of pions because of the much larger mass. At high transverse momentum this can be understood in terms of perturbative quantum chromodynamics (pQCD) . The fragmentation functions $`D_{ah}(z)`$ describe the probability for a parton $`a`$ with momentum $`p`$ to turn into a hadron with momentum $`zp`$, $`0<z<1`$. These fragmentation functions were measured for pions and protons, mainly in $`e^+e^{}`$ collisions, and give a ratio $`p/\pi ^0<0.2`$ for $`P_T>2`$ GeV/$`c`$ when used in $`p+p`$ and $`N+N`$ collisions. The energy loss of partons in a medium can be taken into account by a rescaling of the parton momentum . However this should affect baryons and mesons in the same way.

The lack of nuclear suppression for baryons is a challenge for our understanding of hadron production. The currently accepted picture assumes that a parton with large transverse momentum is produced in a hard scattering reaction between initial partons, propagates through the surrounding hot matter and loses energy by interactions, and finally hadronizes. Apparently the creation and interaction of a parton will happen independently of its later fate during hadronization. Therefore any unusual behavior that is different between baryons and mesons can only be attributed to hadronization itself. We propose to use recombination of quarks from the surface of the hot fireball as an alternative hadronization mechanism .

In the fragmentation process a parton with transverse momentum $`p_T`$ is leaving the interaction zone while still being connected with other partons by a color string. The breaking of the string creates quark antiquark pairs which finally turn into hadrons. The distribution of one of these hadrons, which is bound to have less transverse momentum $`P_T=zp_T`$, is then described by a fragmentation function. The average value $`z`$ is about 0.5 for pions in $`p+p`$ collisions. In other words the production of a, say, 5 GeV/$`c`$ pion has to start with a 10 GeV/$`c`$ parton in average, which are rare to find due to the steeply falling parton spectrum. Jet quenching even enhances the lack of high $`p_T`$ partons. On the other hand, the 5 GeV/$`c`$ pion could be produced by the recombination of a quark and an antiquark with about 2.5 GeV/$`c`$ each in average. 2.5 GeV/$`c`$ and 10 GeV/$`c`$ are separated by orders of magnitude in the parton spectrum. The price to pay is of course that two of these partons have to be found close to each other in phase space. However we do have a densely populated phase space in central heavy ion collisions at RHIC where we even expect the existence of a thermalized quark gluon plasma.

Recombination of quarks has been considered before in hadron collisions  and was also applied to heavy ion collisions . In QCD the leading particle effect in the forward region of a hadron collision is well known. This is the phenomenon that the production of hadrons that share valence partons with the beam hadrons are favored in forward direction. It has been realized that this can only be explained by recombination. This has fueled a series of theoretical work, see e.g.  and references therein. Recently the recombination idea for heavy ion collisions, stimulated by the RHIC results, has been revived for elliptic flow  and hadron spectra and ratios .

The formalism of recombination has already been developed in a covariant setup for the process of baryons coalescing into light nuclei and clusters in nuclear collisions . We give a brief derivation for the case of mesons. By introducing the density matrix $`\widehat{\rho }`$ for the system of partons, the number of quark-antiquark states that will be measured as mesons is given by

$$N_M=\underset{ab}{}\frac{d^3P}{(2\pi )^3}M;𝐏|\widehat{\rho }_{ab}|M;𝐏$$
(1)
Here $`|M;𝐏`$ is a meson state with momentum $`𝐏`$ and the sum is over all combinations of quantum numbers — flavor, helicity and color — of valence partons that contribute to the given meson $`M`$. This can be cast in covariant form using a hypersurface $`\mathrm{\Sigma }`$ for hadronization

$$E\frac{dN_M}{d^3P}=C_M\underset{\mathrm{\Sigma }}{}\frac{d\sigma Pu(\sigma )}{(2\pi )^3}\frac{d^3q}{(2\pi )^3}w_a(\sigma ;\frac{𝐏}{2}𝐪)\mathrm{\Phi }_M^W(𝐪)w_b(\sigma ;\frac{𝐏}{2}+𝐪).$$
(2)
Here $`E`$ is the energy of the four vector $`P`$, $`d\sigma `$ a measure on $`\mathrm{\Sigma }`$ and $`u(\sigma )`$ is the future oriented unit vector orthogonal to the hypersurface $`\mathrm{\Sigma }`$. $`w_a`$ and $`w_b`$ are the phase space densities for the two partons $`a`$ and $`b`$, $`C_M`$ is a degeneracy factor and

$$\mathrm{\Phi }_M^W(𝐪)=d^3r\mathrm{\Phi }_M^W(𝐫,𝐪)$$
(3)
where $`\mathrm{\Phi }_M^W(𝐫,𝐪)`$ is the Wigner function of the meson .

Since the hadron structure is best known in a light cone frame, we write the integral over $`𝐪`$ in terms of light cone coordinates in a frame where the hadron has no transverse momentum but a large light cone component $`P^+`$. This can be achieved by a simple rotation from the lab frame. Introducing the momentum $`k=P/2q`$ of parton $`a`$ in this frame we have $`k^+=xP^+`$ with $`0<x<1`$. We make an ansatz for the Wigner function of the meson in terms of light cone wave functions $`\varphi _M(x)`$. The final result can be written as

$$E\frac{N_M}{d^3P}=C_M\underset{\mathrm{\Sigma }}{}𝑑\sigma \frac{Pu(\sigma )}{(2\pi )^3}\underset{0}{\overset{1}{}}𝑑xw_a(\sigma ;xP^+)\left|\varphi _M(x)\right|^2w_b(\sigma ;(1x)P^+).$$
(4)
For a baryon with valence partons $`a`$, $`b`$ and $`c`$ we obtain

$`E{\displaystyle \frac{N_B}{d^3P}}`$ $`=C_B{\displaystyle \underset{\mathrm{\Sigma }}{}}𝑑\sigma {\displaystyle \frac{Pu(\sigma )}{(2\pi )^3}}{\displaystyle \underset{0}{\overset{1}{}}}𝑑x_1𝑑x_2𝑑x_3\delta (x_1+x_2+x_31)`$
$`\times w_a(\sigma ;x_1P^+)w_b(\sigma ;x_2P^+)w_c(\sigma ;x_3P^+)\left|\varphi _B(x_1,x_2,x_3)\right|^2.`$
$`\varphi _B(x_1,x_2,x_3)`$ is the effective wave function of the baryon in light cone coordinates.

A priori these wave functions are not equal to the light cone wave functions used in exclusive processes. We are recombining effective quarks in a thermal medium and not perturbative partons in an exclusive process. Nevertheless, as an ansatz for a realistic wave function one can adopt the asymptotic form of the light cone distribution amplitudes

$`\varphi _M(x)`$ $`=\sqrt{30}x(1x),`$ (6)
$`\varphi _B(x_1,x_2,x_3)`$ $`=12\sqrt{35}x_1x_2x_3`$
as a model. However, it turns out that for a purely exponential spectrum the shape of the wave function does not matter. In that case the dependence on $`x`$ drops out of the product of phase space densities

$$w_a(\sigma ;xP^+)w_b(\sigma ;(1x)P^+)e^{xP^+/T}e^{(1x)P^+/T}=e^{P^+/T}.$$
(7)
One can show that a narrow width approximation, using $`\delta `$ peaked wave functions that distribute the momentum of the hadron equally among the valence quarks — 1/2 in the case of a meson and 1/3 in the case of a baryon — deviates by less then 20% from a calculation using the wave functions in (6). That there is a small deviation can be attributed to the violation of energy conservation. Since recombination is a $`21`$ or $`31`$ process, energy will generally not be conserved if we enforce momentum conservation and a mass shell condition for all particles. However, one can show that violations of energy conservation are of the order of the effective quark masses or $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ and can therefore be neglected for transverse hadron momenta larger than 2 GeV/$`c`$.

Fragmentation of partons is given by

$$E\frac{d\sigma _h}{d^3P}=\underset{a}{}\underset{0}{\overset{1}{}}\frac{dz}{z^2}D_{ah}(z)E_a\frac{d\sigma _a}{d^3P_a}.$$
(8)
The sum runs over all parton species $`a`$ and $`\sigma _a`$ is the cross section for the production of parton $`a`$ with momentum $`P_a=P/z`$. We use a leading order (LO) pQCD calculation of $`\sigma _a`$  together with LO KKP fragmentation functions . Energy loss of the partons is taken into account by a shift of the parton spectrum by

$$\mathrm{\Delta }p_T=\sqrt{\lambda p_T}.$$
(9)
We summarize that the transverse momentum dependent yield of mesons from recombination can be written as $`C_Mw^2(P_T/2)`$ in the simple narrow width approximation, whereas from fragmentation we expect $`D(z)w(P_T/z)`$. For an exponential parton spectrum $`w=e^{P_T/T}`$ the ratio of recombination to fragmentation is

$$\frac{R}{F}=\frac{C_M}{D}e^{\frac{P_T}{T}\left(1\frac{1}{z}\right)}$$
(10)
where $`D<1`$ and $`z<1`$ are average values of the fragmentation function and the scaling variable. Therefore $`R/F>1`$. In fact, recombination always wins over fragmentation from an exponential spectrum (as long as the exponential is not suppressed by small fugacity factors). The same is true in the case of baryons.

Now let us consider a power law spectrum $`w=A(P_T/\mu )^\alpha `$ with a scale $`\mu `$ and $`\alpha >0`$. Then the ratio of recombination over fragmentation is

$$\frac{R}{F}=\frac{C_MA}{D}\left(\frac{4}{z}\right)^\alpha \left(\frac{P_T}{\mu }\right)^\alpha $$
(11)
and fragmentation ultimately has to win at high $`P_T`$. We also note that we can expect a constant baryon/meson ratio from recombination, when the parton spectrum is exponential. In this case the ratio is only determined by the degeneracy factors

$$\frac{dN_B}{dN_M}=\frac{C_B}{C_M}.$$
(12)
For our numerical studies we consider fragmentation of perturbative partons and recombination from a thermal phase

$$w_{\mathrm{th}}=e^{p_T\mathrm{cosh}(\eta y)/Te^{y^2/2\mathrm{\Delta }^2}}$$
(13)
with an effective temperature $`T`$. $`\eta `$ is the space-time rapidity and $`\mathrm{\Delta }2`$ the width of the rapidity distribution. We fix the hadronization hypersurface $`\mathrm{\Sigma }`$ by the condition $`\tau _f=\sqrt{t^2z^2}=\mathrm{const}.`$ . We set $`\tau _f=5`$ fm/$`c`$. A two phase parton spectrum with a perturbative power law tail and an exponential part at low transverse momentum is also predicted by parton cascades like VNI/BMS .

The parameters, $`\lambda `$ for the average energy loss and $`T`$ for the slope of the exponential parton spectrum, are determined by a fit to the inclusive charged hadron spectrum measured by PHENIX . We obtain $`\lambda 1`$ GeV and $`T350`$ MeV. The temperature contains the effect of a blue-shift because of the strong radial flow. An additional normalization factor of about $`1/30`$ for the recombination part is necessary to describe the data. This is due to the use of an effective temperature which gives a too large particle number compared to the physical temperature that we expect to be around 175 MeV.

In Figure 1 we show the inclusive charged hadron spectrum using recombination and fragmentation for pions, kaons, protons and antiprotons. We note that the crossover between the recombination dominated and fragmentation dominated regimes is around 5 GeV/$`c`$. It is in fact earlier for pions ($``$4 GeV/$`c`$) than for protons ($``$6 GeV/$`c`$). In Figure 1 we also give the spectrum of up quarks as an example for the partonic input of our calculation. We note that the crossover between the perturbative and the thermal domain is here around 3 to 3.5 GeV/$`c`$. It is characteristic for recombination that features of the parton spectrum are pushed to higher transverse momentum in the hadron spectrum.

In Figure 2 we give the proton/pion ratio and the nuclear modification factor $`R_{AA}`$ for pions, protons and charged hadrons. The proton/pion ratio shows a plateau between 2 and 4 GeV/$`c`$ in accordance with experiment and a steep decrease beyond that. This decrease has not yet been seen and is a prediction of our work.

In the quantity $`R_{AA}`$ the effect of jet quenching is manifest for all hadrons at large transverse momentum. For pions $`R_{AA}`$ grows only moderately at low $`P_T`$ in accordance with PHENIX data . On the other hand recombination is much more important for protons, because they are suppressed in the fragmentation process. Therefore $`R_{AA}`$ nearly reaches a value of 1 below 4 GeV/$`c`$ for protons.

In summary, we have discussed recombination as a possible hadronization mechanism in heavy ion collisions. We have shown that recombination can dominate over fragmentation up to transverse momenta of 5 GeV/$`c`$. Recombination provides a natural explanation for the baryon/meson ratio and the nuclear suppression factors observed at RHIC.

Acknowledgments. This work was supported by RIKEN, Brookhaven National Laboratory, DOE grants DE-FG02-96ER40945 and DE-AC02-98CH10886, and by the Alexander von Humboldt Foundation.


text:
# Chapter 1 Limitations on Space-based Air Fluorescence Detector Apertures obtained from IR Cloud Measurements

## Chapter 1 Limitations on Space-based Air Fluorescence Detector Apertures obtained from IR Cloud Measurements

### Abstract

The presence of clouds between an airshower and a space-based detector can dramatically alter the measured signal characteristics due to absorption and scattering of the photonic signals. Furthermore, knowledge of the cloud cover in the observed atmosphere is needed to determine the instantaneous aperture of such a detector. Before exploring the complex nature of cloud-airshower interactions, we examine a simpler issue. We investigate the fraction of ultra-high energy cosmic ray events that may be expected to occur in volumes of the viewed atmosphere non-obscured by clouds. To this end, we use space-based IR data in concert with Monte Carlo simulated $`10^{20}`$ eV airshowers to determine the acceptable event fractions. Earth-observing instruments, such as MODIS, measure detailed cloud configurations via a CO2-slicing technique that can be used to determine cloud-top altitudes over large areas. Thus, events can be accepted if their observed 3-dimensional endpoints occur above low clouds as well as from areas of cloud-free atmosphere. An initial analysis has determined that by accepting airshowers that occur above low clouds, the non-obscured acceptance can be increased by approximately a factor of 3 over that obtained using a cloud-free criterion.

### 1. Introduction

Using the air fluorescence technique, the space-based experiments EUSO  and OWL  will image the UV nitrogen fluorescence signal from extended airshowers with a spatial resolution of $`1`$ km2 (on the ground) while instantaneously monitoring areas approaching $`10^6`$ km2. The effects of clouds must be well understood as clouds can obscure or modify the signal strengths and profiles of the airshowers. The incorporation of a LIDAR system has been proposed to be used in conjunction with the EUSO and OWL detectors, but the scanning requirements imposed by the nearly mega-pixel count and the rapidly moving footprint(s) lead to severe operational requirements. An alternate approach to determine the cloud properties is the use of meteorological measurements. The MODIS instruments  on the Terra and Aqua satellites provide derived cloud property measurements from IR measurements with a spatial resolution of $`1`$ km2 on the ground. Using the technique of CO2 slicing , cloud-top altitudes can be determined and lead to 3-dimensional cloud profiles. The availability of these and future meteorological measurements could be used by space-based air fluorescence experiments to define clear apertures and relax LIDAR operational constraints.

In a previous paper , we used simulated airshowers superimposed over MODIS measured scenes and obtained the fraction of tracks that occur in cloud-free portions of the atmosphere. This analysis used a cloud-mask product that provides a probabilistic determination of whether a particular 1 km2 pixel contains a cloud. Using a conservative selection criteria defined by requiring no cloudy pixels within 3 km of the observed portion of a simulated airshower, we determined that the mean of the distribution was 6.5%. Assuming a 15% duty cycle imposed by moonless night along with other operational constraints, a 6.5% clear fraction leads to a disconcerting 1% effective duty cycle.

In this paper, we present the first results on the potential of recovering airshowers with observed endpoints that occur above low clouds and thus are unobscured. This is accomplished by employing another MODIS data product which yields cloud-top pressures, i.e. cloud-top altitudes, but on a courser $`5\times 5`$ km2 pixel scale. The unobscured event fraction is then determined by combining the cloud-free fraction with the low-cloud fraction.

### 2. Method of Analysis

A simulated, airshower event sample was generated using the OWL Monte Carlo  assuming an isotropic flux of $`10^{20}`$ eV protons. The two OWL instruments  were configured to be in 1000 km orbits with a 500 km satellite separation. The 3-dimensional airshower tracks were reconstructed using the stereo reconstruction technique with the track endpoints determined by requiring observation by both instruments. The approximately 1700 tracks were then superimposed upon 78 different MODIS scenes defined by the MODIS Collection4 cloud-top product. Each pixel was coded with either being clear, having a cloud top with an altitude of 3 km or less, or a cloud top with an altitude higher than 3 km. The spatial location of the xy-projected track lengths were then compared to the MODIS cloud pixel information to determine whether the tracks occurred in a clear, low-cloud, or high-cloud portion of the MODIS scene. The extended nature of the airshowers will lead to a reduced acceptance as compared to that expected by determining the cloud-free or low-cloud pixel fractions.

### 3. Results

Figure 1 shows the distribution of the cloud-free track fraction using the 1 km2 cloud-mask product for the 78 MODIS scenes located near the equator used in this study. The analysis follows that employed in our previous paper : a track is considered in a cloud-free portion of the atmosphere if no cloudy pixels are within 3 km of the observed track. The results presented in Figure 1 employed a newer version, Collection4, of the cloud-mask determining algorithm. The mean of 7% matches well to the 6.5% we had obtained for a slightly larger MODIS data set by using a previous version of the MODIS cloud-mask algorithm. The 78 scenes used in this study are a subset of the 85 scenes used in the original study. Figure 2 shows the distribution of cloud-free track fractions for the 78 MODIS scenes using the MODIS cloud-top data product. The cloud-free fractions were determined by requiring no cloudy pixels to be within 3.5 km of the track. The slightly larger value of 3.5 km is imposed by the courser $`5\times 5`$ km2 pixel size used in the cloud-top product. Note that the cloud-free fraction of 10.7% obtained by this analysis is in relatively good agreement with the 7% obtained from the 1 km2 cloud-mask analysis.

Figure 3 shows the distribution of the fraction of tracks that occur over either cloud-free areas or have a portion of the xy-projected track over clouds with cloud-top heights of less than 3 km. A track over low clouds was accepted if it was at least 3.5 km from a pixel with a $`>3`$ km altitude cloud and has an observed track endpoint at an altitude 4 km or higher. The simulated event sample has 53% of the tracks with observed endpoints above 4 km in altitude. The distribution in Figure 3 exhibits a mean of over 30% and is flatter than that for the cloud-free fraction (Figure 2). The 30% value is a factor of 3 larger than that for the cloud-free fraction. If the track endpoint requirement is increased to 5 km for low-cloud acceptance, the mean of the combined cloud-free and low-cloud distribution is reduced to 26%.

### 4. Discussion

The ability to determine the altitudes of airshowers and cloud tops allows for a significant increase in the acceptance of events by including tracks whose observable endpoints occur above low clouds. The potential for the increase is echoed in the 16% cloud-free versus 53% cloud-free+low-cloud pixel fractions obtained from the MODIS cloud-top data used in this study.

We would like to thank Bill Ridgway for invaluable discussions and for the generation of the MODIS data files.

### 5. References

1. http://www.euso-mission.org/

2. Krizmanic, J. for the OWL Collaboration 2001, Proc. 27th ICRC, Vol. 2, 861

3. http://modis-atmos.gsfc.nasa.gov/

4. http://owl.gsfc.nasa.gov/

5. Sokolsky P., Krizmanic J. 2003, astro-ph/0302501, submitted to Astroparticle Physics

6. Streitmatter R. et al. 2002, http://owl.gsfc.nasa.gov/WP10c.PDF

7. Stumpf R., Pennock J. 1989, Jour. Geophys. Res. 94(C10), 14,363


text:
# Fluctuating Dimension in a Discrete Model for Quantum Gravity Based on the Spectral Principle

## Abstract

The spectral principle of Connes and Chamseddine is used as a starting point to define a discrete model for Euclidean quantum gravity. Instead of summing over ordinary geometries, we consider the sum over generalized geometries where topology, metric and dimension can fluctuate. The model describes the geometry of spaces with a countable number $`n`$ of points, and is related to the Gaussian unitary ensemble of Hermitian matrices. We show that this simple model has two phases. The expectation value $`n`$, the average number of points in the universe, is finite in one phase and diverges in the other. We compute the critical point as well as the critical exponent of $`n`$. Moreover, the space-time dimension $`\delta `$ is a dynamical observable in our model, and plays the role of an order parameter. The computation of $`\delta `$ is discussed and an upper bound is found, $`\delta <\mathrm{\hspace{0.17em}2}`$.

The quest for a quantum theory of gravity is one of the major goals in contemporary research. Any attempt in this direction involves the understanding of the space-time structure at very short distance. It is generally believed that at the Planck scale space-time may not be described by a manifold Kempf . The conventional geometrical setting of general relativity seems to be inadequate to describe the non-manifold micro-structures of space-time. Since the manifold structure has to appear at some macroscopic limit, it is natural to expect that one needs a generalization of ordinary geometry, such as noncommutative geometry (NCG), rather than a completely new formalism.

The basic idea coming from NCG NCG  is that one can describe a Riemannian manifold $`(M,g_{\mu \nu })`$ in a purely algebraic way. There is no loss of information if, instead of the data $`(M,g_{\mu \nu })`$, one is given a triple $`(𝒜,,D)`$, where $`𝒜`$ is the C\*-algebra $`C^0(M)`$ of smooth functions on $`M`$, $``$ is the Hilbert space of $`L^2`$-spinors on $`M`$, and $`𝒟`$ is the Dirac operator acting on $``$. From the Gelfand-Naimark theorem it is known that the topological space $`M`$ can be reconstructed from the set $`\widehat{𝒜}`$ of irreducible representations of $`C^0(M)`$. Metric is also encoded, and the geodesic distance can be computed from $`𝒟`$. Here we will consider only commutative spectral triples - this is enough to go much beyond ordinary geometry. In particular one can treat all Hausdorff topological spaces in this way. Given a pair $`(M,g_{\mu \nu })`$, one can promptly construct the corresponding triple $`(C^0(M),L^2(M),𝒟)`$. However, not all commutative spectral triples, or generalized geometries, come from a pair $`(M,g_{\mu \nu })`$. Nevertheless one can always associate a Hausdorff space $`M=\widehat{𝒜}`$ to a commutative spectral triple, where $`\widehat{𝒜}`$ denotes the set of irreducible representations of $`𝒜`$. However, the space $`M`$ may not be a manifold. Once we trade the original Riemannian geometry for its corresponding commutative triple we need a replacement for the Einstein-Hilbert action $`S_{EH}`$. The so-called spectral action of Chamseddine and Connes CC  is one possible candidate. It depends only on the eigenvalues of $`𝒟`$ and contains $`S_{EH}`$ as a dominant term. In this paper however we shall use another spectral action.

The spectral action can be written for any triple, regardless of whether it comes from a manifold $`(M,g_{\mu \nu })`$ or not. In the spectral geometry approach it is conceivable to write the partition function

$$Z=\underset{x𝒳}{}e^{S[x]},$$
(1)
were the “sum” is over the set $`𝒳`$ of all possible commutative spectral triples and $`S`$ depends on the spectrum of $`𝒟`$. It includes all Hausdorff spaces and therefore all manifolds of all dimensions.

The framework of spectral NCG is appropriate to explore some difficult questions in quantum gravity. For example, Rovelli used a simple dynamical model based on a finite dimensional spectral triple to construct a quantum theory with a quantized physical distance between points Rovelli . In our model we look at the dimension as a dynamical quantity and study its expectation value.

The exact computation of (1) is a formidable task, not yet accomplished. However, the algebraic approach suggests ways of defining discrete approximations to the full theory. For instance one may replace the algebra $`𝒜`$ by a finite dimensional algebra $`A_n`$. In this approach to discretization there is no need to introduce a lattice or simplicial decomposition of the underlying space. The approximation of $`𝒜`$ by a finite algebra works even if the spectral triple does not come from a manifold. Thus, it gives us a generalization of ordinary discretizations discrete .

In this paper we discretize (1) by sampling the set $`𝒳`$ with finite commutative spectral triples. We will think of it as a useful toy model, which we believe captures some of the main features of the full one, Eq. (1).

The key role played by the eigenvalues of the Dirac (or Laplace) operator in the spectral action approach was emphasized in Landi . In our model they are also the natural dynamical variables due to the connection with random matrix theory (RMT). An important ingredient of the model is that the number of points can fluctuate. Moreover, in our simple model the space-time dimension is a dynamical observable and its expectation value can be computable from first principles.

Let us describe the ensemble $`X𝒳`$ of geometries we will consider. A point of $`xX`$ is a commutative spectral triple $`x=(𝒜,,D)`$ where the commutative C\*-algebra $`𝒜`$ has a countable spectrum $`\widehat{𝒜}`$. We divide $`X`$ into subspaces $`X_n`$ consisting of triples $`(𝒜_n,_n,D)`$ such that $`\widehat{𝒜}_n`$ has a fixed number $`n`$ of points. From the Gelfand-Naimark theorem it follows that elements of $`𝒜_n`$ are the (possibly infinite) sequences $`a=(a_1,a_2,\mathrm{},a_n)`$, $`a_j`$. The Hilbert space $`_n`$ is given by vectors $`v=(v_1,\mathrm{},v_n)`$ with norm $`v^2_{i=1}^nv_i^2<\mathrm{}`$. The elements of $`𝒜`$ are represented by diagonal matrices $`\widehat{a}=\text{diag}(a_1,\mathrm{},a_n)`$ acting on $`_n`$. Finally, the operator $`D`$ is a $`n\times n`$ self-adjoint matrix. We will sample the space $`X`$ by $`X_1,X_2,\mathrm{},X_N`$ and eventually take the limit $`N\mathrm{}`$.

Let $`L`$ be a length scale such that the operator $`𝒟`$ given by $`𝒟=D/L`$ will be the analogue of the Dirac operator. The Chamseddine–Connes action depends on a cutoff function of the eigenvalues of $`D/L`$. The cutoff function is zero for eigenvalues of $`D`$ greater than $`L`$ and one otherwise CC ; Landi . In other words, the Boltzmann weight in Eq.(1) would be one outside a compact region in the eigenvalue space, leading to a divergent partition function (see Eq. (11). Let us consider a quadratic action instead:

$$S[x]=\mathrm{Tr}\left(\frac{𝒟}{\mathrm{\Lambda }}\right)^2\beta \mathrm{Tr}(D^2),$$
(2)
where $`\mathrm{\Lambda }`$ is the inverse of Planck’s length $`l_p`$, and $`\beta =(l_p/L)^2`$. Finally, we define the partition function $`Z_N(\beta )=_{n=0}^Nz_n(\beta )`$ where

$$z_n(\beta )=[dD]e^{\beta \mathrm{Tr}(D^2)}$$
(3)
is the partition function restricted to $`X_n`$, in other words, an integral over all independent matrix elements $`D_{ij}`$, where $`[dD]`$ is the usual measure for $`n\times n`$ Hermitian matrices Meh . The partition function $`z_n(\beta )`$ defines the one-matrix Gaussian Unitary Ensemble (GUE) Meh . A straightforward computation gives $`z_n(\beta )=2^{\frac{n}{2}}(\frac{\pi }{2\beta })^{\frac{n^2}{2}}.`$

The expectation values of an observable $`𝒪(D_{ij})`$ restricted to $`X_n`$ and for the entire ensemble are

$`𝒪_{n,\beta }`$ $``$ $`{\displaystyle [dD]𝒪e^{\beta \mathrm{Tr}(D^2)}/z_n(\beta )\text{and}}`$ (4)
$`𝒪(\beta )`$ $``$ $`{\displaystyle \underset{n=1}{\overset{N}{}}}P(n,\beta )𝒪_{n,\beta },`$ (5)
respectively, where the function $`P(n,\beta )=\frac{z_n(\beta )}{_nz_n(\beta )}`$ is interpreted as the probability of having a universe with $`n`$ points. The simplest observable in our model is $`n`$, the number of points in $`\widehat{𝒜}`$. By definition, $`n`$ is constant in $`X_n`$, therefore $`n_{n,\beta }=n`$. Thus we get

$$n(\beta )=\frac{_nn2^{\frac{n}{2}}(\frac{\pi }{2\beta })^{\frac{n^2}{2}}}{_n2^{\frac{n}{2}}(\frac{\pi }{2\beta })^{\frac{n^2}{2}}}.$$
(6)
The mean $`n`$ (“average number of points in the universe”) is not a continuous function of $`\beta `$ at $`\beta _c=\pi /2`$, signaling the onset of a phase transition. Besides straitforward numerical calculation, there are other ways to show that the sum (6) converges for $`\beta >\beta _c`$ and diverges for $`\beta <\beta _c`$. Consider, for instance, the approximation of $`Z_N(\beta )`$ in the limit $`N\mathrm{}`$ by the first term of a Euler-Maclaurin expansion,

$$Z(\beta )=\underset{N\mathrm{}}{lim}\underset{n=0}{\overset{N}{}}z_n(\beta )=_0^{\mathrm{}}𝑑x\mathrm{e}^{\frac{1}{2}a_\beta x^2+bx}+(\beta ),$$
(7)
where $`a_\beta =\mathrm{ln}(\beta /\beta _c)`$, $`b=\mathrm{ln}2/2`$. It is easily seen that

$$n(\beta )=\frac{}{b}\mathrm{ln}Z(\beta ),$$
(8)
neglecting the remainder $``$ in (7) footnote . Eq.(7) suggests a nice interpretation. For $`\beta >\beta _c`$ ($`\beta <\beta _c`$) the quadratic term is positive (negative) and the integral converges (diverges). The phase transition at $`\beta =\beta _c`$ is triggered by the change of signal of a bilinear term in the “fields”. The integral in (7) can be solved in the region $`\beta >\beta _c`$. After expanding $`a_\beta `$ around the critical point, $`a_\beta =(\beta \beta _c)/\beta _c+\mathrm{}`$, we compute $`n`$ by means of (8). The result is $`n(\beta )\left(\frac{\beta \beta _c}{\beta _c}\right)^1\tau ^1.`$ Hence, for $`\beta >\beta _c`$ the system is in a “finite” phase, characterized by a finite value of $`n`$. As $`\beta \beta _c^+`$ $`n`$ diverges with a (mean field) critical exponent $`\nu =1`$. The r.m.s. deviation of $`n`$ may be computed, with the result $`\mathrm{\Delta }n\sqrt{n^2n^2}\tau ^{1/2}`$. However, the relative width of the distribution, $`\mathrm{\Delta }n/n`$, decreases like $`\tau ^{1/2}`$.

For $`\beta <\beta _c`$ the relevant universes have $`n=\mathrm{}`$ and $`\mathrm{\Delta }n/n=0`$. For a $`\mathrm{}`$-dimensional $`D`$ one can define the dimension $`\delta `$ of the space $`\widehat{𝒜}`$ from the eigenvalues of $`D`$. Let $`\{\mu _0(D),\mu _1(D),\mathrm{}\}`$ be the modules of the eigenvalues (i.e. the singular values) of $`D`$ organized in an increasing order. By the Weyl formula NCG , the dimension $`\delta `$ is related to the asymptotic behavior of the eigenvalues for large $`k`$: $`\mu _k(D)k^{\frac{1}{\delta }}`$. By definition $`\delta =0`$ for finite dimensional spectral triples. We can argue that $`\delta `$ is of the form

$$\delta (\beta )=\{\begin{array}{cc}f(\beta )\hfill & \text{if }\beta <\beta _c,\hfill \\ 0\hfill & \text{if }\beta >\beta _c.\hfill \end{array}$$
(9)
This follows from the fact that for $`\beta >\beta _c`$ the probability $`P(n,\beta )`$ is localized around some finite $`n`$. Hence $`\delta `$ works as an order parameter. The value $`\beta _c=\pi /2`$ separates $`\delta =0`$ from the rest.

In order to study the dimension we need to consider the spectral $`\zeta `$-function

$$\zeta (z)=\underset{n\mathrm{}}{lim}\underset{k=0}{\overset{n}{}}\mu _k^z=\mathrm{Tr}(|D|^z),$$
(10)
where $`D`$ is an $`\mathrm{}`$-dimensional matrix ($`\mu _0>0`$). The relation between the dimension and $`\zeta (z)`$ has been discussed in Connes . For large enough values of $`\alpha =\mathrm{Re}(z)`$, $`\mathrm{Tr}(|D|^z)`$ is well defined. One says that $`D`$ has dimension spectrum $`Sd`$ if a discrete subset $`Sd=\{s_1,s_2,\mathrm{}\}`$ exists, such that $`\zeta (z)`$ can be holomorphically extended to $`/Sd`$. This definition is consistent with the Weyl formula. The set $`Sd`$ has more than a single point when for example the geometry is the union of pieces of different dimensions Connes . In what follows we will look at an upper bound for the dimension: It may happens that $`\mathrm{Tr}(|D|^\alpha )=0`$ for large enough $`\alpha `$, whereas for small values of $`\alpha `$, $`\mathrm{Tr}(|D|^\alpha )=\mathrm{}`$. Eventually, there is a value of $`\alpha `$ (say, $`\alpha _c`$) for which $`\mathrm{Tr}(|D|^{\alpha _c})`$ is finite and non-zero. The upper bound for the dimension will be $`\delta =\alpha _c`$.

In order to estimate $`\delta `$ by means of (10), we rewrite (3) and (4) as integrals over the eigenvalues $`\lambda _k`$ of $`D`$. The procedure is well-known Meh , and leads to ($`C_n\pi ^{\frac{n(n1)}{2}}/_{k=1}^nk!`$)

$$z_n(\beta )=C_n_{\mathrm{}}^{\mathrm{}}[d^n\lambda ]\mathrm{\Delta }^2(\lambda _k)e^{\beta _{i=1}^n\lambda _i^2}C_n\mathrm{\Psi }_{n,\beta },$$
(11)
$$𝒪(\lambda _i)_{n,\beta }=_{\mathrm{}}^{\mathrm{}}[d^n\lambda ]𝒪(\lambda _i)\left\{\frac{2^{\frac{n(n1)}{2}}\beta ^{\frac{n^2}{2}}}{\pi ^{\frac{n}{2}}_{k=1}^nk!}\mathrm{\Delta }^2(\lambda _k)e^{\beta _{i=1}^n\lambda _i^2}\right\}_{\mathrm{}}^{\mathrm{}}[d^n\lambda ]𝒪(\lambda _i)𝒫_{n,\beta }(\lambda _k),$$
(12)
where $`\mathrm{\Delta }(\lambda _k)=_{i<j}(\lambda _j\lambda _i)`$ is the Vandermonde determinant (Jastrow factor), and $`[d^n\lambda ]_{i=1}^nd\lambda _i`$.

In RMT, $`\mathrm{\Psi }_{n,\beta }(\gamma )`$ is interpreted as the positional partition function of an ensemble of equal charged particles (with positions given by $`\lambda _i`$) in two dimensions, moving along an infinite line, in thermodynamic equilibrium at temperature $`\gamma `$ \- the so-called “Dyson gas” Dy . Then, $`𝒫_{n,\beta }(\lambda _1,\mathrm{},\lambda _n)`$ defined in (12) is the probability of finding one particle at $`\lambda _1`$, one at $`\lambda _2`$, etc. The value of $`\mathrm{\Psi }_{n,\beta }(\gamma )`$ is known from the Selberg’s integral.

In the region $`\beta \beta _c`$ the partition function $`Z(\beta )`$ is dominated by $`\mathrm{}`$dimensional matrices. Thus, one may try to select the $`\mathrm{}`$dimensional matrices out of the whole ensemble, and then compute the mean of the $`\zeta `$-functions following (10). However, from the standpoint of our statistical approach this procedure does not seems natural since the sum over $`n`$ is a key ingredient in the whole construction. Hence, we look for a quantity related to the $`\zeta `$-function that captures the statistical nature of our model. Let us compute the mean value

$$\mathrm{Tr}_\kappa |D|^\alpha _{n,\beta }=\underset{k=1}{\overset{n}{}}|\lambda _k|^\alpha \theta (|\lambda _k|\kappa )_{n,\beta }.$$
(13)
Bearing in mind (10), we want to study $`\mathrm{Tr}_\kappa |D|^\alpha _{n,\beta }`$ only in the region of eigenvalues where it is a decreasing function of $`\alpha `$. Hence, we introduced an “infrared” cutoff $`\kappa `$, since there is a non-zero probability to find a configuration in the volume $`|\lambda _k|\kappa `$ around the origin. Using the symmetry of the integrand in (13) under permutations of the position indexes, we arrive at

$$\mathrm{Tr}_\kappa |D|^\alpha _{n,\beta }=\sqrt{\beta }_{\mathrm{out},\kappa }𝑑\lambda _n|\lambda _n|^\alpha \sigma _n(\sqrt{\beta }\lambda _n),$$
(14)
where we have used the definition of the spectral density $`\sigma _n(\lambda )=\delta (\lambda \lambda _k)_n`$. Besides, $`_{\mathrm{out},\kappa }_{\mathrm{}}^\kappa +_\kappa ^{\mathrm{}}`$. The computation of (13) was reduced to a one-particle problem. It is known that $`\sigma _n(\lambda )=_k\varphi _k^2(\lambda )`$, $`\varphi _k(x)`$ being the Weber-Hermite functions Meh . In the large $`n`$-limit, $`\sigma _n(\lambda )`$ converges to a nonrandom function, Wigner’s “semi-circle law”,

$$\sigma _n(\sqrt{\beta }\lambda )=\frac{\sqrt{2n}}{\pi }\sqrt{1\frac{\beta }{2n}\lambda ^2},$$
(15)
for $`|\lambda |<_{n,\beta }\sqrt{2n/\beta }`$, and zero otherwise. The average eigenvalue density for Gaussian ensembles of Hermitian matrices with different values of $`n`$ differ only by a change of scale, for large enough $`n`$ (since $`_{n,\beta }\sqrt{n}`$). We propose the following quantity in order to extract information on the $`\zeta `$-function (10):

$$\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha (\beta )=\underset{N\mathrm{}}{lim}\underset{n=1}{\overset{N}{}}P(n,\beta )\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha _{n,\beta }.$$
(16)
Hence, we are sampling the partial traces of the ensemble of $`\mathrm{}`$dimensional Hermitian matrices by the total trace of finite dimensional matrices in the complete Gaussian ensemble of Hermitian matrices. Now, we study the asymptotic properties of $`\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha _{n,\beta }`$ as a function of $`\alpha `$ and $`n`$ for $`nN`$ large. In this case, one may use the semi-circle law in (14),

$$\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha _{n,\beta }\frac{2}{\pi }\sqrt{2n\beta }_{\kappa _{n,\beta }}^{_{n,\beta }}𝑑\lambda _n|\lambda _n|^\alpha \sigma _n(\sqrt{\beta }\lambda _n).$$
(17)
In (17) we are neglecting the contribution from the “tail” of $`\sigma _n(\sqrt{\beta }\lambda )`$, outside of the semi-circle radius $`_{n,\beta }`$. There is an exponential decrease in the tail, so that the total number of particles (or eigenvalues) in this region is of order 2 Bronk . This finite-size correction is not relevant to our asymptotic analysis. Besides, the cutoff $`\kappa `$ is a function of $`n`$ and $`\beta `$: $`\kappa _{n,\beta }=ϵ\sqrt{2n/\beta }`$. The mean particle (or level) spacing is $`\overline{s}(\lambda ,\beta )\sigma _n^1(\lambda ,\beta )`$, so that the mean spacing is almost uniform near the origin. Near the endpoints (edges) the spacing is highly non-uniform. This is the asymptotic region we are interested in. Our procedure is to select a slice of size $`\mathrm{\Delta }_{n,\beta ,\kappa }=_{n,\beta }\kappa _{n,\beta }=(1ϵ)_{n,\beta }`$ near the edge, so that the relative size of the slice, $`\mathrm{\Delta }_{n,\beta ,\kappa }/_{n,\beta }`$, does not depend on $`n`$ and $`\beta `$. This choice ensures that we are treating matrices of different sizes $`n`$ on the same footing. From (17) and (14) we obtain, after some manipulations,

$$\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha _{n,\beta }\frac{2}{\pi }(2n)^{1\frac{\alpha }{2}}\beta ^{\frac{\alpha }{2}}_ϵ^1\frac{dy}{y^\alpha }\sqrt{1y^2}.$$
(18)
The asymptotic behavior of $`\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha _{n,\beta }`$ does not depend in an essential way on the particular choice of $`ϵ`$, as long as we keep $`ϵ0`$.

Now we use the asymptotic formula (18) in (16) and search for the value $`\alpha _c`$ for which, as $`N\mathrm{}`$ and $`\beta \beta _c`$, $`\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha `$ diverges (converges to zero) if $`\alpha <\alpha _c`$ ($`\alpha >\alpha _c`$), with $`\mathrm{Tr}_{\kappa _{n,\beta }}|D|^{\alpha _c}`$ finite and non-zero. This gives an upper bound for the dimension of the “condensed” manifold in the infinite phase ( $`\beta \beta _c`$), which is $`\delta <\alpha _c`$. We obtain:

$$\mathrm{Tr}_{\kappa _{n,\beta }}|D|^\alpha (\beta )\underset{N\mathrm{}}{lim}\underset{n=1}{\overset{N}{}}P(n,\beta )n^{1\frac{\alpha }{2}}.$$
(19)
In the finite phase ($`\beta >\beta _c`$) the sum in (19) converges for $`\alpha 0`$. We conclude that $`\alpha _c=0`$ (i.e. $`\delta =0`$) for $`\beta >\beta _c`$ , as expected. From the behavior of $`P(n,\beta )`$ in the infinite phase it follows that the convergence of the sum in (19) is dictated by the behavior of $`\mathrm{\Gamma }_{n,\alpha }=n^{1\frac{\alpha }{2}}`$ in the limit $`nN\mathrm{}`$. For $`\beta \beta _c`$ we get $`\mathrm{\Gamma }_{n,\alpha }\mathrm{}`$ if $`\alpha <2`$, and $`\mathrm{\Gamma }_{n,\alpha }0`$ if $`\alpha >2`$. For $`\alpha =2`$ it turns out that $`\mathrm{Tr}_{\kappa _{n,\beta }}|D|^2(\beta )1`$. Therefore, we obtain the upper bound $`\delta <2`$.

Notice that we do not have a definition of the dimension as an operator (as suggested by (9)), from which it will be possible to compute its average in the ensemble of all Hermitian matrices, like we have do for $`n`$ in (6). Instead, our guide is the operational definition encapsulated in (16). In order to go beyond the upper bound computed above we need to find a suitable observable which reduces to the classical dimension (as given by the Weyl formula) in the limit $`n\mathrm{}`$. This eventually may lead to a numerical computation of $`\delta `$. Besides, one sees that the reason for the upper bound $`\delta <2`$ lies in the semi-circle law and its leading $`\sqrt{\lambda }`$ behavior near the edge. It is known from $`2D`$ models of discretized pure quantum gravity Bre ; Kan  that some special matrix polynomial potentials $`V(D)=a_k\mathrm{Tr}D^{2k}`$ may lead, in a suitable scaling limit, to a behavior near the edges different from the square root. Thus, a possible way to obtain a bound of higher dimension would be to include higher polynomial interactions, or modify the quadratic one in (2) including more (internal) symmetries besides the unitary one. These and other questions are under study, and will be reported elsewhere.

To conclude, in this Letter we proposed a discrete model for quantum gravity based on the spectral principle. The model is connected with the GUE of Hermitian matrices, and contains the mean number of points, $`n`$, and the dimension of the space-time, $`\delta `$, as dynamical observables. We have shown that the model has two phases: a finite phase with a finite value of $`n`$ and $`\delta =0`$, and an infinite phase with a diverging $`n`$ and a finite $`\delta 0`$. The critical point was computed, $`\beta _c=\pi /2`$, as well as the critical exponent of $`n`$. Moreover, an upper bound for the order parameter $`\delta `$ was found, $`\delta <2`$.

We thank S. Vaidya for suggestions and valuable comments, A. P. Balachandran, D. Dalmazi, J. C. A. Barata and A. Bertuola for helpful conversations. L.C.A. would like to thank the Mathematical Physics Department of USP at São Paulo for their kind hospitality. This work was partially supported by FAPESP, grant 00/03277-3 (L.C.A.), and CNPq (P.T.S.).


text:
# The particle-in-cell model for ab initio thermodynamics: implications for the elastic anisotropy of the Earth’s inner core.

## 1 Introduction

Ab initio calculations based on density functional theory (DFT) (Hohenberg and Kohn (1964); Kohn and Sham (1965); Jones and Gunnarsson (1989)) are widely used to calculate the properties of materials at the extreme pressures found in the interior of the Earth and other planets (Stixrude et al. (1998)), and are known to be capable of high accuracy. For many years, such calculations did not generally include thermal effects, and explicitly treated only zero-temperature materials. However, the past few years have seen an increasing effort to apply DFT calculations to high-temperature solids and liquids of geological interest, using ab initio lattice dynamics or molecular-dynamics simulation (Alfè et al. (2001, 2002); Belonoshko et al. (2000); Laio et al. (2000); Stixrude et al. (1997); Wasserman et al. (1996); Brodholt et al. (2002); Oganov et al. (2001); Steinle-Neumann et al. (2002)). Because dynamical DFT calculations demand large computer resources, a simplifying approximation known as the ‘particle-in-cell’ (PIC) method (Hirshfelder et al. (1954); Holt and Ross (1970); Holt et al. (1970); Ree and Holt (1973); Westra and Cowley (1975); Cowley et al. (1990)) has been used in some work on high-$`T`$, high-$`p`$ solids, an important example being the very recent use of the PIC approximation to re-interpret the elastic anisotropy of the Earth’s solid inner core (Steinle-Neumann et al. (2001, 2002)). However, the reliability of the PIC method has not gone unchallenged (Alfè et al. (2001)) and conclusions based on it are not necessarily secure. To shed light on this question, we present here the results of our own PIC calculations of a range of thermodynamic properties of solid iron at Earth’s core conditions, which we compare with more exact calculations that avoid the PIC approximation.

The purpose of the PIC approximation is to provide a way of calculating the free energy of a vibrating crystal. The essence of the approximation is that correlations between the vibrational displacements of different atoms are neglected, so that each atom is treated as vibrating independently of every other, as in the Einstein model for a vibrating solid. Although the approximation seems at first sight rather crude, empirically it has been shown to yield satisfactory predictions for the thermodynamic properties of a number of solids (Holt and Ross (1970); Holt et al. (1970); Ree and Holt (1973); Westra and Cowley (1975)), at least for temperatures above the Debye temperature. In the geological context, it was shown by Wasserman et al. (1996) that, when implemented using DFT methods, it gives predictions for the $`p(V)`$ relation on the shock Hugoniot in excellent agreement with experiment for Fe up to the melting curve. The PIC approach based on DFT has also been applied successfully to the high-$`p`$, high-$`T`$ properties of other transition metals (Gülseren and Cohen (2001); Cohen and Gülseren (2002)).

In the work on the elastic anisotropy of the Earth’s inner core (Steinle-Neumann et al. (2001, 2002)) mentioned above, PIC was used to calculate the elastic constants of hexagonal-close-packed (hcp) Fe over a range of pressures and temperatures. A key prediction was that the $`c/a`$ ratio increases strongly with $`T`$, and that this leads in turn to a strong $`T`$-dependence of the elastic constants. It emerged from this that the observed anisotropy of seismic velocities in the inner core (Creager (1992); Song and Helmberger (1993); Tromp (1993)) could be interpreted in terms of a partial alignment of hcp crystallites. Remarkably, the calculations required the hcp basal plane to be preferentially aligned along the Earth’s rotational axis, which is the opposite of what had been proposed earlier (Stixrude and Cohen (1995); Steinle-Neumann et al. (1999)). Since the correctness of this interpretation depends heavily on the predicted variation of $`c/a`$ with $`T`$, it is clearly essential to be confident that the PIC approximation can be relied on to give this variation correctly, and one of our aims in this paper is to test this point.

An objective assessment of the errors incurred by the PIC approximation is made possible by the fact that ab initio free energies and other thermodynamic functions of high-temperature solids can now be calculated with statistical-mechanical errors that can be made arbitrarily small. The new methods, described in detail in previous papers (Alfè et al. (2001, 2002)), are based on the DFT calculation of phonon frequencies in the harmonic approximation, supplemented by the ‘thermodynamic integration’ technique (Frenkel and Smit (1996)) for calculating anharmonic contributions. For a given material, and with a given DFT technique for calculating the electronic total energy as a function of ionic positions, we can therefore compute thermodynamic functions either with the PIC approximation or almost exactly. The differences between the two sets of results represent the errors caused by PIC. This is what we do in the present work. Since the main difference between PIC and the newer, more exact methods is that the latter fully include vibrational correlations, we shall refer to these in the following as ‘vibrationally correlated’ calculations.

Two main claims have been made for the PIC method (Wasserman et al. (1996)): First, that it provides a simple and reasonably accurate way of including the effect of the anharmonicity of lattice vibrations on thermodynamic properties; and second, that, even in the absence of anharmonicity, it is computationally much less demanding than more precise methods. Our assessment of these claims for the case of high-$`p`$/high-$`T`$ Fe will suggest that neither is necessarily true, but that nevertheless the PIC method does yield surprisingly accurate predictions for many thermodynamic properties. We shall elucidate the reasons why PIC is often accurate. As we shall see, however, the present work does not support the prediction that the $`c/a`$ ratio of hcp Fe increases strongly with temperature, and this casts doubt on the proposed re-interpretation of the elastic anisotropy of the Earth’s inner core.

The remainder of the paper is organised as follows. In Sec. 2, we describe how we have applied PIC to hcp Fe and how we have separated the various contributions to the free energy and other thermodynamic functions; at the end of Sec. 2, we summarise the DFT techniques. In Sec. 3, we describe the details of our free-energy calculations and highlight their implications for the temperature dependence of the $`c/a`$ ratio. We then present results for several thermodynamic quantities calculated in the PIC approximation, which we compare with the earlier PIC results of Cohen and co-workers (Stixrude et al. (1997); Steinle-Neumann et al. (2001, 2002); Wasserman et al. (1996)) and with the vibrationally correlated results of Alfè et al. (2001). Discussion and conclusions follow in Sec. 4.

## 2 Techniques

### 2.1 Particle-in-cell approximation

In classical statistical mechanics, the ab initio Helmholtz free energy of a vibrating solid containing $`N`$ ions is:

$$F_{\mathrm{AI}}=k_\mathrm{B}T\mathrm{ln}\left\{\frac{1}{\mathrm{\Lambda }^{3N}}𝑑𝐫_1\mathrm{}𝑑𝐫_N\mathrm{exp}\left[\beta U_{\mathrm{AI}}(𝐫_1\mathrm{}𝐫_N)\right]\right\},$$
(1)
where $`U_{\mathrm{AI}}(𝐫_1\mathrm{}𝐫_N)`$ is the ab initio total electronic free energy of the system when the ionic positions are $`𝐫_1,\mathrm{}𝐫_N`$ and $`\mathrm{\Lambda }`$ is the thermal wavelength. It should be noted that $`F_{\mathrm{AI}}`$ depends on the volume $`V`$ and temperature $`T`$ of the system ($`\beta =1/k_\mathrm{B}T`$), and for an hcp crystal it also depends on the axial $`c/a`$ ratio, denoted here by $`q`$; for the moment we do not indicate explicitly the dependence on $`V`$, $`T`$ and $`q`$. Two other points should also be recalled. First, the quantity $`U_{\mathrm{AI}}(𝐫_1,\mathrm{},𝐫_N)`$ is a free energy, because the electrons are treated as being in thermal equilibrium for each set of ionic positions, at a temperature equal to the temperature $`T`$ of the system as a whole. Second, the use of classical statistical mechanics, i.e. the neglect of quantum effects in the nuclear motions, is fully justified at Earth’s core temperatures (Alfè et al. (2001)).

The PIC approach consists in replacing $`U_{\mathrm{AI}}`$ in the above formula by the approximate form $`U_{\mathrm{AI}}^{\mathrm{PIC}}`$, given by:

$$U_{\mathrm{AI}}^{\mathrm{PIC}}(𝐫_1,\mathrm{}𝐫_N)=U_{\mathrm{perf}}+\underset{i=1}{\overset{N}{}}\varphi (𝐮_i).$$
(2)
Here, $`U_{\mathrm{perf}}U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_N)`$ is the ab initio free energy of the system when all ions are at their perfect-lattice positions $`𝐑_i`$, and $`\varphi (𝐮_i)`$ is defined to be the change of energy $`U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_i+𝐮_i,\mathrm{}𝐑_N)U_{\mathrm{perf}}`$ when ion $`i`$ is displaced from its perfect-lattice position to the position $`𝐑_i+𝐮_i`$, all other ions being held fixed at their perfect-lattice positions. With $`U_{\mathrm{AI}}`$ replaced by $`U_{\mathrm{AI}}^{\mathrm{PIC}}`$ in Eqn. (1), the PIC approximation for the free energy is:

$$F_{\mathrm{AI}}^{\mathrm{PIC}}=U_{\mathrm{perf}}+Nf_{\mathrm{vib}}^{\mathrm{PIC}},$$
(3)
where the PIC vibrational free energy per atom is:

$$f_{\mathrm{vib}}^{\mathrm{PIC}}=k_\mathrm{B}T\mathrm{ln}\left[\mathrm{\Lambda }^3𝑑𝐮\mathrm{exp}(\beta \varphi (𝐮))\right].$$
(4)
The problem is thus reduced to the calculation of a three-dimensional integral, which can be performed numerically.

Our later analysis depends on a clear separation of harmonic and anharmonic contributions to $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$, which can be accomplished by considering the series expansion of $`\varphi (𝐮)`$ in powers of ionic displacement $`𝐮`$:

$$\varphi (𝐮)=\frac{1}{2!}\underset{\alpha \beta }{}\mathrm{\Phi }_{\alpha \beta }^{(2)}u_\alpha u_\beta +\frac{1}{3!}\underset{\alpha \beta \gamma }{}\mathrm{\Phi }_{\alpha \beta \gamma }^{(3)}u_\alpha u_\beta u_\gamma +\frac{1}{4!}\underset{\alpha \beta \gamma \delta }{}\mathrm{\Phi }_{\alpha \beta \gamma \delta }^{(4)}u_\alpha u_\beta u_\gamma u_\delta +\mathrm{},$$
(5)
where the Greek indices $`\alpha `$, $`\beta `$… indicate Cartesian components. The value of $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$ when we retain only the quadratic part of $`\varphi (𝐮)`$ is the harmonic vibrational free energy in the PIC approximation, denoted by $`f_{\mathrm{harm}}^{\mathrm{PIC}}`$. This can be expressed as:

$$f_{\mathrm{harm}}^{\mathrm{PIC}}=k_\mathrm{B}T\underset{\nu =1}{\overset{3}{}}\mathrm{ln}(\mathrm{}\omega _\nu /k_\mathrm{B}T)=3k_\mathrm{B}T\mathrm{ln}(\mathrm{}\overline{\omega }/k_\mathrm{B}T).$$
(6)
Here $`\omega _\nu `$ are the three PIC vibrational frequencies, given by $`det|M\omega _\nu ^2\delta _{\alpha \beta }\mathrm{\Phi }_{\alpha \beta }^{(2)}|=0`$, with $`M`$ the ionic mass. The geometric-mean frequency $`\overline{\omega }`$ is defined by Eqn. (6). The anharmonic contribution $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$ to the free energy is the part of $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$ not accounted for by $`f_{\mathrm{harm}}^{\mathrm{PIC}}`$, so that $`f_{\mathrm{anharm}}^{\mathrm{PIC}}f_{\mathrm{vib}}^{\mathrm{PIC}}f_{\mathrm{harm}}^{\mathrm{PIC}}`$.

The PIC results reported in Sec. 3 were obtained by computing $`\varphi (𝐮)`$ for a set of vector displacements $`𝐮`$, and fitting the results using the power series expansion of Eqn. (5); in practice, for the range of displacements that occur with appreciable probability at temperatures below the melting temperature, we find that an extremely accurate fit is obtained if we retain only terms up to quartic order in $`𝐮`$. To show what parameters enter this fit, we write this quartic polynomial explicitly for the case of hcp symmetry:

$$\begin{array}{cc}\hfill \varphi (𝐮)=& \frac{1}{2!}\left[M\omega _a^2(u_x^2+u_y^2)+M\omega _c^2u_z^2\right]+\frac{1}{3!}K^{(3)}(u_y^33u_x^2u_y)\hfill \\ & +\frac{1}{4!}\left[K_a^{(4)}(u_x^2+u_y^2)^2+K_{\mathrm{mix}}^{(4)}u_z^2(u_x^2+u_y^2)+K_c^{(4)}u_z^4\right].\hfill \end{array}$$
(7)
Here, $`u_x`$ and $`u_y`$ are perpendicular Cartesian components in the basal plane, the $`x`$-axis being oriented towards a nearest-neighbour, and $`u_z`$ is the displacement along the hexagonal axis. The frequencies of harmonic vibration in the basal plane and along the hexagonal axis are $`\omega _a`$ and $`\omega _c`$ respectively. This symmetrised polynomial expression may be obtained to a given order by writing down all possible terms in a polynomial of that order, and retaining only those terms which leave the expression invariant under all the point symmetry operations of the chosen cell. In the case of hcp symmetry, these include reflection in the hexagonal ($`z=0`$) plane, rotations of $`\frac{2\pi }{3}`$ about the hexagonal axis, and reflection in the $`y`$$`z`$ plane. With the above form of $`\varphi `$, the geometric-mean frequency $`\overline{\omega }`$ given by

$$3\mathrm{ln}\overline{\omega }=2\mathrm{ln}\omega _a+\mathrm{ln}\omega _c.$$
(8)
The vibrational free energy $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$ is obtained essentially exactly by numerical evaluation of the integral in Eqn. (4) on a regular grid. In analysing the anharmonic contributions, it is useful to note that $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$ can be expanded as a power series in temperature:

$$f_{\mathrm{anharm}}^{\mathrm{PIC}}=dT^2+O(T^3).$$
(9)
For temperatures below the melting point, we find that only the term in $`T^2`$ is significant, so that in practice the anharmonic contributions are completely specified by the coefficient $`d`$ as a function of volume $`V`$ and $`c/a`$ ratio $`q`$. The following exact expression for the coefficient is readily obtained:

$$d=K_a^{(4)}/3(M\omega _a^2)^2+K_{\mathrm{mix}}^{(4)}/12(M\omega _a\omega _c)^2+K_c^{(4)}/8(M\omega _c^2)^2\left(K^{(3)}\right)^2/3(M\omega _a^2)^3.$$
(10)
The above methods are used to obtain contributions to the free energy at a reference value of $`q`$, chosen to be close to the $`T=0`$ equilibrium value. We then perform further calculations over a range of values of $`q`$, to obtain corrections to the free energy components, which are parameterised in terms of $`V`$ and $`T`$. The details of this parameterisation are given in Sec. 3.1 and 3.2. Equilibrium values of $`q`$ are found by direct minimisation of the total free energy.

### 2.2 Ab initio methods

The DFT electronic-structure techniques used to perform the PIC calculations are essentially identical to those described by (Alfè et al. (2001, 2002)). The exchange-correlation functional is the generalised gradient approximation (GGA) of Perdew and Wang (Wang and Perdew (1991); Perdew et al. (1992)). The implementation of DFT is the projector augmented wave (PAW) scheme (Blöchl (1994); Kresse and Joubert (1999)), with core radii, augmentation charge radii etc set to the values reported in Alfè et al. (2000). As before, all atomic states up to and including $`3p`$-states are treated as core states, but the high-pressure response of $`3s`$\- and $`3p`$-states, known to be important at Earth’s-core conditions, is included via an effective pair potential; the accuracy of this procedure has been demonstrated earlier (Alfè et al. (2001)). Thermal excitation of electrons, also important at core conditions, is treated with the usual finite-temperature formulation of DFT (Mermin (1965); Gillan (1989); Wentzcovitch et al. (1992)). We used a plane-wave cut-off of 300 eV, as in our previous work. All calculations were performed with the VASP code (Kresse and Furthmüller (1996a, b)).

The PIC vibrational potential $`\varphi (𝐮)`$ (see Eqn. (2)) should in principle be calculated by displacing a single atom in an infinite crystal. Because PAW calculations require periodic boundary conditions, we must instead use supercell geometry, so that the displaced atom has periodic images. To adhere to the PIC scheme, we must therefore ensure convergence of all results with respect to the size of the supercell, as described in the following Section.

## 3 Results

### 3.1 Free energy of perfect lattice

DFT results for the free energy of the hcp Fe perfect lattice $`U_{\mathrm{perf}}(V,q,T)`$ for the fixed $`c/a`$ ratio $`q`$ equal to 1.60 were reported earlier (Alfè et al. (2001)) for $`V`$ values from 5.2 to 11.4 Å3/atom at temperatures from 200 to $`10^4`$ K. At each $`T`$ value, the results were fitted to a Burch-Murnaghan equation of state (Poirier (2000)), the parameters of which were then fitted as polynomial functions of $`T`$. We have repeated these calculations for the present work, and as expected the results are virtually identical to those reported earlier.

To allow for variation of $`q`$, we have performed additional calculations of $`U_{\mathrm{perf}}(V,q,T)`$ for $`q`$ in the range from 1.48 to 1.72, with $`V`$ going from 5.5 to 10.5 Å3/atom, and $`T`$ going from 2000 to 8000 K. All technical parameters, such as the Monkhorst-Pack sampling set, were kept the same as in the $`q=1.60`$ calculations. To represent the results, we define the deviation $`\mathrm{\Delta }U_{\mathrm{perf}}(V,q,T)`$ of the perfect-lattice free energy from its value at a chosen of $`q`$, denoted by $`q_0`$: $`\mathrm{\Delta }U_{\mathrm{perf}}(V,q,T)U_{\mathrm{perf}}(V,q,T)U_{\mathrm{perf}}(V,q_0,T)`$. In the present case, $`q_0`$ has the fixed value 1.60 throughout. We find that $`\mathrm{\Delta }U_{\mathrm{perf}}(V,q,T)`$ can be very accurately represented at all $`(V,T)`$ by the quadratic form:

$$\mathrm{\Delta }U_{\mathrm{perf}}(q)=\alpha (qq_0)(qq_1).$$
(11)
At each $`V`$ value, $`\alpha `$ and $`q_1`$ can be accurately fitted as linear functions of $`T`$, the coefficients of which are in turn fitted linearly in $`V`$. These fits give a virtually perfect representation of the $`\mathrm{\Delta }U_{\mathrm{perf}}(V,q,T)`$ results. Writing $`\alpha =a+bT`$, $`q_1=r+sT`$ and with $`a=a^{(0)}+a^{(1)}V`$ and similarly for $`b`$, $`r`$ and $`s`$, the numerical values of the parameters are $`a^{(0)}=18.39`$, $`a^{(1)}=1.532`$, $`b^{(0)}=3.070\times 10^4`$, $`b^{(1)}=1.64908\times 10^5`$, $`r^{(0)}=1.63`$, $`r^{(1)}=7.909\times 10^3`$, $`s^{(0)}=7.902\times 10^6`$ and $`s^{(1)}=1.817\times 10^6`$. Units of coefficients are such that $`\alpha `$ is in eV, $`T`$ in K and $`V`$ in Å3.

The equilibrium value $`q`$ of the perfect lattice, denoted by $`q_{\mathrm{eq}}`$, is obtained by minimising $`\mathrm{\Delta }U_{\mathrm{perf}}(V,q,T)`$ with respect to $`q`$. At $`T=0`$, we find a very weak dependence of $`q_{\mathrm{eq}}`$ on $`V`$, going from 1.590 at 7.0 Å3/atom to 1.578 at 10.0 Å3/atom. The $`T`$-dependence of $`\mathrm{\Delta }U_{\mathrm{perf}}(q)`$ is also very weak, with $`\alpha `$ decreasing linearly from 7.4 to 6.2 eV, and $`q_1`$ increasing from 1.59 and 1.62 as temperature varies from 2000 to 8000 K. If we now use the results to predict $`q_{\mathrm{eq}}`$ for the (hypothetical) high-$`T`$ perfect lattice, we find a variation of at most 0.04 at a volume of 10.0 A3. The insignificant $`V`$\- and $`T`$-dependence of $`q_{\mathrm{eq}}`$ is noteworthy, because it implies that any significant variation of $`q_{\mathrm{eq}}`$ at high temperatures can come only from lattice vibrations, to which we now turn.

### 3.2 Vibrational free energy

The basal-plane and axial frequencies $`\omega _a`$ and $`\omega _c`$, and the four anharmonic vibrational coefficients $`K^{(3)}`$, $`K_a^{(4)}`$, $`K_{\mathrm{mix}}^{(4)}`$ and $`K_c^{(4)}`$ were calculated as follows. For given values of $`q`$ and $`V`$, a supercell was constructed in which one of the atoms, the ‘walker’ is displaced from its equilibrium position. The ‘walker’ is given a series of equally spaced displacements in the three directions $`\widehat{𝐮}_y`$, $`\widehat{𝐮}_z`$ and $`\frac{1}{\sqrt{2}}(\widehat{𝐮}_x+\widehat{𝐮}_𝐳)`$, the maximum displacements $`𝐫_{\mathrm{max}}`$ in each direction being chosen so that the Boltzmann factor $`\mathrm{exp}\left[\beta \varphi (𝐫_{\mathrm{max}})\right]0.1`$ for the maximum temperatures of interest. Directions were selected to allow the terms in $`\varphi (𝐫)`$ (Eqn. (7)) to be determined as simply as possible by least squares fitting. Size convergence in terms of supercells is essential to ensure that the walker cannot ‘see’ its image. We discuss the results of a simple convergence test below. All the calculations described here are performed at a fixed electronic temperature of 6000 K, since statistical mechanical considerations will dominate the temperature dependence of $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$. At this temperature, vibrational coefficients were fully converged with respect to $`k`$-point sampling at each cell size. Monkhorst-Pack (Monkhorst and Pack (1976)) sampling was used, with $`9\times 9\times 5`$ $`k`$-points for the 8 atom cell, $`7\times 7\times 5`$ $`k`$-points for 16 atoms and $`3\times 3\times 3`$ for 36. Calculations were performed over a range of volumes from 5.5 Å3/atom to 11.5 Å3/atom with a fixed axial ratio of $`q_0=1.60`$ as before. Using 5th order polynomials to obtain a high quality fit, the quantities $`K^{(3)}`$, $`K_a^{(4)}`$, $`K_{\mathrm{mix}}^{(4)}`$, $`K_c^{(4)}`$ and $`\mathrm{ln}\overline{\omega }`$ were then parameterised in terms of volume. Finally, a correction due to relaxation of $`q`$ is added to the harmonic free energy, as for that of the perfect lattice, by performing a parameterisation of $`\overline{\omega }`$ in terms of $`q`$, which we describe below. Parameterisation of higher order $`K`$ coefficients in terms of $`q`$ yielded only negligible corrections to the $`p`$$`V`$ curve.

As noted in Sec. 2.1, the harmonic free energy is completely determined by the geometric-mean frequency $`\overline{\omega }`$ (Eqns. (6) and (8)), and we consider first our PIC results for $`\overline{\omega }`$ as a function of $`V`$ for the fixed $`c/a`$ ratio $`q_0=1.6`$. We report in Fig. 1 our $`\overline{\omega }`$ results from calculations using supercells of 8, 16 and 36 atoms, compared with the $`\overline{\omega }`$ results from our earlier calculations of the full phonon spectrum (Alfè et al. (2001)). We note two important points. First, the PIC $`\overline{\omega }`$ results are almost independent of supercell size, so that the calculations appear to be fully converged with respect to size effects. Second, the PIC $`\overline{\omega }`$ differs from the full-phonon $`\overline{\omega }`$ only by a volume-independent shift of $`\mathrm{ln}\overline{\omega }`$. Since such a $`V`$-independent shift cannot affect most thermodynamic quantities, we expect the PIC calculations to agree rather closely with full-phonon calculations.

We now consider the dependence of $`\overline{\omega }`$ on $`c/a`$ ratio $`q`$, defining the logarithmic deviation $`\mathrm{\Delta }_\mathrm{L}\overline{\omega }(V,q)`$ from its value at $`q_0`$ by $`\mathrm{ln}\overline{\omega }(V,q)=\mathrm{ln}\overline{\omega }(V,q_0)+\mathrm{\Delta }_\mathrm{L}\overline{\omega }(V,q)`$. In terms of this, the harmonic free energy can be written as:

$$f_{\mathrm{harm}}^{\mathrm{PIC}}(V,q,T)=3k_\mathrm{B}T\mathrm{ln}(\mathrm{}\overline{\omega }/k_\mathrm{B}T)=f_{\mathrm{harm}}^{\mathrm{PIC}}(V,q_0,T)+3k_\mathrm{B}T\mathrm{\Delta }_\mathrm{L}\overline{\omega }.$$
(12)
Calculations of $`\overline{\omega }(q)`$ were made for $`q`$-values going from 1.48 to 1.72 at a series of volumes. We find that $`\overline{\omega }`$ decreases with increasing $`q`$. Not surprisingly, increase of $`q`$ at constant $`V`$ yields an increase of the basal-plane frequency $`\omega _a`$ and a decrease of axial frequency $`\omega _c`$. The decrease of $`\omega _c`$ succeeds in outweighing the increase of $`\omega _a`$ in the formula for $`\overline{\omega }`$ (Eqn. (8)), but nevertheless the resulting decrease of $`\overline{\omega }`$ is a fairly marginal effect. Over the $`q`$-range studied, we found that the quadratic form $`\mathrm{\Delta }_\mathrm{L}\overline{\omega }(q)=\beta (qq_0)(qq_2)`$ gives an accurate fit, and this fit was performed at volumes of 5.5, 7.0 and 10.5 Å3/atom. The $`V`$-dependence of the resulting $`\beta `$ and $`q_2`$ coefficients was obtained by a quadratic fit. Writing $`\beta =_{i=0}^2\beta ^{(i)}V^i`$ and similarly for $`q_2`$, the values of the parameters are $`\beta ^{(0)}=2.083`$, $`\beta ^{(1)}=0.375`$, $`\beta ^{(2)}=0.0264`$, $`q_2^{(0)}=1.990`$, $`q_2^{(1)}=0.141`$ and $`q_2^{(2)}=8.79\times 10^3`$. Units of parameters are such that $`\overline{\omega }`$ is in rad s-1 and $`V`$ is in Å3.

From the values of the four anharmonic coefficients at each $`V`$ and $`q`$, we calculate also the anharmonic vibrational free energy $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$. As explained in Sec. 2.1, the most accurate way of doing this is by explicit calculation of $`f_{\mathrm{vib}}^{\mathrm{PIC}}`$ (Eqn. (4)) by numerical evaluation of the integral over displacement $`𝐮`$, from which $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$ is obtained as the difference $`f_{\mathrm{vib}}^{\mathrm{PIC}}f_{\mathrm{harm}}^{\mathrm{PIC}}`$. Alternatively, we can use Eqn. (10) to calculate the $`d`$ coefficient of the leading term in the temperature expansion of $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$ (see Eqn. (9)). In practice, we find that the two methods yield almost indistinguishable results, and everything that follows is based on evaluation of the $`d`$ coefficient.

We report in Fig. 2 the volume dependence of $`d`$ for $`q=1.60`$, compared with the corresponding anharmonicity coefficient obtained from the vibrationally correlated results of Alfè et al. (2001). The coefficients of this fit are as follows: $`d^{(0)}=6.556\times 10^8`$, $`d^{(1)}=3.399\times 10^8`$, $`d^{(2)}=6.487\times 10^9`$, $`d^{(3)}=5.400\times 10^{10}`$, $`d^{(4)}=1.639\times 10^{11}`$, where $`d(v)=_{i=0}^4d^{(i)}V^i`$, $`d`$ being in eV K-2 and $`V`$ in Å3. Corresponding data is given in Alfè et al. (2001). The striking and important feature of this comparison is that the PIC anharmonic free energy is positive, whereas the results of Alfè et al. (2001) show that in reality it is negative. This means that for high-$`p`$ hcp Fe, the PIC approximation gives a completely incorrect account of anharmonicity, the electronic structure techniques in both works being identical. However, we note that $`f_{\mathrm{anharm}}^{\mathrm{PIC}}`$ is in any case very small. For example, at $`V=7`$ Å3/atom, $`T=6000`$ K, $`f_{\mathrm{anharm}}^{\mathrm{PIC}}15`$ meV/atom. Because it is so small, we ignore its dependence on $`q`$ in the following, and use a polynomial fit to the $`d(V)`$ results for $`q=1.60`$.

### 3.3 Temperature dependence of $`c/a`$ ratio

The equilibrium value $`q_{\mathrm{eq}}`$ of the $`c/a`$ ratio $`q`$ at given $`V`$ and $`T`$ is obtained by minimising the total free energy $`F_{\mathrm{AI}}^{\mathrm{PIC}}`$ with respect to $`q`$. The variation of $`q`$ with $`T`$ comes from the $`T`$\- and $`q`$-dependence of the perfect-lattice free energy $`U_{\mathrm{perf}}^{\mathrm{PIC}}`$, for which we presented results in Sec. 3.1, and from the $`q`$-dependence of the mean vibrational frequency $`\overline{\omega }`$ (Eqn. (6)). (Since we neglect the $`q`$-dependence of the anharmonic component of free energy, anharmonicity does not contribute to the $`T`$-dependence of $`q_{\mathrm{eq}}`$ here.)

Our results for $`q_{\mathrm{eq}}`$ at three different volumes are reported in Fig. 3, where we compare them with the earlier PIC predictions of Steinle-Neumann et al. (2001, 2002). The present results are very different from the earlier ones. At all volumes, we find only a very weak increase of $`q`$ with $`T`$, which is between 5 and 10 times smaller than the variation reported in Steinle-Neumann et al. (2001). We note that this gross discrepancy can come only from a difference in the $`q`$-dependence of the harmonic mean frequency $`\overline{\omega }`$. The reason for this is that the roughly linear dependence of $`q_{\mathrm{eq}}`$ on $`T`$ seen in both our results and those of Steinle-Neumann et al. (2001) cannot originate either from thermal electronic excitations or from anharmonicity, since both of these free-energy components vary as $`T^2`$. The cause of the discrepancy, and its implications for understanding the elastic anisotropy of the inner core, are discussed further in Sec. 4.

### 3.4 Thermodynamic functions

All thermodynamic functions can be calculated by taking appropriate derivatives of the PIC free energy $`F_{\mathrm{AI}}^{\mathrm{PIC}}`$, with its perfect-lattice, harmonic and anharmonic components parameterised as described above. All the results to be presented include the dependence of the equilibrium $`c/a`$ ratio $`q_{\mathrm{eq}}`$ on $`V`$ and $`T`$. Our PIC predictions will be compared with the earlier PIC results of Wasserman et al. (1996) and with the vibrationally correlated results of Alfè et al. (2001).

We begin by considering the thermal pressure $`\mathrm{\Delta }p`$, which is one of the most basic quantities in interpreting the properties of the inner core. This is defined as the difference between the pressure in the material at a given $`V`$ and $`T`$ and the pressure at the same volume but at zero temperature:

$$\mathrm{\Delta }p=(F_{\mathrm{AI}}^{\mathrm{PIC}}/V)_T+(F_{\mathrm{AI}}^{\mathrm{PIC}}/V)_{T=0}.$$
(13)
Fig. 4 shows a comparison of our PIC results for $`\mathrm{\Delta }p`$ with the vibrationally correlated results of Alfè et al. (2001) on isotherms at $`T=2000`$, 4000 and 6000 K. The rather close agreement shows that the PIC approximation gives an accurate account of the thermal pressure. This is expected from our results for $`\overline{\omega }`$, since in the harmonic approximation $`\mathrm{\Delta }p`$ is given by the electronic thermal pressure, which is exactly the same in both calculations, plus a vibrational contribution equal to $`3k_\mathrm{B}T\mathrm{ln}\overline{\omega }/V`$. We have seen that for hcp Fe, the PIC value of $`\mathrm{ln}\overline{\omega }`$ differs from the exact value by an almost constant offset, which has no effect on the derivative of $`\mathrm{ln}\overline{\omega }`$ with respect to $`V`$. Small differences between the PIC and vibrationally correlated thermal pressure may be due to small variations in the offset between the values of $`\mathrm{ln}\overline{\omega }`$.

The situation is similar for the thermal energy $`\mathrm{\Delta }E`$, defined as the difference between the internal energy $`E`$ at a given $`V`$ and $`T`$ and the internal energy at the same $`V`$ but zero temperature: $`\mathrm{\Delta }E=E(V,T)E(V,0)`$. The electronic part of the thermal energy is identical in both PIC and vibrationally correlated calculations, and the harmonic vibrational energy is exactly $`3k_\mathrm{B}T`$/atom in both cases, so that any difference arises only from the small anharmonic contribution. Results for the thermal energy are not presented here.

The good accuracy of PIC for thermal pressure and thermal energy explains why it has also been found to give a good account of shock measurements (Wasserman et al. (1996)). A conventional sequence of shock experiments on samples in the same initial state generates a path through thermodynamic state-space known as the Hugoniot. On this path, the pressure $`p`$, volume $`V`$ and internal energy $`E`$ are related by the Rankine-Hugoniot formula (Poirier (2000)):

$$\frac{1}{2}(p+p_0)(V_0V)=EE_0,$$
(14)
where $`p_0`$, $`V_0`$ and $`E_0`$ refer to the initial state, which usually corresponds to ambient conditions. For given $`V`$, the Hugoniot $`p`$ and temperature $`T`$ can be determined by going to the calculated $`p(V,T)`$ and $`E(V,T)`$ functions and seeking the value of $`T`$ for which $`p`$ and $`E`$ satisfy eqn (14). Fig. 5 shows our calculated PIC $`p(V)`$ on the Hugoniot compared with the experimental results and with the fully correlated results of Alfè et al. (2001). The almost exact agreement between the three sets of results confirms the excellence of PIC for this particular purpose. Similar comparisons are shown for the $`T(p)`$ relation on the Hugoniot in Fig. 6. Here there is a discrepancy, which corresponds to a difference in internal energy between the two theoretical models. This difference is fully accounted for by the neglect of temperature dependence of $`\overline{\omega }`$ in the present work.

We conclude the presentation of results by examining three thermodynamic quantities for which our comparisons between vibrationally correlated calculations and the earlier PIC results of Wasserman et al. (1996) and Stixrude et al. (1997) revealed significant differences (Alfè et al. (2001)).

Fig. 7 compares the present PIC results for thermal expansivity $`\alpha `$ with earlier calculations. At 2000 K, the present results are in very close agreement with the vibrationally correlated results, and in strong disagreement with the earlier PIC results at low pressures. At higher $`T`$, the present results fall somewhat below both previous sets of results. The product $`\alpha K_T`$ (Fig 8) of expansivity and isothermal bulk modulus is important in high-pressure work, because it can sometimes be assumed to be independent of $`p`$ and $`T`$ over a wide range of conditions. The vibrationally correlated calculations showed that for high-$`p`$/high-$`T`$ hcp Fe constancy with $`p`$ is a good approximation, but constancy with $`T`$ is not. The present PIC results show a reduction of $`\alpha K_T`$ with respect to vibrationally correlated and earlier PIC results of, at most, around 15% at 6000K.

The thermodynamic Grüneisen parameter $`\gamma =V(p/E)_V`$ is particularly important, because it relates thermal pressure to thermal energy, and assumptions about $`\gamma `$ are often used in reducing shock data from Hugoniot to isotherm. The large differences between the earlier PIC results and the vibrationally-correlated results for $`\gamma `$ are therefore a cause for concern (Alfè et al. (2001)). The present PIC results (Fig. 9) agree quite closely with the vibrationally-correlated results, and this suggests that the cause of the earlier disagreement was not the PIC approximation itself.

## 4 Discussion

Our comparisons with more exact calculations have shown that the PIC approximation gives very good results for a range of important thermodynamic properties of hcp Fe at Earth’s core conditions. A noteworthy example of this is that $`p(V)`$ on the Hugoniot agrees almost perfectly with the more exact results. Our analysis of the free energy into different components makes clear why PIC is so good. The perfect-lattice component is exactly the same in the two approaches. For the harmonic component, the sole requirement for good results is that the logarithmic derivative $`d\mathrm{ln}\overline{\omega }/d\mathrm{ln}V`$ of the geometric mean frequency $`\overline{\omega }`$ be correct. But we have seen that for hcp Fe the PIC $`\overline{\omega }`$ differs from the $`\overline{\omega }`$ given by calculation of the full phonon spectrum by an almost constant factor over a wide range of volumes, so that this requirement is satisfied. The basic reason for this is that the phonon dispersion relations of hcp Fe scale by a uniform factor with changing volume (see Fig. 3 of Alfè et al. (2001)). Finally, the anharmonic component of free energy is very small, and has only a very minor effect on most thermodynamic functions. The reliability of PIC actually requires that anharmonic effects be small, since we have shown that the PIC approximation misrepresents these effects in predicting the wrong sign of the anharmonic free energy.

Surprisingly, even though PIC appears to be so good, we find important discrepancies with the earlier PIC results of Stixrude et al. (1997) and Wasserman et al. (1996). In particular, our calculations of the thermodynamic Grüneisen parameter agree much more closely with the calculations of Alfè et al. (2001). These discrepancies are clearly not due to PIC itself, but must come from other technical differences. We note that in the work of Wasserman et al. (1996), the PIC calculations actually employed a tight-binding representation of the total energy function, the parameters in the tight-binding model being fitted to ab initio calculations. Conceivably, the tight-binding fit might have led to errors.

Even more surprising is that the strong increase with temperature of the axial $`c/a`$ ratio predicted by recent PIC calculations is not reproduced at all by the present PIC work. According to Steinle-Neumann et al. (2001), at the atomic volume of 7.11 Å3, $`c/a`$ increases from 1.63 to 1.75 as $`T`$ goes from 2000 to 8000 K, whereas in the present PIC calculations at the similar volume of 7.0 Å3, $`c/a`$ increases only from 1.594 to 1.610 over the same temperature range. The correctness of the weak $`c/a`$ variation found here is supported by preliminary calculations (Alfé (2002)) using the techniques of Alfè et al. (2001, 2002). The reasons for this discrepancy are completely unclear at present. The discrepancy has major implications for our understanding of the Earth’s inner core, because the recently proposed re-interpretation of the elastic anisotropy of the inner core appears to depend crucially on a strong $`T`$ variation of $`c/a`$. We believe it is highly desirable that this question be re-examined by other research groups.

Our work sheds light on the usefulness of the PIC approximation. Since we have seen that PIC cannot be relied on for anharmonic contributions, it should be regarded as a way of calculating the geometric-mean harmonic frequency $`\overline{\omega }`$. But PIC requires ab initio calculations on a periodic system in which a single atom is displaced from its perfect-lattice site. These are precisely the same calculations that are performed in order to obtain the force-constant matrix used to compute the full phonon spectrum (Alfè et al. (2001, 2002)). For an ab initio method that yields forces on all ions – which includes the pseudopotential and PAW implementations of DFT, among others – the net result of the PIC procedure is to discard all the information contained in the ionic forces, retaining only the variation of total energy with displacement. This suggests that it may be better to use the force information to compute the force-constant matrix and hence the full phonon spectrum, rather than adopting the PIC approach.