Compositional effects in the liquid Fe–Ni–C system at high pressure

Abstract

We performed molecular dynamics simulations based on density functional theory to systematically investigate the Fe–Ni–C system including (1) pure Fe and Ni; (2) binary Fe–Ni, Fe–C, and Ni–C; and (3) ternary Fe–Ni–C liquid compositions at 3000 K and three simulation volumes corresponding to pressure (P) up to 83 GPa. Liquid structural properties, including coordination numbers, are analyzed using partial radial distribution functions. Self-diffusion coefficients are determined based on the atomic trajectories and the asymptotic slope of the time-dependent mean-square displacement. The results indicate that the average interatomic distance between two Fe atoms (r_Fe–Fe) decreases with P and is sensitive to Ni (X_Ni) and C (X_C) concentration, although the effects are opposite: r_Fe–Fe decreases with increasing X_Ni, but increases with increasing X_C. Average r_Fe–C and r_Ni–C values also decrease with increasing X_Ni and generally remain constant between the two lowest P points, corresponding to a coordination change of carbon from ~ 6.8 to ~ 8.0, and then decrease with additional P once the coordination change is complete. Carbon clustering occurs in both binary (especially Ni–C) and ternary compositions with short-range r_C-C values (~ 1.29 to ~ 1.57 Å), typical for r_C-C in diamond and graphite. The self-diffusion results are generally consistent with high-P diffusion data extrapolated from experiments conducted at lower temperature (T). A subset of additional simulations was conducted at 1675 and 2350 K to estimate the effect of T on diffusion, yielding an activation enthalpy of ~ 53 kJ/mol and activation volume of ~ 0.5 cm³/mol.

Introduction

The cores of the Earth and terrestrial planets (Mercury, Venus, Mars, and the Moon) are composed of multi-component Fe–Ni alloys with substantial concentrations of at least one light element (LE) (e.g., S, O, Si, C, N, H) (e.g., Hirose et al. 2021). In this context, the partitioning behavior of an element between Fe-rich metal and silicate melt is a critical parameter for constraining the chemical and physical differentiation of core and mantle during planetary formation. Nickel has a larger metal-silicate partition coefficient than Fe (e.g., Corgne et al. 2008; Fischer et al. 2015) and is thus consistently found in iron meteorites, its concentration ranging from approximately 4.3 wt.% to as high as 34 wt.% (e.g., Vander Voort 2001) with the proportion increasing with decreasing oxygen fugacity (e.g., Benedix et al. 2014; Corgne et al. 2008).

The geoscience community has traditionally considered Ni to be the “geochemical twin” of Fe owing largely to their similar atomic mass and radius (e.g., Wells 1984), and Ni has therefore been frequently excluded from studies in mineral physics to minimize the complexity of both experiments and computational work on core analog materials. However, the appropriateness of this simplifying assumption is debatable: (1) the steelmaking industry has extensively documented the non-negligible influence of Ni on the properties of commercial steel (e.g., Knowles 1987; Keehan et al. 2006), particularly in terms of yield strength, hardness, and hardenability; (2) even in the binary liquid at ambient pressure (P), Fe and Ni have been shown to not mix ideally, with some excess volume measured (Watanabe et al. 2016); and (3) several geochemical studies have shown that the presence of Ni reduces the solubility of several elements (e.g., C, Re, Ir) in liquid Fe, and enhances that of others (e.g., S, Sn) (e.g., Tsymbulov and Tsemekhman 2001; Fonseca et al. 2011; Capobiano et al., 1999). Furthermore, Ni has been suggested to potentially affect the diffusion coefficient (D) of carbon in liquid Fe alloys based on a comparison of carbon diffusion data determined from high-P experiments in liquid Fe–C (Rebaza et al. 2021; Dobson and Wiedenbeck 2002) and simulations on liquid Fe–Ni–C (Wang et al. 2019). While previous studies have directly verified the effect of Ni on the material and geochemical equilibrium properties of Fe alloys relative to Ni-free compositions, as mentioned above, the influence of Ni on the transport properties of liquid Fe–LE alloys has not been investigated. In fact, Ni has been excluded from all experimental diffusion studies involving liquid iron alloys at high-P, which presently remain limited to pure liquid Fe (Dobson 2002) and binary alloys, Fe–S (Dobson 2000), Fe–O (Posner et al. 2017a), Fe–C (Rebaza et al. 2021; Dobson and Wiedenbeck 2002), Fe–Si, and Fe–Cr (Posner et al. 2017b).

Because Ni is expected to be present in the core of the Earth and other terrestrial planets, a comprehensive understanding of its influence on the properties of liquid Fe alloys is of great significance. Liquid metal transport properties (e.g., diffusion, viscosity) are crucial parameters for constraining the kinetics of mass exchange and extent of chemical equilibrium between core-forming metals and silicate melt during core formation in a magma ocean (e.g., Rubie et al. 2003, 2015), as well as ongoing reactions across the core-mantle boundary (e.g., Jeanloz 1990; Knittle and Jeanloz 1991) including the formation and stability of a light-element enriched layer at the top of the core (e.g., Gubbins and Davies 2013). A direct comparison of the experimental diffusion data in liquid Fe–C by Rebaza et al. (2021) and Dobson and Wiedenbeck (2002) with the computational results on Fe–Ni–C by Wang et al. (2019) is not possible largely because (1) the computational study was conducted at a single temperature (T), several hundred K lower than the experiments and possibly below the melting point for P > 10 GPa (e.g., Mashino et al. 2019); (2) the T dependence of D_C in liquid Fe–Ni–C required for extrapolation is unknown; and (3) the C concentrations vary between the different studies (~ 10–25 at.%), which has been shown to influence D_C in carbon steel (e.g., Liu et al. 1991). Furthermore, Dobson and Wiedenbeck (2002) reported a relatively strong effect of P on D_C, whereas the C diffusion data of Rebaza et al. (2021) showed a practically negligible P effect, which they attributed to differing C contents but did not test.

To systematically address this question, we investigated the Fe–Ni–C system at 3000 K using density functional theory (DFT) based molecular dynamics (MD) simulations on (1) pure liquid Fe and Ni; (2) binary Fe–Ni, Fe–C, and Ni–C; and (3) ternary Fe–Ni–C compositions. The Ni (X_Ni) and C (X_C) contents were varied between 1 and 90 at% and 1 and 20 at%, respectively. Three simulation volumes were set up to determine the effect of P on each composition, and a subset of simulations were run at 1675 and 2350 K to determine the effect of T.

Methods

Molecular dynamics simulations were performed in the canonical (N‐V‐T) ensemble within Kohn‐Sham DFT to compute the electronic structure of the liquid. The Vienna ab initio simulation package (Kresse and Furthmüller 1996) was used to calculate total energy and Hellmann‐Feynman forces. The electronic density at each simulation step was determined by the projector augmented wave formalism (Kresse and Joubert 1999) with the generalized gradient approximation to exchange and correlation (Perdew et al. 1996).

Following our previous studies (Posner et al. 2017a, c; Posner and Steinle-Neumann 2019), non-spin-polarized simulations were performed for cells with 150 atoms, reciprocal space was sampled only at the Γ‐point, and wave functions were expanded into plane waves with a cutoff energy of 550 eV. Cells were initially set up based on pure Fe configurations that were overheated and Fe atoms were then replaced by Ni or C in proportions listed in Table SM1 of the Supplemental Material. We consider three cell volumes V_x (7.121 cm³/mol = 1173 ų/cell), 0.9V_x, and 0.8V_x for all compositions. For the high-carbon (20 at%) compositions, runs at V_x were replaced with 0.7V_x because the V_x simulations yielded negative P values from the evaluation of the Hellmann–Feynman stress tensor (i.e., < 0 GPa). Simulations were conducted at T = 3000 K using a Nosé‐Hoover thermostat (Anderson 1980). Additional simulations were conducted at T = 1675 and 2350 K using the same composition as that reported by Wang et al. (2019) (Fe–7%Ni–20%C) for the following reasons: (1) to explore whether structures accessed at 1675 K represent a liquid despite being below the melting curve; (2) to see whether we reproduce the spin-polarized results by Wang et al. (2019) by our non-spin-polarized simulations; and (3) to determine the T dependence of diffusivities. To further test whether spin-polarized and spin-degenerate simulations yield equivalent results, we performed and analyzed one run for Fe–50%Ni–20%C at low P, and found that results do not differ in any significant way, despite the existence of magnetic moments on the Fe and Ni atoms. All simulations were run for at least 29 ps with a time step of 1 fs, and the first 2 ps were discarded to allow for equilibration.

Structural properties were investigated by analyzing the partial radial distribution functions (RDFs), namely, g_Fe–Fe(r), g_Fe–Ni(r), g_Fe–C(r), g_Ni–Ni(r), g_Ni–C(r), and g_C–C(r) and the corresponding partial and total coordination numbers. Considering an atom of species a, the probability of finding an atom of species b in a spherical shell (r, r + dr) is defined as:

$${g}_{a-b}\left(r\right)=4\pi {r}^{2}{\rho }_{b}{g}_{a-b}\left(r\right)dr,$$

(1)

where ρ_b = X_bV is the number density of species b with a mole fraction X_b and V is the volume per atom. The coordination number N_a-b represents the average number of nearest neighbors of species b surrounding an atom of species a and is calculated as

$$N_{a - b} = 4\pi \rho _b \int_0^{r_{a - b} } {r^2 g_{a - b} \left( r \right)dr} ,$$

(2)

where r_a–b is the position of the minimum after the first peak of g_a–b. In this paper, we refer to N_a as the total coordination of species a (i.e., ${N}_{a}={\sum }_{b}{N}_{a-b}$).

The atomic trajectories and asymptotic slope of the time-dependent mean-square displacement (MSD) in the simulation cells were used to check that cells are in the liquid state and to calculate the self-diffusion coefficient, D_a, for each species a following the Einstein relation (Allen and Tildesley 1991),

$${D}_{a}=\underset{t\to \infty }{\mathrm{lim}}\frac{1}{{N}_{a}^{*}}\sum_{i=1}^{{N}_{a}^{*}}\frac{\langle \left({r}_{i}\left(t+{t}_{0}\right)-{r}_{i}{\left({t}_{0}\right)}^{2}\right)\rangle }{6t},$$

(3)

where ${N}_{a}^{*}$ is the total number of atoms of species a, r_i(t) is the position of the ith atom at time t, and the angular brackets represent the ensemble average computed over different origin times (t₀) along the atomic trajectories. The Arrhenius relation was used to determine the effect of P and T according to

$$D={D}_{0}\mathrm{exp}\left(-\frac{Q}{RT}\right),$$

(4)

where D₀ is the pre-exponential diffusion coefficient, R is the universal gas constant, and Q the activation energy,

$$Q=\Delta H+P\Delta V.$$

(5)

ΔH and ΔV are the activation enthalpy and activation volume (i.e., P- and T-dependence terms), respectively, given as:

$$\frac{\partial \mathrm{ln}D}{\partial \left(1/T\right)}=-\frac{\Delta H}{R},$$

(6)

$$\frac{\partial \mathrm{ln}D}{\partial P}=-\frac{\Delta V}{RT}.$$

(7)

Results and discussion

Structural properties and compression mechanisms

The partial radial distribution functions (RDFs) of pure liquid Fe and Ni and Fe–Ni alloys at the three volumes for the simulations at 3000 K are shown in Fig. 1, and those of Fe–C, Ni–C, and Fe–Ni–C in Fig. 2 for simulations with 4 at% C. Equivalent RDF plots of the Fe–C, Ni–C, and Fe–Ni–C simulations with 1 and 20 at% C, in addition to those at lower T, are provided in Figs. SM1–SM3 of the Supplemental Material. The radial position (r) and intensity (I) of the first RDF peaks are summarized for each atomic pair in Figs. 3 and 4, respectively, as a function of P, and the total coordination number of each species is shown in Fig. 5. The results indicate that atomic pairs can generally be classified into one of three categories: (1) Fe/Ni–Fe/Ni framework with r = 2.18–2.45 Å; (2) shorter scale Fe/Ni–C framework–quasi-interstitial interactions with r = 1.83–1.97 Å; and (3) short-scale C–C interstitial–interstitial interactions with r = 1.29–1.57 Å.

Fig. 1

Partial radial distribution functions (RDFs) for Fe–Ni alloys and end-member liquids at 3000 K for the large (7.121 cm³/mol; left column), intermediate (6.409 cm³/mol; middle column), and small (5.697 cm³/mol; right column) simulation cell volumes. The radial positions (r) of the first RDF peaks, which represent the average nearest neighbor distances between two atom pairs, decrease with P (i.e., cell volume), as summarized in Fig. 3. The radial positions of the first r_Fe–Fe peaks gradually decrease with increasing X_Ni. The first RDF peak intensities (I) of all atomic pairs increase with P, as summarized in Fig. 4. I_Fe–Fe and I_Ni–Ni also increase with decreasing X_Fe and X_Ni, respectively, which implies stronger Fe–Fe and Ni–Ni ordering when present in dilute form

Fig. 2

Partial radial distribution functions (RDFs) for (a–c) Fe–4%C, (d–o) Fe–Ni–4%C, and (p–r) Ni–4%C liquid alloys at 3000 K at the same simulation cell volumes as in Fig. 1. Short-range C–C peaks (green curves, referred to as C–C′) indicate tight carbon clustering in most cells, with peak intensities (I) greater than 1 in some Ni-rich compositions (m, p, q), whereas most other compositions show more broad, low-I_C-C′ peaks, as shown in Fig. 4f as a function of P and Fig. SM4 as a function of X_Ni. The r_Fe–C and r_Ni–C values are approximately 16%–22% shorter than r_Fe–Fe in the respective compositions, and the difference generally increases with P. The positions of the first g_Fe–Fe(r) and g_Fe–C(r) peaks in the most Ni-rich ternary alloy (Fe–90%Ni–4%C) (m–o), respectively, occur at 2%–5% and 2%–3% shorter r than in the other Fe–C and Fe–Ni–C alloys, as summarized in Fig. 3. The I_Fe–C values are consistently higher than those of I_Ni–C, as summarized in Fig. 4

Fig. 3

Average radial positions (i.e., interatomic distances) of the first g(r) peaks for all investigated compositions as a function of P. The symbol color represents nickel content (X_Ni) and follows a basic rainbow scheme with increasing X_Ni: red, orange, green, blue, and purple symbols represent X_Ni =  ~ 0.01, ~ 0.05, 0.25, 0.5, and 0.9, respectively. The symbol shape represents carbon content (X_C): crosses, circles, squares, and triangles represent X_C = 0, 0.01, 0.04, and 0.2, respectively. The symbol filling distinguishes if Fe AND Ni are included in the liquid composition (unfilled) versus if Fe OR Ni are present (filled). The dotted trend lines in a–c show the linear fit to r_Fe–Fe in pure Fe (referred to as r_Fe–Fe*), and the dashed trend lines in (d) and (e) show the step-wise fit to r_Fe–C for Fe–4%C. Average 2σ error bars are shown in the lower corners of each panel. The results show that (a–c) r_Fe–Fe, r_Fe–Ni, and r_Ni–Ni generally decrease with P, whereas (d, e) r_Fe–C and r_Ni–C values, which are 16%–22% shorter than r_Fe–Fe and r_Ni–Ni, respectively, remain mostly constant between the two lowest P points, coinciding with a coordination change of carbon from ~ 6.8 to ~ 8.0 (Fig. 5), and then decrease upon further P increase once the coordination change is complete. In contrast, the r_C-C′ values in (f), which represent the position of the short-range carbon clustering peak (i.e., not the broad second peak near ~ 3.5 Å), generally increase with P, which reflects the larger average void size for the higher coordination carbon. The average distances between all atomic pairs except for C–C′ generally decrease with X_Ni, particularly for r_Fe–Fe, r_Fe–C, and r_Ni–C. The r_Fe–Fe and r_Fe–Ni values also increase with X_C. Some Ni-poor compositions show no Ni–Ni peak or only very large-r peaks (~ 6.2 Å; Fig. 1f, inset) and are therefore not included here

Fig. 4

Intensities of the first g(r) peaks (I), which generally reflect the degree of chemical ordering. The dotted trend lines in (a–c) show the linear fit to I_Fe–Fe values in pure liquid Fe (i.e., I_Fe–Fe*). The dashed trend lines in (d) and (e) show a linear fit to I_Fe–C values in liquid Fe–4%C. (a) For most compositions, I_Fe–Fe are similar and increase at a similar rate with P, with the notable exceptions of the alloys with X_Ni = 0.9, which show ~ 21%–64% higher I_Fe–Fe than I_Fe–Fe* (inset) and the difference generally increases with P. (b) I_Ni–Ni are generally greater than I_Fe–Fe* and increase with decreasing X_Ni, whereas (c) I_Fe–Ni vary both above and below I_Fe–Fe*. Some Ni-poor (X_Ni ≤ 0.05) simulations show exceptionally high I_Ni–Ni values, similar to the very high I_Fe–Fe values in the Fe-poor compositions in (a). As in Fig. 3b, I_Ni–Ni values are not shown for Ni-poor simulations with peaks occurring at very large r (e.g., Fig. 2f, inset). (d) I_Fe–C values increase with X_Ni, which implies higher Fe–C ordering with decreasing Fe concentration. In contrast, (e) I_Ni–C are generally insensitive to composition and lower than I_Fe–C. (f) I_C–C′ values of the short-range g_C-C′(r) peak are mostly very low, except for the Ni-rich compositions, and generally decrease with P and increase with X_Ni (Fig. SM4)

Fig. 5

Total coordination numbers for (a) Fe, (b), Ni, and (c) C. Fe and Ni are both generally close packed over the full range of study, whereas C coordination increases from ~ 6.2–7.4 (quasi-octahedral) to ~ 8.1–8.6 (B2 packing) with increasing P. Fe-rich compositions tend to show slightly higher N_Fe than Fe-poor compositions, whereas a compositional trend is less notable for N_Ni and N_C, although N_Ni in pure liquid Ni is generally higher than N_Ni in most other compositions

The compression behavior of liquids is dominated by a reduction of interatomic distances, increased coordination numbers, or a combination of the two. Although only three simulation volumes were investigated here owing to this study’s primary focus on compositional effects, important insight on the effect of P can be obtained by a careful inspection of the response of the liquid structural properties to isothermal compression. We therefore discuss the results in terms of the interatomic spacings (r), compression rates (dr/dP), intensities (I) of the first g(r) peaks, and coordination changes (dN_a/dP), with r and N_a values listed in Table SM1.

Fe/Ni substitutional atoms

The average interatomic distances (r) between all substitutional atom pairs (i.e., r_Fe–Fe, r_Fe–Ni, r_Ni–Ni) generally decrease with P (Fig. 3a–c). The r_Fe–Fe values are sensitive to both X_Ni and X_C, but the effects are opposite: r_Fe–Fe decreases with X_Ni but increases with X_C (Fig. 3a). The compression rates (dr_Fe–Fe/dP) are considerably lower in the alloys (− 1.51 ± 0.26 mÅ/GPa) than in pure liquid iron (labeled r_Fe–Fe* with dr_Fe–Fe*/dP =  − 2.25 mÅ/GPa) and decrease linearly with X_Ni (R² = 0.72). The I_Fe–Fe values (Fig. 4a) are mostly insensitive to composition, with the notable exception of the highly Ni-rich simulations, for which binary Fe–90%Ni and ternary Fe–90%Ni–4%C show ~ 46% and ~ 21%–64% higher I_Fe–Fe than I_Fe–Fe*, respectively, and the differences in the latter composition increase with P. This demonstrates that Fe–Fe pairs are highly ordered—and closer to one another, as mentioned above—when present in a relatively dilute concentration.

Average r_Ni–Ni values in pure liquid Ni (hereafter r_Ni–Ni*) are approximately 2% (~ 0.05 Å) longer than r_Fe–Fe* with a slightly higher compression rate (dr_Ni–Ni*/dP =  − 2.35 mÅ/GPa) (Fig. 3b). The addition of Ni therefore slightly modifies the structure of the metallic liquid and average quasi-interstitial void volumes, in contrast to Si, which directly substitutes for Fe in binary and ternary liquid alloys (e.g., Pozzo et al. 2013; Posner et al. 2017b; Posner and Steinle-Neumann 2019) with essentially zero size mismatch (i.e., r_Fe–Si = r_Fe–Fe = r_Fe–Fe*). The I_Ni–Ni values are consistently higher than those of I_Fe–Fe* (Fig. 4b) and generally decrease with X_Ni. Some low-Ni (X_Ni ≤ 0.05) compositions show exceptionally high I_Ni–Ni values, similar to the very high I_Fe–Fe values in low-Fe compositions (Fig. 4a), which demonstrates that Ni–Ni pairs are also somewhat more correlated when present in dilute form. On the other hand, some of the g_Ni–Ni(r) curves for compositions with X_Ni ≤ 0.02 do not present a well-defined Ni–Ni peak within the typical range for substitutional pairs (e.g., 2.1–2.5 Å) and are therefore excluded from the summary plots in Figs. 3b and 4b.

Average r_Fe–Ni (Fig. 3c) are generally slightly larger (up to ~ 5%) than r_Fe–Fe* and r_Fe–Fe in the respective compositions, which is similar with the results of our previous computational study involving Fe–4%Ni over a wider T range (Posner and Steinle-Neumann 2019), with the exception of Fe–90%Ni, for which r_Fe–Ni < r_Fe–Fe*. Similar to the r_Fe–Fe results, r_Fe–Ni values are also found to be sensitive to composition: r_Fe–Ni decreases with X_Ni and increases with X_C, which further demonstrates that Ni and C alloy components impose opposite structural effects on liquid Fe. The I_Fe–Ni values generally decrease with X_Ni and scatter both above and below the I_Fe–Fe* curve (Fig. 4c).

The coordination numbers of Fe (Fig. 5a) and Ni (Fig. 5b) range between ~ 11.6 and ~ 13.7, which indicates that both species are essentially close-packed over the full investigated P range. The coordination numbers of Fe and Ni in the pure end-member liquids generally define the upper limits of both, and compositions with X_Ni = 0.9 yield slightly lower N_Fe than the others.

Relation between Fe/Ni substitutional and C quasi-interstitial atoms

Average r_Fe–C (Fig. 3d) and r_Ni–C (Fig. 3e) are both approximately 16%–22% shorter than r_Fe–Fe and r_Ni–Ni, respectively; these differences decrease with P, which is similar to previous results for liquid Fe–4%C, Fe–4%O, and Fe–4%N by Posner and Steinle-Neumann (2019). Both r_Fe–C and r_Ni–C values are sensitive to Ni content, decreasing essentially linearly with X_Ni by up to 4% for X_Ni = 0.9 compared with values for X_Ni ≤ 0.02. This is also demonstrated in a comparison of Fig. 3d and e, which shows that Ni-free Fe–C compositions have the largest r_Fe–C and Fe-free Ni–C compositions have the smallest r_Ni–C. In contrast with the r_Fe–Fe and r_Fe–Ni results presented in the previous section, X_C has a negligible influence on r_Fe–C and r_Ni–C.

The r_Fe–C and r_Ni–C compression trends also differ from those of the Fe/Ni atomic pairs. In nearly all compositions, r_Fe–C values (Fig. 3d) remain constant over the two lowest-P simulations (i.e., dr_Fe–C/dP = 0), which is accompanied by a coordination increase of ~ 6.8 to ~ 8.0 (Fig. 5c), and then decrease with increasing P at a rate of approximately dr_Fe–C/dP =  − 0.5 mÅ/GPa (Fig. 3d) with a continuing increase in coordination number to ~ 8.6. A similar observation was previously reported for oxygen in liquid Fe–O over a larger P range (e.g., Ichikawa and Tsuchiya 2015; Posner et al. 2017c), with similar structural properties to carbon in liquid iron alloys (i.e., r_Fe–C ~ r_Fe-O, N_C ~ N_O; Posner and Steinle-Neumann 2019): r_Fe-O remained essentially constant with increasing P while N_O increased to ~ 8, then decreased monotonically upon a further P increase to ~ 300 GPa, over which range N_O stayed essentially fixed and did not exceed ~ 9.

Figure 4d shows that I_Fe–C values strongly depend on X_Ni, ranging from ~ 2.8–3.3 for Ni-poor compositions to 3.8–4.0 for X_Ni = 0.5 and as high as 5.0–5.2 for X_Ni = 0.9. This indicates that Fe–C pairs become increasingly ordered (and shorter) with decreasing Fe concentration, which is consistent with the Fe–Fe results presented in the previous section. The I_Ni–C values in Fig. 4e are consistently lower than I_Fe–C and insensitive to composition.

The r_Ni–C values are shown in Fig. 3e as a function of P alongside the r_Fe–C trend line of Fe–4%C for comparison. Especially at low P, the range of r_Ni–C values is wide and shows a clear decrease of r_Ni–C with increasing X_Ni, ranging 1.93–1.97 Å for simulations with X_Ni = 0.01–0.02 to 1.86–1.87 Å for simulations with X_Ni > 0.9. The results generally show that r_Ni–C > r_{Fe–C,Fe–4%C} for X_Ni < 0.5 and r_Ni–C < r_{Fe–C,Fe–4%C} for X_Ni > 0.5. The r_Ni–C compression trend is less straightforward than that of dr_Fe–C/dP: approximately half of the studied compositions show the same behavior as dr_Fe–C/dP (i.e., no change in r between the lowest two P points, followed by a decrease upon further increasing P), whereas most of the others decrease continuously with P. There is no systematic trend related to composition to explain this difference; both patterns occur in compositions with high and low Ni and C contents. In contrast, 11 of 15 compositions show a “plateau-decrease” trend of r_Fe–C with P. This difference might be explained by preferential packing of Fe around C.

Carbon quasi-interstitial clustering

As also reported in a previous study on binary liquid Fe–C (Ohmura et al. 2020), most compositions show short-range C–C clustering (Fig. 2) with peak positions (referred to as r_C-C′) on the order of ~ 1.29–1.57 Å, which is within the range of r_C-C in graphite (1.42 Å) and diamond (1.54 Å) (e.g., Harrison 1980). Differing from the dr/dP behavior of the Fe/Ni–Fe/Ni and Fe/Ni–C pairs, most r_C-C′ values are observed to continuously increase with P, with the exception of Fe–4%C, for which r_C-C′ decreases with P, and some carbon-rich (X_C = 0.2) compositions where r_C-C′ values remain largely unchanged over the full P range.

It is important to note that when present, most C–C′ peaks exhibit low g(r) intensities (I < 1; Fig. 4f), whereas others—particularly in the Fe-poor or Fe-free compositions—exhibit I > 1, reaching as high as ~ 2.1 in liquid Ni–20%C at low P, as summarized in Fig. SM4. We interpret the general increase of I_C-C′ with X_Ni to result from Ni–C repulsion that likely reduces the stability (i.e., solubility) of carbon in the liquid. The quantification of chemical potentials in the liquid (e.g., using techniques such as thermodynamic integration) is required to better resolve this question.

The overall positive dr_C-C′/dP trend correlates with generally decreasing I_C-C′ values with P (i.e., dI_C-C′/dP < 0; Fig. 4f), which is opposite to the behavior for the Fe/Ni–Fe/Ni and Fe/Ni–C pairs, where dI/dP > 0. These structural features are related to the gradually increasing occupation of carbon atoms from ~ 6.2 to 7.4 at low P (representing quasi-octahedral voids in the quasi close-packed metal framework) to ~ 8.1–8.6 at high P (representing an approximate B2 packing structure). The occurrence and intensity of carbon clustering generally decreases once the coordination change is completed.

Repulsion forces and preferential coordination environments

If Fe and Ni are assumed to be essentially equivalent species (i.e., “geochemical twins”) in liquid Fe–Ni and Fe–Ni–C alloys, their proportional contribution to a given total coordination number should be equal to their atomic fraction. If true, the following deviation f (in %) would be expected to be approximately zero,

$$f= \frac{A-C}{C}\times 100,$$

(8)

with

$$A=\frac{{N}_{a-\mathrm{Fe}}}{{N}_{a-\mathrm{Fe}}+{N}_{a-\mathrm{Ni}}},$$

(9)

and

$$C=\frac{{X}_{\mathrm{Fe}}}{{X}_{\mathrm{Fe}}+{X}_{\mathrm{Ni}}}.$$

(10)

We tested this hypothesis and found that it generally does not hold in liquid Fe–Ni and Fe–Ni–C, especially for compositions with X_Ni > 0.25 (Fig. 6). The results for carbon were most striking, showing a very strong tendency to coordinate with Fe over Ni, with a strong increase of f with Ni content. This can be rationalized by the metal–carbon phase relations: while Fe₃C is stable at ambient P as a eutectic phase (Wood 1993) and Fe₇C₃ forms at high P (Nakajima et al. 2009), Ni₃C has only recently been synthesized at much higher P (Fedotenko et al. 2021). More surprisingly, the f values for Fe and Ni differ from zero and the difference grows with X_Ni, increasing for Fe and decreasing for Ni. This implies that Fe and Ni more strongly concentrate around themselves than around each other, and a direct comparison of coordination ratios, (N_Fe–Fe/[N_Fe–Fe + N_Fe–Ni]) > X_Fe, (N_Ni–Ni/[N_Ni–Ni + N_Ni–Fe]) > X_Ni, (N_Ni–Fe/[N_Ni–Ni + N_Ni–Fe]) < X_Fe and (N_Fe–Ni/[N_Fe–Fe + N_Fe–Ni]) < X_Ni, confirms this observation, including for the binary Fe–Ni. This implies non-ideal mixing of Fe and Ni both in the C-free and C-bearing compositions, a result that is consistent with the observation of an excess volume of mixing in the Fe–Ni binary liquid at ambient P (Watanabe et al. 2016). Combined, these two observations negate the assumption that Fe and Ni are “geochemical twins”.

Fig. 6

Percent deviation f of the proportion of Fe atoms contributing to the total number of framework atoms (Fe + Ni) comprising the average coordination number of (a) Fe, (b) Ni, and (c) C, as defined in Eqs. (8)–(10). Dashed vertical lines indicate zero deviation, for which Fe and Ni are essentially identical substitutional species. Despite some scatter in (a) for high X_Ni, the results show that Fe atoms generally tend to preferentially coordinate with other Fe atoms and C atoms compared with Ni atoms, i.e., positive deviations in (a) and (c), and less preferentially around Ni atoms, i.e., negative deviations in (b)

Transport properties

Beyond an equilibration time of < 1 ps, the slopes of the MSD curves for all simulations (Figs. SM5 and SM6), as well as those at the lowest T (1675 K) considered for a comparison with the results of Wang et al. (2019), show quasi-linear behavior, which indicates that these are in the liquid state and that the MSD can be used reliably to determine a diffusion constant based on the Einstein relation (Eq. (4)). However, MD simulations with the number of atoms we use here are well known to require significant undercooling before a solid forms from the liquid (e.g., Luo et al. 2003; Braithwaite and Stixrude 2019).

Effect of pressure

The self-diffusion coefficients of Fe, Ni, and C at 3000 K in all compositions are shown as a function of P in Fig. 7a–c, respectively, and listed in Table SM2. The dotted lines shown in Fig. 7a, b indicate the linear trend of D_Fe with P in pure liquid. The results indicate that D_Fe values are generally equal to D_Ni and approximately 2–3 times slower than D_C for all investigated alloy compositions, which agrees with previous studies of both binary Fe–Ni and Fe–C (Posner and Steinle-Neumann 2019) and ternary Fe–Ni–C liquids conducted over a considerably smaller compositional range (Wang et al. 2019).

Fig. 7

Self-diffusion coefficients D of (a) Fe, (b) Ni, and (c) C for the investigated compositions at 3000 K. The dotted lines in (a) and (b) show the linear fit to D_Fe in pure liquid Fe. The dashed lines in (c) show curves extrapolated using the Arrhenian diffusion parameters reported from high-P experiments in liquid Fe–C by Dobson and Wiedenbeck (2002) (DW02) and Rebaza et al. (2021) (R21). The simulation results yield activation volumes mostly on the order of 0.4–0.6 cm³/mol for all species, which is consistent with previous DFT-MD studies and falls in between the experimental values reported by Dobson and Wiedenbeck (2002) and Rebaza et al. (2021). Average 2σ error bars are shown in the left corners of each panel

The pressure dependence (ΔV) at 3000 K was determined using a least-squares best fitting to Eq. (4), yielding 0.52 ± 0.06, 0.53 ± 0.08, and 0.46 ± 0.10 cm³/mol for Fe, Ni, and C, respectively, which is consistent with previous experiments and computations on Fe–Si, Fe–O, and Fe–Cr liquid compositions (Posner et al. 2017a–c). The highest obtained ΔV value was 0.77 ± 0.06 cm³/mol for Ni diffusion in pure liquid Ni; the lowest obtained ΔV value was 0.24 ± 0.06 cm³/mol for C diffusion in liquid Ni–4%C (Table SM2).

Figure 7c compares our D_C results with the diffusion data from high-P experiments involving binary Fe–C liquid alloys by Rebaza et al. (2021) and Dobson and Wiedenbeck (2002) extrapolated to 3000 K and the relevant P range (labeled dashed lines). The extrapolation curve of Rebaza et al. (2021) is comparable with the low-P D_C values for liquid Fe–C (filled black symbols) obtained here, but the mismatch increases with P owing to their very small reported ΔV. Dobson and Wiedenbeck (2002) reported a much larger activation volume and their extrapolation curve is substantially above and below the D_C obtained here at low- and high-P, respectively, with an overlap limited to values at approximately 30 GPa. Our findings indicate that an intermediate value between the two studies of ΔV ~ 0.5 cm³/mol best represents the diffusion of Fe, Ni, and C over the investigated P–T range.

Effect of temperature

The T dependence of diffusion was determined for Fe–7%Ni–20%C, which was also chosen for direct comparison with the results by Wang et al. (2019) at 1675 K and P ≥ 5 GPa. To complement our 3000-K results for this composition, we performed five additional simulations: two at 1675 K using cell volumes of V_0.8 (16 ± 2 GPa) and V_0.7 (61 ± 2 GPa), and three at 2350 K using cell volumes of V_0.9 (0.1 ± 2 GPa), V_0.8 (23 ± 2 GPa), and V_0.7 (69 ± 3 GPa). Our 1675-K diffusion results (Fig. 8) are in good agreement with those of Wang et al. (2019) despite the different treatment on spin-polarization, with ΔV values of 0.65, 0.75, and 0.61 cm³/mol for Fe, Ni, and C, respectively. Wang et al. (2019) performed a detailed isothermal study (1675 K) on liquid Fe–7%Ni–20%C at P up to 67 GPa and reported a substantial reduction of the ΔV of Fe/Ni and C diffusion at ~ 5 GPa, from 1.4 ± 0.3 and 1.2 ± 0.3 cm³/mol for P < 5 GPa to 0.77 ± 0.06 and 0.61 ± 0.09 cm³/mol for P > 5 GPa, respectively, the latter of which are consistent with our values. This change in diffusivity may be related to the onset of the loss of magnetic moments in the spin-polarized simulations of Wang et al. (2019), which can explain why we recover their diffusivities in the non-magnetic simulations, although Wang et al. (2019) report magnetic moments to persist throughout the compression range they consider. With the rapid decrease of effective magnetic moments with T shown in a DFT-MD study on liquid Fe (Korell et al. 2019), we do not expect that magnetism would influence the diffusion results at 2350 and 3000 K. The ΔV values at 2350 K of 0.54 ± 0.03, 0.54 ± 0.11, and 0.50 ± 0.01 cm³/mol for Fe, Ni, and C, respectively, are more consistent with the results at 3000 K (Table SM2). A slight reduction of ΔV (i.e., $\partial D/\partial P$) with T implies that the differences of D values between the metal atoms and carbon at 1675 and 3000 K increase with P, e.g., from ~ 0.7 log unit at 5 GPa to ~ 1.5 log units at 67 GPa.

Fig. 8

Self-diffusion coefficients D of (a) Fe/Ni and (b) C in liquid Fe–7%Ni–20%C at 1675, 2350, and 3000 K. Values reported by Wang et al. (2019) at 1675 K and ≥ 5 GPa are shown as open symbols and those obtained here are shown as filled symbols. The dashed curves show the global model fit (including points at all three T) according to Eq. (4) and the dotted curves show the individual P trends obtained at each isotherm. The activation volumes (i.e., slope of log D versus P) at 2350 and 3000 K are generally similar (∆V_Fe/Ni ~ 0.51 cm³/mol, ∆V_C ~ 0.47 cm³/mol) and the global model fit well reflects the P trend of the results, whereas those at 1675 K are slightly larger (∆V_Fe/Ni ~ 0.70 cm³/mol, ∆V_C ~ 0.61 cm³/mol) and thus slightly mismatch with the global fit

Combining our results with those of Wang et al. (2019) to fit Eq. (4) yields D₀ = (1.00 ± 0.08) × 10⁻⁷ m²/s, ΔH = 56 ± 2 kJ/mol, and ΔV = 0.49 ± 0.02 cm³/mol for Fe/Ni and D₀ = (1.4 ± 0.1) × 10⁻⁷ m²/s, ΔH = 50 ± 2 kJ/mol, and ΔV = 0.47 ± 0.02 cm³/mol for C. The results of this global fit are shown as dashed curves in Fig. 8 alongside dotted curves for the fits to each isotherm. The fitting results provide a good match with the diffusion coefficients at 2350 and 3000 K, but the P-slope of the global model at 1675 K is approximately a factor of two less steep than the isothermal fit for both Fe/Ni and C. A possible explanation for this discrepancy is that the 1675-K runs may actually involve a metastable liquid, as discussed above. The low-T results should therefore be applied with caution. Nevertheless, the model fits to Eq. (4) are comparable regardless of whether the 1675-K results are in- or excluded.

The obtained magnitudes of ΔH and ΔV are consistent with those previously reported for self-diffusion in binary Fe–Si and Fe–Cr liquids (Posner et al. 2017b). Therefore, these parameters are expected to be good approximations for the full range of pure, binary, and ternary liquids investigated here owing to the negligible effect of composition on D, although the effect of T was only tested for one composition.

Negligible effect of composition

The diffusive transport properties are found to be insensitive to composition over the P–T–X range of study, which does not support the hypothesis proposed by Rebaza et al. (2021) that the addition of Ni to liquid Fe–C affects carbon diffusivity. This discrepancy can be rationalized by two main factors: (1) The negligible pressure effect on D_C reported by the experiments by Rebaza et al. (2021) is not reproduced here or elsewhere (e.g., Dobson and Wiedenbeck 2002; Wang et al. 2019), which implies that the extrapolation of D_C using the parameters reported in Rebaza et al. (2021) to higher P than their experimental range (> 15 GPa) most likely result in an overestimate; (2) as mentioned in the previous section, the relatively low D values reported by Wang et al. (2019) are well reproduced here and may reflect liquid metastability at 1675 K, further increasing the disparity between their results with those of Rebaza et al. (2021). Nevertheless, the question of coupled compositional effects on transport properties in ternary or higher order systems will be interesting to test with respect to planetary core analogs, particularly if they contain sulfur or hydrogen.

Conclusions

Nickel has been excluded from all previous high-P diffusion experiments involving liquid iron alloys, as well as the majority of computational studies. The results presented here indicate that Ni and C additions to liquid Fe affect the liquid alloy structural properties. One might intuitively expect that higher concentrations of slightly larger Ni atoms (or smaller C atoms) would gently expand (contract) the Fe framework, i.e., increase (reduce) r_Fe–Fe; however, the opposite is observed: r_Fe–Fe decreases (increases) with increasing X_Ni (X_C). Fe and Ni also show preferential packing behaviors: they each prefer to pack around themselves than each other, and Fe more strongly coordinates around C than Ni. The effect of composition on the transport properties of the investigated alloys is essentially negligible, by contrast. The combined effect of Ni and other LEs (e.g., S, O, Si, H, N) on the structural and transport properties of multi-component liquid Fe alloys should be explored in the future, both computationally and in the laboratory, to more systematically understand the influence of alloy composition on the metal-silicate partitioning behavior of relevant LEs during core formation.