Equation of state for-iron at high pressures and temperatures

The equation of state for hexagonal close packed (hcp or ) phase of Fe at high pressure is created by employing molecular dynamics (MD) simulations in conjunction with the embedded atom method based on the full potential linear muffin tin orbital (FPLMTO) method. Comparison between the existing experimental data and our calculations suggests that the obtained equation of state can be reliably used for calculating iron volumetric properties under conditions appropriate for the Earth’s core. We demonstrate that some experimental data on iron might be subjected to a systematic error. I suggest a model which describes the temperature dependence of the volume better than the Mie-Grüneisen equation.


Introduction
An equation of state (EOS) for -Fe is required to provide constraints on the phase and chemical composition of the Earth core as well as its temperature.This is because the Earth core is likely to consist of iron.Iron is the only element which can satisfy the density derived from seismic measurements and that is also sufficiently abundant in the Earth [1,2].This subject has attracted considerable attention.Significant efforts, both experimental [3][4][5][6][7][8][9][10][11][12][13][14][15][16] and theoretical [2,[16][17][18][19][20][21][22], have been undertaken to determine the properties of iron under extreme conditions.Several structures have been suggested as relevant for iron at the conditions of the Earth core, namely hcp ( ), bcc, dhcp and orthorombic structures.Recently, equation of state for the bcc iron phase was created [23].We demonstrated earlier that the bcc phase fits a number of properties for the Earth inner core [24][25][26].Surprisingly, no robust study of the temperature impact on the -iron was performed.The existing studies are subjected to a number of approximations, such as a quasi-harmonic or particle-in-cell.Therefore, it is not known what the impact of temperature on the properties of -iron really is.The direct measurements of the iron volume at pressures above the pressure of the triple point (about 60 GPa [8]) and elevated temperatures are technically difficult.While the 300 K isotherm for iron is rather reliably determined [3], the highest pressure and temperature achieved until recently in volumetric experiments were limited to about 35 GPa and 1500 K.The comparably high temperature measurements at the highest pressure [15] manifested a major breakthrough in this field, because for the first time iron volumes were measured at pressures up to about 300 GPa and at temperatures up to 1500 K.Even though this temperature is rather modest compared to the temperatures relevant for the Earth core (up to 7000-8000 K), the experiment provides direct measurements of thermal expansion at the pressure of the Earth core.
However, a close examination of the equation of state and the experimental data to which the EOS has been fitted revealed certain problems both with the EOS and the data itself.
Due to these problems and the importance of the subject, we decided to study the EOS foriron theoretically.Of course, theoretical studies cannot disprove experimental data.In our analysis we concentrate mostly on the laws to which the data should obey.The data obtained theoretically allow us to verify the validity of various models of the EOS, because molecular dynamics treats the models of solids in an exact manner unlike other phenomenological descriptions.This is why the MD method is also known as a computer experiment.If the experimental data do not obey the well established theories of solids, there is something wrong either with the data or with the theory.Either way, a close examination of the situation is required.
The paper is composed as follows.First, we describe how we calculated the volumes (V) of -Fe at a number of PT points.Second, we show what kind of problems we encountered when analysing experimental data and the EOS fitted to those data.Third, we show that the Mie-Grüneisen equation makes it possible to describe both experimental and our data.However, while the behaviour of γ, the Grüneisen parameter, obtained from our data, exhibits the expected volume dependence, the γ obtained from the experimental data shows a quite unusual dependence on volume.Fourth, we show that a simple model can describe our data better than Mie-Grüneisen equation.Fifth, we demonstrate that the calculated thermal EOS for -Fe is in very good agreement with the existing experimental data and can be reliably used for determining the iron density under conditions of the Earth core.Sixth, we discuss possible deficiencies of our simulated data and the experimental data.

Method and results
The Fe phase diagram is relatively complicated.At ambient pressure, Fe is stable in the bcc structure due to its magnetic nature.Therefore, any attempt to create a model of the interatomic interaction in Fe, which is parameterized using low-pressure experimental data, should somehow account for this magnetic nature of the bcc phase.This creates certain problems with the general concept of parameterization.Reliable high-pressure data on iron are rather scarce, therefore, a parameterization using the properties of a high-pressure hcp Fe phase is also problematic.
The solution to this problem is to parameterize some model to fit the results of first principles calculations.However, considerable computer time is required for a calculation with even a modest number of atoms prohibiting long runs (essential for high precision calculations of the volume) and restricts the number of points which can be calculated.While bcc phase is comparatively simple case [23], the hcp phase has to be optimized for the c/a ratio and this represents a considerable difficulty.On the other hand, if a parameterized model provides a correct energy (correct means here the results from first principles calculations) for any configuration, then the use of such a model is equivalent to ab initio molecular dynamics.We should mention, though, that at high temperature the electronic entropy might be important.However, the effect of the electronic term on the volume of iron was estimated to be within a few percent of the thermal change of volume [2].Since the precision of first principles calculations is less than that, we can safely use a parameterized model.
Such a model was created in our earlier paper [40].Here we elaborate on how we did that.In order to study the electronic structure of Fe we have used the full-potential linear muffin-tin-orbital (FPLMTO) method [29].The calculations were based on the local-density approximation and we used the Hedin-Lundqvist [30] parametrization for the exchange and correlation potential.Basis functions, electron densities, and potentials were calculated without any geometrical approximation [29].These quantities were expanded in combinations of spherical harmonic functions (with a cut-off max = 8) inside non-overlapping spheres surrounding the atomic sites (muffin-tin spheres) and in a Fourier series in the interstitial region.The muffin-tin sphere occupied approximately 50% of the unit cell.The radial basis functions within the muffin-tin spheres are linear combinations of radial wave functions and their energy derivatives, computed at energies appropriate to their site and principal as well as orbital atomic quantum numbers, whereas outside the spheres the basis functions are combinations of Neuman or Hankel functions [31,32].In the calculations reported here, we made use of pseudo-core 3p and valence band 4s, 4p and 3d basis functions with corresponding two sets of energy parameters, one appropriate for the semi-core 3p states, and the other appropriate for the valence states.The resulting basis formed a single, fully hybridizing basis set.To sample the irreducible wedge of the Brillouin-zone we used the special k-point method [33].In order to speed up the convergence we have associated each calculated eigenvalue with a Gaussian broadening of 20 mRy width.
The FPLMTO energy-volume data for hcp and liquid iron were fitted with an EAM (embeddedatom method [34]) potential, because the EAM method has been shown [34] to give a thermal expansion of metals in good agreement with experiment.The particular form of the applied potential is as follows: where with Here E conf is the potential energy of a system of N atoms, E i is energy of atom i, φ is the pairwise interaction between atoms i and j, r ij is the distance between them, F (ρ) is the embedding function and ρ is another pairwise interaction leading to the density term ρ i .The functions φ, ρ i , and F (ρ) are defined as follows As a result of the fit to first principle calculated energies, the adjustable parameters were calculated to be n = 8.137, m = 4.788, = 0.0173 eV, a = 3.4714 Å and C = 24.939.
This potential was employed in our MD simulations to calculate the volumes of hcp Fe from 60 GPa to 1000 GPa and from 300 K to 12000 K.
A description of the molecular dynamic method can be found elsewhere [35].Most of the simulations were performed using the package DL POLY version 3.10 [36].To ensure the reliability of our results, some of the simulations were duplicated using our own MD code and no relevant difference was found.Simulations in the NTP ensemble were performed.The results of MD simulations at constant P and T (NTP ensemble) with the chosen model of interatomic interaction depend, apart from the initial arrangement of atoms, on the number of time steps (n timesteps ), size of timestep (∆t), number of atoms (N ), cut-off (r cut-off ) of the interatomic potential, the specified time constants for temperature (τ T ) and pressure (τ P ) fluctuations.Therefore, the effect of these parameters was carefully studied by carrying out a number of test runs at various T and P .It was found that the correct results can normally be obtained with N > 384, n timesteps = 16000, ∆t = 0.002 psec, r cut-off = 6 Å, τ T = 0.5 psec and τ P = 0.8 psec.These values were normally used unless it was specifically intended to study the behaviour of a small system.At very high pressures and temperatures the time step was smaller.The assumption of a mean-field distribution of the density was applied for the calculations of the energy and forces at r > r cut-off = 6 Å.The initial configuration for all runs was an ideal hcp lattice.

Analysis of experimental data and high-temperature Birch-Murnaghan equation of state
When we wanted to compare the results of our MD calculations with the experimental EOS [15] in the HTBM (high-temperature Birch-Murnaghan) form we noticed an unphysical behaviour of the isochores calculated with the experimental EOS (figure 1 (a)).For example, the volume at 500 GPa first increased with temperature, then decreased and subsequently increased again.This behaviour is quite evident at 500 GPa.Because of the gradual transition of isochores to the (a) (b) Volume change (a) as a function of temperature at a number of pressures according to the recently published equation of state [15].Pressure (GPa) on isobars is indicated in the boxes.The pressure step from one isobar to the next is equal to 25 GPa, except between the first and the second isobar, where the pressure step is equal to 24 GPa.∂P/∂TV calculated with the same EOS (b) at V equal to 3, 4, and 5 cm 3 /mole.Note the negative values of the derivative (unphysical behaviour).
evidently erroneous behaviour, it is likely that there is an error already at lower pressures and temperatures.As can be seen in figure 1 (b), the temperature derivative of the pressure becomes negative starting at some values of PT.Such a dependence is clearly unphysical (except for the case of molecular crystals in a very narrow PT range, e. g. ice at ambient pressure and temperature, pressure always increases when the temperature increases at constant volume).Let us consider the HTBM EOS which was used for fitting to the experimental data [15].According to this EOS we can calculate the pressure from the following expression where B T,0 and B T,0 are the bulk modulus and its pressure derivative.The relative volume η is defined as and where α is thermal expansion coefficient.The change of the bulk modulus at zero pressure with temperature is described by a polynomial The idea behind HTBM EOS is simple.The Birch-Murnaghan EOS is an isothermal equation.If at each temperature we can find a corresponding V 0 , B, and B , it will be sufficient to calculate the pressure at a given V and temperature.Presumably, the simultaneous fit should provide these values.However, if the temperature is higher than the melting temperature at zero pressure these values are fictious.Clearly, the polynomial expression for the temperature dependence of the bulk modulus does not have physical meaning and an extrapolation is problematic.Of course, if the HTBM is used in the PT range where it was fitted, then the error of the calculation is less or equal to the error of the fit.Since the form of the HTBM EOS does not intrinsically guarantee a correct physical behaviour at extrapolation it is natural to expect that at some PT range the behavior might become unphysical.Therefore, the use of the HTBM EOS for calculating iron properties at the conditions of the Earth core [15] is doubtful.
Therefore, we decided to fit the experimental data to other kinds of thermal equations of state.A widely used thermal EOS is the Mie-Grüneisen equation which was analysed in detail by Anderson [2], in particular with regard to its application to iron.According to this equation the pressure can be calculated from the following expression P (V, T ) = P (V, 300) + P thermal (V, T ). ( The term P (V, 300), usually referred to as the cold pressure, can be calculated using the Birch-Murnaghan EOS, actually the same as the HTBM EOS (see above) at the temperature of 300 K.This term for iron is comparatively well known [3].The second term can be calculated using the Mie-Grüneisen expression Anderson [2] discussed the reliability of this expression and suggested the following expression for the calculation of γ We fitted the experimental data by means of the Mie-Grüneisen equation (which has 5 fitting parameters, i.e. one less than what was used in the HTBM EOS [15]) and obtained a better fit than using the HTBM EOS.The fitting errors are shown in figure 2 as a function of pressure (figure 2 (a)) and temperature (figure 2 (b)).One can see that the maximum error is smaller when using the Mie-Grüneisen EOS.The mean square deviation from experiment is also smaller with the "gamma" (Mie-Grüneisen) EOS.Since the "gamma" equation has less parameters, provides a better fit to the experimental data and is physically justified, it should be considered preferable over the HTBM EOS.The calculated γ 0 and q, obtained from the fit are equal to 2.123 and -0.978, correspondingly.The value for γ 0 is normal, but the value of q 0 is anomalous.Normally, in all descriptions of iron PVT data, including the Hugoniot shockwave adiabat, the value of q 0 ranges from 0.5 to slightly above 1.The negative value of q 0 indicates that there is something very unusual about the experimental data.This can also be seen in the figure 3 where values for γ fitted to the experimental data are compared with the γ, calculated from the "cold" EOS of iron [37] with the expression Equation ( 14) is quite general and, not being particularly precise, correctly predicts the general change of a thermal expansion with volume.The opposite behavior of γ exp and the theoretical γ (equation ( 14)) prevents us from using the experimental data for deriving an EOS for iron under extreme conditions.Therefore, we decided to attempt to use the MD generated data for volumes of -Fe at high pressure and at a series of temperatures.[15] and MD data as a function of volume.These values of γ are compared with the γ calculated from 300 K isotherm for iron using the method by Dugdale and MacDonald [37].While γ fitted to MD data and the theoretical γ shows similar behaviour, the γ fitted to experimental data [15] shows an anomalous behaviour.

Equation of state of -iron and comparison with available data
We calculated a set of MD volumes at a number of pressures and temperatures.The largest difference between experimental [3] and calculated volumes was less than 3%.This is a normal precision when metal volumes are calculated using an embedded-atom method [34].However, if the purpose is to create an EOS which could be useful, it is desirable to minimize the errors.Therefore, to derive an EOS for iron we used the following procedure.The cold isotherm was accepted exactly as the isotherm with equation ( 11) at the temperature of 300 K, fitted to experimental data [15], which, in turn, is close to the previously determined isotherm [3] at 300 K.The isotherm is defined by three parameters, V 0 = 6.695 cm 3 /mole (volume at the zero pressure and temperature 300 K), B = 173.98GPa (bulk modulus), and B = 5.297 (bulk modulus derivative).Using the volumes calculated with these parameters and the set of δV = V P,T − V P,300 , calculated using MD (figure 4, symbols), we composed the set of P − V − T points, which was fitted using the equations ( 11)- (13).As a result of the fit, in addition to the three above mentioned parameters, we obtained the parameters γ 0 = 2.434 and q = 0.4894.This set of parameters completely defines the volumetric properties of -Fe at high pressures and temperatures.
Figure 2 shows that the errors of fitting MD data using equations ( 11)-( 13) are comparably large (although smaller than fitting the experimental data) at lower P and T .While P and T increase, the errors decrease.This is what one should expect because the Mie-Grüneisen EOS in the used approximation becomes more precise at high P and T .Despite comparatively large absolute errors, the maximum relative error is less than 0.5%.
The volume dependence of γ fitted to the MD data (figure 3) is similar to that calculated according to the Dugdale and MacDonald [37] expression (equation ( 14)).This is remarkable because the Dugdale and MacDonald's γ is calculated using experimental iron isotherm at T = 300 K while our γ was extracted from MD data at T > 300 K.
Summarizing, the final expression for the -Fe EOS is where and The parameters B, B , V 0 are respectively the bulk modulus, its pressure derivative, volume at zero pressure and room temperature, and γ 0 and q are parameters for calculating Grüneisen parameter.
The calculated optimal values are given in table 1.

Discussion
A judgement of the results from a simulation requires comparison with experimental data and an analysis of the sources for possible errors.Since our simulation mainly concerned the temperature volume change of -Fe at high pressure, it is natural to compare our EOS with high pressure-high temperature volumetric experiments.We should mention that thermal expansion of metals calculated using EAM method agrees well with experimental values [34].As it follows from the consideration given above, not all of the experimental data are reliable.As regards the measurements of thermal expansion, the most reliable data come from the experiments that used comparatively large samples.A use of large samples allows us to significantly decrease the errors related to non-hydrostatic and thermal stresses and to increase the precision of X-ray structural resolution [5].Thermal expansion can also be extracted from shockwave data [6].In this case, unless the temperature in a shockwave experiment was measured, thermal expansion is merely an estimate, because thermal expansion and temperature are not independent parameters.Even if the temperature was measured, one should remember that those measurements are a subject of debates because the measurements themselves involve a number of poorly determined parameters.The experiment conducted by Funamori and co-authors [5] should, in our opinion, be considered as the most reliable determination of iron volume at high temperature and pressure.Figure 5 shows a comparison between thermal expansion as it is defined from experiment and our equation.The nearly exact match of the experimental data [5] is possibly fortuitous, because our MD simulations do involve certain approximations.Also, some errors were introduced due to the fitting procedure.Nevertheless, a very close match between our equation and experiment suggests that our equation based on a clear physical basis, might be useful in calculating the EOS for iron under the conditions of the Earth core.It is worth mentioning that some years ago I predicted the thermal expansion for MgSiO 3 -perovskite [38] which later was measured [39] to be in perfect agreement with our prediction.15)- (17) with the parameters listed in table s1 (solid curves).The expansions calculated for a number of temperatures starting with 1000 K and up to 9000 K with a step of 1000 K between the curves.The symbols indicate experimental measurements of the thermal expansion.The experimental data are for temperatures between 1500 K and 2000 K [10], 1000 K [5], and 5200 K [6].
In our MD calculations we neglected the effect of the electronic entropy on the volume change of iron.In terms of our expression, which we used for fitting the MD calculated volumes (equation ( 15)- (17), this means that the parameter γ should be somewhat larger.However, as it follows from general considerations [2], the change of γ due to the electronic degrees of freedom is comparably small, about 5% at the Earth's core conditions.This introduces, therefore, about 5% correction to the value V P,T − V P,300 and about 0.5% of the V P,T .This is a precision which is not essential from a practical point of view.
The obtained EOS is in good agreement with shock-wave data [13].The agreement becomes nearly perfect if instead of the parameters quoted in table 1 for B, B , and V 0 one uses the values 178.2 GPa, 5.15, and 6.74 cm 3 /mole as obtained from the reduced shock-wave data proper.This confirms that the temperature volume change as it comes out from MD simulations is in good agreement with shock-wave data.
The errors of fitting MD data can be substantially decreased if instead of using the Mie-Grüneisen expression one uses the following dependence where A = 7.416×10 −3 and B = 2.087×10 −7 if the pressure is given in GPa units and δV in cm 3 /mole.The quality of the fit using equation ( 18) can be seen in figure 4. The physical meaning of this equation is rather transparent.The elastic energy (P δV ) at constant P is equal to the thermal energy, i.e. to c P T .The heat capacity at constant pressure (c P ) is equal to This explains why a second degree polynomial is necessary and sufficient to describe the volume change with temperature (equation ( 18)).We do not here provide a rigorous justification of this equation and it should be treated rather as an empirical observation based on "computerexperimental" data.It would be interesting to check how well this equation works for other materials.
Figure 6.Comparison between the density in the Earth inner core [41] and the density of -Fe.The iron density was calculated at the melting temperatures of iron corresponding to these pressures, i.e. from T = 7080 K at P = 330 GPa and up to T = 7520 K at P = 360 GPa.
Anyhow, the difference between using equations ( 15)-( 17) and equation ( 18) is very small.Therefore, we have chosen to use the Mie-Grüneisen equation because it is widely accepted.Using the obtained EOS we calculated the density of iron at the pressures of the Earth's inner core.At each pressure the temperature was chosen as the temperature of iron melting at the same pressure, as it comes out from our recent work [40] (from 7100 K at 330 GPa and up to 7580 K at 360 GPa).One can see (figure 6) that the calculated iron density is very close to the density of the Earth's inner core as it comes from the seismic measurements generalized by the Preliminary Earth Model (PREM [41]).The comparison shows that the hcp iron is about 1.5% denser than the Earth's inner core.This is close to the precision with which the 300 K isotherm of iron is measured.Assuming that the introduction of impurities in an iron matrix does not change its volume, it would be sufficient to substitute 2%-5% of iron atoms by impurities, such as S or C.This is significantly less than what was suggested previously [42].In our recent study on EOS for the bcc iron [43] we found that to match the core density, one needs about 6 to 7% of Si added to pure iron.This means that, contrary to common beliefs, bcc is more dense than hcp under core P and T conditions.At room temperature and high pressure, where bcc phase is dynamically unstable, bcc phase is less dense than the hcp phase.This clearly shows that the attempts to use the data on the low temperature bcc phase to make conclusions on the iron in the Earth' core should be abandoned.

Conclusions
The present study was motivated by a close examination of some iron equations of state based on the experimental data at high pressures and temperatures [15].The analysis revealed certain problems with the EOS and the data with which the EOS was fitted to.First principles calculations for hcp and liquid iron combined with molecular dynamic simulations gave us the set of temperature volume change points which, in combination with comparatively well known experimental 300 K iron isotherm allowed us to derive a thermal EOS for iron in the Mie-Grüneisen form.The EOS is in very good agreement with experimental data.This equation might be useful when applied to the development of various models of the Earth.We suggest an expression for describing thermal energy which is alternative to the Mie-Grüneisen expression.The expression describes "computer experimental" data better than the Mie-Grüneisen one.A test for other materials is required to make a decisive conclusion on the usefulness of the suggested alternative.It follows from the comparison of our EOS and the density of the inner core that the inner core can be nearly pure iron, with a lower than previously believed concentration of impurities.

Figure 2 .
Figure 2. Difference between experimental and calculated volumes of iron.Sources for the experimental data and the equations used for the fitting are indicated in the legend (MD data are considered as a result of computer experiment).For a detailed description of the equations used for the fitting see the text.The errors of the fit are plotted against pressure (a) and temperature (b).

Figure 3 .
Figure 3. Parameter γ, fitted to experimental[15] and MD data as a function of volume.These values of γ are compared with the γ calculated from 300 K isotherm for iron using the method by Dugdale and MacDonald[37].While γ fitted to MD data and the theoretical γ shows similar behaviour, the γ fitted to experimental data[15] shows an anomalous behaviour.

Figure 4 .
Figure 4.The calculated (MD) temperature change of the -Fe volume for a number of pressures indicated on the legend (symbols).The solid lines represent the fit provided by equation (18).

Figure 5 .
Figure 5. Thermal expansion, calculated using equation (15)-(17) with the parameters listed in table s1 (solid curves).The expansions calculated for a number of temperatures starting with 1000 K and up to 9000 K with a step of 1000 K between the curves.The symbols indicate experimental measurements of the thermal expansion.The experimental data are for temperatures between 1500 K and 2000 K[10], 1000 K[5], and 5200 K[6].