First-principles calculations of high pressure and temperature properties of Fe7C3

Li-Li Fan(范莉莉), Xun Liu(刘勋), Chang Gao(高畅), Zhong-Li Liu(刘中利),Yan-Li Li(李艳丽), and Hai-Jun Huang(黄海军),†

1School of Science,Wuhan University of Technology,Wuhan 430070,China

2School of Materials Science and Engineering,Harbin Institute of Technology,Harbin 150001,China

Keywords: iron carbide,phase stability,thermoelastic properties,sound velocities,inner core

1.Introduction

As is well known, the Earth’s core is composed of Fe–Ni alloy with a small quantity of light elements.Based on geochemical and meteorological studies, light elements most likely originate from carbon (C), hydrogen (H), oxygen (O),sulfur (S), and silicon (Si).[1]Until now, the light elemental composition of the Earth’s core has been poorly understood.Among these light elements,C has been suggested as a major light elemental component of the Earth’s core,as C is enriched in the solar system[2]and has a high affinity for liquid iron under core–mantle differentiation conditions.[3]In addition, the shear wave velocity and Poisson’s ratio of Fe7C3are similar to those of the inner core[4,5]and Fe7C3melting experiments indicate that C might be a potential light element in the inner core.[6]

However, whether Fe7C3is stable under core condition and the nature of its relatively stable structure remain unclear.Based on previous experimental results, two potential stable structures can be reasonably hypothesized at high temperature and pressure: hexagonal structure (space groupP63mc; h-Fe7C3) and orthorhombic structure (space groupPnma,Pmn21mc, orPbca; o-Fe7C3).[5,7–9]Nakajimaet al.found that h-Fe7C3is stable up to 71.5 GPa and 1973 K by usingin situx-ray diffraction.[10]Prescheret al.suggested that o-Fe7C3is the more stable phase up to 205 GPa and 3700 K by using single-crystal and powder x-ray diffraction.[5]Theoretical calculations have been similarly inconclusive and often inconsistent.Fanget al.found that the o-Fe7C3phase is more stable than h-Fe7C3under ambient conditions.[11]Razaet al.found that o-Fe7C3(space groupPbca)is stable below 100 GPa,while h-Fe7C3becomes more stable above 150 GPa and requires higher temperatures for stability.[12]

The thermal equation of state(EOS)of Fe7C3was measured experimentally only at low pressures and temperatures.Litasovet al.[13]and Nakajimaet al.[10]measured the EOS below 30 GPa; however, their results deviate widely when extrapolated to the high temperature and pressure conditions of the inner core.Chen and Zhang acquired compressional(VP) and shear wave velocity (VS) data for Fe7C3up to 154 GPa at 300 K by using x-ray diffraction.[4]Subsequently,Prescheret al.investigated theVPandVSfor Fe7C3up to 158 GPa and 300 K by using the same method.[5]However,these sound velocity measurements of Fe7C3exhibited large discrepancies.[4,5]Therefore,it is necessary to investigate the crystal structure and thermoelastic properties of Fe7C3under core condition to ascertain the amount of carbon present in the inner core.

Here, using first-principles calculations, we estimate the free energy,EOS,and sound velocity of Fe7C3.The EOS and sound velocity data for Fe7C3are calculated at 0 K–5000 K,which is close to the inner core boundary temperature obtained by phonon spectrum calculations.We also use our findings to discuss whether carbon is the main light element present in the inner core.

2.Methods

All calculations were performed by using the Viennaab initiosimulation package (VASP) based on density functional theory (DFT).[14]The conjugate gradient (CG) algorithm was adopted for the geometric optimization of the crystal structures.Exchange–correlation interactions were considered as the Perdew–Burke–Ernzerhof (PBE) parametrization of the gradient approximation (GGA).[15]The electron–ion interactions were represented by projected augmented-wave(PAW) pseudopotentials.[16,17]The simulation cell contained 40 atoms for o-Fe7C3.Prior to formal calculation, we calculated the free energy values of 40 and 20 atoms for h-Fe7C3at some pressure points, obtaining highly similar results.To balance the computational efficiency and precision, we selected 20 atoms for h-Fe7C3calculation at all pressure points.The cutoff energy was 500 eV.The Monkhorst–Pack mesh was 6×4×7 for o-Fe7C3and 6×6×9 for h-Fe7C3.[18]The phonon spectra and free energy of o-Fe7C3and h-Fe7C3were calculated and analyzed by using the Phonopy[19]and Phasego programs.[20]The phonon spectra and Helmholtz free energy were calculated up to 360 GPa in the quasi-harmonic approximation (QHA).During the thermodynamic simulations, we constructed 1×1×2 and 1×2×2 supercells for o- and h-Fe7C3,and set denserk-point meshes to 7×5×3 and 6×6×9 grids,respectively.[21–23]

Previous studies have reported that Fe7C3undergoes a magnetic transition from a ferromagnetic (FM) to nonmagnetic (NM) phase at high pressure.[5]Therefore, we performed spin-limited and unlimited calculations to investigate the effect of spin on the EOS.

The elastic constant of Fe7C3was estimated by using the stress-strain method.[24]Force constant was calculated in real space by using density functional perturbation theory.The bulk modulus and shear modulus were calculated from the elastic tensor by using the Voight–Hill–Reuss method.[24,25]Sound velocities were obtained from the following formulas:

whereVPis the compressional wave velocity,VSis the shear wave velocity,Bis the bulk modulus,Gis the shear modulus,andρis the density.Considering the symmetry of Fe7C3,BandGcan be expressed as

whereCx(x=11,22,33,44,55,66,12,13,or 23)values are all elastic constants.

3.Results

3.1.Relative stabilities of o-and h-phase Fe7C3

To determine which structure is stable, we calculate the Gibbs free energy values of o-Fe7C3and h-Fe7C3at high pressure and temperature (Fig.1).The difference in Gibbs free energy in o-Fe7C3and h-Fe7C3is a function of temperature.At 120 GPa,the free energy value of o-Fe7C3is 14 meV/atom at 0 K and 72 meV/atom at 4000 K, less than the counterparts of h-Fe7C3.At 360 GPa, the free energy of o-Fe7C3is 17 meV/atom more than that of h-Fe7C3at 0 K, but is 29 meV/atom less than that of h-Fe7C3at 4000 K.Thus, the free energy of o-Fe7C3is less than that of h-Fe7C3at high temperatures.In this study,the energy value of self-consistent optimization convergence is set to 10−5eV, which is much smaller than the difference between the free energy values.Therefore, o-Fe7C3appears to be more stable than h-Fe7C3under inner core condition.Hereafter,we investigate the thermal properties of o-Fe7C3(instead of h-Fe7C3[26,27]) at high temperature and high pressure.

Fig.1.Free energy values of o-Fe7C3 and h-Fe7C3 under simultaneous action of high pressure and high temperature,showing temperature-dependent(a)Gibbs free energy at 120 GPa,(b)Gibbs free energy at 360 GPa,and(c)at 120 GPa temperature-dependent difference in Gibbs free energy between o-Fe7C3 and h-Fe7C3; (d)at 360 GPa temperature-dependent difference in Gibbs free energy between o-Fe7C3 and h-Fe7C3.

3.2.Thermal EOS of o-Fe7C3

The Fe7C3is magnetic at low pressures and undergoes magnetic collapse at high pressures.[5]Consequently,we calculate the EOSs for o-Fe7C3with and without spin at 0 K.Figure 2 shows a comparison between the results of this study and previous results from the static compression method at room temperature[5,10,28,29]andab initiosimulations,[13,26,30]whose detailed comparisons can be found in supplementary Table 3.Below 50 GPa, the EOSs with spin consideration are consistent with static measurements at 300 K[5,10,28,29]andab initiosimulations,[13,30]all of which are higher than the simulated results without spin consideration.[13,26]With the increase of pressure, the EOSs with spin restrictions began to converge with the EOSs without spin restrictions due to demagnetization at∼70 GPa–100 GPa, implying the presence of a spin phase transition.The transition pressure is consistent with that observed by using the static compression method.[5]The discontinuity in the EOS of Fe7C3near 53 GPa (or 70 GPa) is presumably attributed to a transition from a high-to low-spin state.[5,29]At 300 K, the extrapolations of the experimental measurements are different from each other considerably at inner core pressures(Fig.2),and the results obtained herein accord with the Prescheret al.’s results.[5]The differences from previous experimental results may be related to their samples, some of which were based on powder x-ray diffraction data,[29]whereas the structure reported by Prescheret al.was based on single-crystal x-ray diffraction data.[5]The pressure and density data of the non-magnetic phase fit the third-order Birch–Murnaghan(BM)equation[31]

wherePis the pressure,ρ0is the initial density,B0is the bulk modulus, andis the pressure derivative.Thus, the bulk modulusB0=322±1 GPa and its pressure derivative=4.38±0.02 can be obtained.

Figure 3 shows the linear compression ratio of Fe7C3at high pressure.The largest compression ratio is along theaaxis,while the smallest one is along thecaxis,indicating that Fe7C3turns anisotropic under compression.The compressibility along theaaxis is greater than that along thebaxis and thecaxis because the Fe–Fe bond is along theaaxis,which is stronger and tighter than the Fe–C bond along thebaxis andcaxis.Although the compression ratios are different from those observed in static experiments, their variation trends along all three axes are consistent with previous experimental results.[28]

To ascertain the carbon content of the inner core,we calculate the thermal EOS of o-Fe7C3at high temperature and pressure by using QHA[22](Fig.4).

Fig.3.Compression ratios varying with pressure in different directions.

Fig.4.Equation of states of Fe7C3 at specific temperatures from 0 GPa to 360 GPa(dashed lines indicate Fe isotherms).

In general,the results obtained herein are consistent with the experimental measurements of Lai.[28]In addition,we calculate the EOSs for h-Fe7C3in this study,but they almost coincide with those obtained for o-Fe7C3at high temperature and pressure.The data for pure Fe are also shown here for comparison.[32]At inner core pressures,the density of Fe7C3is 9.47%lower than that of pure Fe at a temperature of 5000 K.

3.3.Sound velocities

To ascertain the carbon content of the inner core,it is necessary to determine the elastic properties of Fe7C3at high temperatures and pressures.However, there are clear differences among the experimental measurements of sound velocities in Fe7C3at high pressures (Fig.5), even at room temperature.In this study,the sound velocity of Fe7C3at 0 K is calculated by using the method described in Section 2.Because we are more interested in the sound velocity of Fe7C3at core pressure, we only calculate the sound velocity for Fe7C3without spin restriction(Fig.5).

Fig.5.Sound velocities of Fe7C3 obtained in this and previous studies,where Xs are cited from the preliminary reference Earth model[PREM],[33] showing density-dependent(a)compressional waves,(b)Debye velocities,and(c)shear wave velocities.

The calculatedVPof Fe7C3increases linearly with density, complying with Birch’s law:VP=−1.72(±0.12)+1.15(±0.01)ρ.At 0 K, the calculatedVSis also linearly related to the density,VS=0.66(±0.05)+0.34(±0.00)ρ.The Debye velocity(VD)of Fe7C3can be calculated from the following formula:

The calculatedVP,VD, andVSat 0 K are broadly consistent with the experimental results at 300 K obtained by Prescheret al.,[5]and higher than those obtained by Chen and Zhang,[4]particularly at high pressures.Figure 5 also shows the ab initio values.[27,30]The results obtained herein are similar to those obtained by Liet al.[27]and Daset al.,[34]but they are all less than those calculated by Ghoshet al.[30]Ghoshetal.also suggested that this deviation might be attributed to the nonlinear mixing of phonon modes.[30]

4.Geophysical implications for inner core composition

According to the results of this study, the effect of temperature on the sound velocity of Fe7C3is unknown.Our previous experimental results for Fe and Fe–Si alloys suggest that Birch’s law holds true at high temperatures and pressures.[35]Assuming that theVPof Fe7C3follows Birch’s law at high temperatures,we obtain theVPof Fe7C3after determining its density at inner core.The bulk wave velocity(VB)of Fe7C3at high temperatures is calculated according to its EOS given as follows:

whereKSandKTare the bulk modulus along the isentrope and the isotherm,respectively;ρis the density;αis the thermal expansion coefficient; andγis the Grüneisen parameter.The value ofα,γ, andKTeachare obtained from the phonon spectrum, andVSat high temperatures is obtained from the following formula:

and Poisson’s ratio is calculated from the following equation:

The temperature profile of the inner core is calculated by usingT=TICB(ρ/ρICB)γ.The preferred value ofTICBandγare 5400 K and 1.5.[1]However,in this study,we only calculate the data forTICB=5000 K,as QHA often introduces large errors when the temperature is sufficiently high or close to the melting temperature of Fe7C3.[36]The reliability of QHA at high temperatures has been discussed previously.[37]Figure 6 shows the calculated value ofρ,VP,VS, andσof Fe7C3and the comparison with the values of preliminary reference Earth model (PREM).[33]Except for the density, the simulatedVP,VS,andσof Fe7C3obtained in this study deviated from those obtained by Daset al.[34]

The density of Fe and Fe7C3are about 5.87%higher and 6.34%lower than the PREM values respectively in the whole range of Earth’s inner core pressure.The Fe7C3content of 44.71 wt% can explain the density deficit of the inner core(Fig.6(a))if carbon is the only light element in the core and is present as Fe7C3.Under the inner core condition, the value ofVSandσof Fe7C3are very close to the PREM values compared with those of Fe (Fig.6(c)).TheVPof Fe7C3is slightly higher than that of Fe,and both are much larger than the PREM values.Thus, carbon might be ruled out as the major light element in the Earth’s inner core,which is consistent with findings of Fe–C phase,[38]Fe–C solid solution,[39]and geochemical[40]investigations.By analyzing the Fe–C phase diagram at high temperatures and pressures, Feiet al.suggested that the carbon content of the inner core is about 2.24%.[38]More recently, Huanget al.experimentally determined Fe–C solid solution densities up to core pressures,and concluded that interstitial carbon can lower the density of iron and be present in the inner core.[39]Using geochemical methods,McDonough concluded that the carbon content of the inner core is less than 0.2%.[40]Therefore, if carbon is present in the inner core, other light elements must also be present,which has inspired many studies of the Fe–C–Si ternary system.Pamatoet al.proposed an Fe–C–Si model for the inner core composition based on the EOS of Fe–C–Si at high pressures and temperatures.[41]Daset al.found that doping a small quantity of Si impurities at carbon sites in Fe7C3carbide phases can reduce theVPandVS,which are closer to the PREM data.[34]According to theVPandVSvelocities of Fe and Fe-8.6Si,and combining theVPandVSresults of Fe7C3obtained by Chen and Zhang,[4]Huanget al.[35]suggested that the density and sound velocity for an Fe–C–Si inner core both match the PREM values.However, theVPandVSdata obtained by Chenet al.are less than our calculated results and theab initiovalues obtained by Liet al.,[27]as well as the experimental data obtained by Prescheret al.[5]Therefore, further experiments are required to investigate theVPandVSof Fe7C3,particularly at high temperatures.

Fig.6.(a)Pressure-dependent calculated densities of Fe(black line)[35] and Fe7C3 (red line), assuming a temperature of 5000 K,compared with the PREM[33] density profile.Fe7C3 (pink line) and Fe7(C,Si)3 (blue dotted line) are cited from Ref.[34].Fe–1Si–5C (green dotted line)is cited from Ref.[35].(b)Pressure-dependent compressional wave velocities and their comparisons.(c)Pressure-dependent shear wave celocities and their comparisons.(d)Pressure-dependent Poisson’s ratios and their comparisons.

5.Conclusions

In the present study, the thermal equation of state(EOS)and sound velocities of Fe7C3at high pressures are calculated by using first-principles methods.The simulated results are generally in agreement with the previous experimental results.Although the derivedVSandσof Fe7C3can reproduce the inner core values,the calculatedVPof Fe7C3is higher than that of Fe and the inner core.Therefore,it is less likely for carbon to be present only as Fe7C3in the Earth’s core.

Acknowledgement

Project supported by the National Natural Science Foundation of China (Grant Nos.41904085, 41874103, and 42274124).