Ab Initio Thermoelasticity of Liquid Iron-Nickel-Light Element Alloys

Abstract

The earth’s core is thought to be composed of Fe-Ni alloy including substantially large amounts of light elements. Although oxygen, silicon, carbon, nitrogen, sulfur, and hydrogen have been proposed as candidates for the light elements, little is known about the amount and the species so far, primarily because of the difficulties in measurements of liquid properties under the outer core pressure and temperature condition. Here, we carry out massive ab initio computations of liquid Fe-Ni-light element alloys with various compositions under the whole outer core P, T condition in order to quantitatively evaluate their thermoelasticity. Calculated results indicate that Si and S have larger effects on the density of liquid iron than O and H, but the seismological reference values of the outer core can be reproduced simultaneously by any light elements except for C. In order to place further constraints on the outer core chemistry, other information, in particular melting phase relations of iron light elements alloys at the inner core-outer core boundary, are necessary. The optimized best-fit compositions demonstrate that the major element composition of the bulk earth is expected to be CI chondritic for the Si-rich core with the pyrolytic mantle or for the Si-poor core and the (Mg,Fe)SiO₃-dominant mantle. But the H-rich core likely causes a distinct Fe depletion for the bulk Earth composition.

1. Introduction

The earth’s core is thought to be composed of Fe-Ni alloy including substantially large amounts of light elements. These light elements account for observed density deficits of ~10% for the liquid outer core and ~5% for the solid inner core [1,2,3,4,5,6,7]. Determination of the light element (LE) composition of the outer core (OC) has long been one of the central research topics in the deep earth sciences. The density (ρ) and adiabatic bulk (K_T) and shear (K_S) moduli of iron and iron-LE alloys are key to interpreting seismological observations and then constructing a compositional model of the core [5,8]. However, those of the liquid states at the OC pressure (P) and temperature (T) (from ~136 to ~329 GPa and from ~4000 to ~6000 K) are still limitedly clarified experimentally. So far, static experiments have been performed up to less than 100 GPa [9,10,11]. Higher-P behavior of liquid iron was on the other hand investigated by shock wave experiments in multi-Mbar condition [2,12,13,14,15,16]. The temperature, however, changes along the principal Hugoniot and dramatically increases with increasing pressure to more than 8000 K at the pressure of the inner core (IC)-OC boundary (P_ICB) of ~329 GPa, which is far higher than the expected actual ICB temperature (T_ICB) of ~5000–6000 K [17,18,19,20,21,22]. Experimental determination of thermoelasticity of liquid iron alloys in the whole P, T condition of the earth’s OC thus remains technically impractical.

In contrast, ab initio molecular dynamics (AIMD) simulations have been widely applied to clarify ρ and P-wave velocity (V_P) of liquid iron and iron-LE alloys at the OC conditions in order to constrain the OC composition by interpreting seismological observations. These parameters for the Fe-O, Fe-Si, Fe-S, Fe-C, Fe-Ni, and Fe-Si-O system were calculated [23,24,25]. However, the data points in these studies were limited; two particular compositions of Fe_0.82Si_0.10O_0.08 and Fe_0.79Si_0.08O_0.13 only were considered [24] and two particular pressures of the core-mantle boundary (CMB) and ICB only were considered [23]. In particular, in the latter, empirical pressure corrections of 10 GPa and 8 GPa were adopted at the CMB and ICB respectively, though the optimized OC compositions are essentially sensitive to these corrections. Meanwhile, some studies have been performed throughout the whole OC P, T conditions for pure Fe [7,26], Fe-S [27], and Fe-H [28]. However, different formulations were employed to model their thermal equations of state, making a quantitative comparison of the reported thermoelasticity not easy.

In this study, ab initio MD simulations are performed for binary and ternary Fe-Ni-LE alloys with several different LE and Ni fractions from ~100 to ~450 GPa and from 4000 to 8000 K. Equations of state (EoS) and thermoelasticity are then analyzed for each alloy through the same internally consistent way [7]. Using modeled thermoelasticity, we optimize light element compositions for each alloy as a function of the T_ICB and discuss the possible OC composition.

2. Results and Discussion

2.1. Effects of LE on the Thermoelasticity of Liquid Iron

Calculations with several different LE concentrations clarify systematic trends on the effects of LEs on the thermoelasticity of liquid Fe. Incorporations of LEs always decrease ρ but increase V_P (

Table 1. Effects of light element (LE) incorporation on $V_{\mathrm{P}}$ and $\rho$ of liquid Fe calculated at the P_CMB and 4000 K and at the P_ICB and 5300 K. X_LE represents the fraction of LE in atom%.

P, T Condition		O	Si	S	C	H
P = PCMB T = 4000 K	$\frac{\partial \ln V_{\mathrm{P}}}{\partial X_{\mathrm{LE}}}$	0.05(1)	0.13(1)	0.06(1)	0.16(1)	0.02(1)
	$\frac{\partial \ln \rho }{\partial X_{\mathrm{LE}}}$	−0.34(1)	−0.51(1)	−0.41(1)	−0.30(1)	−0.24(1)
	$\frac{\partial \ln V_{\mathrm{P}}}{\partial \ln \rho }$	−0.14(1)	−0.26(1)	−0.16(1)	−0.54(1)	−0.10(1)
P = PICB T = 5300 K	$\frac{\partial \ln V_{\mathrm{P}}}{\partial X_{LE}}$	0.09(1)	0.21(1)	0.16(1)	0.20(1)	0.07(1)
	$\frac{\partial \ln \rho }{\partial X_{\mathrm{LE}}}$	−0.31(1)	−0.48(1)	−0.38(1)	−0.31(1)	−0.21(1)
	$\frac{\partial \ln V_{\mathrm{P}}}{\partial \ln \rho }$	−0.29(1)	−0.44(1)	−0.42(1)	−0.63(1)	−0.35(1)

), but trends are different depending on the type of LE. It is found that incorporations of larger Si and S atoms have only marginal effects on the volume (volume per atom), then the EoS is nearly unchanged (Figure S1). In contrast, incorporations of smaller O, C, and in particular H atoms reduce the volume considerably in the whole OC P range. These are related to the fact that the Fe-Si and S alloys are so-called substitutional-type, while the Fe-O, C, and H alloys are interstitial-type as recognized generally in lower P condition.

Because of these volume reductions, ρ variations associated with the O and H incorporations are smaller than those expected from the small masses. As a result, the effects of Si and S incorporations on ρ are larger than those of O and H incorporations (Table 1). These behaviors are consistent with a recent study reporting structural and dynamical properties of Fe-LE alloy liquids [29]. A similar tendency is seen in V_P, but the systematics is less pronounced since the effects of LEs on ρ and K_S are partially cancelled. Perturbation ratios ($\partial \ln V_{\mathrm{P}}/\partial \ln \rho$) are sometimes referred to discuss the chemical heterogeneity in Earth’s deep interior [30,31]. In the present cases, absolute values of this ratio are always smaller than 1, indicating that the effects of LE incorporations are always much larger in ρ than in V_P.

2.2. Optimized Compositions

Misfits in ρ and V_P between the Fe-Ni-X liquid alloys and the preliminary reference earth model (PREM) [32] are then evaluated as ${\displaystyle \sum }[{\left(\frac{\rho -{\rho }_{PREM}}{{\rho }_{PREM}}\right)}^2+{\left(\frac{V_{\mathrm{P}}-V_{{\mathrm{P}}_{PREM}}}{V_{{\mathrm{P}}_{PREM}}}\right)}^2]$ (Figure 1), where the summation is taken over the whole OC pressure range. It is clearly demonstrated that the misfits are sensitive to the LE concentration and temperature but not so sensitive to the Ni concentration. The best-fit compositions along two adiabats (T_ICB = 5000 K and 6500 K) with three different Ni/(Fe + Ni) ratios, which can be defined by the minima of misfits, are listed in

Table 2. Best-fit compositions of binary and ternary alloys at T_ICB = 5000 K and 6500 K. Misfit, Mg/Si, and Mg/Fe represent a misfit in V_P and ρ from the PREM, Mg/Si, and Mg/Fe ratios expected to the bulk earth with the pyrolytic mantle, respectively. Errors from the fitting procedures are represented in parentheses.

TICB	Best-Fit Composition	$\mathbf{Misfit}(\times {10}^{-2})$	Mg/Si	Mg/Fe
$5000$ K	Fe0.78O0.22	1.8(1)	1.25(1)	1.03(1)
	Fe0.85Si0.15	2.7(1)	1.06(1)	1.04(1)
	Fe0.81S0.19	1.6(1)	1.25(1)	1.08(1)
	Fe0.80C0.20	11.2(1)	1.25(1)	1.01(1)
	Fe0.70H0.30	1.9(1)	1.25(1)	0.97(1)
	Fe0.73Ni0.05O0.22	0.8(1)	1.25(1)	1.10(1)
	Fe0.80Ni0.05Si0.15	1.3(1)	1.06(1)	1.10(1)
	Fe0.76Ni0.05S0.19	0.6(1)	1.25(1)	1.14(1)
	Fe0.75Ni0.05C0.20	11.8(1)	1.25(1)	1.07(1)
	Fe0.64Ni0.04H0.32	1.1(1)	1.25(1)	1.03(1)
	Fe0.69Ni0.09O0.22	0.7(1)	1.25(1)	1.15(1)
	Fe0.74Ni0.10Si0.16	1.4(1)	1.05(1)	1.17(1)
	Fe0.71Ni0.10S0.19	1.2(1)	1.25(1)	1.21(1)
	Fe0.7Ni0.09C0.21	10.4(1)	1.25(1)	1.13(1)
	Fe0.6Ni0.08H0.32	0.9(1)	1.25(1)	1.09(1)
$6500$ K	Fe0.82O0.18	4.8(1)	1.25(1)	1.02(1)
	Fe0.88Si0.12	1.7(1)	1.09(1)	1.02(1)
	Fe0.85S0.15	0.8(1)	1.25(1)	1.05(1)
	Fe0.84C0.16	7.7(1)	1.25(1)	1.00(1)
	Fe0.74H0.26	0.1(1)	1.25(1)	0.97(1)
	Fe0.77Ni0.05O0.18	0.9(1)	1.25(1)	1.08(1)
	Fe0.82Ni0.05Si0.13	1.8(1)	1.08(1)	1.08(1)
	Fe0.79Ni0.05S0.16	0.2(1)	1.25(1)	1.12(1)
	Fe0.78Ni0.05C0.17	7.5(1)	1.25(1)	1.06(1)
	Fe0.69Ni0.04H0.27	2.3(1)	1.25(1)	1.02(1)
	Fe0.71Ni0.1O0.19	2.2(1)	1.25(1)	1.15(1)
	Fe0.77Ni0.1Si0.13	1.3(1)	1.08(1)	1.15(1)
	Fe0.74Ni0.1S0.16	0.6(1)	1.25(1)	1.18(1)
	Fe0.72Ni0.1C0.18	7.7(1)	1.25(1)	1.13(1)
	Fe0.64Ni0.09H0.27	1.0(1)	1.25(1)	1.09(1)

with misfits and the ρ and V_P of best-fit compositions along two adiabats are shown in Figure 2.

Figure 1 and Figure 2 and Table 2 indicate that among the best-fit compositions, the misfit of the Fe(-Ni)-C system is distinctly large. This suggests that carbon could be eliminated from the major LE in the OC, though there is possibility that the situation could change in ternary or higher-order multicomponent systems. In contrast, the misfits of the best-fit compositions of the other LEs have marginal differences, which are almost indistinguishable from each other within the computational uncertainty (shaded regions in Figure 2). There have been some previous studies which constrained the OC composition within the similar manner, taking two LEs into account, suggesting an oxygen-depleted OC [12] or oxygen-rich OC (3.7 wt. % O, 1.9 wt. % Si) [23]. However, according to the results of the present study, the difference between the misfits of best-fit composition models (Figure 1) is very small except for C, indicating that the information of ρ and V_P are insufficient to determine the OC composition uniquely. Therefore, some other information, e.g., melting phase relations, partitioning behavior between solids and liquids, the bulk earth (BE) compositional property and so on, are quite helpful to place further constraints on the LE composition, but all of these are not well understood at the moment.

The best-fit compositions vary depending on the setting of T_ICB since the amount of LEs required to reproduce the PREM decreases for higher T (Table 2). However, the misfit values are insensitive to temperature without any systematic variations, meaning that it is difficult to constrain T_ICB through this optimization. Based on the calculated results at two different T_ICB, the best-fit compositions are represented as a function of T_ICB within the first-order as follows:

$$X_{\mathrm{O}} \left(\mathrm{atom}\%\right)=-2.60\times {10}^{-3} T_{\mathrm{ICB}}\left(\mathrm{K}\right)+35.0\mathrm{for}\mathrm{O},$$

(1)

$$X_{\mathrm{Si}} \left(\mathrm{atom}\%\right)=-1.80\times {10}^{-3} T_{\mathrm{ICB}}\left(\mathrm{K}\right)+24.4\mathrm{for}\mathrm{Si},$$

(2)

$$X_{\mathrm{S}} \left(\mathrm{atom}\%\right)=-2.27\times {10}^{-3} T_{\mathrm{ICB}}\left(\mathrm{K}\right)+30.6\mathrm{for}\mathrm{S},$$

(3)

and

$$X_{\mathrm{H}} \left(\mathrm{atom}\%\right)=-3.33\times {10}^{-3} T_{\mathrm{ICB}}\left(\mathrm{K}\right)+48.3\mathrm{for}\mathrm{H}.$$

(4)

These are found to change only marginally when Ni is incorporated. Here, the Fe(-Ni)-C systems are eliminated since they show larger misfits than others.

Since the T_ICB should correspond to the freezing temperature of the OC liquid, melting phase relations of the Fe-LEs systems at the P_ICB are quite important to place further constraints on the OC composition. There are however almost no available data with enough quality at the moment. Some experiments, though all were conducted at substantially lower pressures than P_ICB, suggest that the eutectic temperature of Fe-FeS system is more than 1000 K lower than the melting temperature of pure Fe [19,34], while solidus or eutectic temperatures are not quite different (within a few 100 K) in the Fe-FeSi [35] and Fe-O systems [36]. A large drop in the melting temperature might also be expected in the Fe-H system [37]. The T_CMB is usually thought to be ~4000 K [38,39] and its adiabatic extrapolation leads to ~5200 K for the T_ICB [7]. This might be ~1000 K lower than the T_M of pure Fe at the P_ICB, suggesting S and H as the potential LEs in the OC. But even larger temperature drops could exist in the ternary or quaternary systems, so it is hard to exclude Si and O from the LE candidates based on this discussion.

Another point is the density jump across the ICB ($\mathsf{\Delta }{\rho }_{\mathrm{ICB}}$), which is reported seismologically to reach ~4.7% [32]. This observed $\mathsf{\Delta }{\rho }_{\mathrm{ICB}}$ is however too large to be reconciled simply by the solid–liquid transition of pure Fe⁴¹. Partitioning of LEs between solid and liquid phases is therefore thought to be required, namely LEs dissolving in the OC should strongly prefer liquid to solid. Again, nothing can be conclusive before the melting phase relations of the Fe-LEs systems are clarified at P_ICB, but extrapolations of experimental knowledge obtained at lower pressures suggest that the strong partitioning occurs in the Fe-O system [36] but not in the Fe-S [19,34], Fe-Si [35,40], and Fe-H systems [37].

In this study, we select the PREM model as a reference P-wave velocity and density of the OC. It is however well-known that the velocity structure of the earth’s interior depends on the reference model. For example, the AK135 model [41,42] has the P wave velocity different from the PREM model in particular at the uppermost and lowermost outer core. Although contrasts between the AK135 and PREM reach ~0.11 km/s and ~0.065 km/s at the uppermost and lowermost part respectively, these make no significant changes in the insights obtained from our analyses.

2.3. Bulk Earth Composition

The BE composition can be affected by the LE compositions of the OC. We next examine what major element composition of the BE are lead from each best-fit composition for the OC (Table 2). In this modeling, the mantle composition is assumed to be the pyrolytic [43], and the IC and crust are ignored because of its negligibly small masses. The chemistry of the mantle, in particular of the lower mantle, is still under debate, but some recent studies similarly suggested the pyrolytic one might be reasonable [44,45].

Table 2 shows the Mg/Si and Mg/Fe ratios of the BE expected from the best-fit compositions for the core combined with the pyrolytic composition for the mantle, which are calculated using the total mass of the earth, weight and atomic % of major elements (Mg, Fe, Si, and O) in the pyrolytic model for the mantle and the optimized composition models for the core. CI chondrite, one of the major candidates of earth’s building block, is known to have the Mg/Si and Mg/Fe ratios of ~1.05 and ~1.23, respectively [43]. Table 2 shows that this Mg/Si ratio is achieved only when Si is the major LE in the OC, but no case can explain the Mg/Fe ratio. Enstatite (EH) chondrite is another candidate of earth’s building block, which is known to have the Mg/Si and Mg/Fe ratios of ~0.79 and ~0.91. Very Si-rich core and mantle are required to explain this small Mg/Si and Table 2 shows that such composition is incompatible with the observations of the OC. An Mg/Si value similar to CI chondrite, ~1.03, is proposed in the OCCAM model [46] with an Mg/Fe ratio of ~1.12. These ratios are close to the values expected for our Si-bearing best-fit compositions. In summary, Si is a geochemically plausible candidate for the major LEs in the core. But if the lower mantle is assumed to be MgSiO₃-dominant, the Mg/Si ratios of the BE expected with the best-fit composition for the core decreases by ~0.2. Then, the Mg/Si and Mg/Fe ratios of all the best-fit compositions except for Si and H-bearing cases match the ratios of CI chondrite and OCCAM model. In this case, Si-rich and H-rich OC with MgSiO₃-dominant lower mantle lead to a too Si-rich and Fe-poor BE composition, respectively.

3. Conclusions

Ab initio thermoelasticity of Fe-Ni-LEs alloy liquids in the whole OC P, T condition indicates that all the LEs have the effects to decrease ρ and increase V_P of pure Fe, but the effects are counterintuitively larger for the Si and S incorporations than for the O, C, and H incorporations. Any best-fit alloy composition except the C-rich case can reproduce the ρ and V_P of the actual OC in the comparable level, so that the information of ρ and V_P only are insufficient to determine the OC composition uniquely. Melting phase relations and LE partitioning in the Fe-Ni-LE systems at the P_ICB are therefore essential to place a tighter constraint on the OC chemistry. The Si-rich best-fit composition for the core with an assumption of the pyrolytic mantle predicts the CI chondritic BE composition, but the O and S-rich best-fit compositions for the core with the MgSiO₃-dominant mantle also lead to the similar chemistry for the BE. The H-rich best-fit composition however causes a distinct deficit of Fe for the BE. In future studies, it might be important to investigate correlations between LEs in higher-order multicomponent systems, which are ignored in this study.

4. Methods

4.1. Ab Initio Molecular Dynamics Simulations

To determine the P-V-T equation of state (EoS) of liquid iron-light element alloys, total internal energy (E) and total pressure (P) are calculated by means of the AIMD technique within the canonical (NVT) ensemble in the same manner as our previous study [7] using a PWSCF code [47] for electronic structure with an original implementation of the constant temperature molecular dynamics (MD) module [48]. The simulations are performed on binaries and ternaries, (Fe-Ni)_1−XO_X; (Fe-Ni)_1−XSi_X; (Fe-Ni)_1−XS_X; (Fe-Ni)_1−XC_X; (Fe-Ni)_1−XH_X, at different atomic fractions ($X_{\mathrm{O}}\le 0.5,X_{\mathrm{Si}}\le 0.3,X_{\mathrm{S}}\le 0.3,X_{\mathrm{C}}\le 0.3,X_{\mathrm{H}}\le 0.4$). Three Ni/(Fe + Ni) ratios of 0, 0.06 (consistent with the geochemically modeled value) [43], and 0.12 are examined.

The Newton’s equation of motion is numerically integrated by using the velocity Verlet algorithm with time steps of 1 fs (10⁻¹⁵ s) for the Fe-LE systems, which is the same for previous studies [7,17,33,49], and 0.5 fs for the Fe-H system. Some results (pure Fe and Fe_1−XO_X systems) are, in part, already reported in the previous studies [7,49]. MD cells basically contain 50 atoms as in our previous study [7] but 100 atoms for the optimized compositions, and T is controlled by the kinetic energy scaling method. The validity of the cell size with 50–100 atoms for liquid iron can be seen in previous calculations [17,33], where a minor variation of the melting temperature of iron (~100 K) was found with changing the cell size from 67 to 980 atoms. Thermodynamic properties of liquid iron were also found to be sufficiently converged for this cell size.

For electronic structure calculations, we apply the generalized gradient approximation (GGA) [50] to the exchange correlation functional instead of the local density approximation (LDA). This is essential since many previous studies reported that GGA shows significant improvements over LDA when it comes to correctly describing ground-state properties and compression behaviors for iron [51,52]. The ultrasoft pseudopotential and plane-wave basis set are used to describe electronic structures. Here, an electronic configuration of 3s²3p⁶3d^6.54s¹4p⁰ is pseudized with a sufficiently small core radius of 2.0 a.u. for Fe; 2s²2p⁴, with a core radius of 1.5 a.u. for O; 3s²3p⁴, with a core radius of 1.7 a.u. for S; 3s²3p²3d⁰, with a core radius of 1.4 a.u. for Si; 2s²2p², with a core radius of 1.1 a.u. for C; 1s¹, with a core radius of 0.8 a.u. for H; and 3s²3p⁶3d⁸4s²4p⁰ with a core radius of 2.0 a.u. for Ni by the Vanderbilt scheme [53] with non-linear core corrections [54]. We apply a kinetic energy cutoff of 50 Ry and spin polarization is not taken into account. These conditions are already well tested in our previous calculations [7,49,55] and are fairly similar to those in calculations by other groups [24]. Liquids in principle have no periodic structure, thus the Γ point only is sampled in our simulations. All the MD simulations are conducted in P,T condition from ~80 to ~500 GPa and from 4000 to 8000 K (Figure S1), which covers the whole P,T range of the core. Standard deviations in calculated T and P are found to be typically ~50 K and ~3 GPa at 5000 K and ~130 GPa and ~120 K and ~6 GPa at 8000 K and ~400 GPa, respectively.

4.2. EoS Analysis of Liquid Iron Alloys

The calculated E-P-V-T relations of liquid iron alloys are analyzed using a single EoS model basically identical to the one in our previous study [7]. For the isothermal part at a reference temperature T₀, we used the Vinet (Morse-Rydberg) Equation,

$$P_{T_0}\left(V\right)=3K_{T_0}{\left(\frac{V}{V_0}\right)}^{-\frac{2}{3}}\left[1-{\left(\frac{V}{V_0}\right)}^{\frac{1}{3}}\right]\exp \left\{\frac{3}{2}\left(K_{T_0}^{’}-1\right)\left[1-{\left(\frac{V}{V_0}\right)}^{\frac{1}{3}}\right]\right\}$$

(5)

Here, $K_{T_0}$ and $K_{T_0}^{’}$ are the isothermal bulk modulus and its pressure derivative at zero pressure at T₀. The internal thermal energy is represented by a second-order polynomial of temperature with a volume dependent second-order coefficient,

$$E_{th}\left(V,T\right)=3n{k}_B\left[T+e_0{\left(\frac{V}{V_0}\right)}^g{T}^2\right]$$

(6)

where $k_{\mathrm{B}}$ is the Boltzmann constant and n is the number of atoms per formula unit. The first term corresponds to the phonon energy (atomic contribution), while the second term represents the electronic contribution.

The thermal pressure is linked with the internal thermal energy by Grüneisen parameter $\gamma$ as in the following Equation,

$$P_{th}\left(V,T\right)=\frac{1}{V}{\displaystyle \int }\gamma \left(V,T\right)d{E}_{th}\left(V,T\right).$$

(7)

We employ the highest temperature of 8000 K in the present calculations as a reference temperature T₀ in order to constrain the reference isotherm as tightly as possible within the broad pressure range. Although in the previous study [7], the following functional form for $\gamma \left(V\right)$

$$\gamma \left(V\right)={\gamma }_0^{}\left\{1+a\left[{\left(V/V_0\right)}^b-1\right]\right\}$$

(8)

was employed, we realized that three additional parameters ${\gamma }_0^{},$ a, and b make the least-square analyses less stable and less systematic. Instead, in this study, we assume the $\gamma$ to be constant for each composition. The previous study [7] reported that the variation of $\gamma$ of pure Fe is 0.2 only from 100 to 400 GPa and from 4000 and 7000 K and we confirmed that this small variation of $\gamma$ does not affect the results of analyzed thermoelasticity. Consequently, the present EoS model requires six parameters in total (V₀, $K_{T_0}$, $K_{T_0}^{’}$, $\gamma$, e₀, and g) to calculate P at a given V, T. These parameters are determined by least squares analyses on the datasets obtained from the AIMD calculations. The derived EoS parameters for best-fit compositions are summarized in Table S1.

Derivative quantities of EoS such as K_T and α are obtained based on the thermodynamic definitions as (∂P/∂V)_T = −K_T/V and (∂P/∂T)_V = αK_T, respectively. K_T is then converted to K_S as

$$\frac{K_{\mathrm{S}}}{K_{\mathrm{T}}}=1+\alpha \gamma T,$$

(9)

and adiabatic temperature gradient is computed using the relationship,

$${\left(\frac{\partial T}{\partial P}\right)}_S=\frac{\alpha VT}{C_P}=\frac{\gamma T}{K_{\mathrm{S}}}.$$

(10)

V_P is then calculated for each composition (Figure S1) as $V_{\mathrm{P}}=\sqrt{\frac{K_{\mathrm{S}}}{\rho }}$ along two different adiabats explained below.

4.3. Adiabats

The adiabatic temperature profile is calculated numerically by integrating Equation (7) from the ICB pressure. ρ and V_P are calculated along the adiabats anchored by two possible ICB temperatures: T_ICB = 5000 K and T_ICB = 6500 K. The former T_ICB is found to give ~3700 K at 136 GPa, which is close to a proposed core-mantle boundary temperature [38,39]. The latter corresponds to the melting temperature (T_M) of pure iron at 329 GPa [17,18], which would be close to the upper bound of ICB temperature since T_M of iron-LE alloys are expected in general to be lower than the T_M of pure Fe.