Elastic properties of body-centered cubic iron in Earth’s inner core

The solid Earth’s inner core (IC) is a sphere with a radius of about 1300 km in the center of the Earth. The information about the IC comes mainly from seismic studies. The composition of the IC is obtained by matching the seismic data and properties of candidate phases subjected to high pressure (P) and temperature (T). The close match between the density of the IC and iron suggests that the main constituent of the IC is iron. However, the stable phase of iron is still a subject of debate. One such iron phase, the body-centered cubic phase ( bcc ), is dynamically unstable at pressures of the IC (330–364 GPa) and low T but gets stabilized at high T characteristic of the IC (5000–7000 K). So far, ab initio molecular dynamics (AIMD) studies attempted to compute the bcc elastic properties for a small (order of 10 2 ) number of atoms. The mechanism of the bcc stabilization cannot be enabled in such cells and that has led to erroneous results. Here we apply AIMD to compute elastic moduli and sound velocities of the Fe bcc phase for a 2000 Fe atom computational cell, which is a cell of unprecedented size for ab initio calculations of iron. Unlike in previous ab initio calculations, both the longitudinal and the shear sound velocities of the Fe bcc phase closely match the properties of the IC material at P = 360 GPa and T = 6600 K, likely the PT conditions in the IC. The calculated density of the bcc iron at these PT conditions is just 3% higher than the density of the IC material according to the Preliminary Earth Model. This suggests that the widely assumed amount of light elements in the IC may need a reconsideration. The anisotropy of the bcc phase is an exact match to the most recent seismic studies.


I. INTRODUCTION AND MOTIVATION
The Earth's solid inner core (IC) [1] is made of almost pure iron, as follows from seismic and other geophysical data [2].This was well established even before the experimental data on iron at extreme pressure temperature (PT), approaching those in the IC, became available [3,4].However, the iron phase and the amount of light impurities (pure iron is heavier than the IC material [5]) remain debated.While the experimental data are controversial as concerns the phase state of iron under PT conditions of the IC (pressures from 330 to 364 GPa and temperatures-from various sources-between 5000 and 7000 K), the latest seismic data [6] and the bodycentered cubic (bcc) paradigm [7][8][9][10][11][12][13] are a perfect match.The seismic data tell us there is an innermost inner core (IMIC) with a radius about 650 km in the Earth's center.This is the inner part of the IC for which robust seismic data have been recently collected.The directions of the fast and slow Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.Funded by Bibsam.
propagation of the seismic signal in the IMIC form an angle of 54 degrees.This is an exact match to the angle formed by the [111] and [100] directions in the bcc lattice, the directions of the fast and slow sound propagation in the bcc lattice [9,14].No such angle can be found in the hexagonal close-packed (hcp) structure.This match suggests the possibility that there is a single crystal of bcc Fe in the center of the Earth.Theoretically, the bcc phase was demonstrated to be dynamically stable in the IC [7,[11][12][13].Numerous ab initio studies [15][16][17][18][19] addressed the issue of the bcc stability and came to the conclusion that the bcc phase is thermodynamically unstable in the IC-the hcp phase is more, even though slightly, stable than the bcc phase.However, these AIMD simulations [15][16][17][18][19] were considering insufficiently large computational cells, and that was a critical flaw.The most recent AIMD studies [11][12][13] demonstrated that the bcc phase becomes thermodynamically stable under IC conditions if the size of the computational cell is large enough ( 1000 atoms).The mechanism of stabilization is described in detail in Ref. [11].AIMD calculations of elastic properties of the bcc iron under PT conditions of the IC were undertaken in the past, however, for way too small cells where the size was insufficient to establish and study the dynamically stable bcc phase [20].Thus at present there are no ab initio data for elastic properties of the Earth IC material.There are data for the so-called quasi-ab initio embedded atom model (EAM); however, it is highly desirable to compute the elastic properties with as few assumptions as possible using first-principles methods.One of the reasons is that atomistic models, for example, the EAM [21,22], do not take into account the impact of electronic entropy, which plays an essential role at the IC temperature [13,23].Another reason is that the EAM approach uses a prescribed functional dependence of interactions on interatomic distance while it is not clear whether such a preselected function is precise enough at extreme PT conditions in a wide range of compressions for all iron phases.

II. COMPUTATIONAL DETAILS
We calculated the energies and forces within the framework of the frozen-core all-electron Projector Augmented Wave (PAW) method [24], as implemented in the Vienna Ab initio Simulation Package (VASP) [25].In our AIMD runs we used the Fe potential with 3p 6 3d 6 4s 2 valence electrons.The energy cutoff was set to 400 eV.This is the cutoff applied in a number of studies.The cutoff convergence was carefully checked using cutoff up to 700 eV [26].Exchange and correlation potentials were treated within the Perdew-Burke-Erzernhof flavor of the generalized gradient approximation (GGA) [27].The finite temperatures for the electronic structure and force calculations were implemented within the Fermi-Dirac smearing approach [28].All necessary convergence tests were performed; the electronic steps converged within 0.01 meV/atom.The AIMD runs have been performed in the canonical ensemble for a given volume and temperature using the Nosé-Hoover thermostat.We used the cubic supercell obtained as the 10 × 10 × 10 translations of the standard two atom bcc cell, thus the supercell contained 2000 atoms.The integration over the Brillouin zone was done using the point.The time step was set to 1.0 fs.The temperature in all our calculations was set to 6600 K.This temperature was used before in other papers (see Table I) to compute elastic properties of the hcp Fe.Our choice is justified by the knowledge that the IC is close to melting temperature of the IC material, probably lowered by a certain amount of impurities.The temperature of 6600 K is below the melting temperatures of all iron phases; for example, the melting temperature of the bcc Fe at the pressure of 360 GPa is predicted to be equal to 7190 K [13].Notice that the latest assessment of the hcp Fe melting temperature gives about 6800 K [13] indicating its lower stability with respect to the bcc phase.Therefore, a temperature of 6600 K for the bcc Fe does not seem to be unreasonably high though there are arguments in favor of 6000 K and even lower temperatures in the IC [29].The lattice constant of the bcc cell was obtained by several short AIMD runs and the cell constant equal to 2.395 Å 2 produced pressure very close to the pressure of 360 GPa.The calculations were unique in what concerns the size of the simulation cell in the ab initio treatment of iron.First we equilibrated the initial ideal bcc structure for 5000 time steps, then applied strains as discussed below, and ran additional 5000 time steps for each strain.The runs were performed on 16 so-called fat nodes with a memory of 384 Gbytes at the Tetralith supercomputer cluster in Linköping, Sweden.The use of the nodes with small memory (96 Gbytes) leads to very inefficient parallelization.One week of computer time with 16 nodes resulted in about 1000 time steps for one strain.The single-crystal elastic constants were calculated via the generalized form of Hooke's law, where σ i stands for Cauchi stresses, c i j stands for elastic moduli, and j stands for strains, all in the Voigt notation (i, j = 1, 2, ...6) and assuming the Einstein summation convention.The calculated time-averaged stresses associated with strains applied to our 2000 atom supercell were used.To obtain the three nonzero c i j (namely, c 11 , c 12 , and c 44 ) of the cubic cell, just one strain is needed whose six components in the Voigt notation are given as = (δ, 0, 0, 0, 0, δ). ( The applied δ were ±1% and ±2%.For each δ (four in total), the MD run with a duration of 5000 time steps was used.The strains were applied to the unstrained cell thermally equilibrated for 5000 time steps; therefore calculation of thermal averages on the strained cells could be started right away.The convergence of the calculated averages over time was monitored via calculation of intermediate averages every 1000 MD time steps.The calculated dependencies of the time-averaged σ i on δ were subject to linear least-squares fit to directly estimate c i j from the relations (3) Isothermal elastic constants were transferred into adiabatic ones using the same procedure as in Ref. [20] assuming thermal expansion 10 −5 K −1 and heat capacity of 5k B .The sound velocity was computed in the Voigt approximation [30].The shear modulus G and the bulk modulus B of a cubic crystal can be written as and Then the longitudinal sound velocity C L and the bulk sound velocity C B can be calculated as where ρ is the density.The method was tested previously [31] in calculations of elastic properties of bcc Mo along the Hugoniot adiabat and was found to reproduce experimental data [32].

III. RESULTS
The calculated data are summarized in Table I and Table II compared with previously published data and the Preliminary Earth Model [5] (PREM).We immediately notice that the computed density of iron is higher than that from PREM by about 3%.This is less than the difference between the density of the hcp phase and PREM simply because the hcp phase is a close-packed one.While previous calculations performed for the densities matching that of PREM might look as closely matching the PREM velocities, this match is illusory because the calculated pressure in such cases is quite different from the pressure in the IC.Namely, for the hcp Fe with the density from the PREM shall produce pressures around 300 GPa; therefore the PT conditions are outside of the PT conditions of the IC [33].The second remarkable result is that the c 11 and c 12 are quite close to each other, unlike previously obtained results when calculations were done for much smaller computational cells, i.e. actually unstable but forced to stay intact by not allowing for shape relaxation.Note that since the Fe bcc is stable the c 11 -c 12 is positive.The error of c 11 -c 12 calculations is much smaller than the c 11 -c 12 magnitude that is also confirmed by analysis of their running averages.The low magnitude of c 11 -c 12 leads to the lowest shear modulus ever obtained for iron under the IC PT conditions (Table II).The IC material is known to possess remarkably low shear modulus [5,34] and this is also the feature obtained in this work for the bcc iron.The low shear modulus was also obtained [35] for Fe hcp at T = 7200 K.That result, however, should be treated with caution.First, the temperature of 7200 K is way too high for the IC.Second, and most importantly, T = 7200 K is about 600 K higher than the melting temperature of the hcp Fe computed with AIMD in the same approximation at the same pressure of 360 GPa.Forsblom and Grimvall [36] showed that when a single crystal is heated above the melting temperature the shear modulus rapidly decreases.Therefore the low shear modulus obtained by Martorell et al. [35] is an artifact and cannot be observed in reality.The Fe hcp phase, if heated to 7200 K, will melt.Therefore the comparison of our data with the data by Martorell et al. [35] is meaningless.While calculated in this study the longitudinal and shear velocities are slightly higher than those from PREM, there are several factors that might affect such a discourse even without using alloying with light elements as a possible cause.We discuss these factors below.We see also that the elastic constants computed in this work are reasonably close to those obtained using the embedded-atom method [8,21,22,37] (EAM).This suggests that the EAM model [21] is reasonably precise and provides further support to the stability of the bcc iron in the IC as communicated in the recent paper [13].It has to be emphasized that the AIMD provides the ultimate estimate of elastic properties while the validity of EAM and similar semiempirical models always has to be tested against an ab initio approach.
Figure 1 shows the P-wave and S-wave velocities calculated from the elastic constants of bcc-Fe at P = 360 GPa and T = 6600 K as a function of propagation direction.Notice the fast and slow sound propagation in the [111] and [100] directions, respectively, and the higher overall anisotropy of the P-wave, which is in contrast to the results for the simulations of small cells [38], where anisotropy of just 3% was obtained.Such an underestimation is likely due to the assumption that the directions of fast and slow propagation of seismic signal coincide with the axis of rotation and equatorial plane, correspondingly.In fact [6], while the fast direction of the propagation of the signal is indeed aligned with the axis of the Earth's rotation, the slow direction forms the angle of 54 degrees with the fast one.For such an angle the anistropy of the hcp single crystal is zero, while the anisotropy of the bcc single crystal is about 13%.We note that the seismic observations [6] are consistent with the existence of the bcc single crystal in the center of the IC with the [111] big diagonal oriented along the rotational axis and inconsistent with the hcp single crystal.There is the mechanism that allows  recrystallization of Fe polycrystal into a single crystal for the bcc phase via self-diffusion, unique for the bcc phase [11], while in the hcp phase such mechanism is absent.

IV. DISCUSSION AND CONCLUSIONS
The computed density of bcc Fe at the P = 360 GPa and T = 6600 K is 3% higher than the density according to PREM (Table I).The core density deficit, i.e., the difference between the density of the Earth's core and the density of the iron stable phase (bcc) in the IC, is considerably smaller than has been previously believed.In the not so distant past the amount of light elements in the IC was often estimated at the level of 10 weight percent and sometimes even higher.Recent publications estimate the amount of light elements based on ab initio calculations of the hcp Fe in the IC.Since the most recent studies suggest stability of the bcc phase, the amount of light elements remains largely unknown.However, we know now that it is small and it is plausible to assume that the small amount of light elements leads to small changes in elasticity.The longitudinal velocity (Table I) calculated in this work is a very good match to PREM.The density of the bcc iron is slightly higher than the density of the IC and this is only logical that the velocities are also slightly higher.The shear velocity is closer to PREM than for any Fe phase computed before.This is explained by the very small c 11 -c 12 in comparison with that previously computed.The reason for this difference most probably originates from the insufficiently large size of the bcc cells computed in the past leading to the bcc instability.Then, instead of the stable bcc structure, the cell oscillates between close-packed-like states and that leads to very high c 11 -c 12 .The relatively high c 44 leads, despite the small c 11 -c 12 , to a relatively high shear modulus as compared with PREM.We note, however, that this is not critical.The high-PT bcc Fe possess uniquely high self-diffusion [11] that results in dynamic disorder [39].This self-diffusion is similar when the bcc Fe is simulated using AIMD and the EAM (Fig. 2).The self-diffusion coefficient is known to converge at very large sizes of computational cells that are out of reach of AIMD (see Fig. 2).However, we know that the shear resistance of the bcc Fe is decreasing with time where the stress is relieved through self-diffusion.For times, characteristic of sound wave propagating through the compu-tational cell, the decrease of the shear modulus is sufficient to get very good agreement with PREM.Therefore even the already small difference between the G's computed in this work (Table II) and from PREM would become smaller for longer simulation times.We underline that even though the timescale of AIMD runs is small in comparison with those where the EAM model was used, the state is not transient.The EAM and AIMD dynamically stabilize the bcc sizes at exactly the same size of computational cell at the same temperature [13].Figure 2 shows that the MSD is quantitatively the same for EAM and AIMD.Therefore the bcc is not a transient state as some might think but a stable or, at the very least, metastable state.The detailed analysis of the shear relaxation [40] in bcc Fe demonstrates that the shear resistance can be split in a pure elastic response and a comparably small hydrodynamic response.If estimating the pure elastic response we subtract from the computed shear stress the hydrodynamic part (as computed using multimillion cells simulated for nanoseconds) then the remaining elastic response would fit the shear mod- FIG. 2. Size dependence of mean-square displacement (MSD) in bcc Fe.The package LAMMPS [45] was used for these calculations.The MSD is computed as in Ref. [13] using the EAM model [21] at the T = 6600 K and constant volume corresponding to the pressure 360 GPa for the number of particles 2000, 16 million, and 128 million atoms as indicated in the legend (1 time step = 1 fs).The AIMD computed MSD is equal to 0.28 Å 2 in close agreement to calculated using the EAM model.ulus of the IC very precisely.Another source of bringing the AIMD data and PREM in a better agreement comes from the spin fluctuations that might provide extra pressure term (up to 50 GPa [41]).That means that at the pressure of 360 GPa the density of the bcc Fe would be smaller and fit the seismic observation even better.Such extra pressure could erase the 3% density difference between the AIMD bcc Fe density and that from PREM.It is interesting that this effect is not relevant for the hcp phase.A proper calculation of spin fluctuations in high-PT iron might, however, require a better assessment of electron correlations [42], what is currently unfeasible in the simulations of large cells.Therefore the presented results here suggest that the standard attempts to provide rather narrow estimates for particular light elements operating with the equation of state of a likely irrelevant iron phase (hcp) [43] may not be the optimal way.We want to reiterate that the elastic constants provided in this study are, to our knowledge, correctly calculated on a DFT basis for the first time ever.

5 FIG. 1 .
FIG.1.Adiabatic single crystal P-wave and S-wave velocities for bcc Fe at IC conditions as a function of propagation direction: (a) V P , (b) V S1 , and (c) V S2 .This figure was generated by Unicef Careware[44].

TABLE I .
Calculated isothermal (adiabatic) elastic moduli of iron compared to previously published data and PREM.

TABLE II .
Isothermal (adiabatic) primary and secondary sound velocities and bulk and shear moduli of bcc iron compared with previously published data and PREM.