zotero/storage/9CE39IWY/.zotero-ft-cache

Mechanism of the body-centered cubic iron stabilization under the Earth core conditions
Anatoly B. Belonoshko1*, Timofei Lukinov1, Jie Fu1,2, Jijun Zhao2, Sergio Davis3, Sergei I. Simak4
1Condensed Matter Theory, Department of Theoretical Physics, AlbaNova University Center, Royal Institute of Technology (KTH), 106 91 Stockholm, Sweden.
2Key Laboratory of Materials Modification by Laser, Ion and Electron Beams, Dalian University of Technology, Ministry of Education, 116024, Dalian, China.
3Chilean Committee of Nuclear Energy (CCHEN), Casilla 188-D, Santiago, Chile.
4Department of Physics, Chemistry and Biology (IFM), Linköping University, SE-58183 Linköping, Sweden.
*Correspondence to: anatoly@kth.se
The Earth solid inner core is mostly iron1,2, therefore, the question – what is the structure of iron in the Earth inner core – is central to our understanding of the Core. However, the stable phase of iron in the Core is still unknown. Currently, two major candidates are considered – hexagonal close-packed (hcp) and body centered cubic (bcc) structures. Neither of these structures received unanimous support3,4. Here we demonstrate stability of the bcc phase under conditions in the center of the Core by performing novel constant pressure-temperature ab initio molecular dynamics simulations with varying shape and volume of the computational cell. The bcc phase is stabilized by the discovered unique diffusion mechanism that originates in the low temperature dynamical instability of the bcc phase. It appears that the bcc phase has already been observed in the recent experiments, however, the experimental data was misinterpreted. The diffusion of iron atoms in solid state is quite unique and might allow us to explain both the anisotropy and the low shear modulus of the inner Core.
The Earth solid inner and liquid outer cores consist mainly of iron1. The pressure range of the inner core is from 330 to 365 GPa2. The temperature is known less precisely, since the melting temperatures of iron have not been measured at the core pressures. When extrapolated, the data from shockwave (SW)3,4 and diamond anvil cell (DAC)5 experiments place melting temperature of iron somewhat above 6000 K at the pressure of the inner core outer core (ICOC) boundary (330 GPa). The melting temperature of iron in the center of the Earth then could be close to 7000 K. Some DAC studies6 place this number considerably lower, around 5000 K. Theory supports higher melting temperatures7-9.

Structure and properties of a crystal are intimately related. Therefore, elasticity of the core, its seismic response, density, heat capacity and heat conductivity, to name a few, critically depend on the structure of iron in the core. The interpretation of seismic data, our major source of information about deep Earth, is completely dependent on the phase diagram of iron at the core conditions10-12.
At high pressure and low temperature the stable phase of iron is hcp13. There are two experimental studies that support transition of iron to another, yet unknown phase at pressures above 150-200 GPa and temperatures higher than 3000 K3,14. The SW data3 was interpreted as a transition from the hcp to the bcc phase. The DAC data14 is consistent with stabilization of the bcc phase15. Later SW studies4 claim that the Brown and McQueen data3 is not an evidence for a phase transition. We note, however, that a similar situation was recently resolved for the case of Mo. It was suggested16 that there is no phase transition in Mo contrary to the discovery made by Hixson et al.17 earlier. Recently, it was demonstrated18 that the ‘no-transition’ scenario16 is not consistent with theory and likely the original data by Hixson et al.17 is indeed an evidence for a transition. Recent DAC studies19 claim to measure iron at the core conditions. While the pressure is probably in the core pressure range, the temperature, as can be seen from the width of the peaks of their X-ray spectra, does not probably exceed 4000 to 5000 K20. Besides the changes in the X-ray spectra19 can be explained by the iron carbide formation or the bcc phase stabilization20. Anzellini et al.5 did not see the bcc phase in their experiment at pressures and temperatures that are somewhat lower than the PT conditions of the shock-wave experiment3 where the solid-solid transition was suggested, yet these PT conditions are overlapping with the conditions in Ref. 14.
The dynamic stability21,22 and also the thermodynamic22 stability of the high-temperature high-pressure bcc iron was suggested on the basis of molecular dynamics simulations (the stability of the bcc phase of Fe was suggested23 as an explanation of the results of SW3 experiments). However, recently it was proposed24 that the bcc phase is not even dynamically stable at the pressure-temperature (PT) of the Earth inner core. The authors of Ref. 24 ruled out the dynamical stabilization21,22 of the bcc phase. They showed that the stabilization observed earlier is due to the constant volume and shape of the simulated cell. We note, however, that the size of all previous ab initio simulated cells was small and sometimes very small.
All of the above leads to confusion regarding the stable phase in the core. Theory and experiment provide arguments both supporting and rejecting the bcc stability. While properties of the bcc phase allow us to explain a number of enigmatic properties of the core10-12, the thermodynamic stability of the bcc phase remains a controversial topic. Here we explore the mechanism of the dynamical instability of the bcc phase. The insight into the mechanism allows us to understand what stabilizes the bcc phase and then to answer whether and - if yes - at which conditions the thermodynamic stability of the bcc phase sets in.
It is well known that the bcc phase at high P and low T is unstable and its phonon spectrum includes imaginary frequencies. The instability is due to decrease of the energy (for

calculation of energies and description of molecular dynamics see Methods and Supplementary Materials) when (110) planes are shifted along each other with some wave length (Fig. 1).
The depth of the minimum is related to the kinetic energy corresponding to the temperature of dynamical stabilization of the bcc phase21,22 at constant volume. Accordingly, the
mechanism of instability can be explained by the (110) planes sliding along the [1 1 0 ] direction
to reach the local energy minimum. This is happening at any temperature that is insufficient to escape the minimum. However, if the temperature is high enough, the planes can oscillate around the equilibrium bcc positions, thus stabilizing the bcc phase dynamically. This is possible if the shape of the cell and its volume are maintained constant. Such a motion leads to the appearance of deviatoric stresses24 that would destroy the bcc structure if the volume/shape restrictions are removed. However, as we show below, if the cell is large enough, it may accommodate the strains due to the planes motion without destroying the initial crystal structure. Besides, due to the large amplitude of thermal motion at high temperature the planes are diffused in space therefore the plane shift is dispersed with a finite dispersion length. Therefore, one needs to take good care of the convergence of the simulated cell sizes to prove whether the bcc phase becomes dynamically stable or not.
Currently there are no simulations that would provide properties of the bcc phase at constant pressure and temperature without restrictions on the shape of the computational cell. All of the existing simulations of the bcc phase provide its thermodynamics and/or dynamics compromised by the restrictions and small size of the computational cell. This might lead to wrong conclusions regarding the bcc stability. We performed ab initio MD simulations at constant P and T allowing for deformation of the computational cell with the size up to 1024 atoms. Pressure in all simulations was set to 360 GPa, close to the pressure in the Earth center (364 GPa). We see that (Fig.2) if the size of the cell is 432 atoms (6×6×6 unit bcc cells) then the bcc phase is unstable even at T as high as 7000 K, in full agreement with recent simulations24 where cells with 150 atoms have been simulated at a constant volume. However, if the size is increased to 1024 atoms (8×8×8) the bcc phase remains intact (Fig. 2). Second, if temperature is lowered to 6000 K, the bcc phase (1024 atom supercell) spontaneously transforms into the hcp phase (Fig. 3).
Thus, the size of the computational cell and temperature are critical for stabilization of the bcc phase. The enthalpy difference (Hbcc-Hhcp) that also includes electronic entropy is 19.3 kJ/mole in favor of the hcp phase (Fig. 3). That means that to be stable the vibrational entropy difference (Sbcc-Shcp) should be larger than 2.757 J/mole/K or 0.332 kB (where kB is the Boltzmann constant). This is comparable to the entropy difference of similar phase transitions in, e.g., Ti and Zr25,26,27. Direct calculation of the free energy difference for systems of size 1024 atoms and higher (we do not know at which size the properties of the bcc phase will converge) are not currently feasible. Recently, the difference of vibrational free energies was computed28 using the SCAILD (self-consistent ab initio lattice dynamics) method29. This difference is 27.1

kJ/mole at 360 GPa and 7000 K. This means that the bcc phase is more stable than the hcp phase by 7.3 kJ/mole. The authors of Ref. 29, however, considered enthalpies for the ideal bcc and hcp lattices and that caused confusion. In our study we computed the enthalpies without any simplifying assumptions. If the enthalpy and entropy differences do not change with temperature, the bcc would stabilize at 5200 K, which is considerably lower than nearly all current estimates of the inner core temperatures that are essentially temperatures of melting of the core material. Other estimates of entropy difference (Suppl. Materials), based solely on our results, range from 2.87 J/K/mole (low limit) to 8.0 J/K/mole (upper limit) (Suppl. Materials). The entropy difference calculated using SCAILD falls in this range and actually is closer to the lower limit.
What is the mechanism behind the stabilization of the bcc phase? If one looks at the snapshots in Fig. 4 (the whole movie is included in Supplementary materials) one can see that atoms in the bcc phase diffuse in avalanches along the <110> directions. These are the directions of the soft mode that make the bcc structure unstable at T=0 K. At high temperature, the sliding of the planes is transformed into sliding of the pieces of the planes simply because at high temperature the planes do not exist as such – atoms are substantially out of the plane. Such a possibility makes the bcc phase stable dynamically while the increase in entropy (see Fig. 4S in Supplementary Materials) makes the bcc phase stable thermodynamically.
Such a motion in small cells destroys the structure. Bigger cells are apparently capable of relaxing the stress associated with the strain caused by the motion however such strains lead to quite a remarkable mechanism of diffusion, to our knowledge, never observed before. Such diffusion might have profound consequences for the formation of lattice preferential orientation in the inner core30 and explanation of seismic data including anisotropy11,12,30,31, extremely low shear resistance10,30 and high attenuation of seismic waves.
A careful examination of the structure of the bcc phase (note, that this is the first time we are looking at the structure of bcc iron calculated ab initio without constraints on shape and volume) tells us that the atoms in the bcc structure are effectively 11 coordinated (Fig. 4). The first coordination sphere contains the 8 atoms however the 6 atoms in the second shell contribute about half of them to the first coordination sphere due to considerable amplitude of thermal motion. This is the coordination that was measured in the recent shock wave experiments32. However, the authors of the experimental paper32 assumed that 11-coordinated atoms should rather be a feature of close packed structures (formally 12 coordinated). The authors32 explained the coordination 11 by thermal motion in the hcp phase. It is thermal motion indeed – but in the bcc phase that makes iron atoms 11 coordinated. Thus, we can suggest that the stability of the bcc phase at the core conditions has already been observed but misinterpreted.
To conclude, for the first time we demonstrated the dynamical stability of the bcc phase at high pressure and temperature. The bcc phase becomes thermodynamically stable according to first-principles calculations. The stability increases on increasing the size of the computational sample and is likely to converge at larger sizes rendering full convergence difficult to achieve.

However, already at the size of 1024 atoms with periodic boundary conditions the bcc phase at 7000 K and 360 GPa is thermodynamically stable. The discovered mechanism of stabilization – diffusion along the planes with the soft phonon modes at 0 K- is novel but could likely be the same in other metallic refractory systems exhibiting phonon softening. The mechanism might allow one to explain a number of seismic observations related to anisotropy and low shear modulus in the Earth inner core. Indeed, the possibility to efficiently move atoms within the same crystalline structure allows one to efficiently remove the stress – and that makes the bcc phase very soft with a very low resistance to shear. The bcc phase is very anisotropic11,12. The diffusion allows easy texturing of iron in response to any stress. Most of the earlier ab initio molecular dynamics studies of the bcc phase with small (less than 1024 atoms) size are likely inadequate. The unique feature of the Fe bcc phase such as high-temperature self-diffusion even in an ideal lattice might be responsible for the formation of large-scale anisotropic structures needed to explain the Earth inner core anisotropy11,30,31. It is also important to mention that the diffusion might interfere in the X-ray measurements and should be considered when choosing the appropriate experimental technique.
Methods: The energies and forces were calculated within the framework of the frozen-core allelectron Projector Augmented Wave (PAW) method33, as implemented in the Vienna Ab initio Simulation Package (VASP)34,35,36. The energy cut-off was set to 400 eV. Exchange and correlation potentials were treated within the generalized gradient approximation (GGA)37,38. Fourteen electrons were considered as valence, therefore any overlapping of core states at high pressures and temperatures was avoided. The finite temperatures for the electronic structure and force calculations were implemented within the Fermi-Dirac smearing approach39. All necessary convergence tests were performed, the electronic steps converged within 0.0001 meV/atom. The molecular dynamics runs have been performed in the NPT (N – number of atoms, P – pressure, T –temperature) ensemble for a given pressure and temperature. We used supercells up to 1024 atoms. Tests have shown that Γ-point is sufficient at these sizes. The timestep was set to 1 fs and the runs continued for up to 18,000 timesteps. The averages were computed over last 3,000 to 5,000 timesteps depending on the run and the errors were estimated by the blocking technique.
References and Notes:
1. Birch, F. Elasticity and constitution of the Earth's interior. J. Geophys. Res. 57, 227-286 (1952).
2. Dziewonski, A. M. & Anderson, D. L. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297-356 (1981).
3. Brown, J. M. & McQueen, R. G. Phase transitions, Grüneisen parameter, and elasticity for shocked iron between 77 GPa and 400 GPa. J. Geophys. Res. 91, 7480-7494 (1986).
4. Nguyen, J. H. & Holmes, N. C. Melting of iron at the physical conditions of the Earth's core. Nature 427, 339-342 (2004).

5. Anzellini, S., Dewaele, A., Mezouar, M., Loubeyre, P. & Morard, G. Melting of Iron at Earth’s Inner Core Boundary Based on Fast X-ray Diffraction. Science 340, 464-466 (2013).
6. Boehler, R. Temperatures in the Earth's core from melting-point measurements of iron at highstatic pressures. Nature 363, 534-536(1993).
7. Belonoshko, A. B. & Ahuja, R. Embedded-atom molecular dynamic study of iron melting. Phys. Earth. Planet. Inter. 102, 171-190 (1997).
8. Alfé, D., Gillan, M. J. & Price, G. D. The melting curve of iron at the pressures of the Earth's core from ab initio calculations. Nature 401, 462-464 (1999).
9. Belonoshko, A. B., Ahuja, R. & Johansson, B. Quasi-ab initio molecular dynamic study of Fe melting. Phys. Rev. Lett. 84, 3638-3641 (2000).
10. Belonoshko, A. B., Skorodumova, N. V., Davis, S., Osiptsov, A. N., Rosengren, A. & Johansson, B. Origin of the low rigidity of the Earth's inner core. Science 316, 1603-1605 (2007).
11. Belonoshko, A. B., Skorodumova, N. V., Rosengren, A. & Johansson, B. Elastic anisotropy of the Earth's inner core. Science 319, 797-799 (2008).
12. Mattesini, M., Belonoshko, A. B., Buforn, E., Ramirez, M., Simak, S. I., Udias, A., Mao, H.K. & Ahuja, R. Hemispheric anisotropic patterns of the Earth's inner core. PNAS 107, 9507 (2010).
13. Hemley, R. J. & Mao, H.-K. In situ studies of iron under pressure: New windows on the earth's core. Int. Geol. Rev. 43, 1-30 (2001).
14. Mikhaylushkin, A. S., Simak, S. I., Dubrovinsky, L., Dubrovinskaia, N., Johansson, B. & I. A. Abrikosov, Pure iron compressed and heated to extreme conditions. Phys. Rev. Lett. 99, 165505 (2007).
15. Belonoshko, A. B., Derlet, P. M., Mikhaylushkin, A. S., Simak, S. I., Hellman, O., Burakovsky, L., Swift, D. C. & Johansson, B. Quenching of bcc-Fe from high to room temperature at high-pressure conditions: a molecular dynamics simulation. New J. Phys. 11, 093039 (2009).
16. Nguyen, J. H., Akin, M. C., Chau, R., Fratanduono, D. E., Ambrose, W. P., Fat'yanov, O. V., Asimow, P. D. & Holmes, N. C. Molybdenum sound velocity and shear modulus softening under shock compression. Phys. Rev. B 89, 174109 (2014).
17. Hixson, R. S., Boness, D. A., Shaner, J. W. & Moriarty, J. A. Acoustic velocities and phase transitions in molybdenum under strong shock wave compression. Phys. Rev. Lett. 62, 637 (1989).
18. Lukinov, T., Simak, S. I. & Belonoshko, A. B. Sound velocity in shock compressed

molybdenum obtained by ab initio molecular dynamics. Phys. Rev. B 92, 060101(R) (2015).
19. Tateno, S., Hirose, K., Ohishi, Y. & Tatsumi, Y. The structure of iron in Earth’s Inner Core. Science 330, 359 –362 (2010).
20. Dubrovinsky, L., Dubrovinskaia, N. & Prakapenka, V. Is iron at the Earth’s core conditions hcp-structured? Fisica de la Terra 23, 73-82 (2011).
21. Belonoshko, A. B., Ahuja, R. & Johansson, B. Stability of the body-centered-cubic phase of iron in the Earth’s inner core. Nature 424, 1032-1034 (2003).
22. Vocadlo, L., Alfé, D. & Price, G. D. Possible thermal and chemical stabilization of bodycentered-cubic iron in the Earth’s core. Nature 424, 536-539 (2003).
23. Ross, M., Young, D. A. & Grover, R. Theory of the iron phase diagram at earth core conditions. J. Geophys. Res. 95, 21713-21716 (1990).
24. Godwal, B. K., Gonzales-Cataldo, F., Verma, A. K., Stixrude, L. & Jeanloz, R. Stability of iron crystal structures at 0.3-1.5 TPa, Earth and Planetary Science Letters 409, 299-306, (2015).
25. Nishitani, S. R., Kawabe, H. & Aoki, M. First-principles calculations on bcc-hcp transition of titanium. Mater. Sci. Eng. A312, 77-83 (2001).
26. Petry, W., Heiming, A., Trampenau, J., Alba, M., Herzig, C., Schober, H. R. & Vogl, G. Phonon dispersion of the bcc phase of group-IV metals. I. bcc titanium. Phys. Rev. B 43, 1093310947 (1991).
27. Petry, W., Heiming, A., Trampenau, J., Alba, M., Herzig, C., Schober, H. R. & Vogl, G. Phonon dispersion of the bcc phase of group-IV metals. II. bcc zirconium, a model case of dynamical precursors of martensitic transitions. Phys. Rev. B 43, 10948-10961(1991).
28. Niu, Z. W., Zeng, Z. Y., Cai, L. C. & Chen, X. R. Study of the thermodynamic stability of iron at inner core from first-principles theory combined with lattice dynamics. Phys. Earth. Planet. Inter. 248, 12-19 (2015).
29. Souvatzis, P., Eriksson, O., Katsnelson, M. I. & Rudin, S. P. The self-consistent ab initio lattice dynamical method. Comp. Mater. Sci. 44, 888-894 (2009).
30. Deguen, R. Structure and dynamics of Earth’s inner core. Earth Planet. Sci. Lett. 333-334, 211-225 (2012).
31. Romanowicz, B., Cao, A., Godwal, B., Wenk, R., Ventosa, S. & Jeanloz, R. Seismic anisotropy in the Earth’s innermost inner core: testing structural models against mineral physics predictions. Geoph. Res. Lett. 43, 93-100 (2016).
32. Ping, Y. et al., Solid Iron Compressed Up to 560 GPa. Phys. Rev. Lett. 111, 065501 (2013).

33. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953-17979 (1994).
34. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758-1775 (1999).
35. Kresse, G. & Hafner, J. Ab-initio molecular-dynamics for open-shell transition-metals. Phys. Rev. B 48, 13115-13118 (1993).
36. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15-50 (1996).
37. Wang, Y. & Perdew, J. P. Correlation hole of the spin-polarized electron-gas, with exact small-wave-vector and high-density scaling. Phys. Rev. B 44, 13298-13307 (1991).
38. Perdew, J. P., Chevary, J. A., Vosko, S. H., Jackson, K. A., Pederson, M. R., Singh, D. J. & Fiolhais, C. Atoms, molecules, solids, and surfaces- applications of the generalized gradient approximation for exchange and correlation. Phys. Rev. B 46, 6671-6687 (1992).
39. Mermin, N. D. Thermal Properties of the Inhomogeneous Electron Gas. Phys. Rev. 137, A1441 (1965).
40. Sherman, D. M. Electronic, entropy and the high-pressure stability of bcc iron. AIP Conf. Proc. 309, 895 (1994).
41. Hellman, O., Steneteg, P., Abrikosov, I. A. & Simak, S. I. Temperature dependent effective potential method for accurate free energy calculations of solids. Phys. Rev. B 87, 104111 (2013).
42. Antolin, N., Restrepo, O. D. & Windl, W. Fast free-energy calculations for unstable hightemperature phases. Phys. Rev. B 86, 054119 (2012).
43. Frenkel, D. & Ladd, A. J. C. New Monte Carlo method to compute the free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres. J. Chem. Phys. 81, 3188 3193 (1984).
44. Bouchet, J., Mazevet, S., Morard, G., Guyot, F. & Musella, R. Ab initio equation of state of iron up to 1500 GPa. Phys. Rev. B 87, 094102 (2013).
45. Aquilanti, G. et al. Melting of iron determined by X-ray absorption spectroscopy to 100 GPa PNAS 112, 12042-12045 (2015).
46. Belonoshko, A. B. & Dubrovinsky, L. S. A simulation study of induced failure and recrystallization of a perfect MgO crystal under non-hydrostatic compression: Application to melting in the diamond-anvil cell. Amer. Miner. 82, 441-451 (1997).
47. Sutton, A. P. & Chen, J. Long-range Finnis–Sinclair potentials. Phil. Mag. Lett. 61, 139–146 (1990).

48. González, D. & Davis, S. Fitting of interatomic potentials without forces: A parallel particle swarm optimization algorithm. Comput. Phys. Commun. 185, 3090–3093 (2014).
Acknowledgments: Computations were performed using the facilities provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputing Center in Linköping (Sweden). We also wish to thank the Swedish Research Council (VR) for financial support (grants 2013-5767 and 2014-4750). A.B.B., J.Z. and J.F. acknowledge funding from the National Magnetic Confinement Fusion Program of China (2015GB118000) and the China Scholarship Council. S.I.S. acknowledges Linköping Linnaeus Initiative for Novel Functional Materials (LiLi-NFM) and the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linköping University (Faculty Grant SFO-Mat-LiU No. 2009 00971).
Figure Captions.
Fig. 1. Modeling N-point transverse acoustic phonon mode in ideal bcc Fe at 360 GPa. The mode
is the shortest possible wave of atomic displacements in the <110> direction with the [1 1 0 ]
polarization. Left panel: (001) projection of the equilibrium bcc atomic positions at different electron temperatures. Atomic positions in the ideal bcc structure are marked with crosses. Red squares, green rhombuses, blue triangles, and magenta circles indicate the equilibrium atomic positions at 2000, 6000, 8120 and 10000 K, respectively. Their distances from the ideal bcc positions (crosses) correspond to the amplitude of the N-point mode minimizing the electron free energy. This amplitude is equivalent to the atomic displacement from the ideal bcc positions and it is indicated by arrows of the same color as for the equilibrium atomic positions at corresponding electron temperatures. The size of the arrows is proportional to the N-wave amplitude minimizing the electron free energy. Accordingly, one may see that the amplitude decreases with increasing temperature. Right panel: Difference between the electron free energy at a particular electron temperature and the electron free energy of the ideal bcc structure at the same temperature. The difference is given as a function of the amplitude of the atomic displacement according to the N-point transverse acoustic phonon mode in the <110> direction
with the [1 1 0 ] polarization (N-wave shift). The data are presented for the same temperatures as
in the left panel and they are given with respect to the electron free energy of the ideal bcc structure. One may see that the atomic displacements of the equilibrium atomic positions from those in the ideal bcc structure decrease the electron free energy but the energy gain decreases as the temperature increases.
Fig. 2. Time evolution of the initially cubic box sizes (a-red, b-green, and c-black; in the initial cubic structure a=b=c) during molecular dynamics runs. The evolution is computed at pressure 360 GPa and temperature 7000 K. The upper set of three curves shows the evolution for the cell composed by 8×8×8 bcc unit cells (each cell with 2 atoms) that is 1024 atoms in total. The lower set shows the results of the same simulation for the 6×6×6 bcc unit cell with 432 atoms. The 432

atom cell demonstrates instability developing between 7,000 and 8,000 timesteps. After that, the cell shape is deformed and the bcc structure transforms into the hcp one. On the contrary, the 1024 atoms cell shows no such instability. The computations for the 1024 atoms cell were continued longer than the similar runs that demonstrated the instability for the cell with 150 atoms24.
Fig. 3. Results of molecular dynamics simulations in the NPT ensemble at 360 GPa and high temperatures for different sizes of computational cells starting from the bcc and hcp structures. The upper panel shows time evolution of enthalpy that also includes the electronic entropy term and the lower panel time evolution of the electronic entropy. The bcc structure readily transforms into hcp structure at 6000 K even at 1024 atom cell (squares). Both enthalpy and electronic entropy demonstrate drop at the transition from the bcc to the hcp phase. The hcp to bcc enthalpy and entropy differences are computed as the difference between the enthalpy and entropy averages for the hcp and bcc structures at 7000 K and 1024 atom cells.
Fig. 4. Structural data on bcc iron from molecular dynamic NPT simulation at 7000 K and 360 GPa. The A and B parts show 3-dimensional distribution of probability density of atoms as a function of distance from the central atom. The density is averaged over the atoms and the last 8,000 timesteps of the run of 18,000 timesteps duration. Note, that A does not show the distribution of atoms around their equilibrium positions but rather the density of the vectors that connect this atom and the central atom. This density enters the radial distribution function normally shown as a one-dimensional plot. The first neighbor clouds in A have the shape squeezed in the direction to the central atom (shown as a small sphere for clarity), while the second neighbors are elongated in the same direction. The part B shows the projection of A on the 110 plane that is the plane that cuts the bcc unit cell along the face diagonal. One can see that half of the second neighbors density contributes to the coordination number that is the total coordination number is 8 + 6/2 = 11. Parts C and D show two projections of five snapshots of the computational cell. The five snapshots are shown in the picture with atoms colored according to the red-white-blue sequence corresponding to time progression. The snapshots are separated by 200 timesteps. The positions of the atoms are averaged using the 1000 timesteps running window to average out the high-frequency vibrations. The projections along X (part C) and Z (part D) axes are shown. One can clearly see the diffusion along the (110) planes. For the full movie of the simulation and diffusion see Supplementary Materials video.
Author contributions:
A.B.B. designed the study, performed calculations and wrote the paper; T.L. performed calculations and wrote the paper; J.F. and S. D. performed calculations and discussed the paper; J.Z. discussed the paper; S.I.S. performed the calculations and wrote the paper; all authors discussed the results and contributed to paper writing.

Competing financial interests: The authors declare no competing financial interests. Correspondence: anatoly@kth.se Supplementary Materials: Supplementary Text Figures S1-S5 with captions References (33-48) Caption for Movie S1 Movie S1

Supplementary Information for
Mechanism of the body-centered cubic iron stabilization under the Earth core conditions
Anatoly B. Belonoshko1*, Timofei Lukinov1, Jie Fu1,2, Jijun Zhao2, Sergio Davis3, Sergei I. Simak4
correspondence to: anatoly@kth.se This file includes:
Supplementary Text Captions for Figs. S1-S5 Figs. S1-S5 Caption for Movie S1 Other Supplementary Materials for this manuscript includes the following: References 33-48 (included in the manuscript file) Movie S1
1

Supplementary Text Entropy change on transition from hcp to bcc phase

It has been established that bcc iron is dynamically unstable at high pressure and low temperature40. There are several methods that allow calculation of entropy of a dynamically unstable phase29,41,42,43. We note, that thermodynamic integration43 is not
technically possible in the present case because we see that the dynamical instability sets
in only at the supercell size of 1024 atoms which is very time-consuming to simulate. The
bcc phase stability increases on size, as we see from comparison of 432 and 1024 atom
simulations, and most likely is not converged at the size of 1024 atoms. Therefore, one
needs an alternative to the thermodynamic integration method. Recently, the entropy difference (∆Shcp-bcc) has been calculated28 using the self-consistent ab initio lattice dynamics (SCAILD) approach29 (with the very same PAW-potential with 14 valence
electrons that is used in this study). We do not think that the SCAILD is the best possible
method. However, it is better than MD in case of small dynamically unstable systems
because the molecular dynamics of small bcc structures provide us with incorrect data.
That is the bcc structure does not possess the correct dynamics if the size of the structure
is small. At the same time, SCAILD provides the entropy in reasonable agreement with experiment29. The difference of vibrational entropies28 T⋅∆Svibbcc-hcp = 27.1 kJ⋅mole-1, according to SCAILD. This makes the bcc phase by -7.8 kJ more stable than the hcp
phase. This suggests that the stability of the bcc phase sets in at around 5200 K at the pressure of 360 GPa. The ∆Svibbcc-hcp = 3.9 J⋅mole-1⋅K-1. Let us see whether we can verify this number independently. Using the Clausius-Clapeyron equation for melting of both
the hcp and bcc phases we can write

dT/dPliq-hcp=(Vliq-Vhcp)/(Sliq-Shcp)

dT/dP liq-bcc =(V liq -V bcc )/(S liq -S bcc )
Considering that melting temperatures of the bcc and hcp are close44 and, correspondingly, the gradients of their melting curves are close too, we can obtain that

where

Sbcc-Shcp = (Vbcc-Vhcp)/A, A=dT/dP melting

on melting at pressure 360 GPa. The most recent melting curve of iron44 computed up to
the pressure 15 Mbar gives the dT/dP melting = 8.9 K/GPa at the pressure 360 GPa. The
volume difference can be calculated from our MD runs for the bcc and hcp phases at 360
GPa and 7000 K. This is equal to 0.032 cm3/mol. Therefore, we are getting ∆Sbcc-hcp = 3.6 J/mole/K.
This entropy includes both electronic and vibrational terms. Our MD runs give the electronic part Sbcc-Shcp = 0.73 J/mole/K. Therefore, the ∆Svibbcc-hcp = 2.87 J/mole/K which is somewhat less than the SCAILD result (3.9 J/mole/K) yet sufficient to stabilize
the bcc phase. Note, that this is a lower bound for the entropy.

2

It is of interest to compare the entropy change for iron to the entropy change for similar transitions where experimental data exist. The measured phonon spectra were used to calculate the vibrational entropy in quasiharmonic approximation for Zr26 and Ti27 at the temperatures 1156 K and 1188 K, respectively, for the hcp and bcc phases. For example, for Zr ∆Sbcc-hcp = 2.17 J/mole/K and for Ti ∆Sbcc-hcp = 2.41 J/mole/K. One can expect that ∆Svibbcc-hcp increases with temperature.
Therefore, all existing estimates provide consistent range of vibrational entropy changes. To stabilize the bcc phase over the hcp phase one needs ∆Svibbcc-hcp >19.3 kJ/mole / 7000 K = 2.75 J/mole/K. Even if we assume that the SCAILD calculations overestimate the ∆Svibbcc-hcp by more than 30 percent (which is highly unlikely) the entropy difference will be still sufficient to stabilize the bcc phase. Moreover, it is likely that the SCAILD approach underestimates the entropy difference. The discovered by us the unique mechanism of diffusion in the bcc phase (see Fig. 4 and Movie S1) considerably increases the number of available microstates and, therefore, increases the entropy.
Thus, we conclude that the increase of entropy in the bcc iron phase is sufficient to stabilize the bcc phase thermodynamically at the pressure of 360 GPa, likely starting from the temperature around 5000 K. This temperature is estimated using the most reliable based on the same method as used in this paper estimate of entropy in bcc and hcp iron28. Remarkably, this number is very close to the extrapolated temperature of melting if the extrapolation is done using the measurements providing the so-called low melting curve (e.g. Ref. 45). It was suggested, that these experiments do not measure melting temperatures but rather the temperature of the onset of dynamic recrystallization7,46. Indeed, if the possibility of the recrystallization was considered the melting temperatures from experiments increased significantly5 and became consistent with the high melting temperature as suggested by theory7. The recent measurements that are very sensitive to appearance of liquid-like structures insist on the low melting curve45. We want to emphasize that during re-crystallization such liquid-like regions are abundant because the grain boundaries at high T are essentially liquid-like10. Besides, the discovered in this study self-diffusion might be interpreted as the signature of the liquid phase. Also, the diffusion in the bcc phase makes its shear resistance likely very low. There is the evidence that the inner core material has very low shear resistance. Stabilization of the bcc phase might allow us to explain it.
Impact of system size and duration of runs
In order to evaluate the impact of increase of the size of the simulated system well beyond 1024 atoms and continuation of the run well beyond dozens of thousands of time steps one needs to turn to a ‘quasi ab initio’ approach9. In such an approach, a semiempirical potential is fitted to results of ab initio calculations. The choice of the data set which parameters are fitted to depends on what is actually required from the potential. In our case we require that the potential would reasonably well describe the energy of iron phases at the pressure of 360 GPa. Therefore, we have chosen to fit the potential to the run with 1024 atoms at T=6000 K and P=360 GPa starting from the bcc phase (Fig. 3). In that run the bcc phase transforms into the hcp phase. Selecting configurations from those parts of the run were the thermal equilibrium is established and both the bcc and hcp
3

configurations are presented we fitted the potential in the Sutton-Chen form47. This potential was used earlier being fitted to a different set of data9 in a wide range of volumes. Here we restrict our task to match as close as possible the results of ab initio simulations in the narrow range of volumes.
The potential is described as follows.

𝐸𝐸𝑖𝑖

=

1 2

∑𝑁𝑁𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎
𝑗𝑗=1,𝑗𝑗≠𝑖𝑖

𝜑𝜑

<EFBFBD>𝑟𝑟𝑖𝑖𝑖𝑖<EFBFBD>

+

𝐹𝐹(𝜌𝜌𝑖𝑖),

(1)

Where

𝜑𝜑<EFBFBD>𝑟𝑟𝑖𝑖𝑖𝑖<EFBFBD> = 𝜀𝜀 <20>𝑟𝑟𝑎𝑖𝑎𝑖𝑖𝑖<F09D9196>𝑛𝑛,

(2)

𝐹𝐹(𝜌𝜌𝑖𝑖) = −𝜀𝜀𝜀𝜀<F09D9C80>𝜌𝜌𝑖𝑖,

𝜌𝜌𝑖𝑖

=

∑𝑁𝑁𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎
𝑗𝑗=1,𝑗𝑗≠𝑖𝑖

<EFBFBD>𝑟𝑟𝑎𝑖𝑎𝑖𝑖𝑖<EFBFBD>𝑚𝑚

.

(3)

Ei is the potential energy of atom i, Natoms is the number of atoms in the system, ϕ is the pair-wise potential between atoms i and j, F(ρi) is the embedded energy function and rij is the distance between atom i and atom j. The total potential energy of the system is then

𝐸𝐸𝑡𝑡𝑡𝑡𝑡𝑡

=

∑𝑁𝑁𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎
𝑖𝑖=1

𝐸𝐸𝑖𝑖.

(4)

In the present work, the Sutton-Chen potential is fitted to the energies calculated using ab initio molecular dynamics. The configurations are randomly selected from the configurations produced by ab initio MD simulations under the condition of 360 GPa and 6000 K except the first 1200 configurations (see Fig. S1). Half of the input configurations are bcc structured and half are hcp structured. The fit is performed determining the minimum of the following function

𝑓𝑓 = <20>∑𝑁𝑖𝑖𝑁=𝑐𝑐1𝑐𝑐𝑐𝑐𝑐𝑐<F09D9190>𝐸𝐸𝑖𝑆𝑖𝑆𝑆𝑆 − 𝐸𝐸𝑖𝐷𝑖 𝐷𝐷𝐷𝐷𝐷<F09D90B7>2

(5)

where ESC is the potential energy of configurations calculated by the Sutton-Chen potential, EDFT is energy of the same configuration calculated from ab initio methods,
Nconf is the number of configurations considered in this fitting procedure. The method used for minimization is the particle swarm optimization algorithm48. The obtained
parameters are n = 8.7932 , m = 8.14475, ε = 0.0220225 eV , a = 3.48501 Å , and C = 28.8474. The root mean square error of the fit is 6.2 meV/atom. Fig. S1 shows the
comparison of the energies calculated by the fitted Sutton-Chen potential and ab initio
method.

4

Fig. S1. Energies of configurations as computed in ab intio MD and by the embeddedatom Sutton-Chen potential with fitted parameters. The configurations marked by filled triangles are used for fit. Open circles show the energies computed by the SC potential. The rest of configurations has not been used for the fit and provided for the comparison. The decrease of energies at around 3000th configuration is related to the transition from the bcc to the hcp phase. Using this potential we first performed molecular dynamics runs on the system of the same size (1024 atoms). Fig. S2 shows the results.
5

Fig. S2. Volume of the computational cell during the MD run. Simulations are performed in NPT ensemble at N=8×8×8 bcc unit cells (1024 atoms), pressure of 360 GPa and temperature as indicated on the panel. The timestep is equal 1 femtosecond. If the system is stable, volume remains nearly constant for the most of the run as is the case at 7000 K (upper panel). If the volume remains constant for some period of time and then converge to a different value and again remains constant we have instability of the initial configuration (in this case bcc) and a new stable (in this case hcp) configuration. This is what happens at T=6000 K (lower panel) and T=6800 K (middle panel). The bcc is stable at 7000 K and transforms to hcp at temperature 6800 K and lower at this size in full agreement with ab initio MD simulations.
Quite remarkable, the potential shows the same behavior as the ab initio system (Fig. 3). First, the bcc phase remains stable for 200 picoseconds (200,000 timesteps, see the upper panel in Fig. S2). In ab initio MD we could not perform such long runs and we saw that ab initio system is stable for 18 picoseconds. Second, the bcc phase transforms to the hcp phase at 6000 K, exactly as does the ab initio system. Moreover, even the kinetics of the transition is reasonably close (the time from the start of the simulation until the transition is 3 picoseconds in the ab initio run and 8 picoseconds in the EAM run). Third, we can now put more tight constraints on the temperature of the instability. The EAM bcc is unstable at 6800 K (the panel in the middle of Fig. S2). Therefore, the EAM description allows us to reasonably close (just slightly short of ‘perfectly’) reproduce the behavior of the ab initio system.
6

Fig. S3. Volume of the computational cell during the MD run. Simulations are performed in NPT ensemble at N=32×32×32 bcc unit cells (65536 atoms), pressure of 360 GPa and temperature as indicated on the panel. The timestep is equal 1 femtosecond. In the very beginning of each run (first 1000 timesteps) volume increases as the system gets equilibrated. If the system is stable, volume remains nearly constant for the rest of the run as is the case at 7000 K (upper panel) and 5500 K (middle panel). If the volume remains constant for some period of time and then converge to a different value and again remains constant we have instability of the initial configuration and a new stable (in this case hcp) configuration. This is what happens at T=5000 K (lower panel).
We see that increase of the cell from 8×8×8 bcc unit cells (1024) atoms to 32×32×32 bcc unit cells (65536 atoms) makes the bcc phase stable already at the T=5500 K and possibly even lower (but not lower than 5000 K). Thus, increase of the size leads to the increase of the temperature range of stability by about 1500 K.
Now entropies can be computed for large sizes using the Sutton-Chen potential. First, we compute the velocity autocorrelation function (VACF) in the NVE ensemble where volume is chosen to correspond to the pressure of 360 GPa and temperature of 7000 K. The temperature is scaled for the first 10,000 timesteps to equilibrate the system and then truly NVE calculations performed without any interference in the dynamics of the system. The Fourier transform of the VACF gives the density of vibrational states (phonon density of states or PDOS). From PDOS calculation of entropy is performed
7

assuming the model of independent oscillators. This is rather standard procedure and is implemented in a number of molecular dynamics software packages for example in LAMMPS. We computed entropies for the bcc and hcp phases as function of size. The results are presented in Fig. S4.
Fig. S4. Size dependence of entropy of the hcp and bcc phases computed at the temperature of 7000 K and pressure of 360 GPa. The entropy is computed as described in the text using 500,000 timesteps runs. While at small supercell sizes the entropy difference is small (the self-diffusion is not considerable enough to contribute to the low frequencies of the calculated PDOS) at large sizes the low frequencies states become highly populated (as we see it from the PDOS plot in Fig. S5) and that leads to considerable increase of entropy in the bcc phase. The magnitude of the entropy change is similar to that of a solid-liquid transition. This is rather natural since in both the liquid and bcc phase the diffusion is activated. The system size is provided in the lattice constant, i.e. N means the size N×N×N unit cells. The unit cell contains 2 atoms both for the hcp and bcc phase. For the system size N=40 the number of atoms is equal to 128,000.
Some of the calculated PDOS are shown in Fig. S5. Please note that the PDOS in bcc phase are slightly different even in such large systems as 128,000 and 16,000,000 atoms. The low frequencies in the 8×8×8 bcc phase are somewhat more populated than in the hcp phase but this size is not yet sufficient to provide converged PDOS in the bcc
8

phase. The calculations of PDOS with big sizes become rather challenging even for the Sutton-Chen potential. Thus, the memory required for the 200×200×200 was about 5 Terabytes. This number increases as cube of the linear size.
Fig. S5. Phonon density of states (PDOS) computed for hcp and bcc phases of different sizes. The systems have been equilibrated first at P=360 GPa and T=7000 K. After that the simulations have been performed in NVE ensemble and the trajectories were accumulated for 500,000 timesteps. The number of atoms in bcc system were 1024, 128,000, and 16,000,000. The number of atoms in hcp system was 85,184 (here hcp was computed in orthogonal installation with 4 atoms in unit cell; the number of translations was chosen to get a close to cube shape). Simulations with smaller and larger number of atoms in the hcp phase demonstrated that this number is sufficient for convergence. Convergence in the bcc phase requires much larger sizes. Much higher population of the states with a low frequency in the bcc phase is responsible for the larger entropy of the bcc phase (Fig. S4) as compared with the hcp phase. From these calculations we see that ∆Svibbcc-hcp ≈ 8.0 J/mole/K. This is likely the upper limit of the entropy change on this transition since the oscillators are not really independent and, besides, the diffusion is not exactly an oscillation with low frequencies. Nevertheless, it is obvious that the diffusion greatly increases the entropy of the bcc
9

phase. The lower limit for the ∆Svibbcc-hcp = 2.87 J/mole/K. This range of entropy difference renders the bcc phase at 7000 K and 360 GPa more stable than the hcp phase. Moreover, the increase of the system size leads to the increase of the bcc stability range by at least 1300 K. Therefore, it is clear that the ab initio bcc phase would stabilize at about 6000 K and possibly even lower. This is definitely within the temperature range of the Earth Inner Core. Movie S1 Motion in the 1024 atom bcc cell at 360 GPa and 7000 K obtained by molecular dynamics. The coordinates are averaged over the moving window of 1000 timesteps to remove the high-frequency oscillations. Spontaneous diffusion events along (110) directions are clearly visible.
10