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

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 (rFe–Fe) decreases with P and is sensitive to Ni (XNi) and C (XC) concentration, although the effects are opposite: rFe–Fe decreases with increasing XNi, but increases with increasing XC. Average rFe–C and rNi–C values also decrease with increasing XNi 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 rC-C values (~ 1.29 to ~ 1.57 Å), typical for rC-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 cm3/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 1 3 43 Page 2 of 14 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. 2003Rubie et al. , 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 lightelement 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 3 /mol = 1173 Å 3 /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 spinpolarized results by Wang et al. (2019) by our non-spinpolarized 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: where ρ b = X b V 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 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 = ∑ b N a−b ). The atomic trajectories and asymptotic slope of the timedependent 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), 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 0 ) along the atomic trajectories. The Arrhenius relation was used to determine the effect of P and T according to where D 0 is the pre-exponential diffusion coefficient, R is the universal gas constant, and Q the activation energy, ΔH and ΔV are the activation enthalpy and activation volume (i.e., P-and T-dependence terms), respectively, given as:

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) 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 Å. 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)

Fig. 1
Partial radial distribution functions (RDFs) for Fe-Ni alloys and end-member liquids at 3000 K for the large (7.121 cm 3 /mol; left column), intermediate (6.409 cm 3 /mol; middle column), and small (5.697 cm 3 /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 and decrease linearly with X Ni (R 2 = 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 endmember 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 Fefree 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 Page 7 of 14 43 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 compositionsexhibit 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 quasioctahedral 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, with and 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 3 C is stable at ambient P as a eutectic phase (Wood 1993) and Fe 7 C 3 forms at high P (Nakajima et al. 2009), Ni 3 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".

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 quasilinear 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).  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. 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

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  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) and Steinle-Neumann 2019) and ternary Fe-Ni-C liquids conducted over a considerably smaller compositional range (Wang et al. 2019). 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 3 /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 3 /mol for Ni diffusion in pure liquid Ni; the lowest obtained ΔV value was 0.24 ± 0.06 cm 3 /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 3 /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 3 /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 3 /mol for P < 5 GPa to 0.77 ± 0.06 and 0.61 ± 0.09 cm 3 /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 3 /mol for Fe, Ni, 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 and C, respectively, are more consistent with the results at 3000 K (Table SM2). A slight reduction of ΔV (i.e., D∕ 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.
Combining our results with those of Wang et al. (2019) to fit Eq. (4) yields D 0 = (1.00 ± 0.08) × 10 −7 m 2 /s, ΔH = 56 ± 2 kJ/mol, and ΔV = 0.49 ± 0.02 cm 3 /mol for Fe/ Ni and D 0 = (1.4 ± 0.1) × 10 −7 m 2 /s, ΔH = 50 ± 2 kJ/mol, and ΔV = 0.47 ± 0.02 cm 3 /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 The simulation results yield activation volumes mostly on the order of 0.4-0.6 cm 3 /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 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 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 3 /mol, ∆V C ~ 0.47 cm 3 /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 3 /mol, ∆V C ~ 0.61 cm 3 /mol) and thus slightly mismatch with the global fit