Temperature of the inner-core boundary of the Earth: Melting of iron at high pressure from first-principles coexistence simulations
TTemperature of the inner-core boundary of the Earth: Melting of iron at highpressure from first-principles coexistence simulations
Dario Alf`e ∗ Materials Simulation Laboratory, and London Centre for NanotechnologyDepartment of Earth Sciences, and Department of Physics and Astronomy,UCL, Gower Street, London WC1E 6BT, United Kingdom (Dated: November 1, 2018)The Earth’s core consists of a solid ball with a radius of 1221 Km, surrounded by a liquid shellwhich extends up to 3480 Km from the centre of the planet, roughly half way towards the surface(the mean radius of the Earth is 6373 km). The main constituent of the core is iron, and thereforethe melting temperature of iron at the pressure encountered at the boundary between the solidand the liquid (the ICB) provides an estimate of the temperature of the core. Here I report themelting temperature of Fe at pressures near that of the ICB, obtained with first principles techniquesbased on density functional theory. The calculations have been performed by directly simulatingsolid and liquid iron in coexistence, and show that and at a pressure of ∼
328 GPa iron melts at ∼ ±
100 K. These findings are in good agreement with earlier simulations, which used exactlythe same quantum mechanics techniques, but obtained melting properties from the calculation ofthe free energies of solid and liquid Fe [1, 2, 3].
PACS numbers:
The study of iron under extreme conditions has along hystory. In particular, numerous attempts havebeen made to obtain its high pressure melting proper-ties [4, 5, 6, 7, 8, 9, 10, 11, 12]. Experimentally, Earth’score conditions can only be reproduced by shock wave(SW) experiments, in which a high speed projectile isfired at an iron sample, and upon impact high pressureand high temperature conditions are produced. By vary-ing the speed of the projectile it is possible to investi-gate a characteristic pressure-volume relation known asthe Hugoniot [13], and even infer temperatures, althougha word of caution here is in order, as temperature es-timates are often based on the knowledge of quantitieslike the constant volume specific heat and the Gr¨uneisenparameter, which are only approximately known at therelevant conditions [10]. If the speed of the projectileis high enough, the conditions of pressure and temper-ature are such that the sample melts, and it is there-fore possible to obtain points on the melting curve, ofcourse with the caveat mentioned above about temper-ature measurements. An alternative route to high pres-sure high temperature properties is the use of diamondanvil cells (DAC), in which the sample is surrounded by apressure medium and statically compressed between twodiamond anvils. In DAC experiments pressure and tem-peratures can be directly measured, and therefore thesetechniques should in principle be more reliable to inves-tigate melting proeprties. Unfortunately, in the case ofiron is it not so, and there is a fairly large range of resultsobtained by different groups [4, 5, 6, 7, 8, 9].An alternative approach used for the past ten years orso has been to employ theory –and in particular quan-tum mechanics techniques based on density functionaltheory– to calculate the high pressure melting curve ofiron. A number of groups have used different approaches to the problem. Our own strategy has been to calculatethe Gibbs free energy of solid and liquid iron, and thenobtain the melting curve by imposing their equality forany fixed pressure. We obtained a melting temperatureof ∼ ∼ ∼ N V E ensemble, i.e. constant number of atoms a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b N , constant volume V and constant internal energy E .In the N V E ensemble, for each chosen volume V thereis a whole range of energies E for which solid and liquidcan coexist for long time; the average temperature andpressure along the simulation then provide a point on themelting curve. If the energy E is above(below) the rangefor which coexistence can be maintained, the system willcompletely melt (solidify), and the simulation does notprovide useful melting properties informations. It shouldbe pointed out that any finite system will eventually meltor solidify if simulated for long enough, due to sponta-neous fluctuations. However, melting(solidification) re-sulting from a too high(low) value of E typically appearon much shorter time scales.The present calculations have beed performed withdensity functional theory with the generalised gradientapproximation known as PW91 [22] and the projectoraugmented wave method [24, 25] as implemented in the VASP code [23]. An efficient extrapolation of the chargedensity was employed [27]. Single particle orbitals wereexpanded in plane waves with a cutoff of 300 eV, and Iused the finite temperature implementation of DFT asdeveloped by Mermin [28]. These settings are exactlyequivalent to those used in our previous work [1, 2, 3],so the melting properties obtained here will be directlycomparable to those early ones. The simulations havebeen performed on hexagonal closed packed (hcp) cellscontaining 980 atoms (7 × × × − eV. With these prescriptions the drift in theconstant of motion was ∼ . ∼ E , provided tothe system by assigning different initial velocities to the FIG. 1: (Color online) Snapshot of a DFT molecular dynamicssimulation showing solid and liquid iron in coexistence. Thesimulation cell contains 980 atoms.FIG. 2: Temperature (upper panel) and pressure (lowerpanel) for a simulation of solid and liquid iron in coexistence. atoms. The simulation with the highest value of E com-pletely melted after ∼ E solidified after ∼
11 ps. Among the otherthree, one melted after ∼
14 ps, one after ∼
24 ps, whilethe last one has remained in coexistence for the wholelength of ∼
25 ps. However, most of these simulationswere coexisting for long enough, so that useful meltinginformation from the period of coexistence could actuallybe extracted in almost all cases.A snapshot of a simulation with solid and liquid incoexistence is show in Fig. 1 [29].In Fig. 2 I display the temperatures and the pressurescorresponding to the simulation that remained in coexis-tence for the whole 25 ps length, which provides a meltingpoint ( p, T ) = (328 ± , ±
100 K). It is interest-ing to notice a temperature excursion in the simulationafter ∼
15 ps, which lasts for ∼ ac-cidental melting (or freezing), even if the internal energy E is within the range of coexistence. This problem is FIG. 3: Temperature (upper panel) and pressure (lowerpanel) for a simulation of solid and liquid iron in coexistence.The system eventually completely solidifies, with a drop inpressure and an increase in temperature due to release of vol-ume and latent heat of fusion respectively.FIG. 4: As in Fig. 3, but here the system eventually melts. mitigated by the use of large simulation cells, and there-fore this is one of the reasons why large simulation cellsare needed in conjunction with the coexistence approach.In Fig. 3 I show a simulation that eventually solidi-fied, however, as mentioned above, coexistence was main-tained for a long period, and the information gatheredby the central part of the simulation can still be usedto obtain a point on the melting curve, and the result is( p, T ) = (324 ± , ±
100 K), which is consistentwith the previous point.Similarly, in Fig. 4 I show one of the simulations thateventually melted, and by taking the average tempera-ture and pressure from the central part of the simulationI get ( p, T ) = (331 ± , ±
100 K), which is alsoconsistent with the other previous two points.The c/a value used in the simulations was fixed at1.6, which resulted in slightly non-hydrostatic condi-tions. To study the effect of non-hydrostaticity, andthat of the relatively small size of the simulation cells, I have performed simulations using a classical embed-ded atom model (EAM), adapted to deliver results veryclose to the present ab-inition techniques [16]. I haveperformed simulations on cells containing 7840 atoms(14 × × ∼
100 K. Theab-initio non-hydrostatic conditions with c/a = 1 . c/a = 1 .
65 the simulations are almost exactly un-der hydrostatic conditions [30]. The effect of the non-hydrostaticity with c/a = 1 . ∼
100 K, so that the combined effects ofnon-hydrostaticity and small size cancel each other. Afinal check was performed by repeating the simulationsusing 62720 atom-cells (28 × × ∼ ∼
330 GPa is thereforeconfirmed to be in the region of ∼ − FIG. 5: (Color online) Comparison of melting curve of Fefrom present calculations with previous experimental and abinitio results; blue filled square: DFT coexisting simulationfor over 25 ps; blue open squares: DFT coexisting simulationsending up to complete solids or liquids, results are gatheredfrom the region of coexistence; black heavy solid curve: DFTmelting curve from Ref. [2]; long-dashed blue curve: DFT-OPM results of Ref. [15]; light purple solid curve: DFT-EAMresults of Ref. [14]; filled red circles: DFT-EAM “corrected”results (see text); black chained and maroon dashed curves:DAC measurements of Refs. [4] and [5]; green open diamonds:DAC measurements of Ref. [6]; green plus: DAC measurementof Ref. [9]; green filled triangle: DAC measurement of Ref. [7];black stars, black open circle and pink filled diamond: shockexperiments of Refs. [12], [10] and [11]. Error bars are thosequoted in original references.
The next question that one might ask is how accu-rate is DFT for this problem. We argued in our pre-vious work that DFT should indeed be quite reliable,mainly because solid and liquid iron have very similarstructures, so that possible errors would largely cancelbetween the two phases. However, we also pointed outthat DFT does not seem to reproduce the zero tempera-ture pressure-volume equation of state of hcp iron com-pletely correctly, possibly underestimating the pressureby ∼ . ∼
150 K [2]. This would bring the meltingtemperature of Fe at ICB condition to ∼ ∼ . ∗ Electronic address: [email protected][1] D. Alf`e, M. J. Gillan and G. D. Price, Nature, , 462(1999).[2] D. Alf`e, G. D. Price, and M. J. Gillan, Phys. Rev. B, ,165118, (2002).[3] D. Alf`e, G. D. Price, M. J. Gillan, Phys. Rev. B, ,045123 (2001).[4] Q. Williams, R. Jeanloz, J. D. Bass, B. Svendesen, T. J.Ahrens, Science , 181 (1987).[5] R. Boehler, Nature , 534 (1993).[6] G. Shen, H. Mao, R. J. Hemley, T. S. Duffy and M. L.Rivers, Geophys. Res. Lett. , 373 (1998).[7] A. P. Jephcoat and S. P. Besedin, Phi. Trans. R. Soc.Lond. A , 1333 (1996).[8] S. K. Saxena, G. Shen, and P. Lazor, Science, , 405(1994).[9] Y. Ma, M. Somayazulu, G. Shen, H. K. Mao, J. Shu, R. J.Hemley, Phys. Earth Planet. Int., , 455 (2004).[10] J. M. Brown and R. G. McQueen, J. Geophys. Res. ,7485 (1986).[11] J. H. Nguyen, and N. C. Holmes, Nature
339 (2004).[12] C. S. Yoo, N. C. Holmes, M. Ross, D. J. Webb and C.Pike, Phys. Rev. Lett. , 3931 (1993).[13] J.-P. Poirier, Introduction to the Physics of the Earth’sInterior , Cambridge University Press, Cambridge (1991).[14] A. B. Belonoshko, R. Ahuja, and B. Johansson, Phys.Rev. Lett., , 3638 (2000).[15] A. Laio, S. Bernard, G. L. Chiarotti, S. Scandolo and E.Tosatti, Science , 1027 (2000).[16] D. Alf`e, G. D. Price, M. J. Gillan, J. Chem. Phys, ,7127 (2002).[17] D. Alf`e, Phys. Rev. B, , 064423 (2003).[18] L. Voˇcadlo and D. Alf`e, Phys. Rev. B, , 214105 (2002).[19] T. Ogitsu, F. Schwegler, F. Gygi, G. Galli, Phys. Rev.Lett., , 175502 (2003)[20] S. A. Bonev, F. Schwegler, T. Ogitsu, G. Galli, Nature , 669 (2004).[21] D. Alf`e, Phys. Rev. Lett., , 235701 (2005).[22] Y. Wang and J. Perdew, Phys. Rev. B , 13298 (1991);J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson,M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev.B , 6671 (1992).[23] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996).[24] P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994).[25] G. Kresse and J. Joubert, Phys. Rev. B , 1758 (1999).[26] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[27] D. Alf`e, Comp. Phys. Comm. , 31 (1999).[28] N. D. Mermin, Phys. Rev. , A1441 (1965).[29] Figure realised with the xcrysdens software: A. Kokalj,Comp. Mater. Sci. NpH ensemble (constant pressure p and enthalpy H ), using the algo-rithm developed by E.R. Hernandez, J. Chem. Phys.115