Figure 1: The static structure factor S(q), and radial distribution function g(r) of liquid Fe. The ‘calc. 136 GPa’ label was ab initio simulations data. The ‘exp. 116.1 GPa’ label was structural analyses of liquid Fe via in situ x-ray diffraction measurements (Kuwayama et al., 2020), and was shifted up.
Figure 2: The ion-ion dynamic structure factors S(q,ω) of Fe liquid under Earth’s core conditions. (a) is at 136 GPa and 4000 K, and (b) is at 330 GPa and 5500 K. The ‘qi’ labels are wave numbers with qi=2π√i/a (i=1,2…), where a is the equilibrium lattice parameters at the specified state.
Figure 3: Dispersion curve and adiabatic sound velocity of Fe liquid under the Earth’s core condition. The dispersion curves were directly derived from Figure 2.
Table 1: The fitting parameters of pressure and energy of pure Fe liquid under Earth’s core condition calculated by Eq.(6).
Figure 4: The ab initio calculated (a) isothermal bulk modulus KT and (b) adiabatic sound velocity VP at different temperature and pressure states
Figure 5: The geotherm profiles of liquid Fe under Earth core condition. (a) isentropic temperature Tad, (b) density ρ, (c) adiabatic bulk modulus KS and (d) adiabatic sound velocity VP. The ICB temperatures were selected as 5000 K, 5500 K, and 6000 K, respectively. ‘Experimental data’ (Kuwayama et al., 2020), ‘PREM values’ were data from PREM (Dziewonski and Anderson, 1981). ‘FeNiO by Sii’ was from (Li et al., 2023), and the label ‘Anderson XXX’ was from (Anderson and Ahrens, 1994). The four subfigures shared the lables.

The equations of states (EOS), longitudinal sound velocity (VP), and density are of great importance to in understanding Earth’s and terrestrial planets’ core composition behavior. the VP is a crucial physical quantity and directly correlated with the seismic wave. The dynamic structure factor is also crucial for the x-ray experiment to get an insight into the properties of matter and was recently reported to get the EOS of liquid Fe liquid experimentally [1]. Here, we calculate the dynamic structure factors and EOS by ab initio molecular dynamic simulations separately and compared the numerical accuracy of VP derived from both the calculated dynamic structures factors and EOS methods.

The calculation of VP is usually derived from the EOS or elastic constants. . e elasticity and VP of Earth’s interiors were t established by Birch [2], and the definations of the thermoelastic properties were supplied, Then the Birch-Murnaghan finite strain EOS of Fe were presented based on the ultrasonic, thermal expansion, and enthalpy data at 1 bar and on pulse-heating and shock wave compression and sound speed data up to 10 Mbar [3]. The EOS and thermoelastic properties were calculated by ab initio molecular dynamic simulations and Vinet (Morse-Rydberg) equation combined with the quasi-harmonic Debye model [4]. It showed that the density and adiabatic bulk modulus of pure liquid iron was found to be 8-10% and 3-10% larger than the preliminary reference Earth model (PREM) values [5], and the VP agreed well [4]. The VP and bulk modulus of liquid iron binaries were calculated by ab initio simulations combined with the Birch-Murnaghan EOS [6,7], which showed that Ni addition decreased VP while other light elements addition increased VP. The ab initio molecular dynamic simulations were also conducted on the iron-hydrogen alloys, and showed that both the density and VP of liquid iron containing ~1 wt.% hydrogen match seismological observations [8]. On the other hand, the elastic constants tensor were calculated by the first principle simulations to observe hydrogen effect on solid iron, and showed that hydrogen was an undesirable light element in the Earth’s core to match the seismological observations [9].

The experimental methods of the VP include derive from the phonon density of state [10] and structure factors [1]. The nuclear resonant inelastic x-ray scattering measured Fe to 153 GPa, and then the Debye average phonon velocity was obtained by the phonon density of state [10]. Thethermo elastic properties were determined by nuclear resonant inelastic xray scattering and calculated from ab initio theory, which agreed well with each other [10]. The VP was obtained to 45 GPa and 2700 K based on inelastic x-ray scattering measurement [1], then the Mie-Grüneisen EOS for liquid Fe was determined by results and previous shock-wave data. It indicated that Earth’s outer core exhibited 7.5%-7.6% density deflict 3.7%-4.4% velocity excess, and an almost identical adiabatic bulk modulus [1] The EOS and VP of Fe-S of fluid under Martian core conditions were measured using the ultrasonic pulse-echo overlap method combined with a Kawai-type multi-anvil apparatus up to 20 GPa [11]. The EOS was expressed by the adiabatic third-order Birch-Murnaghan EOS, and the VP was derived by its defination [11]. The traditional EOS was derived from empirical equations, such as Birch-Murnaghan EOS [2,3] and Mie-Grüneisen EOS [1].

The structure factors have also been calculated by orbital-free density-function theory simulation [12], molecular dynamics [13], and ab initio molecular simulations [14-16]. The ab initio calculation of ion-ion dynamic structure factors includes warm dense Al [15], warm dense lithium [12]. and mixture [14]. then the VP is the slope of the dispersion relation for small wave vectors, which is determined by analyzing the position of the side peaks of dynamic structure factors The VPs of Fe-Ni-O liquid were calculated by ion-ion dynamic structure factors and ab initio molecular dynamics simulations, and the results showed that VP of the Fe-Ni-O liquid was lower at the core-mantle boundary (CMB) and high at the inner-core boundary (ICB) than PREM values [17]. The calculated structure factors can be compared with available experimental xray and neutron scattering data. However, the ab initio calculated structure factors under Earth’s core condition have never been compared with experimental results and the calculation accuracy of VP is not Verified

The methods of the VP include EOS or elastic constants from independent DFT calculation and dynamic structure factors. In this paper, the ab initio molecular dynamics simulations were adopted in the calculation of the dynamic structures and EOS of liquid Fe. The static structure factors and radial distribution function were compared with the experimental data. The dynamic structure factors and dispersion curve were calculated, and then the VPs at the spicified states were collected. TheEOS of pure Fe was by multivariate polynomial method via separate ab initio molecular dynamics runs, and then the thermo elastic properties and isentropic profiels of Earth’s outer core were collected by definations . At last, thethermoelastic properties under Earth’s conditions were compared with PREM values, and the calculated VPs by the two methods were analyzed.

where q is the wave number, N is the total number of particles, and rv(t) is the vth ion position at time t.

The dynamic structure factor is calculated via the intermediate scattering function Fii (q,t). The ion-ion dynamic struc

Structure Factors

A detailed illustration of structure factors is shown in [18]. The dynamic structure factors are the spectral function of the density-density correlations in the system. The ion density in Fourier space

where q is the wave number, N is the total number of particles, and rv(t) is the vth ion position at time t.

The intermediate scattering function Fii(q,t) is defined as

The dynamic structure factor is calculated via the intermediate scattering function Fii (q,t). The ion-ion dynamic structure factor Sii (q, ω) is defined as the Fourier transform of the intermediate scattering function

where the ω denotes the frequency The ion-ion dynamic structure factor S (q,ω) characterizes the collective dynamics of flucuation of ionic density over both length and time scales. The S (q,ω) can be measured by inelastic x-ray scattering experiments (Kuwayama et al., 2020). By collecting the peak positions of S (q,ω) at different wave vectors q, the dispersion relation is obtained. The VP P is the slope of the dispersion relation at small wave vectors, On the other hand, the sound velocity can be calculated directly by the EOS (shown in Section 2.2). Therotically the sound velocities calculated by EOS and dynamic structure factors are equal.

The static structure factor is obtained via the intermediate scattering function Sii (q) = Fii (q,0). The isothermal bulk modulus (KT) can be computed from S(q) at q→0 with the form [18] of

where ni=N/V is the ionic density, kB is the Boltzmann constant and T is the temperature. The isothermal bulk modulus can also be determined from the EOS via separate ab initio molecular dynamics runs. To be pointed out, S(q) is the Fourier transformation of radial distribution function g(r). The radial distribution function (gr) is defined as

where V is the cell volume, N is the number of atoms, ri, and rj are atomic coordinates of atoms i and j, and <…> means the time or ensemble average.

Equations of States

Traditionally, the EOS of Fe or its alloy under Earth’s core condition is described by empirical equations (such as BirchMurnaghan EOS and Mie-Gruneisen EOS), and then its thermoelastic properties were collected. In this paper, we adopted a fitting method on EOS based on pressure-volume-temperature (P-V-T data) and energy-volume-temperature (E-V-T) data, then the thermo-elastic properties were calculated by definations directly

Based on the hundreds of P-V-T data and E-V-T data, the EOS is

where Aij and Bij are the fitted parameters of pressure and energy, respectively.

The definition of thermal expansion coefficient α is

The specific heat capacity at constant volume CV is

The isothermal bulk modulus KT is

The thermodynamic Gruneisen parameter γ is defined by

The adiabatic bulk modulus KS is

where α is the t of thermal expansion, and γ is the Grüneisen parameter. To be pointed out, KT and KS depend on both temperature and pressure. The adabatic sound velocity VP of liquid is

The isentropic temperature profile Tad is given by the following thermodynamic relationships as

where CP is isobaric heat capacity. An isentropic temperature Profile referred to as geotherm is obtained by integrating Eq. with the Gruneisen parameter γ for liquid Fe The pressure at ICB is 329 GPa, and the temperatures at the ICB are selected as 5000 K, 5500 K, and 6000 K. Then, the thermoelastic properties along the geotherm are collected by the isentropic pressue and temperature profiles

Calculation Details

The ab initio molecular dynamics calculations are implemented in the plane wave density functional VASP code [19,20]. In this calculation, projector augmented waves (PAWs) [21,22] and generalized gradient approximation (GGA) in the parameterization of Perdew, Burke, and Ernzerhof [23] are adopted The pseudopotential of Fe has p electrons as valence electrons. The plane wave cutoff is 500 eV, which is sufficient to ensure that the pressure converges within 1% accuracy. The time-dependent mean square displacement is employed to check the system in the liquid state. The selected time step is 1 fs. In the calculation of structure factors, we selected 128 atoms as the cell, and Fe atoms were randomly distributed in the cell. First, the NPT ensemble is calculated for 2000 time steps, and then the equilibrium volume is computed under certain pressures and temperatures. Then, the long-term correlation function is collected by the NVT ensemble with the equilibrium volume. When the total time exceeds 20 ps, the convergence ion-ion structure factors are obtained. In the calculation of EOS, we selected 64 atoms as the cell. The MD simulation with NVT ensemble was executed for a total of 6 ps with the last 3 ps considered as the production run, where the pressure and fluctutation were less then 1%

Structure Factors

The static structure factors S(q) and radial distribution function g(r) of Fe under the Earth’s core condition were collected in Figure 1. e shapes of S(q) and g(r) were similar to the reported experimental data, which was obtained from the in situ x-ray n measurements (Kuwayama et al., 2020). From the calculations of S(q) and g(r), it was obvious that was liquid. When the temperature was lower, the liquid Fe may be , and the shapes of S(q) and g(r) should be changed. Note that the smallest q available in S(q) is determined by 2π/a with a being the lattice constant of the simulation supercell, and the structure factor at small q represents the long-ranged structural information of systems. e KT can be calculated from S(q) at q→0 with the form of Eq. . In order to accurately compute KT, larger supercells were needed to obtain S(q) at small q. As the limited size of the supercells in the ab initio molecular dynamics simulations, the KT was not calculated by S(q) and just calculated by EOS.

The ion-ion dynamic structure factors were calculated by Eq.. e ion-ion dynamic structure factors S(q,ω) at t wave vectors qs at 136 GPa and 330 GPa were shown in Figure 2. e shape of the calculated ion-ion dynamic structure factors was similar to the reported high-pressure inelastic xray scattering measurements of liquid Fe [1]. e is a central Rayleigh peak and two ion acoustic modes are observed in S(q,ω). e ion acoustic modes are the Stoke and anti-Stoke peak. e sound velocities were derived by the dispersion relations of the ion acoustic modes. Figure 2 mainly provided the ion acoustic modes information of S(q,ω). As the qs increased, the peaks of S(q,ω) moved to the high ω and then move back to the low ω.

Figure 1

Figure 2

The adiabatic sound velocities along the geotherm curve can also be collected by long-time ab initio molecular dynamics simulations. By collecting the peak positions of S(q,ω) at ent wave vectors, the dispersion curves of Fe liquid at specid temperature and pressure states were obtained and plotted in Figure 3. The adiabatic sound velocity VP is the slope of the linear part of the dispersion curve at long wavelengths. The adiabatic sound velocity of Fe liquid at CMB was 8.38 km/s and at ICB was 10.46 km/s. As q increases, the correlation scale is reduced below that of the repulsion, the dispersion curves diverge. The accuracy of the VP is closely correlated with the energy , number of ions, time step and number of time steps.

Figure 3

Equation of State

From the ab initio molecular dynamics simulations, the pressure and energy at specified volume and temperature states (NVT ensemble) were collected. By using the sklearn package in python programming language, the EOS was fitted by Eq. within 0.5% accuracy. The fitting parameters of EOS by ing Eq. were shown in Table 1. From the d EOS, the isothermal bulk modulus KT and adiabatic sound velocity VP at d temperature and pressure state were calculated by Eqs. and , shown in Figure 4. e numerical results of KT and VP were consistent with the reported seismic [5,24] and calculated data [1-3]. e thermoelastic properties were usually derived by the empirical Birch-Murnaghan EOS [2,3] and Mie-- Gruneisen EOS [1]. Our calculated data about KT and VP verid the feasibility and s of our proposed multivariate polynomial method about the EOS.

Table 1

Figure 4

Properties of Earth’s Core

The geotherm curve was calculated by Eq. (13) and the fitted EOS. The isentropic temperature Tad and density ρ profile were calculated with three different ICB conditions (TICB=5000 K, 5500 K, 6000 K), shown in Figure 5 (a) and (b). For the same ICB temperature (5000 K), the Tad at CMB calculated by the fitted EOS was 6.46% lower than the Mie-Grüneisen EOS by experimental data (Kuwayama et al., 2020) and 5.88% lower than EOS (Anderson and Ahrens, 1994). The density ranged from 11.04 g/cm3 to 13.25 g/cm3 when TICB=5500 K. The density was about 10% higher than the PREM density value (Dziewonski and Anderson, 1981), which implied the existence of light elements. The density difference between our fitted data and the Mie-Grüneisen EOS by experimental data (Kuwayama et al., 2020) was within 3%. The adiabatic sound velocities calculated by EOS were both lower than the PREM values and calculated by the dispersion curve.

After the geotherm curve was confirmed, the isentropic bulk modulus KS and (d) adiabatic sound velocity VP were collected in Figure 5 (c) and (d). The KS ranged from 657.05 GPa to 1259.92 GPa when TICB = 5500 K, which was consistent with the PREM values (Dziewonski and Anderson, 1981). The temperature and pressure profiles in the calculation of ion-ion dynamic structure factors in Section 3.1 were also the geotherm curve with TICB = 5500 K. The VP along the geotherm ranged from 7.71km/s to 9.75 km/s by the EOS method and ranged from 8.38 km/s to 10.46 km/s by dynamic structure factors. The reported VP of Al calculated by dynamic structure factors agreed with experimental data within 6% (Rüter and Redmer, 2014). The sound velocity calculated by fitted EOS were coincided with the Mie-Grüneisen EOS by experimental data (Kuwayama et al., 2020).

Though the VP calculated by the EOS and dynamic structure factors are physically equivalent, the difference between the calculated VPs by these two methods was about 5%. When the VPs compared with the PREM values, the two methods may give different conclusions. The size of the supercell in the dynamic structure factors calculations was only 128 atoms, which would affect the accuracy of VP. The lattice parameter decided the number of q in the dispersion curves, which would directly affect the fitting accuracy of VP at small qs. The equivalent of the sound velocity calculated by dynamic structure factors and EOS are verified, but not applicable to small accurate difference analysis. When a sound velocity at a specified state is needed, ab initio molecular dynamic simulations about the dynamic structure factors with a large supercell is appropriate, which calculated a long-time simulation with only one calculation. When sound velocity in a wide temperature and pressure range are needed, the ab initio molecular dynamic simulations about the EOS with a relatively small supercell (64 atoms for example) were appropriate, which calculated relatively short-time simulations with hundreds of calculations.

Figure 5

The static structure factors, ion-ion dynamic structure factors, and EOS of Fe under the Earth’s core condition were calculated by the ab initio molecular dynamics simulations. Then, the sound velocities of Fe are calculated by both dynamic structure factors and EOS, which are physically equivalent. The VP of pure Fe under Earth’s core condition ranged from 7.71km/s to 9.75 km/s by the EOS method and ranged from 8.38 km/s to 10.46 km/s by dynamic structure factors method. Ab initio molecular dynamics simulations make a supplement for the higher energy x-ray to collect diffuse signals in a wider q range. The calculation of structure factors opened a way that can be compared to the high pressure in situ x-ray diffraction and inelastic x-ray scattering measurements.

This work was support by the National Natural Science Foundation of China (NSFC) [grant numbers 11975058, 11775031 and 11625415], and the fund of Key Laboratory of Computational Physics [grant number 6142A05RW202103]. We thanked for the fund support from Laboratory of Computational Physics in Institute of Applied Physics and Computational Mathematics. The data from this paper are available on reasonable request to Wei-Jie Li (liweijie8680@126.com).

  1. Kuwayama Y, Morard G, Nakajima Y, Hirose K, Baron AQR, et al. (2020) Equation of State of Liquid Iron under Extreme Conditions. Physical Review Letters, 124: 165701.
  2. Birch F, (1952) Elasticity and constitution of the Earth's interior. Journal of Geophysical Research, 57: 227-86.
  3. Anderson WW, Ahrens TJ, (1994) An equation of state for liquid iron and implications for the Earth's core. Journal of Geophysical Research: Solid Earth, 99: 4273-84.
  4. Ichikawa H, Tsuchiya T, Tange Y, (2014) The P V T equation of state and thermodynamic properties of liquid iron. Journal of Geophysical Research: Solid Earth 119: 240-52.
  5. Dziewonski AM, Anderson DL, (1981) Preliminary reference Earth model. Physics of the earth planetary interiors, 25: 297-56
  6. Badro J, Côté AS, Brodholt JP, (2014) A seismologically consistent compositional model of Earth’s core. Proceedings of the National Academy of Sciences, 111: 7542-5.
  7. Brodholt J, Badro J, (2017) Composition of the low seismic velocity E layer at the top of Earth's core. Geophysical Research Letters, 44: 8303-10.
  8. Umemoto K, Hirose K, (2015) Liquid iron hydrogen alloys at outer core conditions by first principles calculations. Geophysical Research Letters, 42: 7513-20.
  9. Caracas R, (2015) The infulence of hydrogen on the seismic properties of solid iron. Geophysical Research Letters, 42: 3780-5.
  10. Mao H, Xu J, Struzhkin V, Shu J, Hemley R, et al. (2001) Phonon density of states of iron up to 153 gigapascals. Science, 292: 914-6.
  11. Nishida K, Shibazaki Y, Terasaki H, Higo Y, Suzuki A, ect al. (2020) (2020) effect of sulfur on sound velocity of liquid iron under Martian core conditions. Nature Communications 11: 1954.
  12. White T, Richardson S, Crowley B, Pattison L, Harris J, et al. (2013) Orbital-free density-functional theory simulations of the dynamic structure factor of warm dense aluminum. Physical Review Letters, 111: 175002.
  13. Liu Q, Lu D, Chen M, (2020) Structure and dynamics of warm dense aluminum: a molecular dynamics study with density functional theory and deep potential. Journal of Physics: Condensed Matter, 32: 144002.
  14. Gill NM, Heinonen RA, Starrett CE, Saumon D, (2015) Ion-ion dynamic structure factor of warm dense mixtures. Physical Review E 91: 063109.
  15. Rüter HR, Redmer R, (2014) Ab initio simulations for the ion-ion structure factor of warm dense aluminum. Physical Review Letters 112: 145007.
  16. Witte B, Shihab M, Glenzer S, Redmer R, (2017) Ab initio simulations of the dynamic ion structure factor of warm dense lithium. Physical Review B 95: 144105.
  17. Li WJ, Li Z, Mo CJ, Ma Z, He XT, et al. (2023) self duffusion coffecent and sound velocity of Fe-Ni-O : Implications for the satification of Earth's outer core. Physics of the Earth and Planetary Interiors, 106983.
  18. Hansen JP, McDonald IR, (2006) Theory of simple liquids. Elsevier.
  19. Kresse G, Furthmüller J, (1996) Efficent iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B, 54: 11169.
  20. Kresse G, Hafner J, (1993) Ab initio molecular dynamics for liquid metals. Physical Review B, 47: 558.
  21. Blöchl PE, (1994) Projector augmented-wave method. Physical Review B, 50: 17953.
  22. Kresse G, Joubert D, (1999) From ultraso pseudopotentials to the projector augmented-wave method. Physical Review B, 59: 1758.
  23. Perdew JP, Burke K, Ernzerhof M, (1996) Generalized gradient approximation made simple. Physical Review Letters, 77: 3865.
  24. Garnero EJ, Helmberger DV, Grand SP, (1993) Constraining outermost core velocity with SmKS waves. Geophysical Research letters, 20: 2463-6.