text:
# Physical Properties of Iron in the Inner Core

## I INTRODUCTION

The presence and slow growth of Earth’s inner core is one of the most significant manifestations of the dynamics in the interior of our planet. As it is inaccessible to direct observation, an understanding of the physical state of the inner core requires an integrative approach combining results from many fields in the geosciences. Seismology, geo- and paleomagnetism, geo- and cosmochemistry, geodynamics, and mineral physics have advanced our knowledge of the structure and processes in the inner core, revealing many surprises.

Foremost among these have been the discoveries of anisotropy and heterogeneity in the inner core. Long assumed to be a featureless spherically symmetric body, a higher number and higher quality of seismic data revealed that the inner core is strongly anisotropic to compressional wave propagation \[Morelli et al., 1986; Woodhouse et al., 1986\]. Generally, seismic waves travel faster along paths parallel to the Earth’s polar axis by 3-4% compared to equatorial ray paths \[Creager, 1992; Song and Helmberger, 1993\]. The presence of anisotropy is significant because it promises to reveal dynamical processes within the inner core.

Anisotropy is usually attributed to lattice preferred orientation, which may develop during inner core growth \[Karato, 1993; Bergman, 1997\], or by solid state deformation \[Buffett, 2000\]. The source of stress that may be responsible for deformation of the inner core is unknown, although several mechanisms have been proposed \[Jeanloz and Wenk, 1988; Yoshida et al., 1996; Buffett 1996; 1997; Karato, 1999; Buffett and Bloxham 2000\].

An understanding of the origin of inner core anisotropy will require further advances in our knowledge of the physical properties of iron at inner core conditions, and may rely critically on further observations of the detailed structure of the inner core. For example, recent observations indicate that the magnitude of the anisotropy may vary with position: heterogeneity has been observed on length scales from 1-1000 km \[Creager, 1997; Tanaka and Hamaguchi, 1997; Vidale and Earle, 2000\]. Inner core structure may change with time as well. Song and Richards  interpreted apparent changes in travel times of inner core sensitive phases in terms of super-rotation of the inner core with respect to the mantle. Some recent studies have argued for a much slower rotation rate than that advocated originally, or questioned the interpretation of time dependent structure \[Souriau, 1998; Laske and Masters, 1999; Vidale and Earle, 2000\].

The inner core also plays an essential role in the dynamics of the overlying outer core. The anisotropy and long magnetic diffusion time of the inner core may alter the frequency and nature of reversals, and influence the form of the time-averaged field \[Hollerbach and Jones, 1993; Clement and Stixrude, 1995\]. Moreover, important energy sources driving the geodynamo process are associated with solidification of the inner core: the density contrast across the inner core boundary is due to the phase transition from the liquid to the solid, and chemical differentiation during theincongruent freezing of the inner core. Both of these processes provide energy for the dynamo through the release of latent heat \[Verhoogen, 1961\] and the generation of chemical buoyancy \[Braginsky, 1963\]. Other energy sources for magnetic field generation are secular cooling of the Earth, gravitational energy from thermal contraction of the core, radioactive heat generation, and precession \[Verhoogen, 1980; Buffett et al., 1996\].

Both thermal and compositional contributions to the buoyancy depend on the thermal state of the core. The more viscous mantle controls the cooling time scale of the Earth and facilitates the formation of a thermal boundary layer at the core mantle boundary. The heat flux out of the core controls the rate of inner core growth and light element partitioning during this process \[Buffett et al., 1996\]. Conversely, a reliable estimate on temperature in the Earth’s core would advance our understanding of the current thermal state and evolution of the Earth \[Jeanloz and Morris, 1986; Yukutake, 2000\] with important implications for the dynamics of the Earth.

Because the inner core is inaccessible, the study of model systems by theory and experiments is essential. Here we consider the ways in which mineral physics may lend deeper insight into inner core processes and to the origin of its structure, extending previous reviews by Jeanloz  and Stixrude and Brown . We begin with geophysical background on the inner core including recent seismological advances, constraints on the composition, thermal state, and dynamics of the inner core. As our subsequent discussions draw on various experimental and theoretical approaches in mineral physics we then give an overview of recent developments in methods in the following section, focusing on computational mineral physics (a recent review focusing on advances in experiments has been given by Hemley and Mao ). To the extent that the inner core is composed of nearly pure iron, physical properties of this element at high pressure and temperature govern the behavior of the inner core; we consequently review advances in our knowledge of the high pressure physical properties of iron, focusing on crystalline structure, equation of state, and elasticity at both static condition and high temperature. In the final section we examine the implications of these results for inner core temperature, and integrate aspects of elasticity with considerations of the dynamics in the inner core to develop a simple speculative model of polycrystalline structure that explains major aspects of its anisotropy.

## II GEOPHYSICAL BACKGROUND

### II.1 Aggregate Seismic Properties

Lehmann  discovered the inner core by recognizing weak arrivals of PKiKP within the P-wave shadow zone of the core. The amplitudes of these arrivals were sufficient to invoke a discontinuous seismic boundary in the Earth’s core. The P-wave contrast across this boundary was soon established; Birch  and Bullen  argued that the inner core must be solid based on this estimate. The best evidence for inner core solidity comes from studies of inner core sensitive normal modes \[Dziewonski and Gilbert, 1971\]: Earth models with finite shear modulus of the inner core provide a significantly better fit to eigenfrequency observations than those with a liquid inner core. Recent observations of body wave phases involving a shear wave in the inner core (PKJKP, SKJKP, and pPKJKP) \[Okal and Cansi, 1998; Deuss et al., 2000\] support solidity, but are still controversial. The inferred shear wave velocity $`v_S`$ of the inner core is remarkable: it is low compared to the compressional wave velocity $`v_P`$, a property which can also be expressed in terms of the Poisson’s ratio $`\sigma `$.The value of $`\sigma `$=0.44 for the inner core is nearly that of a liquid (0.5), leading to speculation that this region may be partially molten \[Singh et al., 2000\].

In principle density $`\rho `$, $`v_P`$, and $`v_S`$ also depend on depth. However, constraints on the depth dependence of $`\rho `$ and $`v_S`$ are weak. Seismic observations are consistent with an inner core in a state of adiabatic self-compression.

### II.2 Anisotropy

First evidence for deviations from a spherically symmetric structure came from the observation that eigenfrequencies of core sensitive normal modes are split much more strongly than predicted by ellipticity and rotation of the Earth alone \[Masters and Gilbert, 1981\]. Anomalies in PKIKP travel times were initially interpreted as topography on the inner core \[Poupinet et al., 1983\]. Morelli et al.  and Woodhouse et al.  interpreted similar observations of eigenfrequencies and travel times as inner core anisotropy. Observation of differential travel times PKIKP-PKP$`_{\text{BC}}`$ \[Creager, 1992; Song and Helmberger, 1993\] and a reanalysis of normal mode data \[Tromp, 1993\] confirmed that the inner core displays a hexagonal (cylindrical) pattern of anisotropy with a magnitude of 3-4% and symmetry axis nearly parallel to Earth’s rotation axis. For example, PKIKP arrives 5-6 s earlier along polar paths than predicted from radially symmetric Earth models such as PREM \[Dziewonski and Anderson, 1981\]. It is worthwhile pointing out that some of the travel time differences could be due to mantle structure not accounted for in the reference model \[Bréger et al., 1999; Ishii et al., 2002a; 2002b\]. In particular, small scale heterogeneity in the lowermost mantle could be sampled preferentially for select body wave core paths \[Bregér et al., 1999; Tromp, 2001; Ishii et al., 2002b\].

In further investigation deviations from first order anisotropy have been put forward, for example lateral variations in $`v_P`$ of the inner core on length scales ranging from hemispherical differences \[Tanaka and Hamaguchi, 1997; Creager, 1999; Niu and Wen, 2001\], to hundreds of kilometers \[Creager, 1997\], down to a few kilometers \[Vidale and Earle, 2000\]. Radial variations may also exist: weak anisotropy may be present in the uppermost inner core (to a depth of 50-100 km) \[Shearer, 1994; Song and Helmberger, 1995; Su and Dziewonski, 1995\] and strong uniform anisotropy in its inner half \[Song and Helmberger, 1995; Creager, 1999\]. Seismological studies of the inner core are discussed in more detail elsewhere in this volume.

### II.3 Composition

Seismically determined properties of the core may be compared to laboratory measurements under high compression. Measurements of the equation of state show that only elements with an atomic number close to that of iron satisfy the seismic constraints \[Birch, 1964\]. Additional arguments are necessary to uniquely implicate iron \[Jeanloz, 1990\]: iron is one of the most abundant elements in stars and meteorites, much more so than in the portions of the Earth that are directly observable \[Brown and Mussett, 1993\]; and a conducting liquid is necessary in the outer core to explain the existence of a long lived dynamo process that creates Earth’s magnetic field \[Merrill et al., 1996\].

To the degree that we are certain about the main constituent of the core we are also sure that the core contains other lighter elements: pure iron can not satisfy the seismological constraints for both portions of the core. Liquid iron is about 10% too dense to satisfy both the density and bulk modulus in the outer core \[Birch., 1964; Jeanloz, 1979; Brown and McQueen, 1986\] and while solid iron can explain the bulk modulus of the inner core for reasonable temperatures it overestimates the density even for very high temperature (8000 K) \[Jephcoat and Olson, 1987; Stixrude et al., 1997\]. The identity and amount of the light element is still uncertain, but based on cosmochemical arguments hydrogen, carbon, oxygen, magnesium, silicon, and sulfur have been proposed \[Poirier, 1994\], with oxygen and sulfur being the most popular. To infer information on the composition of the core from geochemistry, two questions are of central importance: did the core form in chemical equilibrium \[Karato and Murthy, 1997\] and what are the physical conditions of the core forming event, as pressure and temperature critically determine the partition coefficient of various elements between silicate and metallic melt \[Ito et al., 1995; Li and Agee, 1996; Okuchi, 1997\].

Alternatively, one may use the available seismological information on the current physical state of the outer and inner core ($`\rho `$, $`v_P`$, $`v_S`$) and compare to the physical properties of candidate iron alloys at the appropriate pressure and temperature condition. The compositional space of Fe-X with X any light element has been sparsely sampled in shock wave experiments at the conditions relevant for the core. Only binary compounds in the Fe-S system (pyrrhotite Fe0.9S \[Brown et al., 1984, troilite FeS \[Anderson and Ahrens, 1996\], and pyrite Fe2S \[Ahrens and Jeanloz, 1987; Anderson and Ahrens, 1996\]) as well as wüstite FeO \[Jeanloz and Ahrens, 1980; Yagi et al., 1988\] have been exposed to shock. The data has been extrapolated to inner core conditions \[Stixrude et al., 1997\] and compared to the required elastic parameters (Fig. 1). This analysis indicates that small amounts of either S or O (few atomic percent)would be sufficient to match the properties of the inner core.

Alfè et al. \[2000a; 2000b; 2002\] combined the geophysical approach with a chemical argument. They evaluated the liquid-solid partition coefficients of candidate light elements assuming thermodynamic equilibrium at the inner core boundary. The results show that neither S, Si, nor O alone can satisfy the observed density contrast at the inner core boundary and that a ternary or higher mixture of small amounts of S or Si with O is required.

### II.4 Thermal State

Like the composition, the temperature of the inner core cannot be determined by direct observation. Assuming that the inner core is growing in equilibrium from incongruent freezing of the outer core liquid a knowledge of the melting behavior of iron-rich systems at the pressure of the inner core boundary (330 GPa) would yield an important fixed-point temperature for the construction of whole Earth geotherms. Because the core is not a pure system, the temperature at the inner core boundary should differ from the melting point of pure iron. Freezing point depression in an eutectic system with no solid solution is given by the van Laar equation \[Brown and McQueen, 1982\] which yields a value of 800 K for a melting point of iron of 6000 K and 10% mole fraction of the light element. The value for the freezing point depression must be viewed as highly uncertain, however, since solid solution almost certainly exists at the high temperatures of the core. An independent estimate of the temperature of the core may be obtained by comparing the elastic properties of iron with those seismologically determined. We describe this approach as applied to the inner core below.

Seismic observations do provide constraints on some aspects of the thermal stateof the core. In the outer core the compressional wave velocity $`v_P`$ equals the bulk sound velocity $`v_B=\sqrt{K_S/\rho }`$. In a homogeneous, convecting system, $`v_B`$ and $`\rho `$ are related by adiabatic self-compression. Deviations from this state are characterized by the Bullen  inhomogeneity parameter $`\eta `$ being different from one. $`\eta `$ is defined as

$$\eta =\frac{v_B^2}{\rho g}\frac{\rho }{r},$$
(1)
where $`g`$ the gravitational acceleration, and $`r`$ the radius. $`\eta `$ for the outer core is constrained by seismology to be 1$`\pm `$0.05 \[Masters, 1979\]. This is consistent with (but does not uniquely require) a vigorously convecting outer core, and a resulting geotherm close to an adiabat, characterized by the gradient:

$$\frac{T}{r}=\frac{\gamma g}{v_B^2}T,$$
(2)
where $`\gamma `$ is the Grüneisen parameter. Adopting the value measured for liquid iron at core conditions ($`\gamma `$=1.5) \[Brown and McQueen, 1986\] and a temperature at the inner core boundary of 6000 K, one finds a temperature contrast of 1500 K across the outer core.

The temperature contrast in the inner core is likely to be small. We can place an upper bound on it by assuming that the inner core is a perfect thermal insulator. If the inner core grows through freezing of the outer core its temperature profile will follow the solidus temperature. Based on this assumption Stixrude et al.  estimated the total temperature difference across the inner core to be less than 400 K. Conductive or convective heat loss will further reduce this temperature gradient relative to the insulating case. Conduction is likely to be a very effective way to extract heat from the inner core, such that the temperature profile may fall below an adiabat \[Yukutake, 1998; Buffett, 2000\].

### II.5 Dynamics

Song and Richards  found that the differential travel time of PKIKP-PKP$`_{\text{BC}}`$ for earthquakes in the South Sandwich islands recorded in Alaska increased by 0.3 s over a period of three decades, and concluded that the inner core rotates relative to the mantle by 1/year, a finding that was confirmed qualitatively using global data sets \[Su et al., 1996\]. Creager  showed that part of the signal could be explained by lateral heterogeneity in the inner core, and reassessed the rotation rate to a lower value. Recent years have seen body wave \[Souriau, 1998\] and free oscillations studies \[Laske and Masters, 1999\] that can not resolve inner core rotation, and put close bounds on rotation rate.

If it is present, differential inner core rotation would provide one of the few opportunities for direct observations of the dynamics in Earth’s deep interior. Moreover, differential rotation could have a significant effect on the angular momentum budget of the Earth yielding an explanation of decadal fluctuations in the length of day \[Buffett, 1996; Buffett and Creager, 1999\]. Its origin is not fully understood, but geodynamo simulations produce super-rotation by electro-magnetic coupling with the overlying outer core \[Glatzmaier and Roberts, 1996; Kuang and Bloxham, 1997; Aurnou et al., 1998\]. Gravitational stresses, arising from mass anomalies in the mantle, are also expected to act on the inner core \[Buffett, 1996; 1997\]. These tend to work against super-rotation by gravitationally locking the inner core into synchronous rotation with the mantle. In detail, the interplay between forces driving and resisting super-rotation depend on the rheology of the inner core, which is currently unknown. If the viscosity is sufficiently low, super-rotation may take place, with the consequence that the inner core undergoes continuous viscous deformation in response to the gravitational perturbations.

The interaction of super-rotation with gravitational stresses is just one of many proposed sources of internal deformation in the inner core. The subject has received substantial attention because solid state flow in the inner core can result in lattice-preferred orientation, thought to be essential for producing seismically observed anisotropy.

Other proposed sources of stress in the inner core include:
(a) Coupling with the magnetic field generated in the overlying outer core \[Karato, 1999; Buffett and Bloxham, 2000; Buffett and Wenk, 2001\]. Karato  considered the radial component of the Lorentz force ($`F_r`$) at the inner core boundary which is caused by the zonal magnetic field ($`B_\varphi `$). This is typically the strongest contribution to the Lorentz force in geodynamo models \[Glatzmaier and Roberts, 1995; Kuang and Bloxham, 1997\]. Buffett and Bloxham  argue that the inner core adjusts to $`F_r`$ in a way to minimize steady solid state flow: only weak flow in the inner core is induced that is largely confined to the outermost portion. Considering additional terms to the Lorentz force $`\stackrel{}{F}`$ by including radial components of the magnetic field $`B_r`$ they conclude that the azimuthal term $`F_\varphi `$ which is proportional to $`B_rB_\varphi `$ induces a steady shear flow throughout the inner core.
(b) Thermal convection \[Jeanloz and Wenk, 1988; Wenk et al., 2000a\]. As in any proposed model of inner core flow, the viscosity of the inner core remains an important uncertainty, as does the origin and magnitude of heat sources required to drive the convection.
(c) Aspherical growth of the inner core \[Yoshida et al., 1996\]. Fundamental considerations, based on the expected cylindrical symmetry of flow in the outer core, and detailed geodynamo simulations \[Glatzmaier and Roberts, 1995\] indicate that heat is transported more efficiently in the equatorial plane than along the poles, leading to an inhomogeneous growth rate of the inner core, and internal viscous relaxation. A key question is whether the magnitude of the effect is sufficient to produce lattice preferred orientation. In particular, the resulting strain rates are very small, and may not be sufficient to generate significant polycrystalline texture via recrystallization.

It has also been proposed that polycrystalline texture in the inner core may be acquired during solidification \[Karato, 1993; Bergman, 1997\]. However, if the inner core does experience solid state deformation, by one or more of the mechanisms described above, it is unclear to what extent the texture acquired during solidification would be preserved. It is possible that texture in the outermost portions of the inner core is dominated by the solidification process, whereas lattice preferred orientation in the bulk of the inner core is produced by deformation.

Further progress in our understanding of the composition, temperature, dynamics, and origin of anisotropy in the inner core is currently limited by our lack of knowledge of the properties of iron and iron alloys at high pressures and temperatures. A better understanding of elastic and other properties of iron at inner core conditions can provide a way to test hypotheses concerning the state and dynamics of the inner core.

## III MINERAL PHYSICS METHODS

As elasticity plays a central role in deep Earth geophysics we will emphasize aspects of mineral physics that are directly related to the determination of elastic properties. To gain deeper insight into complex elastic behavior, such as anisotropy, we need to know the full elastic constant tensor at the conditions in the Earth’s center. We will focus on methods based on first-principles quantum mechanical theory, but also briefly review experimental progress, as it relates to comparison and validation of theory. A full review of experimental work has been given recently by Hemley and Mao .

### III.1 Experimental Progress

Determination of the elastic constants of metals in the diamond cell remains a particular challenge as now-standard techniques, such as Brillouin spectroscopy, cannot readily be used for opaque materials. A variety of alternative methods have been developed and applied to study iron at high pressure and ambient temperature. In the lattice strain technique \[Singh et al., 1998a; 1998b\] X-ray diffraction is used to study the strain induced in a polycrystal by uniaxial stress. A full determination of the elastic constant requires the measurement of $`d`$-spacing for many $`(h,k,l)`$ lattice planes, or additional assumptions such as homogeneity of the stress field in the sample \[Mao et al., 1998\]. Other efforts have exploited the relationship between phonon dispersion in the long-wavelength limit and elastic wave propagation: measurements of phonon dispersion have been used to estimate the average elastic wave velocity \[Lübbers et al., 2000; Mao et al., 2001\], and the longitudinal wave velocity \[Fiquet et al., 2001\]. An approximate calibration has been explored in which the zone-center Raman active optical mode is related to the $`c_{44}`$ shear elastic constant by a Brillouin-zone folding argument \[Olijnyk and Jephcoat, 2000\]. Finally, Anderson et al.  by analyzing pressure induced changes in the intensity of X-ray diffraction patterns from hcp iron, extracted a Debye temperature $`\mathrm{\Theta }_D`$, and thus average elastic wave velocity, which they equated with $`v_S`$.

Whereas diamond anvil cell experiments most readily measure properties at ambient temperature, shock wave experiments achieve pressure and temperature conditions similar to those of the core through dynamic compression.By varying the speed of the driver impacting the sample, a set of different thermodynamic conditions are accessed, along a curve in pressure-density space called the Hugoniot. Temperature is not determined directly by the Rankine-Hugoniot equations, it must be measured using special techniques, such as optical pyrometry \[Yoo et al., 1993\] or calculated on the basis of a thermodynamic model \[Brown and McQueen, 1986\]. Using temperature and Grüneisen paramter an adiabatic bulk modulus ($`K_S`$) on the Hugoniot can be determined. The impact of the driver plate on the sample not only sets up a shock in the sample but also in the plate itself. When the shock wave reflects off the back of the impactor, pressure is released and a longitudinal (compressional) sound wave is set up traveling forward through the system of impactor and sample. This has been exploited to determine $`v_P`$; combining $`v_P`$ with $`K_S`$ the corresponding $`v_S`$ can be calculated \[Brown and McQueen, 1986\].

### III.2 Computational Mineral Physics

With the sparse probing of thermodynamic conditions relevant for Earth’s inner core by the experimental methods discussed in the previous section, and the difficulty to obtain information on single crystal elasticity, first-principles material physics methods provide an ideal supplement to experimental study, with all of thermodynamic space accessible, and various approaches to determine elasticity at hand. In the following sections we will introduce the basic principles of calculating such properties.

### III.3 Total Energy Methods

Density functional theory \[Hohenberg and Kohn, 1964; Kohn and Sham, 1965\] provides a powerful and in principle exact way to obtain the energetics of a material with $`N`$ nuclei and $`n`$ interacting electrons in the groundstate (for a review see Lundqvist and March ), with the electronic charge density $`\rho _e(\stackrel{}{r})`$ being the fundamental variable. It can be shown \[Hohenberg and Kohn, 1964\] that ground state properties are a unique functional of $`\rho _e(\stackrel{}{r})`$ with the total (internal) energy

$$E[\rho _e(\stackrel{}{r})]=T[\rho _e(\stackrel{}{r})]+U[\rho _e(\stackrel{}{r})]+E_{xc}[\rho _e(\stackrel{}{r})].$$
(3)
Here $`T`$ is the kinetic energy of a system of non-inter/-acting electrons with the same charge density as the interacting system, and $`U`$ is the electrostatic (Coulomb) energy containing terms for the electrostatic interaction between the nuclei, the electrons, and nuclei-electron interactions. The final term $`E_{xc}`$ is the exchange-corre/-lation energy accounting for many body interactions between the electrons. Density functional theory allows one to calculate the exact charge density $`\rho _e(\stackrel{}{r})`$ and hence the many-body total energy from a set of $`n`$ single-particle coupled differential equations \[Kohn and Sham, 1965\]

$$\{^2+V_{KS}[\rho _e(\stackrel{}{r})]\}\psi _i=\epsilon _i\psi _i,$$
(4)
where $`\psi _i`$ is the wave function of a single electronic state, $`\epsilon _i`$ the corresponding eigenvalue, and $`V_{KS}`$ the effective (Kohn-Sham) potential that includes the Coulomb and exchange-correlation terms from (3). The Kohn-Sham equations are solved self-consistently by iteration. Density functional theory has been generalized to spin polarized (magnetic) systems \[Singh, 1994\].

While density functional theory is exact in principle the exact solution of the Kohn-Sham equations requires the knowledge of the universal form of the exchange-correlation potential which is yet unknown. Approximations for $`V_{xc}`$ however have been very successful. The local density approximation (LDA) \[Lundqvist and March, 1983\] replaces $`V_{xc}`$ at every point in the crystal with the value of a homogeneous electron gas with the same local charge density. This lowest order approximation yields excellent agreement with experiment for a wide variety of materials, but fails for some metals. Most prominently for iron LDA wrongly predicts hcp as the ground state structure for iron at ambient pressure \[Stixrude et al., 1994\]. Generalized gradient approximations (GGA) include a dependence on local gradients of the charge density in addition to the charge density itself \[Perdew et al., 1996\]. GGA yields the correct ground state of iron at ambient pressure and predicts the phase transition from bcc to hcp iron at the experimentally determined pressure \[Asada and Terakura, 1992; Stixrude et al., 1994\].

In addition to total energy it is possible to calculate directly first derivatives of the total energy with first-principles methods. This allows one to determine forces acting on the nuclei and stresses acting on the lattice \[Nielsen and Martins, 1985\].

All-electron, or full potential methods make no additional essential approximations to density functional theory. Computational methods such as the Linearized Augmented Plane Wave (LAPW) method provide an important standard of comparison. All-electron methods are very costly (slow), and are currently impractical for many problems of interest. More approximate computational methods have been developed, which, when applied with care, can yield results that are nearly identical to the all-electron limit.

In the pseudopotential approximation the nucleus and core electrons are replaced inside a sphere of radius $`r_c`$ (cut-off radius) with a simpler object that has the same scattering properties (for a review see Pickett ). The pseudopotential is much smoother than the bare Coulomb potential of the nuclei, and the solution sought is only for the pseudo-wavefunctions of the valence electrons that show less rapid spatial fluctuations than the real wavefunction in the core region or those of the core electrons themselves. The construction of the pseudopotential is non-unique and good agreement with all-electron calculations must be demonstrated.

Iron provides a particular challenge. For example, all-electron results show that pressure-induced changes in the $`3p`$ band are important for the equation of state \[Stixrude et al., 1994\], and so should be treated fully as valence electrons in a pseudopotential approach.

For our work on the high temperature elasticity of hcp iron \[Steinle-Neumann et al., 2001\] we have constructed a Troullier-Martins \[Troullier and Martins, 1991\] type pseudopotential for iron in which $`3s`$, $`3p`$, $`3d`$, and higher electronic states are treated fully as valence electrons. Agreement with all-electron calculations of the equation of state and elastic constants is excellent: for hcp iron the pressure at inner core densities is within 1 % of all-electron (LAPW) results (Fig. 2), and the elastic constant tensor at inner core density is within 2% rms.

The predictions of density functional theory can be compared directly with experiment, or with geophysical observations. For example, by computing the total energy as a function of volume, one obtains the equation of state. This equation of state is static, that is the effects of thermal vibrations are absent. This athermal state is one that is not attainable in the lab where zero-point motion cannot be eliminated. Static properties are often directly comparable to experimental measurements taken at ambient temperature since the effect of 300 K is small for properties such as the density, or the elastic constants. However, for comparison with the earth’s core, thermal effects are essential. Calculation of the effects of temperature from first-principles is more involved because one must calculate the energies associated with atomic displacements, including those that break the symmetry of the lattice.

### III.4 High Temperature Methods

Statistical mechanics provides the tools to deal with material properties at high temperature. The thermodynamic behavior of any physical system is uniquely defined by the so-called fundamental relation, which, for a non-magnetic or Pauli-paramagnetic solid in the canonical ensemble (particle number $`N`$, volume $`V`$, and temperature $`T`$ held constant) takes the form

$$F(V,T)=E(V,T)TS_{el}(V,T)+F_{vib}(V,T)$$
(5)
where $`F`$ is the Helmholtz free energy. The total energy $`E(V,T)`$ is now a function of $`T`$ as well, because we explicitly have to account for thermal excitation of electrons according to Fermi-Dirac statistics. $`S_{el}`$ is the entropy associated with this excitation of the electrons \[McMahon and Ross, 1977\], and $`F_{vib}`$ is the vibrational part of the free energy. $`F_{vib}`$ is derived from the partition function for a system with $`N`$ atoms, which in the classical limit, appropriate at high temperature conditions (significantly above $`\mathrm{\Theta }_D`$) is

$`Z_{vib}={\displaystyle \frac{1}{N!\mathrm{\Lambda }^{3N}}}`$ (6)
$`{\displaystyle 𝑑\stackrel{}{R_1}𝑑\stackrel{}{R_2}\mathrm{}𝑑\stackrel{}{R_N}\mathrm{exp}\left[\frac{F_{el}(\stackrel{}{R_1},\stackrel{}{R_2},\mathrm{},\stackrel{}{R_N};T)}{kT}\right]}`$
and

$$F_{vib}=kT\mathrm{ln}Z_{vib}.$$
(7)
$`Z_{vib}`$ is a 3N dimensional integral over the coordinates of the nuclei located at $`\stackrel{}{R_i}`$ with the electronic free energy $`F_{el}=ETS_{el}`$ uniquely defined by the coordinates of the atoms and $`T`$. $`\mathrm{\Lambda }=h/\sqrt{2\pi mkT}`$ is the de Broglie wavelength with $`h`$ the Planck and $`k`$ the Boltzmann constant, and $`m`$ the nuclear mass.

A naïve attempt to evaluate the integral (III.4) fails because of the large dimensionality, and because most configurations contribute little to the integral. What is required is a search of configuration space that is directed towards those configurations that have relatively low energy. In the particle-in-a-cell (PIC) method and the lattice dynamics method described next, atoms are restricted to vibrations about their ideal crystallographic sites, that is diffusion is neglected. This is not a severe approximation to equilibrium thermodynamic properties at temperatures below the premelting region. Molecular dynamics, described last, in principle permits diffusion, although in practice computationally feasible dynamical trajectories are sufficiently short that special techniques are often required to study non-equilibrium processes.

#### III.4.1 Particle-in-a-Cell

Here the basic approximation motivates the division of the lattice into non-overlapping sub-volumes centered on the nuclei (Wigner-Seitz cell $`\mathrm{\Delta }_{\text{ WS}}`$) with the coordinates of each atom restricted to its cell. A second basic assumption in the PIC model is that the motions of the atoms are uncorrelated. We can expect this approximation to become increasingly valid with rising temperature above $`\mathrm{\Theta }_D`$ and below melting.

If the energy change resulting from moving one particle be independent of the vibrations of the other atoms, the partition function factorizes, and the $`3N`$ dimensional integral is replaced by the product of $`N`$ identical 3-dimensional integrals, which reduces the computational burden tremendously

$$Z_{vib}=\frac{1}{\mathrm{\Lambda }^{3N}}\left(_{\mathrm{\Delta }_{\text{WS}}}𝑑\stackrel{}{R_w}\mathrm{exp}\left[\frac{\mathrm{\Delta }F_{el}(\stackrel{}{R_w},T)}{kT}\right]\right)^N$$
(8)
with

$$\mathrm{\Delta }F_{el}(\stackrel{}{R_w},T)=F_{el}(\stackrel{}{R_w},T)F_{el}(\stackrel{}{R_{w0}},T).$$
(9)
($`N!`$ no longer appears in the prefactor of (8) as the atoms are distinguished by their lattice site.) One evaluates the energetics of moving one atom $`w`$ (the so-called wanderer) with the equilibrium lattice position $`\stackrel{}{R_{w0}}`$ in the potential of the otherwise ideal lattice; this is a mean field approach to the vibrational free energy.

Because large displacements of $`w`$ are included in the integral (up to about 1/2 nearest neighbor distance) the PIC model treats anharmonicity of the vibrations explicitly. The method is computationally very efficient, and can be sped up even more by minimizing the total number of calculations involved in evaluating the integral (8). The angular integrations can be performed efficiently by developing a quadrature which requires evaluation of the integrand along a small number of special directions that are determined by the point symmetry of the lattice site \[Wasserman et al., 1996\].

The cell model has been used for calculations on iron for the thermodynamics of both the hcp and fcc phase \[Wasserman et al., 1996\], and for high temperature elasticity of hcp iron \[Steinle-Neumann et al., 2001\]. It has also been successfully applied to thermoelasticity of tantalum \[Cohen and Gülseren, 2001; Gülseren and Cohen, 2002\].

#### III.4.2 Lattice Dynamics

The calculation of forces on the atoms allows one directly to compute the vibrational frequencies of the material. The dynamical matrix may be calculated row by row by displacing one atom by a small amount from its equilibrium site in a supercell and evaluating the resulting forces on the other atoms. $`F_{vib}`$ is then calculated by performing the appropriate summation over wavevector and phonon branches.

A fundamental approximation in lattice dynamics is that the vibrations of the atoms about their equilibrium are harmonic, only terms that are quadratic in atomic displacements are retained in the expression for the total energy. However, anharmonicity might become important for conditions of the inner core where the temperature is just below the melting point of the material, and considerable effort must be put into anharmonic corrections.

Lattice dynamics has been used extensively over the past few years to address the thermodynamics of hcp iron \[Alfè et al., 2001\], studies on stability of various phases \[Voc̆adlo et al., 2000\], melting \[Alfè et al., 1999\], and studies on core composition \[Alfè et al., 2000a; 2000b; 2002\].

#### III.4.3 Molecular Dynamics

While lattice dynamics and the PIC model essentially evaluate ensemble averages for the thermodynamics of a system, molecular dynamics explores the time evolution of a single realization of the system. In this case thermodynamic properties are calculated as time averages by appealing to the ergodic hypothesis. To obtain the time evolution of the system the forces acting on the atoms are coupled with Newton’s second law. The set of $`N`$ coupled differential equations is then integrated. A large supercell (100 atoms) and long time series (thousands of time steps) are required for convergence of equilibrium properties. Molecular dynamics has been applied to the study of high temperature properties of solid iron including melting \[Laio et al., 2000\] by using a clever hybrid scheme for the electronic structure, as well as for liquid iron \[de Wijs et al., 1998\]. Laio et al.,  combined first-principles total energy and force calculations for a limited number of time steps with a semi-empirical potential fit to the first-principle results.

### III.5 Elastic Constants

Seismic wave propagation is governed by the elastic constants, and the density. In most cases, seismic frequencies are sufficiently high that the adiabatic elastic constants, $`c_{ijkl}^S`$, are relevant. First-principles calculations of the type described here yield the isothermal elastic constants, $`c_{ijkl}^T`$, most directly, and a conversion must be applied \[Davies, 1974\]

$$c_{ijkl}^S=c_{ijkl}^T+\frac{T}{\rho C_V}\lambda _{ij}\lambda _{kl},$$
(10)
where $`C_V`$ is the specific heat at constant volume and

$$\lambda _{ij}=\underset{l,kl}{}\alpha _{kl}c_{ijkl}^T.$$
(11)
The thermal expansivity tensor $`\alpha _{ij}`$ for a hexagonal system has only entries in the diagonal with

$$\alpha _{11}=\alpha _{22}=1/a(a/T)_P,\alpha _{33}=1/c(c/T)_P,$$
(12)
the linear thermal expansivity of the $`a`$\- and $`c`$-axis, respectively.

Under conditions of isotropic pre-stress, elastic wave propagation is governed by the so-called stress-strain coefficients which are defined by

$$\sigma _{ij}=c_{ijkl}^T\epsilon _{kl},$$
(13)
with $`\sigma _{ij}`$ the stress and $`\epsilon _{ij}`$ the strain. Other definitions of the elastic constants have appeared in the literature \[Barron and Klein, 1965\]: one may define elastic constants as the second strain derivatives of the free energy, which are not equivalent to (13) in general. If the pre-stress is isotropic, and if the applied strain is isochoric to all orders, then the free energy may be directly related to the stress-strain coefficients

$$F(V,\epsilon _{ij}^{},T)=F(V,0,T)+\frac{1}{2}c_{ijkl}^T(V,T)\epsilon _{ij}^{}\epsilon _{kl}^{}.$$
(14)
where the primes indicate the deviatoric strain. This relationship has been used to calculate the elastic constants from first-principles calculations of the total energy alone \[Cohen et al., 1997; Steinle-Neumann et al., 1999\]. The elastic constants may also be calculated by appealing to the dissipation-fluctuation theorem which relates the $`c_{ijkl}^T`$ to fluctuations in the shape of the crystal at constant stress. This provides one means of calculating the elastic constants in molecular dynamics simulations \[Parrinello and Rahman, 1982; Wentzcovitch, 1991\].

Elastic wave velocities are related to the elastic constants by the Christoffel equations

$$\left(c_{ijkl}^Sn_jn_k\rho v^2\delta _{il}\right)u_i=0,$$
(15)
where $`\stackrel{}{n}`$ is the propagation direction and $`\stackrel{}{u}`$ the polarization of the wave, $`v`$ the phase velocity, and $`\delta _{il}`$ is the Kronecker delta function. Solving (15) for a given propagation direction yields three velocities, one with quasi-longitudinal polarization ($`v_P`$), and two with quasi-transverse polarizations ($`v_S`$). From the full elastic constant tensor, we may also determine the bulk ($`K_S`$) and shear ($`G`$) moduli using Hashin-Shtrikman bounds \[Watt and Peselnick, 1980\] which give tighter bounds on $`G`$ than the usually used Voigt-Reuss-Hill averages\[Hill, 1963\]. The isotropically averaged aggregate velocities $`v_P`$ and $`v_S`$ can then be calculated by

$$v_P=\sqrt{\frac{K_S+\frac{4}{3}G}{\rho }},v_S=\sqrt{\frac{G}{\rho }}$$
(16)
Usually the Voigt notation is used to represent elastic constants, replacing the fourth rank tensor $`c_{ijkl}`$ with a 6$`\times `$6 pseudomatrix in which pairs of indices are replaced by a single index utilizing the symmetry of the stress and strain tensors:

$`111;`$ $`222;333;`$ (17)
$`23,324;`$ $`13,315;12,216.`$
In this notation the five single crystal elastic constants for a hexagonal system are: the longitudinal elastic constants $`c_{11}`$ and $`c_{33}`$, the off-diagonal elastic constants, $`c_{12}`$ and $`c_{13}`$, and a shear constant $`c_{44}`$. In the following discussions we also refer to another linearly dependent shear constant for comparison with $`c_{44}`$: $`c_{66}=1/2(c_{11}c_{12})`$. To calculate the five elastic constants with a first-principles total energy method, we must evaluate the effect of five different strains on the free energy (14). The bulk modulus and the change in the equilibrium $`c/a`$ ratio with compression provide two pieces of information, yielding two independent combinations of elastic constants. The application of three isochoric strains, of hexagonal, orthorhombic, and monoclinic symmetry, yield the other three pieces of information necessary to obtain the full elastic constant tensor \[Steinle-Neumann et al., 1999\].

## IV PHYSICAL PROPERTIES OF DENSE IRON

### IV.1 Phase Diagram

To the extent that the inner core is composed of pure iron (or nearly pure iron), the phase diagram of iron determines the crystalline structure of the inner core. Despite considerable progress in experimental determination of the phase diagram and melting at pressures approaching those of the core, the stable phase of iron at inner core conditions can not yet uniquely be identified.

This issue is of great geophysical and geochemical importance. First, it is central in our understanding of inner core anisotropy. Different phases of iron show a distinctly different single crystal anisotropy both in magnitude and symmetry \[Stixrude and Cohen, 1995b\]. The enthalpy of various phases and hence the amount of latent heat released at the inner core boundary may depend strongly on the crystalline phase. Finally, the ability to incorporate impurities may be determined by the structure.

Three phases of iron have been unambiguously identified (Fig. 3): the ambient condition ferromagnetic body center cubic phase (bcc, $`\alpha `$) is stable to about 13 GPa and up to 1200 K; at higher temperature the cubic close packed (fcc, $`\gamma `$) phase exists, with bcc reappearing in a narrow stability field ($`\delta `$) just below melting. The magnetic ground state of the fcc phase is a complex spin density wave \[Tsunoda et al., 1993; Uhl et al., 1994\]; in the stability field of fcc iron the local moments appear not to be ordered, however. The hexagonal close packed (hcp, $`\epsilon `$) phase is the high pressure phase, stable to at least 300 GPa at room temperature \[Mao et al., 1990\]. Theory predicts this phase to be non-magnetic (Pauli paramagnet) at core pressures \[Söderlind et al., 1996; Steinle-Neumann et al., 1999\]. Magnetism is not observed experimentally in the hcp phase \[Taylor et al., 1991\], but magnetic moments on the atoms are predicted on the basis of first-principles theory \[Steinle-Neumann et al., 1999\] up to 50 GPa. Moreover, magnetism is observed in epitaxially grown overexpanded lattices of hcp iron \[Maurer et al., 1991\], consistent with theoretical predictions. The possible presence of magnetic states in hcp iron is important for understanding the equation of state (see below) and the phase diagram in the sub-megabar range. The competition between magnetic and non-magnetic contributions to the internal energy, differences in vibrational and magnetic entropy, and differences in volume all contribute to phase stability in iron \[Moroni et al., 1996\].

Two experimental lines of evidence suggest additional stable polymorphs of iron at high pressure and temperature. First, in the shock wave experiment by Brown and McQueen  (see also Brown et al.  and Brown ) there are two discontinuities in $`v_P`$, one at 200 GPa, the other 243 GPa (Fig. 4).

While the one at higher pressure was associated with melting of the sample, the lower one was originally attributed to the hcp to fcc phase transition. With a better characterization of the phase diagram at lower temperature and pressure today this is an unlikelyscenario: the fcc-hcp phase transition ends with a triple point at much lower pressure (Fig. 3) \[Shen et al., 1998; Boehler, 2000\]. Re-appearance of the bcc phase of iron has been suggested as an explanation of the apparent solid-solid phase transformation at 200 GPa \[Matsui and Anderson, 1997\]. However, first-principles theory shows that the bcc phase is mechanically unstable at high pressure and is unlikely to exist as a stable phase \[Stixrude and Cohen, 1995a; Söderlind et al., 1996; Voc̆adlo et al., 2000\]. An alternative interpretation attributes the first discontinuity to the onset of melting, with melting completed only at 243 GPa \[Boehler and Ross, 1997\]. A recent repetition of the experiment of Brown and McQueen  has been unable to resolve whether or not a phase transformation in addition to melting occurs on the Hugoniot \[Nguyen and Holmes, 2001\].

Second, additional phases of iron have been proposed on the basis of static highpressure experiments. X-ray diffraction patterns measured in laser heated diamond anvil cell experiments have been argued to be incompatible with any known iron polymorph \[Saxena et al., 1995; Yoo et al., 1995; Andrault et al., 1997\]. The anomalous signal is subtle, and the proposed structures, a double hexagonal dhcp \[Saxena et al., 1995; Yoo et al., 1995\] and an orthorhombic structure \[Andrault et al., 1997\], are closely related to hcp. Using in-situ X-ray diffraction Shen et al.  found only fcc and hcp phases at pressure and temperature, while the dhcp phase was observed in temperature quenched samples. Andrault et al.  used high strength pressure media which could induce non-hydrostatic conditions \[Boehler, 2000\]. Voc̆adlo et al.  examined the relative stability of proposed high pressure phases using computational ab-initio methods, and found that the orthorhombic structure is mechanically unstable, and that the dhcp phase is energetically less favored than hcp.

Most measurements of the melting temperature of iron from static experiments show reasonable agreement up to a pressure of 100 GPa, where the melting point is 2800-3300 K \[Shen and Heinz, 1998; Boehler, 2000; Hemley and Mao, 2001\]; the data of Williams et al.  yield a significantly higher temperature at 100 GPa (4100 K). The highest pressure datum from diamond anvil cell experiments is at 200 GPa where Boehler  finds the melting point at 3800 K (Fig. 4).

As mentioned above, in shock wave experiments the temperature is not determined directly; based on their dynamic compression data Brown and McQueen  calculated the temperature at the Hugoniot melting point (243 GPa) to be in the range of 5000 to 5700 K (Fig. 4), consistent with subsequent theoretical calculations of the Hugoniot temperature \[Wasserman et al., 1996\]. Measurements of the temperature by optical pyrometry yield a melting point between 235 and 300 GPa with temperatures of 6350 and 6720 K, respectively \[Yoo et al., 1993\].

Two ab- initio calculations of the melting curve of iron have been carried out, yielding inconsistent results (Fig. 4). Alfè et al.  found a melting temperature of 6700 $`\pm `$ 600 K at the inner core boundary by comparing Gibbs free energies of solid and liquid. Laio et al.  determined a considerably lower temperature of 5400 $`\pm `$ 400 K. The origin of these discrepancies are not clear, but may be related to the quality of the anharmonic corrections in the study of Alfè et al. , the semi-empirical potential used to augment first-principles calculations in the study of Laio et al. , or different statistical sampling and runtime adopted in these two studies. In this context, it is worth pointing out that theoretical calculations of the melting temperature are extremely demanding as they involve the comparison of two large numbers (Gibbs free energies of solid and liquid) which must both be calculated to high precision.

### IV.2 Static Equation of State

Experimental measurements and theoretical predictions of the equation of state of non-magnetic hcp iron agree well at core pressures (Fig. 5). At relatively low pressures, however, a discrepancy develops that is larger than that for the bcc phase and what is typical for other transition metals \[Körling and Hḧaglund, 1992; Steinle-Neumann et al., 1999\] suggesting that fundamental aspects of the physics of hcp iron may not be well understood to date.

Although experiment has so far not detected magnetism in hcp iron, recent first-principles theoretical calculations have found stable magnetic states \[Steinle-Neumann et al., 1999\]. These are more stable than the non-magnetic state by more than 10 mRy at low pressure. The most stable magnetic arrangement found so far is one of antiferromagnetic ordering (afmII, Fig. 5) which retains a finite moment up to 50 GPa, well into the pressure region where hcp iron is stable \[Steinle-Neumann et al., 1999\]. Because magnetism tends to expand the lattice, the presence of magnetism reducesthe discrepancy between experimental and theoretical equations of state considerably. It is likely that still more stable magnetic states will be found, and that the ground state is a more complex magnetic structure involving spin-glass like disorder, incommensurate spin density waves, non-collinear magnetism, or a combination of these \[Cohen et al., 2002\].

### IV.3 Static Elastic Constants

A comparison shows considerable disagreement of the single crystal elastic constants (Fig. 6) between experiment \[Mao et al., 1998; Singh et al., 1998b\] and theory \[Stixrude and Cohen, 1995b; Söderlind et al., 1996; Cohen et al., 1997; Steinle-Neumann et al., 1999\].The difference in the longitudinal constants $`c_{11}`$ and $`c_{33}`$ decreases from $``$50% at low pressure to a little morethan 10% at high compression. As pronounced is the discrepancy for the shear elastic constants: the difference in $`c_{44}`$ increases with compression to 30%, and $`c_{66}`$ differs as much as 40% at high pressure, and by a factor of two at the low density data point.

Because the full elastic constant tensor of hcp iron has been measured only by the new lattice strain technique, the large discrepancy between experiment and theoretical prediction prompted a stringent comparison of both methods for a well characterized hcp metal. For the $`5d`$ transition metal rhenium, ultrasonic measurements of elastic constants were not only performed at ambient condition, but also up to 0.5 GPa, constraining the initial pressure slope (Fig. 6) \[Manghnani et al., 1974\]. Calculated elastic constants at zero pressure and their compression dependence \[Steinle-Neumann et al., 1999\] show excellent agreement with the ultrasonic data. In contrast, lattice strain experiments \[Duffy et al., 1999\] do not compare favorably with either result. One of the longitudinal elastic constants ($`c_{33}`$) is considerably smaller (20%), and as in the case for iron the shear constant $`c_{44}`$ shows the largest discrepancy, being larger than the ultrasonic value by 50%. This comparison indicates that additional factors other than elasticity influence the measurements in lattice strain experiments. Subsequently it has been found that strong lattice preferred orientation developed in the experiments for hcp iron \[Wenk et al 2000b; Matthies et al., 2001\], which may cause the assumption of stress homogeneity to be violated.

For aggregate $`v_P`$ and $`v_S`$, the discrepancies between theory and experiment are smaller, but still significant (Fig. 7). Part of the discrepancy between theory and experiment may be attributed to the equation of state. Since theory overestimates the density of hcp iron at low pressure, one expects the elastic wave velocities to be overestimated. The reduction of bulk modulus for the afmII magnetic state will yield lower compressional wave velocity, but for a quantitative comparison of $`v_P`$ and $`v_S`$ information on the full elastic constant tensor for this orthorhombic structure (with nine independent elastic constants) will be needed. The LA phonon data ($`v_P`$) by Fiquet et al.  at high compression and the $`v_S`$ from Anderson et al.  at low pressure appear to be anomalous as they fall below the shock wave data \[Brown and McQueen, 1986\] which, due to thermal effects, would be expected to yield smaller sound velocities than room temperature experiments.

### IV.4 Thermal Equation of State

The Hugoniot provides the strongest constraint on the equation of state of iron at core conditions. First-principles theoretical calculations of the Hugoniot \[Wasserman et al., 1996; Laio et al., 2000; Alfè et al., 2001\] have typically found excellent agreement with that experimentally measured \[Brown and McQueen, 1986\].

A direct comparison of solid state properties between shock wave experiments and computations at inner core pressures is unfortunately not possible, as the Hugoniot is in the stability field for the liquid above 250 GPa as discussed above. We have performed a detailed comparison of the properties of iron at inner core conditions obtained from several first-principles theoretical calculations \[Laio et al., 2000; Alfè et al., 2001; Steinle-Neumann et al., 2001\] and from static \[Dubrovinsky et al., 2000\] and dynamic compression experiments (Fig. 8). The properties of iron determined from these sources agree with each other to within 2% in pressure and to within 10% in bulk modulus at 13 Mg/m3 and 6000 K. This comparison supports the conclusion that the inner core is not pure iron, but that it must contain a small amount of lighter elements. Iron is consistently found to be denser than the inner core, even at a temperature of 7000 K. The bulk modulus of iron, while showing considerably more scatter, appears to be consistent with that of the inner core at a temperature near 6000 K.

### IV.5 Structure

The axial ratio $`c/a`$ of the hexagonal unit cell is important for understanding the elastic anisotropy. Among transition metals at ambient conditions, the value of $`c/a`$ is found to be correlated with the longitudinal wave anisotropy, which can be characterized by the ratio $`c_{33}/c_{11}`$. Large values of $`c/a`$ are associated with a relatively small $`c_{33}`$ and slow P-wave propagation along the $`c`$-axis \[Simmons and Wang, 1971\]. A change in axial ratio through compression or temperature could change the single crystal anisotropy considerably.

Similar to most ambient condition hcp transition metals iron crystallizes with an axial ratio $`c/a`$ slightly below the ideal value ($``$ 1.6) \[Jephcoat et al., 1986; Stixrude et al., 1994\]. It changes little with compression, experiments show a slight decrease with pressure, theory a minor increase.

At ambient pressure, hcp transition metals typically show significant but small changes in $`c/a`$ up to high homologous temperature \[Eckerlin and Kandler, 1971\]. The ratio $`c_{33}/c_{11}`$ shows correspondingly small changes \[Simmons and Wang, 1971\]. The largest change in this ratio is exhibited by titanium which shows a change of 13%. Experiments on hcp iron at higher pressure observed a small increase in $`c/a`$ with temperature in the pressure range of 15-30 GPa and temperatures up to 1200 K (Fig. 9) \[Huang et al., 1987; Funamori et al., 1996; Uchida et al., 2001\]. Similarly, Dubrovinsky et al.  report a $`c/a`$=1.623 at 61 GPa and 1550 K. However, these results can be questioned on the grounds that non-hydrostatic stress may have influenced the temperature dependence of $`c/a`$, a contention supported by hysteresis on a heating and cooling cycle seen in the data of Funamori et al. , and data by Dubrovinsky et al.  at 185 GPa and 1115 K with a small $`c/a`$ (1.585).

First-principles theoretical calculations have predicted a significant temperature induced increase in the $`c/a`$ ratio of hcp iron \[Wasserman et al., 1996; Steinle-Neumann et al., 2001\], which is apparently consistent with the majority of the existing experimental data. For an inner core density of 13.04 Mg/m3 the axial ratio increases from the static value close to 1.6 to about 1.7 at 6000 K (Fig. 10). This implies that at constant density the $`c`$-axis grows at the expense of the $`a`$-axis. The linear thermal expansivities $`\alpha _{11}`$ and $`\alpha _{33}`$ at constant pressure provide another way to represent the change in structural properties. Over a wide range of thermodynamic conditions, the $`a`$-axis is predicted to compress slightly with increasing temperature (Fig. 11).

The temperature-induced change in $`c/a`$ appears to depend on the absolute temperature rather than homologous temperature: the reason that $`c/a`$ of hcp iron reaches such large values at inner core conditions appears to be that higher sub-solidus temperatures may be reached. The origin of the temperature induced increase in $`c/a`$ can be traced to the vibrational entropy, other contributions to the vibrational and electronic free energy have little effect (Fig. 12). Our calculations show that the entropy increases substantially with increasing $`c/a`$. The entropic contribution to the free energy depends on the absolute temperature

$$F_{vib}=E_{vib}TS_{vib}$$
(18)
so that large values of $`c/a`$ become increasingly more favorable energetically at high temperature.

### IV.6 High Temperature Elasticity

As a consequence of the increase in $`c/a`$ the longitudinal anisotropy changes radically with temperature, $`c_{11}`$ becomes larger than $`c_{33}`$ (Fig. 13), with compressional wave propagation in the basal plane being faster than along the $`c`$-axis (Fig. 14). As the $`c`$-axis expands it becomes more compressible, and the corresponding longitudinal modulus, $`c_{33}`$, softens. $`c_{11}`$ in turn increases slightly. The off-diagonal elastic constants are also affected by the temperature-induced change in structure: $`c_{12}`$ increases rapidly with temperature because the basal plane shrinks withincreasing temperature at constant density.

The shear constants $`c_{44}`$ and $`c_{66}`$ show a strong temperature dependence and decrease almost by a factor of four at 6000 K and change order as well (Fig. 13). The velocity of shear waves is considerably smaller at high temperature and the sense of shear anisotropy is reversed (Fig. 14), with the propagation of the (001) polarized shear wave becoming faster along the $`c`$-axis than in the basal plane.

Our calculations imply a shear instability in hcp iron at very high temperature, $`c_{66}`$ approaches zero at 7000 K yielding an upper bound on its mechanical stability or melting point. This is in qualitative agreement with the results by Laio et al.  who predict a shear instability at somewhat lower temperature.

## V PROPERTIES OF THE INNER CORE

### V.1 Aggregate Elasticity

Our current understanding of the physical properties of iron shows that the high Poisson’s ratio of the inner core can be explained by solid-phase elasticity \[Laio et al., 2000; Steinle-Neumann et al., 2001\]. The shear elastic constants are predicted to decrease rapidly with increasing temperature - by a factor of four at 6000 K. As a result, $`G`$ becomes rapidly smaller with increasing temperature (Fig. 15), leading to a Poisson’s ratio of hcp iron that is in quantitative agreement with seismic models of the inner core. These results confirm inferences on the basis of shock wave measurements of $`v_P`$ and estimates of $`v_S`$ at core conditions \[Brown and McQueen, 1986\]. It does not seem necessary to invoke additional mechanisms to explain the high Poissons ratio of the inner core such as viscoelastic dispersion \[Jackson et al., 2000\] or the presence of partial melt \[Singh et al., 2000\].

### V.2 Temperature

Knowledge of the elasticity of iron permits an estimate of inner core temperature that is independent of the iron melting curve and the associated uncertainties related to freezing point depression. The temperature of the inner core is estimated to be 5700 ($`\pm `$ 500) K. At this temperature Poisson’s ratio, $`K_S`$, and $`G`$ of iron simultaneously agree with the properties of seismological models of the inner core (Fig. 15). We infer that the additional light elements that must be present in the inner core have a greater effect on the density through a decrease in mean atomic weight, than they do on the elastic moduli. Further investigations of the elasticity of iron-light element alloys will be needed to test this hypothesis.

Based on an estimate of the Earth’s central temperature, we may construct a core geotherm. We assume that the temperature distribution in the inner and outer core are adiabatic and adopt $`\gamma =1.6`$ \[Wasserman et al., 1996; Alfè et al., 2001\] for the inner core, and $`\gamma =1.5`$ for the outer core (Fig. 4). The latter value is consistent with shock wave data \[Brown and McQueen, 1986\]. With this choice the temperature is 5750 K in the center of the Earth, 5500 K at the inner core boundary, and 4000 K at the core mantle boundary (Fig. 4).

Other estimates of the temperature at the inner core boundary fall in the range of 4500-6000 K, depending on melting point estimate and degree of melting point depression. Our temperature is somewhat higher than the geotherm given by Brown and McQueen  with 5000 K at the inner core boundary, and also higher than the extrapolation from static experiments ($`<`$ 5000 K) \[Boehler, 2000\] or that based on the melting point from Laio et al.  (5400 K). The melting point from Alfè et al. , which is similar to that measured by Yoo et al. , combined with an assumed melting point depression of 700 K \[Alfè et al., 2002\] yields a inner core boundary temperature at 6000 K.

### V.3 Simplified Model of Texture and Anisotropy

The sense of anisotropy that is found at high temperature changes our view of the polycrystalline texture of the inner core and the dynamic processes that may produce it. On the basis of first-principles calculations of the elastic constants, we propose a simple model of the polycrystalline texture of the inner core that explains the main features of its anisotropy \[Steinle-Neumann et al., 2001\]. We find that if 1/3 of the basal planes are aligned with Earth’s rotation axis in an otherwise randomly oriented medium, compressional wave travel time anomalies are well explained (Fig. 16). This model is almost certainly over-simplified. The key element is the tendency for the fast crystallographic direction ($`a`$) to be aligned with the observed symmetry axis of inner-core anisotropy (approximately polar). It is probable that the actual direction and degree of crystallographic alignment will vary with geographic location. Such variations may account for seismological observations of heterogeneity \[Creager, 1997; Tanaka and Hamaguchi, 1997; Vidale and Earle., 2000\].

Important remaining questions include the origin of polycrystalline texture in the inner core, which may have been acquired during solidification, or may have developed subsequently as a result of plastic deformation, as discussed above. If plastic deformation is the prevalent texturing mechanism as we have argued, crystal alignment must depend on the dominant microscopic deformation mechanism in hcp iron at inner core conditions or patterns of growth and recrystallization, and the source of the stress field. While these are currently unknown, our simple model is consistent with the available information.

Candidates for the deformation mechanism include diffusion and recrystallization \[Yoshida et al., 1996; Stixrude and Cohen, 1995b\] dislocation glide \[Wenk et al., 1988; Wenk et al., 2000b; Poirier and Price, 1999\], or a combination of both \[Buffett and Wenk, 2001\]. The type of deformation mechanism is determined by the magnitude of stress and the grain size of the material. When stresses are large, dislocation creep dominates because a large dislocation density facilitates the slip along crystal planes; if the stress is small with a low dislocation density diffusion creep prevails as diffusion of point defects dominates deformation. Small grains facilitate diffusion, larger grain sizes favor dislocation glide. Estimates on the grain size of the crystals in the inner core range from a few mm \[Buffett, 1996\] to the km scale \[Bergman, 1998\]. With the critical grain size determining the deformation regime being in the meter range \[Yoshida et al., 1996\] the large uncertainty in grain size does not provide the means to favor one over the other.

Yoshida et al.  argue that diffusion and recrystallization would result in a texture with the elastically stiffest (fastest) axis coinciding with the direction of the flow, minimizing the strain energy; for iron at high temperature this would tend to align basal planes with the dominant pattern of flow. Among the crystallographic slip planes that may participate in dislocation glide, the basal plane is the slip system that is most easily activated at high pressure according to recent theoretical and experimental work \[Poirier and Price, 1999; Wenk et al., 2000b\]. The predicted high axial ratio at inner core conditions would probably further enhance basal slip as it is typical for ambient hcp metals with large $`c/a`$. As crystals deform in an external stress field basal planes would rotate in the direction of maximum shear, yielding a texture with basal planes aligned with the direction of flow as well. Active recrystallization during slip would tend to enhance the resulting fabric.

Above we have discussed possible sources for stress in the inner core, which allyield distinct flow patterns. Many of them share a flow that is dominant in the polar direction \[Yoshida et al., 1996; Karato, 1999\] which would yield a texture of basal planes aligned with Earth’s rotation axis. The shear flow invoked by Buffett and Bloxham  also yields a texture with the general characteristics of basal planes aligned with the polar direction \[Buffett and Wenk, 2001\].

To obtain the seismic characteristics of such a textural model from single crystal results we convert the elastic anisotropy of iron to differential travel time anomalies. We start by considering an aggregate in which all $`a`$-axes of hcp iron are aligned with the polar axis of the inner core but otherwise randomly oriented. We obtain the elastic properties of this aggregate by averaging elasticity over the solid angle about the $`a`$-axis \[Stixrude, 1998\] which again yields an aggregate elasticity of hexagonal (cylindrical) symmetry. Using the Christoffel equations (15) we can calculate $`v_P`$ as function of the angle $`\xi `$ between the P-wave propagation direction and the pole. Defining the amplitude of the anisotropy as

$$\delta v_P(\xi )/v_{P,av}=(v_P(\xi )v_{P,av})/v_{P,av}$$
(19)
($`v_{P,av}`$ is the average $`v_P`$) we obtain an amplitude of 10% at $`\xi =0^{}`$. This is about a factor of three to five larger than global seismic anisotropy models for the inner core \[Song, 1997; Ishii et al., 2002a\]. It is unlikely that crystallographic alignment in the inner core is perfect. Consequently, in our simple model, the degree of alignment is reduced by a factor of three in order to match the gross amplitude of the seismically observed anisotropy.

To compare directly with seismic observations we compute differential PKIKP travel time anomalies due to inner core anisotropy by

$$\delta t(\mathrm{\Delta },\xi )=t(\mathrm{\Delta })\delta v_P(\xi )/v_{P,av},$$
(20)
where $`\mathrm{\Delta }`$ is the angular distance from source to receiver, and $`t`$ the travel time of the PKIKP phase through the inner core. PKIKP and PKP$`_{\text{BC}}`$ are seen together over a narrow range of distances near $`\mathrm{\Delta }=150^{}`$ ($`t`$=124 s) which we use as a reference distance \[Stixrude and Cohen, 1995b\]. The resulting differential travel time differences are in good agreement both in amplitude and angular dependence with seismic data (Fig. 16).

While we refer to plastic deformation explicitly in the development of this textural model, it is also consistent with other classes of structure. The general character of solidification texture from dentritic growth for pure zinc (like iron at high temperature, a transition metal with high $`c/a`$) \[Bergman et al., 2000\] is in agreement with the model we propose.

## VI CONCLUSIONS AND OUTLOOK

The seismically observed complexity in inner core structure and importance of inner core crystallization for geodynamo processes have initiated considerable interest in the physical state and dynamics of this innermost portion of our planet. Advances in experimental and theoretical methods in mineral physics have made it possible to address important questions regarding core composition, temperature, crystalline structure, and elasticity in the inner core.

Combining geophysical information on core structure and chemical constraints for the partitioning of elements between the solid and liquid promising steps have been made to characterize the light element composition in the core \[Alfè et al., 2000a; 2000b; 2002\] with first results indicating that at least a ternary mixture is needed to satisfy all data. The effect of the light element on inner core elastic properties is expected to be minor, except for density, as evidenced by the ability of pure iron to reproduce inner core aggregate elastic properties \[Laio et al., 2000; Steinle-Neumann et al., 2001\].

Discrepancies in static properties of the high pressure phase of iron, hcp, between experimental data and theoretical predictions suggest that some aspect of the physics of this phase is not well understood to date. Theory indicates the possible presence of magnetic moments on the atoms in the hcp phase \[Steinle-Neumann et al., 1999\]. Experimental and theoretical efforts should be targeted towards a better characterization of magnetic properties at pressures below 100 GPa. Theoretical investigation should focus on the characterization of more complex magnetic structures including non-collinear, disordered, and incommensurate states \[Cohen et al., 2002\]. Also, advances in studies of phonon spectra at high pressure will hopefully result in independent estimates of compressional and shear acoustic wave velocities and of the full elastic constant tensor in the near future.

The thermoelastic properties of iron appear to depend critically on the $`c/a`$ ratio, especially the elastic anisotropy. A careful study of structural parameters as a function of pressure and temperature using experimental and independent theoretical methods could be used to test the predicted reversal in elastic anisotropy of hcp iron at high temperature \[Steinle-Neumann et al., 2001\].

An extension of the melting curve from diamond anvil cell experiments to higher pressures could help to lower current uncertainties in the temperature of the Earth’s core. However, present uncertainties in the melting point of iron are probably comparable to uncertainties in the freezing point depression. Further investigation of the phase stability and elastic properties of iron light element alloys will be important.

###### Acknowledgements.

We greatly appreciate helpful communication of results and preprints from D. Alfè, B. Buffett, L. Dubrovinsky, O. Gülseren, M. Ishii, J. Nguyen, and J. Tromp. This work was supported by the National Science Foundation under grants EAR-9980553 (LS) and EAR-998002 (REC), and by DOE ASCI/ASAP subcontract B341492 to Caltech DOE W-7405-ENG-48 (REC). Computations were performed on the Cray SV1 at the Geophysical Laboratory, supported by NSF grant EAR-9975753 and by the W. M. Keck Foundation.