Superionic hcp-Fe alloys and their seismic velocities in Earth’s inner core

Earth’s inner core (IC) is less dense than pure iron, indicating the existence of light elements within it 1 . Si, S, C, O, and H have been suggested to be the candidates 2,3 , and the properties of Fe-light-element alloys were studied to constrain the IC composition 4-19 . Light elements have a signicant inuence on seismic velocities 4-13 , melting temperatures 15-17 , and thermal conductivities of Fe-alloys 18,19 . However, the state of the light elements in the IC is rarely considered. Using ab initio molecular dynamics (AIMD) simulations, we found that H, O, and C in hexagonal close-packed (hcp) Fe transform to a superionic state under IC conditions, showing high diffusion coecients like liquid. It suggests the IC can be in superionic state rather than normal solid state. The liquid-like light elements lead to a signicant reduction in the seismic velocities approaching the seismological observation of the IC 20,21 . The signicant decrease in shear wave velocity (V S ) gives an explanation on the soft IC 21 . In adddtion, the light-element convection in the IC has potential inuence on the IC seismological structure and magnetic eld.


Full Text
Earth's inner core (IC), the central part of our planet, was formed and grown because of the solidi cation of liquid Fe alloys in the outer core. Seismological observations suggest that the structure of the IC is very complicated and di cult to understand. An unresolved problem is that the IC is soft with a low shearwave velocity (V S = ~3.6 km/s) 20,21 , which cannot be matched by sound velocities in Fe and Fe alloys [8][9][10][11][12][13] .
Theoretical prediction suggests that the premelting effect of Fe and Fe alloys can lead to strong nonlinear shear weakening at high temperatures [22][23][24] . However, a temperature close to the melting temperature with a small gradient is required. Alternatively, some fractions of melts need to be considered 25 , but the mechanism of the stable presence of the melting fraction since the formation of the IC is unclear. Another long-standing controversy is the reason for the anisotropic seismic velocity structure. Generally, the compression wave velocities in the polar direction are ~3 % larger than those in the equatorial direction in the IC 26 . This phenomenon is explained by the preferred orientation of hexagonal close-packed (hcp) Fe lattices with <001> directions 27 or body-centre-cubic (bcc) Fe lattices with <111> directions 28 along the spin axis. Moreover, earthquake doublets show that seismic waves present different travel times in the IC over the past few decades. This suggests that the IC structure changes with time, which is attributed to the super-rotation or oscillation of the IC [29][30][31] . Investigations on the seismic properties of Fe and Fe alloys have improved the understanding of the IC structure. However, direct experimental measurements under the IC conditions are extremely challenging, while AIMD simulations have increased the knowledge on the properties of Fe alloys 13,[22][23][24]27,28 .
H, O, and C are, respectively, the rst, second, and fourth abundant elements in our solar system 32 . They are also critical volatiles presenting in Earth's interior since its formation. In addition, they can be brought down to the core through mantle convection. In our simulations, H, O, and C atoms were placed at the interstitial sites of hcp-Fe to verify light-element effects on the properties of Fe alloys. AIMD simulations were carried out at high P-T on FeH 0.25 , FeO 0.0625 , and FeC 0.0625 structures, which contained 0.45, 1.75, and 1.33 wt. % of light elements, respectively. The structures were constructed based on previous experimental studies, and the details are provided in the supplementary information S1. Solid-superionicliquid phase transitions with increasing temperature are obtained and shown in Fig. 1. At low temperatures, light elements in hcp-Fe undergo small mean-square displacements (MSDs) and simply vibrate about their lattice positions demonstrating that the material is an ordinary solid. At temperatures above 3000 K, the light elements become highly diffusive like liquid as the MSDs increase monotonically during the simulation last over 100 ps indicating a superionic state ( Supplementary Fig. 2). A further increase in temperature leads to the melting of the lattice. The melting temperatures of these structures were estimated using the two-phase method to exclude the superheating state in the calculations (Supplementary Information S2) 33  The electric conductivities of Fe alloys are very important in understanding the energy budget of geodynamo and the heat ux of Earth's core 18,19 . For superionic Fe alloys, the electric conductivity is composed of ionic conductivity and electronic conductivity. The ionic conductivities were calculated by Nernst-Einstein equation, and the values are in the range 10 2 -10 3 S m -1 (Supplementary Fig. 7). We also estimated the electronic conductivities using the combination of density functional theory with dynamical mean eld theory (DFT+DMFT) method ( Supplementary Information S4). Although the presence of light elements decreases the electronic conductivity of pure Fe, it is still over 2-3 orders of magnitude higher than the ionic conductivity, and therefore the conductivity contribution due to the ionic diffusion is negligible.
The liquid-like light elements in the Fe lattice also affect the seismic velocities. AIMD simulations with the NVT ensemble were conducted on FeH 0.25 , FeO 0.0625 , and FeC 0.0625 at 2000-6000 K and 360 GPa. In order to exclude the effect of the superheating state, the simulation temperatures are lower than the predicted melting temperatures. The structures and cell parameters were carefully determined to maintain hydrostatic stress. The independent elastic constants for C 11 , C 12 , C 13 , C 33 , and C 44 were calculated by distorting the equilibrium structure and solving the stress-strain relations (Supplementary information S5). The compression (V P ) and shear (V S ) wave velocities of these structures at high P-T were deduced from the elastic constants and are plotted in Fig. 3. The V P of the alloys decrease almost linearly with temperature, while V S show an accelerated decline upon the solid-superionic transition (Fig. 3c) Table 1 and Supplementary Fig. 9). Therefore, the superionic Fe alloys at 5500-6000 K show quite similar seismic characteristics as that of the IC. Previous theoretical simulations show a decrease in V P and V S of Fe and Fe alloys with increasing temperature because of the premelting effect [22][23][24] . However, the calculated velocities can be consistent with the PREM at approximately 7200 K for Fe and 6700 K for Fe-Si-C alloy 22,24 , which are higher than the most computationally and experimentally predicted temperatures of IC [15][16][17][18][19] . Here we provide a new mechanism by which highly diffusive light elements bring strong anharmonic vibration to Fe alloys and lead to signi cant reductions of V P and V S close to the PREM model at 5500-6000 K. Elastic softening caused by its superionic effect was also reported in Li 2 O when it is used as a nuclear fusion material at high temperatures 40 .
Seismologic observations suggest that Earth's core is composed of the liquid outer core with Vs equal to zero, and the solid inner core with Vs equal to ~3.6 km/s. At the ICB, solidi cation of the IC generates latent heat and leads to the separation and buoyancy of light-element phases, which promote the convection ows of the outer core, producing and powering the geomagnetic eld (black dashed arrows in Fig. 4). Based on our simulations, some light elements such as H, O, and C can present in the IC in the superionic state. Owing to the high diffusion coe cients, light elements continue their circulation in the IC (red dashed arrows in Fig. 4). In this case, the IC is not conventionally known normal solid Fe alloy, but a mixture of solid Fe and light-element uids. The convection of the ionic light elements may generate magnetic eld in the IC. However, quantitative modelling is needed to study if this magnetic eld is su cient to in uence Earth's magnetic eld.
The liquid-like light elements also bring signi cant in uence on the elastic properties of the IC, and the long-standing seismic puzzles in the IC, as mentioned above, can be rationalized within the superionic IC model. First, the liquid-like light elements have a profound effect on the seismic velocities, presenting the V S quite close to that of the PREM at 5500-6000 K, which explains why the IC is quite soft. In addition, superionic Fe alloys show 8-11 % lower Vp than that of pure Fe. If the distribution of superionic light elements in the IC is inhomogeneous with a higher concentration in the equatorial direction than in the polar direction, seismic velocities along the equatorial direction will be lower, which gives a simple reason for the anisotropic seismological structure. Also, the enrichment of liquid-like light elements may enhance seismic wave attenuation, which has been widely observed in the IC 21 . Moreover, the light-elements convection changes the distribution of light element in the IC leading to the evolution of the seismological structure with time. This phenomenon may provide an alternative interpretation of the differential travel time of similar earthquakes [29][30][31] . Generally, the seismic properties of superionic Fe alloys are consistent with the seismology characteristics of the IC.
Suggested by some geodynamo models, geomagenetic eld also has in uence on the IC 41,42 . In this case, the diffusion of ionic light element can be affected by the inner-core magnetic eld under Lorentz force, and the convection and distribution of light elements may related with geometry and strength of the magnetic eld. Besides, both seismic velocity in the IC and the geomagnetic eld show observable change in last few dacades [29][30][31] . Then, it is important to study if there is a certain relation between the inner core seismic structure and geomagnetic eld, which might be a clue for understanding the structure and evolution of the IC.

Method
The calculations were performed using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional and projector augmented wave (PAW) method 43

as implemented in the Vienna Ab Inito
Simulaton Package (VASP) 44 . A plane wave representation for the wave function with a cutoff energy of 800 eV was adopted. In AIMD simulations, the cutoff energy was set to 400 eV with Γ point for k-space sampling. AIMD simulations within canonical ensemble (NVT) were conducted to study the solidsuperionic-liquid phase transition and diffusion coe cients of superionic light elements. The simulations lasted ~10-100 ps at temperatures from ~2000 to ~7000 K and pressure from ~250 to ~400 GPa with a time step of 1fs. The two-phase method 33 45 . For each equilibrium structure, a 10,000 time-steps (10 ps) simulation was conducted to check for sure that the stress eld is hydrostatic. The elastic constant C ij was calculated by distorting the equilibrium structure and solving the stress-strain relations. Details for the calculation methods and convergence tests with large supercells and long simulation time are provided in Supplementary Information.

Declarations Data availability
The data that support this study are available from the corresponding author upon reasonable request.