In our ab initio MD simulations, atomic forces were obtained from the electronic states calculated by the projector–augmented plane–wave (PAW) method (Blochl, 1994; Kresse and Joubert, 1999) within the framework of density functional theory (DFT). In the calculations, the generalized gradient approximation formulated by Perdew, Burke, and Ernzerhof (GGA–PBE) was used for the exchange–correlation potential (Perdew et al., 1996). The cutoff energies of the plane wave were 30 and 300 Ry for the pseudo–wavefunctions and pseudo–charge density, respectively. These cutoff energies were the same as those used in a previous study (Ohmura et al., 2020). The energy functional was minimized using an iterative scheme based on the preconditioned conjugate–gradient method (Kresse and Hafner, 1994; Shimojo et al., 2001). The gamma point was used for the k–point sampling. As valence electrons, 3d, 4s, and 4p states of Fe, 1s state of H, 2s and 2p states of C and O, and 3s, 3p, and 3d states of Si and S were used. We used 200 atoms in a cubic supercell under periodic boundary conditions. The thermodynamic states investigated in this study were 140 GPa and 5000 K. Properties of liquid Fe–LE (LE = H, O, C, Si and S) binary systems have been previously investigated at 140 GPa and 5000 K using AIMD (Ohmura et al., 2020). To be able to compare the results of the present work with those of the binary systems, we have used the same pressure and temperature (140 GPa and 5000 K). Four alloy compositions were studied: Fe160H20O20 (Fe–H–O), Fe160C20O20 (Fe–C–O), Fe160Si20O20 (Fe–Si–O), and Fe160S20O20 (Fe–S–O). Since this study does not aim to clarify the properties of liquid Fe alloy under the actual outer core condition, the compositions of the four liquids in this study were also not set to reproduce the density of outer core. To obtain each liquid state, we first prepared pure liquid Fe (liquid Fe₂₀₀), and 40 Fe atoms in the liquid were randomly replaced by LEs. To make the system reach a completely disordered state without the effects of the initial configuration, we carried out MD simulations for approximately 2 ps at a temperature 10,000 K. Then, we decreased the temperature of the systems gradually to a target temperature to 5000 K. After these preprocessing, to determine the density of each liquid system, a constant–pressure MD simulations was performed for 2.4 ps. By averaging over 1.8 ps after equilibration, we determined the densities (volumes) in each liquid. The structural properties were investigated using MD simulations in a canonical ensemble (Nose, 1984; Hoover, 1985). Using the Nose–Hoover thermostat technique, the equations of motion were solved via an explicit reversible integrator (Tuckerman et al., 1992) with a time step of Δt = 1.2 fs for all systems except the system including hydrogen, in which the time step was 0.48 fs. The quantities of interest were obtained by averaging over approximately 4.8–9.6 ps after an initial equilibration taking approximately 2 ps. In all the calculations, spin polarization was not considered.

Figure 1 shows the densities of liquid Fe–H–O, Fe–C–O, Fe–Si–O, and Fe–S–O at 140 GPa and 5000 K obtained from the simulations. The densities of liquid Fe and Fe160O40 (Fe–O) at the same P,T are also shown as solid and dashed lines, respectively. The densities of all ternary liquids are smaller than that of liquid Fe. However, in comparison with liquid Fe–O, the densities of liquid Fe–H–O and Fe–C–O are larger than that of liquid Fe–O, while the densities of liquid Fe–Si–O and Fe–S–O are smaller than that of liquid Fe–O. At a lower temperature condition (T = 3700 K), the same behaviors can be observed (Supplementary Figure S1).

The partial pair distribution functions gαβ(r) for liquid Fe–H–O (A), Fe–C–O (B), Fe–Si–O (C), and Fe–S–O (D) are shown in Figure 2. The g _Fe–Fe(r) has the first peak at approximately 2.2 Å and the second peak at approximately 4 Å in all the liquid systems. For g _O–O(r), the profile is almost the same in each liquid, in which we can see the first peak at approximately 2.2 Å and a small second peak at approximately 3.5 Å. For the interactions between Fe and light elements, the first peak positions of Fe–H, Fe–C, and Fe–O are smaller than that of Fe–Fe, and the first peak positions of Fe–Si and Fe–S are almost the same as that of Fe–Fe. These features are also seen in the Fe–LE binary systems (Ohmura et al., 2020). The results of densities and gαβ(r) indicate that H, C, and O are interstitial, and that Si and S show substitutional properties in the liquid Fe–LE–O ternary systems as well as the liquid Fe–LE binary systems.

Regarding the LE–LE interactions, we can see the different features of each liquid. For homoatomic interactions, H–H does not have clear correlations. In contrast, C–C shows strong interactions in the sense that there is a first peak at approximately 1.5 Å, which is the smallest first peak position in all gαβ(r). The strong C–C interaction is found in liquid Fe–C–H system as well (Belonoshko et al., 2015). The first peak positions of g _O–O(r) are almost the same as those of g _Fe–Fe(r) in all liquids. However, the second–peak positions are smaller than those of g _Fe–Fe(r). For Si–Si interactions, the first–and second–peak positions of Si–Si are almost the same as those of Fe–Fe or Fe–Si. These features align with the results of earlier studies for liquid Fe–O or Fe–Si binary (Alfe et al., 1999; Morard et al., 2014; Posner and Steinle–Neumann, 2019; Ohmura et al., 2020) and for Fe–Si–O ternary systems (Pozzo et al., 2013; Arveson et al., 2019; Huang et al., 2019). For S–S interactions, the first peak is broad, and its position is much farther than those of other gαβ(r) in liquid Fe–S–O. In terms of the heteroatomic interactions, H–O, C–O, and S–O have similar behaviors in the sense that there are broad first peaks. In contrast, the first peak of g _Si–O(r) is located at 1.8 Å and sharp, meaning strong interactions between Si and O.

To clarify the detailed bonding nature in the liquid Fe–LE–O ternary systems, we use the Mulliken population analysis, where the electronic wave functions are expanded in an atomic orbital basis set (Mulliken, 1955; Shimojo et al., 2008). Figure 3 shows the time–averaged distribution p α β ( O ¯ ) of the bond overlap population O i ∈ α , j ∈ β , which gives a semiquantitative estimate of the covalent–like bonding between the α and β atomic types in liquid Fe–H–O (A), Fe–C–O (B), Fe–Si–O (C), and Fe–S–O (D). It should be noted that, since the atomic–orbital basis used in the expansion of the wavefunctions is not unique, the absolute magnitudes of O i j are little physical meaning. However, the trends remain unchanged for any choice of atomic–orbital basis sets, and therefore we can discuss their relative variations meaningfully.

Figure 3A shows p α β ( O ¯ ) in liquid Fe–H–O. All of p α β ( O ¯ ) in liquid Fe–H–O do not show strong overlaps. In contrast, as shown in Figures 3B,C, p C − C ( O ¯ ) , p S i − S i ( O ¯ ) , and p S i − O ( O ¯ ) have values in the regions O ¯ > 1.0, while other p α β ( O ¯ ) are almost zero in the same region. For liquid Fe–S–O, p S − S ( O ¯ ) shows a relatively stronger interaction in the sense that there is a shoulder at approximately O ¯ = 0.5, but p S − S ( O ¯ ) = 0 in the region O ¯ > 1.0. Although the Si–O and Si–Si atomic distances are almost the same as those for Fe–O and Fe–Fe, respectively (Figure 2C), the bonding nature is different. Covalency are clearly found between Si–O and Si–Si, whereas not between Fe–O and Fe–Fe. These results indicate that C–C, Si–Si, and Si–O have covalent–like interactions in liquid Fe mixtures even under high–pressure conditions (>100 GPa).

In order to analyze the covalent–like interactions of C–C, Si–Si, and Si–O in more detail, we investigate the time evolution of the bond overlap population O i j ( t ) . Figure 4 shows O C 1 − C 2 ( t ) (overlap population between C1 and C2 atoms) displayed with the time evolution of the atomic distance between C1 and C2 atoms, d _C1–C2(t), and the snapshots of atomic configurations. At t = 0.00 ps, d _C1–C2(t) is more than 2.0 Å and O C 1 − C 2 ( t )   is zero, meaning no covalent bonding between C1 and C2. O C 1 − C 2 ( t ) increases when d _C1–C2(t) decreases, meaning that a covalent bond forms between C1 and C2 (the snapshot at t = 0.13 ps). After O C 1 − C 2 ( t ) and d _C1–C2(t) oscillate approximately 10 times, O C 1 − C 2 ( t ) becomes almost 0 and d _C1–C2(t) becomes more than 2.5 Å at approximately 0.40 ps, meaning that the covalent bond between C1 and C2 is lost, as shown in the snapshot at t = 0.40 ps. From the result, the lifetime of the C–C covalent bond is estimated to be approximately 0.3 ps at this P,T.

Figure 5 shows the time evolution of the distance and overlap population of focusing atomic pair (Si1 and Si2), dSi1–Si2(t) and the O Si 1 − Si 2 ( t ) , with the snapshots of atomic configuration in liquid Fe–Si–O. At t = 0.00 ps, dSi1–Si2(t) is almost 3.5 Å, and O Si 1 − Si 2 ( t ) is almost 0. Similar to the C–C interactions shown in Figure 4, with time, an Si–Si covalent–like bond forms (the snapshot at 0.70 ps), and the Si–Si bond oscillates several times, then the bond dissolves (at 1.30 ps).

In liquid Fe–Si–O, Si–O covalent bonding can also be observed. OSi–O(t), dSi–O(t), and snapshots of the atomic configuration are shown in Figure 6. At t = 0.00 ps, OSi–O(t) is almost 0, meaning no covalent bond between Si and O (the snapshot at t = 0.00 ps). OSi–O(t) begins to increase at approximately 0.05 ps, meaning the covalent bond formation between Si and O. Similar to C–C and Si–Si covalent bonds, after the Si–O bond oscillates several times, OSi–O(t) decreases drastically at approximately 0.22 ps. At the same time. dSi–O(t) increases, meaning that covalent bond between Si–O dissolves (the snapshot at 0.25 ps). From this result, we can estimate the lifetime of the Si–O covalent bond to be approximately 0.15 ps at this P,T.

Our calculations clearly indicate that in liquid Fe–Si–O and Fe–C–O, covalent–like interactions with finite lifetimes exist between C and C, Si and Si, and also Si and O under high–pressure more than 100 GPa.

Charge distribution of Si atoms in liquid Fe–Si–O.

Figure 7A shows the averaged distribution qα(Q) of the gross charge Q i ∈ α for α type atoms in liquid Fe–Si–O obtained by population analysis as same as O i j . As shown in Figure 7A, the peak positions of q _Fe(Q), q _Si(Q), and q _O(Q) are approximately Q = 0.1, –0.2, and –0.7, respectively. While the Si–O covalent bond affects the atomic charge of the Si atom, the C–C and Si–Si interactions have no effects on the atomic charge of C and Si (Supplementary Figure S2). To analyze q _Si(Q) in more detail, we decompose q _Si(Q) into q _Si–O(Q) and q _Si–nonO(Q). Here, q _Si–O(Q) means the distribution of the atomic charge for Si atoms coordinated to O atoms, which is displayed as a thick dashed line in Figure 7B and q _Si–nonO(Q) means the distribution of the atomic charge for Si atoms unbonded to O atoms, which is displayed as a thin solid line in Figure 7B. As shown in Figure 7B, q _Si–O(Q) is distributed to a relatively large Q region, which means that Si atoms have more positive charges when bonded to the O atom compared to the charges when the Si atoms are unbonded to O.

According to an earlier study (Arveson et al., 2019), liquid Fe–Si–O shows might have immiscibility between liquid Fe–Si–O and Fe–Si (there are two domains with and without oxygen in liquid Fe–Si–O) just above their solidus temperatures at Earth’s outer core pressures. Atomic configurations in our simulation also seem to be suggestive of this phase separation though the actual immiscibility region might exist at some lower temperatures (Supplementary Figure S3). The shift of atomic charge due to the Si–O covalent bond is possibly related to the gathering of oxygen atoms with negative charges around Si, leading to the separation of ionic Fe–Si–O and metallic Fe–Si domains in liquid Fe–Si–O. However, at the moment, the upper temperature bound of the immiscibility region at the outer core pressures is yet to be determined. In addition to these atomic perspectives, energetic information regarding thermodynamic stability of miscible and immiscible liquids has to be determined to understand the immiscibility of liquid Fe–Si–O in more detail.