Elasticity and Viscosity of hcp Iron at Earth's Inner Core Conditions From Machine Learning‐Based Large‐Scale Atomistic Simulations

Although considerable efforts have been made in the last years to examine the physical properties of hexagonal close‐packed (hcp) iron at extreme conditions, it remains challenging to explain many geophysical observations in Earth's inner core. Here we examine the elastic and plastic behavior of hcp iron and the effects of structural defects at inner core conditions using large‐scale atomistic simulations coupled with machine learning‐based interatomic potential. Our results suggest that the seismic anisotropy pattern in the inner core can be ascribed to the elastic anisotropy (6%) of hcp iron. The observed low shear wave velocity is largely produced by viscous grain boundaries in iron polycrystal. We also found highly mobile and abundant vacancies in hcp iron yield a viscous strength (1015±1) that is consistent with the geophysical observations. Therefore, our findings highlight the role played by structural defects and lessen the demand for light elements to explain the observed seismic data.

conditions. By using this approach, the full elastic tensor for hcp iron has been computed (Belonoshko et al., 2008;Li et al., 2018;Mattesini et al., 2010;Sha & Cohen, 2010;Vočadlo, 2007;Vočadlo et al., 2009). However, the obtained results are not always consistent with each other. For example, Belonoshko et al. (2008), Mattesini et al. (2010), and Sha and Cohen (2010) found that hcp iron is isotropic at Earth's inner core conditions, while Vočadlo et al. (2009)   conclude that the elastic anisotropy of hcp iron is around 6%-8%. Ritterbex and Tsuchiya (2020) have also applied ab initio methods coupled with the quasi-harmonic approximation to study vacancy formation and migration in hcp iron. The resulted viscous strength of hcp iron (10 17±1 Pa s) is low enough to induce seismic attenuation but fails to match the estimation from the Earth's nutation study, which is in the range of 2-7 × 10 14 Pa s (Koot & Dumberry, 2011).
We note that previous ab initio simulations were performed with a relatively small simulation cell and a short time scale due to the high computational cost of such methods. These limitations might lead to discrepancies among different theoretical studies and when comparing theoretical results with geophysical observations. In this study, we aim to investigate the physical properties of hcp iron and the effects of structural defects near melting at Earth's inner core conditions using molecular dynamics (MD) simulations. Specifically, we study the full elastic tensor and the effects of grain boundaries (GB). The formation and migration of atomic vacancies are also considered. To overcome the limits on system size and temporal length scale faced by ab initio methods, we have developed a deep-learning interatomic potential for Fe that retains ab initio accuracy but is more cost-effective. This method allows us to simulate systems that are large enough even to study the effects of GB on the shear wave velocity of hcp iron while retaining ab initio accuracy.

Molecular Dynamics Simulations With Deep Learning Potential
We only sketch the principal points for the techniques used in our study and leave all details to Text S1 and S2 in Supporting Information S1. To construct a deep-learning potential (DP) for iron that is applicable for MD simulations at Earth's core conditions, we first employed an iterative strategy to generate the training data set. Then the electronic structure of these configurations was determined using first-principles calculations with Quantum ESPRESSO (Giannozzi et al., 2009) to obtain forces, energies, and pressures, which provides an input for the training of the DP model using the DeePMD package Zhang et al., 2018). We quantify the error relative to DFT using an independent data set for a wide range of pressure-temperature conditions. For all conditions, the maximal root-mean-squared error (RMSE) for the energy, force, and pressure are 6.0 meV/ atom, 0.30 eV/Å, and 0.65 GPa, respectively. A similar performance is attained for configurations with different number of atoms at different pressure and temperature conditions, and for configurations with or without defects. Molecular dynamics simulations were performed using LAMMPS (Plimpton, 1995). In most cases, the Nośe-Hoover thermostat and barostat (Evans & Holian, 1985;Parrinello & Rahman, 1981) was used to control temperature and pressure, respectively.

Viscosity of hcp Iron
We focus on the contribution of dislocation creep to the plasticity of hcp iron (Wenk et al., 2000) as the role played by diffusion creep is very limited (Ritterbex & Tsuchiya, 2020). At inner core conditions with a high temperature and low deviatoric stress, dislocation climb is expected to be the rate-limiting step, and thus we can describe the strain rate (Reaman et al., 2011) as, where A is the dimensionless parameter and can be calculated using the method developed in Ritterbex and Tsuchiya (2020), G is the shear modulus, b is the Burgers vector, σ is stress, k B is Boltzmann constant, T is temperature, and D denotes vacancy diffusivity which can be directly estimated from our MD simulations. The viscosity will be stress dependent and can be derived by η = σ/2 .

The c/a Ratio of hcp Iron
We first calculate the c/a ratio for hcp iron, which is a necessary structural parameter for determining elastic properties. As shown in Figure 1, our results underestimate the c/a values by 9% compared to the theoretical predictions of Mattesini et al. (2010). The discrepancy might be due to the short simulation time (at most 16 ps) and length scales (128 atoms) in their ab initio MD simulations. To check this, we performed a simulation with 108 atoms at 6000 K and the resulted c/a ratio for a total simulation time of 16 ps is 1.625 ± 0.001, very close to the value of 1.627 obtained in Mattesini et al. (2010). Tateno et al. (2010) performed high pressure and temperature experiments up to 377 GPa and 5700 K. Their reported c/a values are slightly lower than our results, but the c/a uncertainties in their study are unclear. Belonoshko et al. (2008), Mattesini et al. (2010), and Sha and Cohen (2010) found that the P-wave anisotropy for hcp iron at 360 GPa and 6000 K is very weak and thus cannot match the observed seismic anisotropy of about 3% (Tkalčić et al., 2022). They attribute the vanishing of the P-wave anisotropy at high temperature to the formation of an isotropic and close-packed local environment when c/a is close to 1.63. To verify whether there is a strong connection between c/a ratio and elastic anisotropy, we calculated c 11 and c 33 at 360 GPa and 6000 K with different c/a ratios since they strongly affect the compressional wave anisotropy. We found the correlation between c 33 − c 11 and c/a up to 1.63 is very weak and the elastic properties of hcp iron remain anisotropic at the ideal value of c/a = 1.63. We conclude that elastic properties in hcp iron are influenced by the long-range structural properties of crystal and only weakly affected by the anisotropy in the first coordination shell.

The Melting Temperature
We computed the thermodynamic melting temperature of hcp iron at 360 GPa based on the solid-liquid coexistence method (Laio et al., 2000), which sets a basis for discussing the behavior of shear modulus near melting. For this, a 28 × 10 × 10 supercell (11,200 atoms) was first set up and equilibrated at 360 GPa and 6600 K. We fixed the first half of the simulation cell along the a direction and heated the other half up to 9000 K to melt it. The melted part was slowly quenched back to 6600 K. At this stage, the constraint was released to put the solid and liquid parts in contact. We then run MD simulations in the constant-pressure and constant-enthalpy ensemble to equilibrate the system. The melting temperature at 360 GPa was found to be 6620 ± 90 K, which agrees very well with previous estimations of 6600 K (Alfè et al., 2002) and 6410 K (Sun et al., 2018).

Elastic Constants for hcp Iron
We calculated the elastic constants of hcp iron at high temperature and 300 and 360 GPa, respectively, using a simulation box of 1,440 atoms ( Figure 2). We first discuss the results at 360 GPa, which is relevant for the inner core. From 4000 to 7,250 K, c 33 is always larger than c 11 , indicating that longitudinal wave velocity traveling along the c direction is faster than that along the basal plane. Our results compare well with  up to 7000 K with a difference of less than 5%. However, a very large disagreement is present at 7250 K for c 11 , where we report a value of 1,646 GPa compared to 1,424 GPa in . In order to understand the cause for this large deviation, we performed simulations with 64 atoms for 20 ps. The size of the simulation box and the short total simulation time were intentionally chosen to stay consistent with . We found a substantial reduction in c 33 compared to the simulation with 1,440 atoms. The difference is comparable to the discrepancy between our values for c 11 calculated with 1,440 atoms and the one reported by , and thus we believe the latter is due to the small simulation size used in .
We also calculated the elastic constants at 300 GPa. We observe a similar trend for the temperature dependence of elastic constants as that at 360 GPa, and c 33 is always larger than c 11 in the temperature range of 4000-6,000 K. Vočadlo et al. (2009) found that there is a crossover in the magnitude of c 33 and c 11 at 5500 K and 300 GPa, which is not present in our study (see Figure SM3 in Supporting Information S1). The discrepancy might also be due to system size since our simulations with 64 atoms yield a c 33 − c 11 that is comparable to that in Vočadlo et al. (2009)  exceeds T m by 700 K. They argue that the drop in elastic constants and shear modulus could occur also below T m due to a rapid increase of atomic defects defined as atoms that have a coordination number in the first coordination shell differing from 12, but do not provide evidence for this in their simulations.
We employ the solid-liquid coexistence method, where the liquid phase plays the role of a source or sink of atomic defects, to examine the fraction of defective atoms in the thermodynamic melting process (see Text SM3 in Supporting Information S1 for more details). Our results show no significant change in the fraction of defective atoms between the solid hcp iron in contact with a liquid phase and the perfect hcp iron at the same temperature ( Figure SM4 in Supporting Information S1). Therefore, there is no structural evidence to support the hypothesis of a strong shear modulus softening near T m . We conclude that the reduction observed in  is a consequence of the fact that the concentration of atomic defects at T k is higher than at T m , which is not surprising considering that T k is the temperature where homogeneous nucleation occurs.

Formation and Migration of Vacancies
The formation and migration of vacancies strongly affect dislocation mobility at high temperature and thus the viscous strength of hcp iron. We use the thermodynamic integration method to determine the Gibbs free energy associated with the formation of a vacancy. As shown in the inset of Figure 3b, the obtained free energy decreases from 9.19 eV at 4000 K to 5.82 eV at 7000 K. The corresponding equilibrium vacancy concentration varies from 2.7 × 10 −12 to 6.5 × 10 −5 . At 360 GPa and 6620 K, the vacancy concentration is about 1.4 × 10 −5 , which compares well with other metals like copper near the thermodynamic melting temperature (Siegel, 1978). However, our reported formation energy values at 360 GPa and 5000 K are considerably lower than that in Ritterbex and Tsuchiya (2020), which is 8.89 eV compared to our value of 8.05 eV. The discrepancy is likely caused by anharmonicity effects which were not included in their study. A recent theoretical study (Cheng & Ceriotti, 2018) also shows that the quasi-harmonic approximation overestimates the vacancy formation energy in bcc iron at high temperature.
We then used the method proposed by Mendelev and Mishin (2009) to determine diffusion coefficient of vacancies. Our results are plotted in Figure 3a. At 4000 K, unlike the case of perfect hcp iron, there is a long-time drift in the mean-squared displacement of atoms caused by the presence of a vacancy. The calculated diffusivity dramatically increases from 9.8 × 10 −22 m 2 /s at 4000 K to 1.6 × 10 −12 m 2 /s at 7,000 K, nearly nine orders of magnitude. Near melting, the diffusion coefficient is around 3.1 × 10 −13 m 2 /s. Our calculated diffusivity is almost three order of magnitude higher than the previous study (Ritterbex & Tsuchiya, 2020). The discrepancy is due to the combined effects of anharmonicity and high-temperature atomic motions missing in their study. Our results agree with a previous experimental work (Reaman et al., 2012), which predicts an inner core solid-state diffusivity of 1 × 10 −14 m 2 /s at T/T m = 0.95 although the result was based on a very large extrapolation from 65 GPa to 2000 K to inner core conditions.

Seismic Anisotropy
With the use of a deep learning interatomic potential, a more reliable estimation of the elastic anisotropy of hcp iron is achieved by effectively enhancing time and length scales in MD simulations. Here we only focus on P-wave anisotropy as the strength of shear wave anisotropy in the inner core still suffers from considerable uncertainties (Tkalčić et al., 2022). Figure 4a shows the maximal seismic anisotropy as a function of temperature at 360 GPa. We found that the maximal anisotropy for the compressional wave is about 6%, and it is rather independent of temperature, especially close to Earth's inner core temperature.
We scale the calculated single-crystal anisotropy of hcp iron to the global seismic observations based on the multiscale model developed by Lincot et al. (2016). In their model, polycrystal plasticity, inner-core formation models, elastic moduli, and crystal alignment mechanism are integrated, and possible combinations that   and the preliminary reference Earth model (PREM) (Dziewonski & Anderson, 1981). The red start represents the effects of viscous grain boundaries (GB) on shear wave velocity. The magenta line denotes the pre-melting effects as observed in . The green dotted lines in a and b are guides to the eye obtained by a linear fitting and fitting second-order polynomials, respectively. (c) Shows the polycrystal iron sample used for the shear wave velocity calculations. The length of the cubic side is about 200 Å and the simulation box contains 1,164,419 atoms. We color each atom based on their mobility. The high mobile atoms always lie at grain boundaries. Even with these extended defects, the structure of the hcp iron grain is still hexagonal closed-packed by comparing its pair distribution function with that of the perfect hcp iron as shown in (d).
match the present-day seismic observations have been determined. They also identified two dimensionless elastic parameters (Δ c−a and Δ45°) that most affect the inner core anisotropy, where Δ c−a defines the difference of P waves velocities along c and a directions, and Δ 45° describes the concavity of the single-crystal P-wave velocities (see Text SM2 in Supporting Information S1 for more details). Our calculated Δ c−a and Δ45° are about 5% and −4%, respectively, and are relatively insensitive to temperature. Based on these values, we found that only one model in Lincot et al. (2016) can produce enough anisotropy to be consistent with the seismic observations, where the dominant deformation mechanism is pyramidal slip, and the inner core preferentially grows in the equatorial plane (Yoshida et al., 1996) with solidification texturing (Bergman, 1997) and density stratifications (Deguen & Cardin, 2009). It should be noted that our study does not consider the effects of light elements such as Si, S, and C.

Shear Wave Velocity
We employ the Voigt average scheme to estimate the shear modulus and shear wave velocity for a polycrystalline aggregate and show the results in Figure 4b. Our reported values are much higher than geophysical observations (Tkalčić & Phạm, 2018). Since we ruled out the possibility of a strong pre-melting shear modulus softening near the thermodynamic melting temperature as proposed in , the low shear velocity is likely caused by the presence of light elements (He et al., 2022) or GB (Belonoshko et al., 2007).
The incorporation of light elements is of particular interest as their presence in Earth's core is supported by geochemical evidence (Hirose et al., 2013). A previous study suggests that for light elements staying at substitutional sites, there is only a mild reduction in the shear wave velocity (Vočadlo, 2007). Therefore, a large concentration of light elements in substitutional sites are needed in order to match the seismological data, something that might violate geochemical constraints. In contrast, a considerable decrease in the shear modulus is observed for light elements like O, C, and H sitting at interstitial positions (He et al., 2022). However, it remains uncertain which site is more thermodynamically stable. He et al. (2022) show that O, C, and H prefer to stay at interstitial sites. But their calculations were estimated using the quasi-harmonic approximation. As shown in our study, this theory might break down at high temperature. Therefore, more research is needed to clarify whether light elements can be considered as the primary reason for the decrease of the shear wave velocity of hcp iron to a level that matches geophysical observations.
Another plausible origin of the observed low shear velocity is the presence of viscous GB, a mechanism that was first proposed by Belonoshko et al. (2007). These extended defects can relax the shear stress and dramatically decrease the shear bulk modulus. To check this hypothesis, we have performed MD simulations with the DP model for a hcp iron polycrystal, which contains 10 grains and more than 1 million atoms, as shown in Figure 4. The polycrystal was first equilibrated at 360 GPa and 6000 K. Then we applied a shear strain ϵ 44 with a magnitude of 2%. Since the polycrystal is isotopic on average, the resulted off-diagonal stress was used to calculate the shear modulus c 44 . We found that the polycrystalline nature of the sample leads to a 17% reduction in shear velocity (Figure 4), which is in a better agreement with the seismic data. The remaining discrepancy might be due to the effects of light elements that are not considered. We stress that the grain size used in our MD simulations is probably unrealistic if compared to the grain size in Earth's inner core, and thus whether our results are applicable to the inner core is subject to further study.

Viscosity
We use the obtained self-diffusivity and shear bulk modulus to evaluate the viscosity of hcp iron at 360 GPa. We found that dislocation creep yields a viscosity of 3.5 × 10 14 -3.8 × 10 16 Pa s at 6000 K for possible strain rates between 10 −14 -10 −18 s −1 , which is typical for the convection process in the inner core (Sumita & Bergman, 2007). The viscosity of hcp iron could be further reduced by grain boundary sliding if present.
Our calculated viscosity is about three orders of magnitude lower than the value reported by Ritterbex and Tsuchiya (2020), which is mainly caused by the increased equilibrium vacancy concentration and lower migration enthalpy. The lower viscosity of 10 15±1 Pa s is more consistent with the estimation from the Earth's nutation study (Koot & Dumberry, 2011), which predicts an inner core viscosity of 2-7 × 10 14 Pa s. It also matches the geodynamic observations according to which the inner core viscosity must be low enough (<10 16 Pa s) to allow its shape to quickly adjust when the inner core moves out of the alignment with the mantle (Buffett, 1997). It should 8 of 10 be noted that our study does not consider the effects of light elements. Recent experiments up to 60 GPa and 1650 K suggest that incorporating Si would harden hcp iron and thus increase its viscosity (Brennan et al., 2021). However, whether this effect persists at the inner core conditions where diffusion is pervasive as shown in our study is unclear and is subject to further study.

Conclusions
In this study, we employ MD simulations to examine the physical properties of hcp iron and the effects of structural defects at Earth's inner core conditions. To overcome the limits on system size and temporal length scale faced by ab initio methods, we have developed a deep-learning interatomic potential for Fe that retains ab initio accuracy but is more cost-effective. Our results suggest that at Earth's inner core conditions, the anisotropy of the compressional wave velocity is about 6%. Therefore, a preferential alignment of the c-axis of iron crystals along the Earth's rotation axis is consistent with the observed seismic anisotropy pattern. In addition, pre-melting effects on the elastic properties of hcp iron at core conditions predicted in earlier calculations are not observed. We also found the viscous GB can dramatically reduce the shear wave velocity of hcp iron, providing a possible mechanism to explain the observed low shear wave velocity in the inner core, with the caveat that grain sizes in the core are likely to be orders of magnitude larger than those simulated here. Both the concentration and mobility of vacancies in hcp iron are very high, which considerably weakens hcp iron and causes a low viscosity of 10 15±1 Pa s that is in a reasonable agreement with geophysical measurements. Our results highlight the role played by structural defects and lessen the demand for light elements to explain the observed seismic data.

Data Availability Statement
The source code for first-principles calculations is available at https://www.quantum-espresso.org/. The DeePMD package has been used to train the deep-learning interatomic potential, which can be found at https://deepmodeling.com/. According to AGU data policy, the raw data supporting the findings of this study are available at https://doi.org/10.5281/zenodo.7331669.