Carlos Pinilla

Aldemar De Moya

Segolene Rabin

Guillaume Morard

Mathieu Roskosz

Marc Blanchard

First-principles investigation of equilibrium iron isotope fractionation in Fe1-xSx alloys at Earth's core formation conditions

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Iron is the major element of solar system planetary cores, in relation with the nucleosynthesis sequence. Iron plays a major role in every stance of planetary formation and differentiation.

Geophysical approaches have confirmed the presence of iron for Earth's core but with almost 10% of light elements required to make up its mass [START_REF] Hirose | Composition and State of the Core[END_REF]. Traditionally the effects due to the presence of light elements in the core are being inferred by geophysical studies where the density and sound velocity propagation are studied as function of light element concentration to match those values reported from geophysical models such as PREM. However, recently it has been suggested that the Fe isotopic composition in the mantle could be used to infer the composition of Earth's core [START_REF] Bourdon | Isotope tracers of core formation[END_REF][START_REF] Craddock | Iron isotopic compositions of geological reference materials and chondrites[END_REF][START_REF] Craddock | Abyssal peridotites reveal the near-chondritic Fe isotopic composition of the Earth[END_REF][START_REF] Lesher | Iron isotope fractionation at the core-mantle boundary by thermodiffusion[END_REF][START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF][START_REF] Sossi | Iron isotope systematics in planetary reservoirs[END_REF]. This is based on the assumption that the presence of elements such as C, O, Si, S should affect the partitioning of iron isotopes between the mantle and the core. Originally, it has been suggested that core-mantle differentiation will leave imprints on the iron isotope signature of Earth's mantle because of the difference of Fe valence state and coordination in the mantle (Fe 2+ ) and in the core (Fe 0 metal) [START_REF] Polyakov | Equilibrium iron isotope fractionation at core-mantle boundary conditions[END_REF]. That study suggested that the mantle should be enriched in heavy isotopes, product of equilibrium isotope fractionation during differentiation. However, subsequent explanations argue that the bulk silicate Earth is chondritic in its iron isotopic composition and that any difference in isotope composition seen in basalts is due to fractionation during partial melting of the rock from which they have been formed [START_REF] Craddock | Abyssal peridotites reveal the near-chondritic Fe isotopic composition of the Earth[END_REF]. Indeed, isotope compositions in natural rocks are scattered in values, with 56 Fe -0.01 ‰ for carbonaceous chondrites, 56 Fe -0.1 to +0.1 ‰ for enstatite chondrites and 56 Fe +0.1 ‰ for basalts [START_REF] Craddock | Iron isotopic composition of geological reference materials and chondrites[END_REF][START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF][START_REF] Poitrasson | On the iron isotope heterogeneity of lithospheric mantle xenoliths: implications for mantle metasomatism, the origin of basalts and the iron isotope composition of the earth[END_REF]. These scattered values are therefore interpreted to suggest that the accessible mantle has a Fe isotope composition that is indistinguishable from chondritic composition (i.e. 56 Fe 0 ‰). Apart from expected evidence of the presence of light elements in the core due to difference in chemical bonds, one must also consider the effects of pressure and temperature. It is well known that equilibrium fractionation effects should vanish at high temperatures. However, in the case of lower mantle and core conditions where pressures are also high, these fractionation effects might still be important. Recently, [START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF] have carried out experiments and computational modelling on the effects of pressure and light elements (such as H, C, O) on the equilibrium isotope fractionation of 57 Fe/ 54 Fe at the core-mantle boundary. Their work has found a significant imprint on isotope fractionation due to light elements such as C and H ( 56 Fe 0.03-0.05 ‰) whilst elements such as O leave almost no imprint ( 56 Fe 0.007-0.01 ‰). Their results support the suggestions by [START_REF] Craddock | Abyssal peridotites reveal the near-chondritic Fe isotopic composition of the Earth[END_REF] that in order to have a mantle with a 56 Fe near 0‰, light elements such as C and H that provide a large imprint in fractionation should not be present in the core at a significant level. Additionally, it favours the presence of O as a light element in the core. However, the effects of other light elements such as S on equilibrium Fe isotope fractionation at core conditions are yet unclear as well as any pressure induced magnetic effects that can be of relevance. There has been recent experimental work by [START_REF] Ni | Heavy iron isotope composition of iron meteorites explained by core cristallization[END_REF]; [START_REF] Shahar | Sulfur-controlled iron isotope fractionation experiments of core formation in planetary bodies[END_REF] on the effects of S on Fe isotope fractionation between solid metal and liquid metal, or between metal and silicate. In both cases, sulfur enters into the metallic phases and preferentially into the liquid metal phase. [START_REF] Shahar | Sulfur-controlled iron isotope fractionation experiments of core formation in planetary bodies[END_REF] found that metal-silicate fractionation increases significantly with the sulfur content in the metal, whereas [START_REF] Ni | Heavy iron isotope composition of iron meteorites explained by core cristallization[END_REF] observed that the solid-liquid fractionation in this system does not depend on the sulfur content of the liquid metal. Obviously, much work needs to be done to understand these effects at temperature and pressure conditions of the terrestrial core, far away from experimental studies performed below 2 GPa. In this context, the difficulty to carry out experiments at high pressure and high temperature leads to the use of crystalline phases such as hcp Iron or I4 Fe3S as proxies of the molten systems leading to large extrapolations and possible errors in the interpretation of final results. This approach is used for instance in [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF] where highpressure data (up to 206 GPa) suggest a minuscule Fe isotope fractionation between metal and silicate that is one order of magnitude lower than the one found by [START_REF] Shahar | Sulfur-controlled iron isotope fractionation experiments of core formation in planetary bodies[END_REF] at 1-2 GPa.

In recent years, computational techniques have become useful tools to estimate isotope fractionation factors in mineral systems. This is particularly true for methods based on quantum mechanics that allow to describe the vibrational properties of light and heavy isotopes at any pressure and temperature conditions [START_REF] Blanchard | Equilibrium fractionation of non-traditional isotopes: A molecular modelling perspective[END_REF]. However, the weakness of these methods typically lies in the use of the quasi-harmonic approximation that falls when temperature increases or when the systems become fully molten. In addition, in dynamical systems such as liquids and melts, it has been shown that configurational disorder needs to be taken into account in order to obtain meaningful equilibrium fractionation factors [START_REF] Blanchard | Equilibrium fractionation of non-traditional isotopes: A molecular modelling perspective[END_REF][START_REF] Pinilla | Equilibrium magnesium isotope fractionation between aqueous Mg2+ and carbonate minerals: insights from path integral moeclular dynamics[END_REF].

In this work, we use ab-initio computational methods to perform a systematic study of the 56 Fe/ 54 Fe equilibrium isotope fractionation in molten and solid Fe1-xSx alloys at conditions of the core formation.

By comparing the fractionation factors obtained from solid and molten systems, we estimate the validity of the experimental approximation consisting in using solid metals as a proxy for molten alloys. Additionally, we comment on the effects of S concentration on Fe isotope fractionation in liquid systems. We discuss our findings on view of latest results on equilibrium Fe isotope fractionation and their relevance for the formation of the Earth's core.

METHODS

Equilibrium isotope fractionation factor from harmonic vibrational modes

The equilibrium isotope fractionation factor of an element displaying two isotopic forms Y and Y* between two phases a and b is related to the ratio of the isotope concentration ratios:

𝛼(𝑎, 𝑏, 𝑌) = ( 𝑛 𝑌 * 𝑛 𝑌 ) 𝑎 ( 𝑛 𝑌 * 𝑛 𝑌 ) 𝑏 ,
where nY is the mole fraction of isotopes Y in phase a. The equilibrium fractionation factor between two phases can be related to the reduced partition function ratio β(a,Y) of each phase by:

ln(α(a,b,Y))=ln(β(a,Y)) -ln(β(b,Y)).
Isotopic reduced partition function ratios are usually presented as 10 3 ln(β(a,Y)) in permil. The reduced partition function ratio β(a,Y) describes the isotopic fractionation properties of a given phase normalized to the properties of a classically behaving system. It can be shown that within the harmonic approximation and using the Teller-Redlich rule for frequencies [START_REF] Wilson | Molecular vibrations[END_REF], the β(a,Y) factor in a solid phase can be written as [START_REF] Meheut | Equilibrium isotope fractionation in the kaolinite, quartz, water system: prediction from first-principles density functional theory[END_REF]:

𝛽(𝑎, 𝑌) = [∏ ∏ 𝜐 𝑞,𝑖 * 𝜐 𝑞,𝑖 𝑞 3𝑁 𝑎𝑡 𝑖=1
𝑒 -ℎ𝜐 𝑞,𝑖 * /(2𝑘 𝐵 𝑇)

1 -𝑒 -ℎ𝜐 𝑞,𝑖 * /(𝑘 𝐵 𝑇)

1 -𝑒 -ℎ𝜐 𝑞,𝑖 /(𝑘 𝐵 𝑇) 𝑒 -ℎ𝜐 𝑞,𝑖 /(2𝑘 𝐵 𝑇) ]

1/𝑁𝑁 𝑞 (2)
where νq,i and ν*q,i are the phonon frequencies of the i-th branch with wave vector q for the two isotopologues. Nat is the number of atoms in the unit-cell and N is the number of sites occupied by the Y atom in the unit-cell. h and kB are the Planck and Boltzmann constants respectively, and T the temperature. The second product in equation 2 is usually performed on a grid containing Nq q-vectors to ensure full convergence of the β(a,Y) value.

Equilibrium isotope fractionation factor from the kinetic energy

To calculate the equilibrium fractionation factor directly from molecular dynamics simulations of liquid systems, we can use the equation relating the β-factor to the atomic kinetic energy (K) [START_REF] Polyakov | On anharmonic and pressure corrections to the equilibrium isotopic consntants for minerals[END_REF][START_REF] Polyakov | Equilibrium iron isotope fractionation at core-mantle boundary conditions[END_REF]:

𝑙𝑛𝛽(𝑎, 𝑌) = ∆𝑚 𝑚 * ( 𝐾 𝑘 𝐵 𝑇 - 3 2 ) , (3)
where Δm=m*-m with m* the mass of the heavier isotope. In this work we have used the definition of kinetic energy (K) as a function of the partial vibrational density of states (pDOS) that can be computed directly from computer simulations or experiments using nuclear resonant inelastic X-ray scattering spectroscopy (NRIXS) [START_REF] Dauphas | Sciphon: a data analysis software for nuclear resonant inelastic X-ray scattering with applications to Fe, Kr, Sn[END_REF][START_REF] Dauphas | A general moment NRIKX approach to the determination of equilibrium Fe isotopic fractionation factors: Application to goethite and jarosite[END_REF][START_REF] Polyakov | Equilibrium iron isotope fractionation at core-mantle boundary conditions[END_REF]. In this way, the kinetic energy of an atomic specie in the phase of interest can be defined as a function of the vibrational density of states of this atomic specie, g(ν), as [START_REF] Vos | Determining the band gap and mean kinetic energy of atoms from reflection electron energy loss spectra[END_REF]:

𝐾 = 3 2 ∫ 𝑔(𝜈)𝛼(𝜈, 𝑇)𝑑𝜈 𝜈 𝑓 𝜈 𝑖 ∫ 𝑔(𝜈)𝑑𝜈 𝜈 𝑓 𝜈 𝑖 , ( 4
)
with α(ν,T) the kinetic energy of a quantum harmonic oscillator. This last term is defined as:

𝛼(𝜈, 𝑇) = ℎ𝜈 ( 1 𝑒 ℎ𝜈/𝑘 𝐵 𝑇 -1 + 1 2 ). (5)
The vibrational density of states of a system of atoms can be retrieved from the Fourier transform of the velocity auto-correlation function (VCF) thermodynamically averaged over all atoms of a given specie. The pDOS for a system of N atoms is given by

𝑔(𝜐) = ∫ 𝑑𝑡 〈𝒗(𝑡) • 𝒗(0)〉 〈𝒗(0) • 𝒗(0)〉 ∞ -∞ 𝑒 𝑖2𝜋𝜈𝑡 (6)
The VCF has been calculated directly from ab-initio molecular dynamics simulations at the temperature and pressure of interest.

Computational details

The calculations in this work are based on density functional theory as implemented in the VASP code developed by [START_REF] Kresse | Ab initio molecular dynamics simulations of the liquid metal-amorphous semiconductor transition in germanium[END_REF] and updated by [START_REF] Kresse | Efficient iterative schemes for ab-initio total energy calculations using a plane-wave basis set[END_REF]. The exchange and correlation energy was represented by the generalized-gradient-approximation GGA approach in the version for solids, PBEsol, developed by Von [START_REF] Barth | Basic Density Functional Theory: An overview[END_REF]. Electronic wave functions were expanded on a plane-wave basis with a maximum cutoff energy of 600 eV, which guarantees convergence of all calculated properties. Pseudopotentials of the projector augmented wave type [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] were used for Fe and S, with an electronic configuration of valence states 3s 2 3p 6 3d 6 4s 2 and 3s 2 3p 4 , respectively.

Ab-initio molecular dynamics (AIMD) calculations were performed using a Nosé-Poincaré thermostat [START_REF] Bond | The Nosé-Poincaré tmethod for constant temperature molecular dynamics[END_REF] and a Parrinello-Raman barostat [START_REF] Hernandez | Metrc-tensor flexible-cell algorithm for isothermal-isobaric molecular dynamics simulations[END_REF] to work at temperatures and pressures of interest. All systems were thermodynamically equilibrated at core formation conditions (60 GPa, 3000 K) by means of a NPT ensemble for a minimum time of 2 ps, which was found sufficient to guarantee fully thermodynamic equilibrium. Once the system was equilibrated to the desired pressure and temperature, calculations were subsequently performed in a NVT ensemble at the same temperature, using a Nosé-Hoover thermostat [START_REF] Hoover | Canonical dynamics: equilibrium phase-space methods[END_REF][START_REF] Nosé | a unified formulation of the constant temperature molecular-dynamics methods[END_REF]. All AIMD calculations were carried out with an integration time of 0.5 fs.

In order to understand the effect of S on the Fe isotope fractionation, we selected the following systems: a pure solid hcp iron structure (54 atoms) and its molten counterpart (0 at%S), a solid I4 Fe3S structure (64 atoms) and its molten counterpart (25 at%S) and two liquid alloys with intermediate compositions of 9.4 at%S and 17.2 at%S (64 atoms each). All calculations besides the one for pure liquid iron have been carried out at 60 GPa and 3000 K. In the case of pure liquid iron we have carried out calculations at a higher temperature of 6000 K due to the fact that at this pressure of 60

GPa and 3000 K, the most stable phase of iron is still the solid hcp structure and not a molten system, as shown by previous phase diagram calculations [START_REF] Alfé | Composition and temperature of the Earth's core constrained by combining ab-initio calculations and seismic data[END_REF]. In order to understand the effect of temperature on force constants and isotope fractionation factors, we have run an extra AIMD calculation of solid hcp Fe at 300 K.

Additional to the AIMD calculations mentioned above, we have also carried out zero-temperature structural optimizations of the two crystalline systems, to estimate the 56/54 Fe equilibrium isotope fractionation factor from vibrational properties using the harmonic approximation through equation ( 2). We estimated vibrational frequencies by using first the perturbation theory method as implemented in VASP to obtain the dynamical matrix. Then, force constants and vibrational frequencies have been post-processed using the PHONOPY code [START_REF] Togo | First principles phonon calculations in material sciences Scr[END_REF]. Phonon frequencies are given by the eigenvalues of the dynamical matrix. The k-point for total energy convergence as well as the q-mesh used to obtain convergence of vibrational properties used in this work can be found in Table 1. Using these parameters, frequencies and -factors are converged to 0.001 cm -1 and 0.0001 ‰ respectively. Additionally, we have produced solid solutions at similar concentration as the two liquid alloys discussed above using the solid solution method reported by D'Arco et al. ( 2013); [START_REF] Mustapha | On the use of symmetry in configurational analysis for the simulation of disordered solids[END_REF]. The structures of these solid solutions were found to be dynamically non-stable and therefore their isotope fractionation factors are not reported in this work.

It should also be noted that we have found magnetic effects to be non-relevant at the studied pressure and temperature conditions and therefore have not been taken into account in this work (See Valencia et al. (2020)).

Keeping in mind these methodological part, we can now outline the four different methods used in this work to estimate -factors and explain the limitations and advantages of each one of them. Firstly, the trajectories obtained from finite-temperature AIMD calculations have been used to obtain the velocity auto-correlation function (VCF) and subsequently, the partial density of states (pDOS) using equation ( 6). One limitation of the pDOS method lies on the length of the trajectory used to obtain meaningful results given that the noise in the velocity auto-correlation function increases with the length of the time interval and could lead to non-converged pDOS spectra. In this work, the reported pDOS are obtained as an average of several 6 ps trajectories (at least 5 trajectories) obtained over a longer MD trajectory. We have estimated the kinetic energy and subsequently the -factor using equations 4 and using the following equation (e.g. [START_REF] Dauphas | A general moment NRIKX approach to the determination of equilibrium Fe isotopic fractionation factors: Application to goethite and jarosite[END_REF]:

〈𝐹〉 = 𝑚 ℏ 2 ∫ 𝑔(𝜐) +∞ 0 𝜐 2 𝑑𝜐 (7)
This method has the advantage to include anharmonic effects as well as configurational disorder, of relevance in the case of liquid systems. We will refer to this method as the "pDOS-melt" results.

Furthermore, we have used the same pDOS methodology applied to Fe-hcp and Fe3S-I4 crystalline systems with AIMD calculations performed at 300 K and 60 GPa. The -factors obtained through this method will be presented as the "pDOS-solid" results.

Using the snapshots from the same AIMD trajectories of the molten systems, we have performed zerotemperature constant-pressure optimisations of the atomic positions in order to minimise the total energy. Following this, we have used the relaxed configurations to calculate vibrational properties and estimate the -factor using equation 2. This method is similar to the one presented in previous works by [START_REF] Pinilla | Equilibrium magnesium isotope fractionation between aqueous Mg2+ and carbonate minerals: insights from path integral moeclular dynamics[END_REF] and results obtained from this method will be referred as the "RELAX" results.

In the "RELAX-melt" method the obtained snapshots can be seen as those resembling a glassy system obtained through rapid quenching and characterised by structural disorder but where bond distance distribution has been evened out due to optimization of atomic positions. Finally, in the case of crystalline Fe-hcp and Fe3S-I4, we have calculated the vibrational properties within the harmonic approximation from the optimised structure at zero temperature and used equation 2 to obtain thefactor. Results from these calculations will be known as "RELAX-solid" results. It should be noted that in the case of solid structures, the RELAX-melt method (i.e., AIMD at finite temperature followed by zero-temperature relaxations) converges to the same results as the RELAX-solid method (i.e., direct zero-temperature structural relaxation). A summary of the methods used in the present study can be found in Table 2.

RESULTS AND DISCUSSION

Iron and sulfur kinetic energies

To determine the kinetic energy for Fe and S atoms we calculated the velocity auto-correlation function (VCF) for all the studied systems as shown in Figure 1. In the case of liquid alloys, the typical behaviour for dense liquids was observed, with an intermediate negative region produced by atomic collision and a subsequent rise before oscillating around a zero correlation [START_REF] Hansen | Theory of simple liquids[END_REF]. The fact that the positive part of the VCF is larger than the negative one describes a system with net diffusion, opposite to what is seen in the case of solids. The corresponding pDOS are also displayed in Figure 1. It is interesting to note that in the case of solids the pDOS has a zero value at zero frequency. At low frequency (long times) the pDOS of solids should follow a quadratic trend just like the one described by the Debye model for the long wavelength limit [START_REF] Hansen | Theory of simple liquids[END_REF] and corresponding to acoustic vibrational modes. This is different from what is seen in the case of liquids where diffusive effects lead to a net pDOS value at zero frequency. Additionally, in the case of solids we have also displayed pDOS from zero-temperature vibrational frequencies to show the effect of temperature producing a broader spectrum.

Using the pDOS we have determined the kinetic energy via equation 4. Figure 2 shows the iron kinetic energy as a function of temperature for some of the studied systems. We note the large error bars reported in these results, where kinetic energies for liquid and solid system are close. However, we observe in details, systematic trends. Fe atoms have a slightly smaller kinetic energy than S atoms with a KS -KFe = 0.8 meV for solid Fe3S at 300 K and 60 GPa. An Fe atom in pure solid Fe hcp lattice has an average kinetic energy KFe-solid = 43.4 meV at 300 K whereas a similar atom in a Fe3S structure has a KFe3S-solid = 44.2 meV. This difference complies with the fact that the Fe-S bond has a more covalent nature thus producing a stronger bond than a pure metallic Fe-Fe bond. This conclusion could also be extended to liquids but with a smaller magnitude. Similar behaviour has been also reported by [START_REF] Morard | In situ structural investigation of Fe-S-Si immiscible liquid system and evolution of Fe-S bond properties with pressure[END_REF] in the case of Fe-S-Si liquid alloys. Overall, the Fe kinetic energy of Fe1-xSx alloys tends to be slightly higher than those based on pure Fe. Furthermore, the kinetic energy of an Fe atom in a solid, liquid and gas phase of these alloys is different due to the bond nature, but the trends are not the same depending on the composition investigated. All systems possess a larger kinetic energy than that for a free Fe atom (3kBT/2) = 38.7 meV/atom at 300 K, but in the case of pure Fe, the solid has a lower kinetic energy than that found for the liquid phase, whereas, this effect is reverted in the case of the Fe3S system. Finally, we can observe that the kinetic energy difference decreases as temperature increases, leading to a constant value due to the increase of diffusive-like contributions.

Iron β-factors from solid and molten alloys

The 56 Fe/ 54 Fe reduced partition function ratios (β-factors) obtained from all methods described in Table 2 are shown in Figure 3. A further insight into our results can be obtained from the iron force constants displayed in Table 3. We need to differentiate the results estimated from calculations at zero temperature (RELAX-melt and RELAX-solid methods) from those at finite temperature (pDOS-melt and pDOS-solid methods). At 3000 K, -factors from the later methods seem to be 0.02 ‰ lower than the -factors from the former one (Fig. 3). Consistently, force constants from pDOS calculations (pDOS-melt and pDOS-solid) are ~ 100 N/m smaller than those from frequency methods (RELAX-melt and RELAX-solid) (Table3). This difference is found to be related to shorter Fe-Fe and Fe-S bonds at zero temperature as well as to the truncation in the kinetic energy integral to a certain maximum cutoff value of 246 meV (See SI). This truncation procedure is much on the lines to what is done for NRIXS experiments [START_REF] Dauphas | A general moment NRIKX approach to the determination of equilibrium Fe isotopic fractionation factors: Application to goethite and jarosite[END_REF].

If we now compare solid, glass and liquid results, we can make several observations. For both Fe and Fe3S compositions, RELAX-melt and RELAX-solid methods present comparable results for the -factors (Fig. 3) or the force constants (Table 3). The RELAX-melt method includes some configurational disorder like in a glass, even if the distribution of bond distances has been diminished by the optimization of atomic positions. Therefore, these results suggest that the isotopic properties of a glass system, as that given by the RELAX-melt method, may be equivalent to that of a perfect solid structure.

When dealing with liquid systems taking full account of thermal effects such as in the case of the pDOS calculations, we see that the -factors extrapolated at conditions of core formation (3000 K, 60 GPa) from a liquid Fe simulation run at 6000 K and from a solid simulation run at 300 K are comparable with β Fe solid-β Fe liquid ≈ 0.02 ‰ at 3000 K. The same conclusion can be drawn for the Fe3S composition.

Overall, our results highlight that although configurational disorder is different in all systems, bond type and local environment are the controlling parameters of this modelling. This shows that β-factors obtained from solid systems may in principle be used in place of molten ones, with differences not larger than 0.02 ‰ at conditions of core formation. Hence, we can conclude that crystalline or glassy iron alloys are good proxies for molten systems for deep Earth studies as has already been done in several experimental works.

Finally, we can compare the values of iron force constants calculated for pure Fe-hcp and Fe3S using the RELAX-melt or RELAX-solid methods to the equivalent NRIXS-derived room-temperature values reported in previous experimental works by [START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF] and [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF]. Calculated and experimental values are all in good agreement with the largest difference between these values of about 10% for the pure Fe system (Table 3). It is noteworthy that there are no data for solid alloys at S concentrations of 5.6 and 10.6 wt% because these solid alloys were found to be dynamically unstable under the pressure used.

Thermal extrapolation of the -factor

The temperature dependence of the -factor for zero-temperature calculations (RELAX-solid and RELAX-melt methods) and finite temperature MD simulations (pDOS-solid and pDOS-melt methods) has been determined using equations 2 and 3 respectively. The procedure of using vibrational frequencies or the partial density of states obtained at a given temperature to extrapolate the functional behaviour of -factors to other temperatures is a widely used approach both in theory and experimental analysis, but can be prone to errors. Here we assess the validity of such procedure by comparing the -factors obtained using the pDOS method from AIMD simulations performed at several temperatures (3000 K, 4000 K, 5000 K, 6000 K and 7000 K) with the values obtained using an extrapolation procedure (Fig. 4). Our results show that the data extrapolation from a simulation at 3000 K toward higher temperatures is perfectly valid. On the other hand, the data extrapolation towards 3000 K from simulations performed at higher temperature (4000-7000 K) lead to small uncertainties. Extrapolated values oscillate from one simulation temperature to the other but globally these oscillations do not exceed the error bar already calculated for the pDOS-melt method (about ±0.01 ‰ at 3000 K, Fig. 4). More critical is the linear extrapolation up to 3000 K of the pDOS-solid simulation run at 300 K. In order to assess the validity of this extrapolation, we ran an overheated Fehcp structure at 60 GPa and 3000 K. This simulation leads to a -factor close to the melt one at the same conditions, suggesting that this extrapolation overestimates the force constant by about 30 N/m (Fig. 4b) and the -factor by about 0.01 ‰ (Fig. 4a). Data for pDOS-solid displayed in Figure 6 as well as Tables 3 and4 have been corrected to account for this overestimation.

Overall, our results suggest that the force constant determined by spectroscopic analysis at room temperature will be overestimated by 30 N/m when discussing implications at conditions near the 3000 K. This overestimation is of the same order of magnitude as error bars reported in spectroscopic methods [START_REF] Dauphas | A general moment NRIKX approach to the determination of equilibrium Fe isotopic fractionation factors: Application to goethite and jarosite[END_REF][START_REF] Roskosz | Spinel-olivinepyroxene equilibrium iron isotopic fractionation and applications to natural peridotites[END_REF]. Therefore, these results warn on the use of solids as proxies for high-temperature liquids and make this approximation acceptable for a range of temperatures between 1000-2000 K. It is important to say that the behaviour observed here for the force constant as a function of temperature is true at the high-pressure regime where this work has been carried out. However, we expect a more significant correction at low pressures where vibrational effects, thermal expansion and even magnetic effects might be more pronounced.

Correcting for fluidicity effects on estimated -factors

The pDOS methodology described above has been widely used for the NRIXS-based and static computation study of isotope fractionation properties in crystalline systems where the harmonic approximation is most of the time a good assumption [START_REF] Bourdon | Isotope tracers of core formation[END_REF][START_REF] Dauphas | Sciphon: a data analysis software for nuclear resonant inelastic X-ray scattering with applications to Fe, Kr, Sn[END_REF].

However, a direct extension to liquids may not be appropriate since entropy and fluidicity effects are dominated by low-frequency modes. In particular, the pDOS of liquids has a non-zero value at zero frequency, which corresponds to the diffusion mode (Fig 1).

To assess the effect of fluidicity on the calculated -factors we have used the two phase approach proposed by [START_REF] Lin | The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids[END_REF] and [START_REF] Lin | Thw-phase thermodynamic model for efficient and accurate absolute entropy of water from molecular dynamics simulations[END_REF] in which the partial density of states (𝑔(𝜈)) of a system with 3N degrees of freedom can be partitioned into a gas-and a solid-like component:

𝑔(𝜈) = 𝑔 𝑔 (𝜈) + 𝑔 𝑆 (𝜈), (8)
where the gas-like diffusive component 𝑔 𝑔 (𝜈) corresponds to 3N g = 3fN degrees of freedom with f being the gas fraction. The remainder, 𝑔 𝑆 (𝜈), describes a non-diffuse system in which 𝑔 𝑆 (0) = 0 (i.e.

no diffusion). Thus, the degrees of freedom from the system that behaves in a solid-like way are 3N S = 3N(1-f).

In order to estimate the solid-like density of states, it is necessary to define the gas-like diffusive component 𝑔 𝑔 (𝜈) and the gas fraction f. In their original paper, [START_REF] Lin | The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids[END_REF] approximate 𝑔 𝑔 (𝜈) to that of a hard sphere fluid defined as:

𝑔 𝑔 (𝜈) = 𝑠 0 1+[ 𝜋𝑠 0 𝜈 6𝑓𝑁 ] 2 , (9)
where 𝑠 0 is the original 𝑔(𝜈) zero-frequency value. Likewise, the gas fraction, f, can be written as

2Δ -9/2 𝑓 15/2 -6Δ -3 𝑓 5 -Δ -3 2 𝑓 7 2 + 6Δ -3 2 𝑓 5 2 + 2f -2 = 0, (10)
where

Δ(T, ρ, m, 𝑠 0 ) = 2𝑠 0 9𝑁 ( 𝜋𝑘 𝐵 𝑇 𝑚 ) 1/2 𝜌 1/3 ( 6 𝜋 ) 2/3 (11)
with 𝑚 the mass of the particles and 𝜌 = N/V, the density. is a function of the studied alloy. Note that through MD calculations, the obtained value s0 determines and subsequently 𝑔 𝑔 (𝜈) and the gas fraction also called fluidicity factor f. In this way, any thermodynamic property can be estimated as a sum of gas and solid contributions.

The respective decompositions at conditions of the core formation can be seen in Figure 5. Results

show that the pDOS for a gas decays rapidly and monotonically as frequency increases. This behaviour is expected since in gases usually the mean free path is larger than the particle diameter. In the case of a solid-like system, and as shown earlier, the velocity auto-correlation function oscillates around zero with amplitudes decreasing as time increases. In a single particle system oscillating under a harmonic potential, this oscillation can be fitted to a cosine function and the 𝑔 𝑆 (𝜈) resembles a Dirac delta. In the case of a solid the 𝑔 𝑆 (𝜈) has several peaks reflecting the ordered nature of crystal structures.

Furthermore, the values of the fluidicity factor f increases as a function of temperature showing how the diffusive character of the alloy increases with temperature (Fig. 5c). We have used these results to estimate the -factor corrected for fluidicity effects in the case of Fe and Fe3S melts. Results are reported in figure 6. The results with or without fluidicity correction are similar with -factors differing by less than 0.01 ‰ among them. Therefore, fluidicity effects present in the pDOS are not relevant when estimating the equilibrium fractionation factor of Fe in liquid Fex-1Sx alloys.

Effect of S concentration on Fe β-factor

In order to evaluate the effect of S on the iron β-factor we have performed calculations on four systems with S concentrations ranging from 0.0 to 25 at% (Fig. 6). As already mentioned in part 3.2, results from zero-temperature methods (RELAX-melt and RELAX-solid) must be discussed separately from pDOS-based methods. Overall, all methods show a small variation of -factors (i.e. not larger than 0.02 ‰) over the whole S concentration range. However, calculations show subtle trends where the Fefactor changes continuously with S concentration (see green and black lines in Fig. 6). In addition, the change seems to be different for the case of crystalline and non-crystalline systems (Figs. 3 and6). In the case of crystalline systems (RELAX-solid and pDOS-solid results) we see the -factor for pure Fehcp to be larger than those of Fe3S-I4. The same trend is observed for the force constants (Table 3).

Conversely, pDOS-melt and RELAX-melt results show the opposite behaviour where -factors and force constants increase with S content. However, we note here that experimental uncertainties likely exceed these differences calculated here for P-T conditions relevant to the Earth's core formation.

The fact that -factors and force constants depend on concentration in the case of Fe1-xSx alloys can be understood from the local structural and chemical environment described through the coordination number and bond lengths [START_REF] Smith | The relationship of force constant and bond length[END_REF]. Therefore, we have derived the partial radial distribution functions and coordination numbers as a function of S content for solid and liquid systems (Fig. 7). In the case of molten alloys, we see that sulphur shows some covalent bonding with iron in the first shell, creating shorter Fe-S bonds and increasing slightly the length of the remaining Fe-Fe bonds. Indeed, Figure 7c clearly shows that the Fe-Fe coordination number (CN Fe-Fe ) goes from 13.4 at 0 at%S to 9.8 at 25 at%S, showing a decrease in the number of Fe atoms surrounding another Fe atom. Likewise, the Fe-S coordination number increases from 0.0 in absence of sulphur to 3.1 at 25 at%S. By adding the Fe-Fe and Fe-S coordination numbers together, we see that the total number of atoms surrounding a Fe atom in molten systems is fairly constant through all the S concentration range, suggesting that S atoms substitute readily for Fe atoms within the melt. Therefore, the nearly constant total coordination number and the increase in the number of shorter Fe-S bonds with sulphur concentration lead to an increase of -factors and force constants for all molten alloys. In contrast, in crystalline systems, atoms are more localized around their equilibrium position and as a result, bond length distributions tend to be narrower than in melts, as seen in figure 7. Additionally, Fe3S-I4 displays an average Fe-Fe distance significantly larger (+6%) than in Fe-hcp (Fig. 7) while the total coordination number remains nearly the same. This observation explains why crystalline Fe3S presents a smaller -factor and force constant than Fe-hcp.

A closer look at Fe3S-I4 shows that this structure has three distinct Fe crystallographic sites: FeI, FeII and FeIII whose -factors are 0.094, 0.103, 0.117 ‰ at 60 GPa and 3000 K respectively. The first coordination shell of the FeI site displays 12 Fe atoms (average distance dFeI-Fe=2.44 Å) and 2 S atoms (average distance dFeI-S=2.08 Å). FeII site is surrounded by 10 Fe atoms (dFeII-Fe=2.40 Å) and 4 S atoms (dFeII-S=2.14 Å) while FeIII site is surrounded by 10 Fe atoms (dFeIII-Fe=2.38 Å) and 3 S atoms (dFeIII-S=2.13 Å). Iron in site III has the highest -factor because its total coordination number is the lowest (13) and is dominated by iron with the shortest Fe-Fe bonds. Between sites I and II, displaying the same total coordination number of 14, iron in site II has the highest -factor in relation to a larger number of S neighbours and a shorter average bond length.

Results thus suggest that the small changes seen in the force constants and -factors when S concentration increases in Fe1-xSx alloys are mainly due to a combination of attenuated structural changes related to atomic coordination and bond stiffness. The geochemical implications and possible reappraisals of existing experimental data related to these small variations are discussed below.

Implications in the context of terrestrial core formation

Our theoretical calculations alleviate a number of controversies and speculations in the field of HP-HT behaviour of Fe isotopes in solid and liquid alloys. To a certain extent, these conclusions likely hold for other isotopic systems.

First of all, it is commonly proposed that solids and glasses could be considered as representative of molten counterparts as a first approximation [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF][START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF][START_REF] Yang | Iron isotopic fractionation in mineral phases from Earth's lower mantle: Did terrestrial magma ocean crystallization fractionate iron isotopes?[END_REF].

This was obviously proposed for lack of better alternatives. For pure iron, the simple close-packed structure of liquids, the similarities between density and compressibility relative to its solid phase under core conditions made this assumption reasonable [START_REF] Laio | Physics of Iron at Earth's Core Conditions[END_REF]. Similar assumptions were also made for silicates for which structural data also indicates strong similarities between glasses and melts at the local scale [START_REF] Dauphas | Magma redox and structural controls on iron isotope variations in Earth's mantle and crust[END_REF][START_REF] Mysen | Silicate Glasses and Melts[END_REF][START_REF] Sanloup | Structural change in molten basalt at deep mantle conditions[END_REF]. Our results demonstrate, in the case of pure iron and iron sulphides that this assumption holds at least in the high pressure and temperature conditions prevailing during the Earth's core formation.

Second, another common assumption made to extrapolate experimental spectroscopic results to planetary conditions [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF][START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF][START_REF] Yang | Iron isotopic fractionation in mineral phases from Earth's lower mantle: Did terrestrial magma ocean crystallization fractionate iron isotopes?[END_REF] or, conversely to try to explain flagrant inconsistencies between spectroscopic and experimental results [START_REF] Shahar | An assessment of iron isotope fractionation during core formation[END_REF] focusses on the validity of harmonic approximations. Anharmonicity must indeed be present at some level as it is the mechanism responsible for thermal expansion, which is present in geological materials. For example, the volume of forsterite expands by 8% when the temperature is increased from 300 to 2300 K [START_REF] Bouhifd | Thermal expansion of forsterite up to the melting point[END_REF] and no expansion would be present if the bonds were perfectly harmonic. If the bonds were perfectly harmonic, then the force constant of those bonds should not vary with temperature. Conversely, anharmonicity in interatomic potentials would be manifested as an apparent variation of the force constant of the iron bonds as temperature changes.

Our results show such a variation (Fig. 4b). The Fe force constant decreases from 288 N/m at 300 K to 262 N/m at 3000 K for solid iron, and the force constant for liquid iron decreases from about 255 N/m at 3000 K to 238 N/m at 7000 K. However, these variations are of the same order of magnitude as error bars reported in spectroscopic studies. Consistent results were recently produced in the case of silicate minerals and glasses (Roskosz et al., in prep.). This is an important result because it strengthens the use of spectroscopic approach to derive isotopic properties of geomaterials at extreme conditions. It also highlights one more time the unexplained inconsistency between fractionation factors derived from experimental petrological experiments (e.g. [START_REF] Shahar | Sulfur-controlled iron isotope fractionation experiments of core formation in planetary bodies[END_REF] and theoretical and spectroscopic determinations [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF][START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF][START_REF] Yang | Iron isotopic fractionation in mineral phases from Earth's lower mantle: Did terrestrial magma ocean crystallization fractionate iron isotopes?[END_REF]. To date this discrepancy remains unexplained but is clearly beyond the effects related to the anharmonic nature of metallic materials.

Third, these data can be used to describe Fe isotopes distribution within core forming materials during its crystallization. At the studied conditions (60 GPa and 3000 K), our calculations suggest that no measurable fractionation should be produced between liquids and solids even if incongruent crystallisation of a Fe-S molten alloy is considered. The limited amount of light alloying element and Ni should not change this conclusion since NRIXS measurements performed on solid analogues confirmed that reasonable amounts of these elements could not significantly modify the force constants of the alloy [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF][START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF]. Moving to milder conditions, experiments were recently performed to mimic core formation in planetesimals where pressure and temperature do not exceed a few GPa and 1600 K. Clearly, extrapolations of our calculations to lower temperature and pressure is limited by phase transitions and no firm conclusion can be drawn in the context of planetesimal differentiation. Nonetheless, as a first attempt, we have compared our 56/54 solid- liquid=10 3 (ln solid-ln liquid) results with those recently reported by [START_REF] Ni | Heavy iron isotope composition of iron meteorites explained by core cristallization[END_REF] concerning planetesimal core crystallization (Table 4). We find our results to be about five times smaller than their reported average value of 56/54 solid-liquid = 0.09 0.04 ‰ at 1573 K for a similar range of S concentration and much lower pressure. If our predicted very small fractionation is confirmed by future calculations performed at lower pressure and temperature, then the relatively heavy iron isotopic signature of iron meteorites could not be explained solely by core crystallization and may also derive from fractionation associated to processes involved during the exhumation of core materials.

CONCLUSIONS

In this work, we have used computational methods to understand the effect of S on Fe equilibrium isotope fractionation in Fe1-xSx alloys as well as the consequences of using solid metal as proxies of molten alloys when performing experimental studies. We have shown that although configurational disorder could be different between these systems, bond type and local environment are the controlling parameters of Fe isotope fractionation and that β-factors obtained from solid systems may in principle be used in place of molten ones, with differences not larger than 0.01 ‰ at the high pressure and temperature conditions investigated here. In addition, while the Fe force constant is found to decrease with increasing temperature, this consequence of the anharmonicity of the interatomic bonds remains smaller than the typical error bars of spectroscopically-derived force constants at high pressure. This result makes valid the extrapolation in temperature of such experimental force constants. Regarding the S concentration effect, we observe changes in Fe -factor and force constant much smaller than the estimated error bars even if the subtle and continuous trends with S concentration can be rationalized by a combination of structural changes related to atomic coordination and bond stiffness. Our calculations have thus suggested that no measurable fractionation should be produced between liquid and solid even if incongruent crystallisation of a Fe-S molten alloy is considered at conditions of the Earth's core formation. However, it is important to note that the results outlined here for Fe1-xSx alloys apply to high pressure and high temperature conditions and that the effect of anharmonicity, S concentration and possibly magnetization on Fe isotopes behaviour at milder conditions could be larger than what has been found in this work. Figure 2. Iron kinetic energy difference as a function of temperature for the studied solid and liquid systems at 60 GPa. The difference has been taken with respect to the kinetic energy per atom of an ideal gas (3kBT/2). The inset is a closeup of the kinetic energy for the region between 3000K and 6000K.

Figure 3. 56/54 Fe -factor in pure Fe and Fe3S liquid and solid systems estimated using the RELAX-solid (blue circles) and RELAX-melt (green diamonds) methods, and the integration of the pDOS for a liquid (black squares) and solid (red triangles) system. Dashed lines correspond to results from NRIXS experiments reported by [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF]. Colour-circles show values from independent AIMD simulations performed at 3000, 4000, 5000, 6000 and 7000 K. Solid lines correspond to temperature dependences estimated from a single simulation performed at the displayed temperatures, and then extrapolated to other temperatures. Violet dashed line corresponds to results from the pDOS-solid method for Fe-hcp run at 300 K and extrapolated at higher temperatures. b) Fe force constant as a function of temperature in liquid iron (pDOS-melt method). Violet-dotted line and bluedashed line correspond to values obtained from pDOS-solid at 300 K and 3000 K respectively.

Figure 5. Partial density of states for a) pure iron at 6000 K and b) Fe3S at 3000 K, both at 60 GPa. Black-solid line represents total pDOS, dashed and dash-dotted lines are the solid and gas contributions respectively. c) fluidicity factor f as a function of temperature for the pure iron system.

Figure 6. 56/54 Fe -factor as a function of S concentration at 60 GPa and 3000 K using different calculation methods. In addition to the four standard methods summarized in Table 2, "pDOS decomp." corresponds to the pDOS-melt result corrected for fluidicity effects (cf. section 3.5). Experimental results from [START_REF] Shahar | Pressure-dependent isotopic composition of iron alloys[END_REF] and [START_REF] Liu | Iron isotopic fractionation between silicate mantle and metallic core at high pressure[END_REF] are shown for comparison. The small panel on the right hand side corresponds to a zoom up of results obtained at 0 and 25 at% and is displayed to facilitate the reading of the reported values.

Crystalline system

Supercell size k-points mesh q-points-mesh Fe -hcp 4x4x4 3x3x3 4x4x4

Fe3S -I4 2x1x2 1x2x1 6x6x6

Table 1. k and q meshes used for the zero-temperature calculation of electronic and vibrational properties of crystalline systems.

Method

Figure captions Figure 1 .

 captions1 Figure captions

Figure 4 .

 4 Figure4. a) Temperature dependence of the Fe -factor estimated using the pDOS-melt method. Colour-circles show values from independent AIMD simulations performed at 3000, 4000, 5000, 6000 and 7000 K. Solid lines correspond to temperature dependences estimated from a single simulation performed at the displayed temperatures, and then extrapolated to other temperatures. Violet dashed line corresponds to results from the pDOS-solid method for Fe-hcp run at 300 K and extrapolated at higher temperatures. b) Fe force constant as a function of temperature in liquid iron (pDOS-melt method). Violet-dotted line and bluedashed line correspond to values obtained from pDOS-solid at 300 K and 3000 K respectively.

Figure 7 .

 7 Figure 7. Structural properties for molten Fe1-xSx alloys. a) and b) panels correspond to Fe-Fe and Fe-S radial distribution functions respectively for various S concentrations. Dotted and Dashed-dotted lines are values for crystalline Fe-hcp and Fe3S-I4 at 300K. c) Fe-Fe (black dots)

  Figure 2Click here to access/download;Figure;fig2.pdf

  Figure 3aClick here to access/download;Figure;fig3a.pdf

  Figure 3bClick here to access/download;Figure;fig3b.pdf

Figure 6

 6 Figure 6Click here to access/download;Figure;fig6.pdf

Figure 7c

 7c Figure 7cClick here to access/download;Figure;fig7c.pdf

ACKNOWLEDGEMENTS

This work was performed using Granado-HPC from the Universidad del Norte, SCARF from the STFC of the UK, and the HPC resources from CALMIP (Grant 2020 -P1037). The authors acknowledge funding from COLCIENCIAS (No. 2015-710-51568; Contract No. 023-2016) and ECOSNORD (C17U01, No. FP44842-143-2017) through research grants.

account.

Frequencies estimated using the harmonic approximation RELAX-solid Optimization of Fe-hcp and Fe3S-I4 crystal structures at T = 0 K. Frequencies estimated using the harmonic approximation.

Calculated using the vibrational frequencies (eq. 2)