## Abstract

The bulk properties of iron at the pressure and temperature conditions of Earth's core were determined by a method that combines first-principles and classical molecular dynamic simulations. The theory indicates that (i) the iron melting temperature at inner-core boundary (ICB) pressure (330 gigapascals) is 5400 (±400) kelvin; (ii) liquid iron at ICB conditions is about 6% denser than Earth's outer core; and (iii) the shear modulus of solid iron close to its melting line is 140 gigapascals, consistent with the seismic value for the inner core. These results reconcile melting temperature estimates based on sound velocity shock wave data with those based on diamond anvil cell experiments.

Iron is thought to be the main constituent of Earth's solid inner core and liquid outer core. However, the bulk properties of Fe at such extreme physical conditions remain uncertain, including (i) the Fe melting temperature *T*
_{m} at the pressure of the ICB (330 GPa) (1–6); (ii) the density of liquid Fe at *T*
_{m}, which is needed to determine whether elements lighter than Fe are present in the outer core (7–10); and (iii) the elastic behavior of solid Fe close to the melting line (11–14).

The Fe melting temperature at high pressures has been determined by diamond anvil cell (DAC) and shock wave experiments (1–6). Recent DAC estimates of the melting line extend to 190 GPa (3) and agree with each other (5,6) within 500 K. Shock wave–based estimates are available at ∼240 GPa, but these result in a wider range of possible melting temperatures of 5800 K (1), 6700 K (2), and 6350 K (4). All of the shock wave–based estimates for *T*
_{m} are inconsistent with the extrapolation of the Fe melting line from static DAC experiments (3), which predict a *T*
_{m} of ∼4000 K at 240 GPa.

We calculated the properties of Fe by a method that combines first-principles and classical molecular dynamics (MD) simulations. A correct account of the electronic structure of Fe at the ab initio level is fundamental for an accurate and reliable description of the properties of Fe at Earth's core conditions (13,15–17). Our calculations were based on a finite-temperature extension of density-functional theory within the gradient-corrected local density approximation (GC LDA) (18) and on a pseudopotential description of the valence electron interaction with the ion core (nucleus plus 1s, 2s, and 2p atomic core states) (19). The calculated low-temperature pressure-density curve for hexagonal close-packed (hcp) Fe agrees with the x-ray data (20) (first-principles densities are ∼1% smaller than experimentally determined densities at all pressures).

First-principles quality information on the high-temperature properties of Fe (melting properties, elasticity, diffusion, and Hugoniot equation of state) was obtained in this work by constructing classical potentials with an explicit dependence on the thermodynamic pressure-temperature (*P-T*) condition, exactly mimicking the first-principles MD at that *P-T* point. The potential, which includes non–two-body terms (21, 22) and angular forces (22), is optimized by matching the classical and first-principles forces and stresses with a self-consistent iterative procedure (23). Thermodynamic properties at that*P-T* condition are then extracted from classical MD simulations. The optimal potential (OP) constructed in this way will not be transferable to a different *P-T* condition, where a different potential must be constructed. Our approach is thus different from previous attempts to estimate the melting temperature through classical MD studies (21, 24), where a single potential was used at all *P-T* conditions. Here, the essential coincidence of forces and stresses with first-principles ones at any given *P-T* point guarantees ab initio quality to the results. Our method also differs from a recent determination of*T*
_{m} by first-principles thermodynamic integration (25) because it allows us to overcome the size and time-scale limitations typical of first-principles simulations.

The overall predictivity of our theoretical approach was tested by calculating the melting properties of Al at ambient pressure (26), where well-constrained experimental data are available. Our approach gives a melting temperature (27) for Al at ambient pressure of 940 ± 30 K [experimental data, *T*
_{m} = 933 K (28)], a liquid density ρ at *T*
_{m} of 2.37 ± 0.02 g/cm^{3} [experimental data, 2.375 g/cm^{3} (28)], a density increase at the melting point of Δρ/ρ = 6.1 ± 0.2% [experimental data, 6.6% (28)], and an entropy of melting ΔS = 1.18 (±0.10)k_{B} [experimental data 1.38*k*
_{B} (28)] (*k*
_{B} is the Boltzmann constant), leading, through the Clausius-Clapeyron law, to a melting slope of 71 (±8) K/GPa [experimental data, 65 K/GPa (29)]. For Al, our agreement with experimental data at melting is similar to that obtained with first-principles thermodynamic integration on the same system (30).

Particular attention has been devoted to the analysis of the error bars on our calculated values. These are of two types: (i) errors introduced by the OP modeling of the first-principles dynamics and (ii) errors intrinsic to the first-principles calculations on which the OP construction is based. In the case of Fe at Earth's core conditions, the second type of errors is larger than the first type of errors.

A measure of how accurately the OP modeling reproduces the first-principles thermodynamic observables at *P-T* is given by the difference between the OP and first-principles free energies δF = F_{OP}(P,T) − F_{FP}(P,T). In particular, the error on *T*
_{m} introduced by the OP modeling is related to the errors δ*F*
^{s,l} in the free energy of the liquid (l) and solid (s) phases as |δT_{m}|ΔS ≃ |δF^{l} − δF^{s}| < |δF^{l}| + |δF^{s}|. Because the first-principles potential energy *U*
_{FP} is similar (Fig. 1) to the OP one (*U*
_{OP}), δ*F* can be calculated as δF ≃ 〈U_{OP} − U_{FP}〉 (31), where the average is calculated on an OP MD trajectory. From Fig. 1 (and a similar calculation for the liquid phase), we calculate, assuming ΔS ∼ 1k_{B}, that δT_{m} ∼ 100 K (32).

A second source of error comes from the finite size of the first-principles simulation cell. This has been estimated by repeating the full procedure with cells of increasing size, for selected*P-T* points. For P = 330 GPa, the melting properties were calculated (27) with OPs generated using first-principles cells containing 64 and 128 atoms. The resulting*T*
_{m}, ρ, Δρ, and Δ*S* are unchanged (within the errors associated with the OP modeling) when going from 64 to 128 atom cells. Such convergence is similar to that observed in the case of Al (26).

A third source of errors may arise because of a different accuracy of the GC LDA in describing the solid and the liquid. As is customary in first-principles calculations, errors due to the GC LDA can only be assessed by comparison with experiments. As for errors in*T*
_{m}, we considered that (i) where*T*
_{m} is known, as in Al, the*T*
_{m} calculated within the LDA is accurate; (ii) free-energy differences on Fe at T = 0 K are accurately described by the GC LDA (16); and (iii) others (25) have estimated the error on*T*
_{m} associated with the GC LDA for Fe to be ±300 K. We followed the error analysis presented in (25) and assigned to the error associated with the GC LDA a value of, at most, ±300 K. In conclusion, our estimate of *T*
_{m}should be affected by an overall error of ±400 K. As for errors in density, our calculated values are ∼1% smaller than the experimental ones in the solid and the liquid phases; this error is systematic and does not affect the density differences.

The calculated melting line of hcp Fe [ɛ phase (6)] from 100 to 330 GPa (Fig. 2) agrees, within the error bars, with recent laser-heated DAC experiments (3, 5, 6), with a slightly better agreement with data from (3) and (5), available up to 190 and 150 GPa, respectively. At the ICB pressure (330 GPa), we find that the hcp Fe melts at 5400 K, higher than Boehler's extrapolation of the melting point at 330 GPa (4900 K) (3). As shown below, our melting line is also consistent with Brown and McQueen's shock wave data (1). The melting temperatures recently calculated by first-principles thermodynamic integration (25), the earlier DAC measurements by Williams*et al*. (2), and that extracted from direct temperature measurements in shock wave experiments (2,4) are higher than our results at any pressure.

To understand the origin of the discrepancy between our results and shock wave experiments, we calculated the thermodynamic properties of Fe along the Hugoniot equation of state, defined as (33)(1)where one assumes that Fe is shocked from a state where pressure is *P*
_{0}, atomic volume is *V*
_{0}, and its internal energy per atom (potential and kinetic) is *U*
_{0}, to a state where these are *P*, *V*, and *U*, respectively. We determined the values of *V* and *T*that fulfill Eq. 1 in a set of pressures ranging from 100 to 400 GPa (34). At a given *P*, *T* is varied until *V* and *U*, as calculated with the OP method (35) satisfy Eq. 1. The Hugoniot relation (Eq. 1) can easily be generalized to include the possibility of the coexistence of the liquid and solid phases (33). In particular, if a system at *P-T* conditions lies in the coexistence line of phases A and B (the melting line in the case of solid and liquid coexistence), the volume and the internal energy of the two phases have to be weighted with their relative molar proportions*x*
_{A} and 1 − x_{A}, so that Eq. 1 reads
(2)If Eq. 2 admits a solution with 0 < x_{A} < 1 at pressure *P*, then the system displays phase coexistence. The shocked Fe remains solid up to 200 (±20) GPa (Fig. 3). By further increasing pressure, the additional energy provided by the impact is partially absorbed by the system in terms of heat of melting, which progressively increases the relative molar fraction of melted Fe. In this regime, in fact, Eq. 2 admits a solution with 0 < x_{A} < 1, indicating that the solid (A) and the liquid (B) phases coexist, and the system lies on the melting line. The melting regime along the Hugoniot equation of state ends at a pressure of 280 (±20) GPa, above which Fe is completely melted. The density and the sound velocity (36) of Fe along the Hugoniot relation are compared with experimental data in Fig. 3, A and B. The agreement with experiments is satisfactory in the solid and in the liquid phase. In the regime of phase coexistence (200 to 280 GPa), the sample is not homogeneous, and sound velocities cannot be extracted from elastic theory.

Our simulations also provide an accurate determination of the temperature along the Hugoniot equation of state at all pressures, in the liquid and in the solid regime (Fig. 3C). Unlike pressure and density, the reliability of temperature determinations in a shock wave experiment on metallic systems is still controversial (37, 38). In the solid portion of the Hugoniot equation of state, our values agree with tight-binding calculations (39) and with Brown and McQueen's estimates based on reasonable thermodynamic parameters (1). Instead, temperatures measured in shock wave experiments are systematically higher than ours in the solid portion of the Hugoniot equation. This may be traced to difficulties in directly measuring the temperature of a shocked sample.

Moreover, our calculations suggest a reinterpretation of Brown and McQueen's data for sound velocities (1). We propose that the first kink (at 200 GPa) in Brown and McQueen's data is related to melting rather than to a solid-solid transition. Brown and McQueen's Hugoniot equation of state would then intercept the melting line at ∼200 GPa and ∼4200 K, in agreement with the DAC melting results (3), thus reconciling Brown and McQueen's shock wave measurements with static DAC data. According to this reinterpretation, the second kink observed by Brown and McQueen is not associated with a phase transition, but may rather be a by-product of the phase coexistence between the solid and the liquid, as suggested by MD simulations performed on Ar by Belonoshko (40), and thus be dependent on experimental conditions. This scenario would be confirmed by the recent repetition of Brown and McQueen's experiment, by Nguyen and Holmes (41), where the second kink was not observed. It should be noted, however, that Brown and McQueen's estimate of*T*
_{m}, based on the second kink in the sound velocity curve, is compatible with our melting line, within the mutual uncertainties.

Our theory also provides information on the thermodynamic, elastic, and diffusion properties of Fe at ICB conditions. At 330 GPa, the calculated density decrease (Δρ/ρ) upon melting of Fe is 1.6%. The enthalpy of melting (Δ*H*
_{m}) and the entropy of melting found at 330 GPa are 0.7 × 10^{6}J/kg and 0.86*k*
_{B}, respectively (42). Δ*H*
_{m} is smaller than previously suggested (10, 37), except for estimates based on dislocation theory (43). With our value for Δ*H*
_{m}, the contribution to the geodynamo energy budget of the freezing of Fe in the liquid outer core should be smaller than suggested (10), or alternatively, the age of the inner core should be shorter (44).

Our calculated room-temperature elastic constants (Table 1) are consistent with full ab initio calculations (13) and with the recently revised DAC data (14). We estimated a shear viscosity of the liquid at *T*
_{m} of 1.3 (±0.1) × 10^{−2}Pa, in agreement with recent calculations (17). For the solid, we calculated the bulk and the shear moduli of hcp Fe at*T*
_{m}, and they matched the compressional and the shear wave seismic data for the inner core (Fig. 4). For the compressional wave velocity, we calculated an anisotropy of <10% [only an upper bound can be provided because of the uncertainty in the calculated elastic constants (45)]. If we were to attribute the observed anisotropy of the inner core extracted from the seismic data to the partial alignment of hcp Fe grains along Earth's rotational axis, we would estimate a degree of alignment >30%, in line with recent suggestions (13). The shear modulus estimated from seismic data in the solid inner core is anomalously low for a close-packed phase, leading to suggestions that additional low shear wave phases may be present there, other than hcp Fe (14, 24). Usually, the shear modulus at the melting point of most solid metals shows, at zero pressure, a reduction of <50% from its low-temperature value (46). The inner-core shear modulus (47) displays a 70% reduction with respect to low-temperature measurements (14) and calculations (13) in pure Fe. Our calculated shear modulus of hcp Fe close to melting is compatible with the seismic data (Fig. 4). This decreased shear modulus fits a Born-Durand model of melting (46). According to this model, the isobaric thermal dilatation within the solid phase and after melting is linearly correlated with the decrease of the shear modulus. This empirical law has been verified in a large class of systems, ranging from molecules to metals (46). The shear modulus of Fe closely follows the Born-Durand law (Fig. 5), even if, in comparison with standard metals, compressed Fe melts much closer to the mechanical instability (vanishing shear modulus), which correlates with the small volume increase at melting.

The density of liquid Fe at 330 GPa and 5400 K is estimated to be 12.80 g/cm^{3}. Its difference with the preliminary reference Earth model value for Earth's outer core [12.166 g/cm^{3}(47)] is ∼5 to 6%, consistent with the values currently assumed in geophysical models (∼7%) (10). Moreover, for solid Fe at 330 GPa and 5400 K, we find a density of 13.0 g/cm^{3}, ∼2 to 3% larger than the density of the inner core at the ICB [12.76 g/cm^{3} (47)], supporting the possible presence of lighter elements also in the inner core (7).

↵* To whom correspondence should be addressed. E-mail: guido{at}sissa.it