Ab initio determination of iron melting at terapascal pressures and Super-Earths core crystallization

,


I. INTRODUCTION
During the last two decades, several thousand exoplanets have been detected [1][2][3].Measurements of their masses and radii have become more and more accurate, which allowed us to place constraints on their composition and better understand their atmospheres, formation, and evolution.Among the detected exoplanets, there are many Super-Earths, which are assumed to have a rocky composition but are larger than Earth.They are a particularly interesting type of exoplanets because they have no analogues in our solar system.Their interior structure has been the subject of numerous studies that have tried to constrain their composition [4][5][6][7][8][9][10][11].Observations combined with modelling suggest that planets larger than 1.6 Earth radii are not purely rocky [12].Recently, it has been proposed that Super-Earth may form from sub-Neptune sized planets that lose their H/He envelopes by irradiation from a supermassive black hole [13].
Developing realistic planetary interior models to infer the composition of Super-Earths requires knowledge of the equation of state (EOS) of the candidate minerals such as iron and silicates at extreme pressure-temperature conditions [14,15], which represent a challenge in planetary and materials science because the conditions of interest often lie outside the reach of laboratory experiments.Models predict that pressures in the interior of Super-Earth planets can exceed 1000 GPa at the core-mantle boundary, and temperatures can exceed 10000 K [4,5,8].More recently Boujibar et al. [16] studied the possible range of temperatures at the core-mantle boundary, for which ca solid inner and a liquid outer core coexist, and showed that it depends on the planet's total mass and its core-mass fraction, but also sensitively depends on knowing the melting temperatures * f gonzalez@berkeley.edu of iron at TPa pressures.The formation of a metallic iron core in these planets is driven by chemical differentiation and gravitational separation of liquid silicate-iron mixtures.The generation of magnetic fields is a direct consequence of the presence of liquid iron in the core.In the Earth, the crystallization of a solid inner core is assumed to be a major driver for the magnetic dynamo.Thus, studying the melting behaviour of iron at pressures of Super-Earth interiors will contribute to a better understanding of their internal structure, core crystallization, and dynamo activity.
Because iron is the main constituent of the Earth's core, the characterization of its high-pressure properties have been of fundamental importance geophysics and condensed-matter physics.Planetary formation models predict it to be the main constituent of the core of the other terrestrial planets as well.The phase diagram of iron and, in particular, its melting line is not well understood for pressures exceeding 300 GPa, but both experiments and first-principles quantum mechanical simulations have shown that iron transforms from the body-centered-cubic (bcc) phase to the hexagonal-closedpacked (hcp) phase under pressure [17,18].
The melting curve of iron at pressures relevant to the Earth innercore boundary (330 GPa) has been explored with different computational methods.Predictions published in the last two decades place the melting temperature of iron between 6300 K and 7300 K for this particular pressure [19][20][21][22].Using thermodynamic integration, Alfè et al. [20] performed density functional theory molecular dynamics (DFT-MD) simulations and obtained a melting temperature of 6350 K, which was later confirmed by first-principles two-phase simulations in the microcanonical ensemble, using a relatively large number (1000) of atoms [21].In a recent DFT study, Bouchet et al. [22] performed two-phase simulations in the canonical ensemble to extend the melting curve of bcc iron up to 1500 GPa.However, since the hcp structure is predicted to be the stable phase at these pressures [18,[23][24][25], two-phase simulations of hcp iron may actually lead to higher melting temperatures.Nevertheless, the results in Ref. [22] are in good agreement with the aforementioned DFT studies [20,21].
A separate study by Sola and Alfè [26] showed that quantum Monte Carlo simulations, which treats correlation effects more accurately than DFT, tend to favor higher temperatures for ICB conditions, placing the melting temperature at 6900 ± 400 K for 330 GPa.Belonoshko et al. [19] fitted an embedded atom potential (EAM) to DFT simulations of hcp iron and performed two-phase simulations that predicted a melting temperature of 7100 K, similar to the quantum Monte Carlo predictions.In a recent work, these authors coupled the EAM simulations to thermodynamic integration to obtain the free energies of bcc and hcp iron and obtained a melting temperature of 7190 K at 360 GPa for bcc and slightly lower, 6800 K, for hcp [27].The most recent experiments on iron at TPa pressures by Kraus et al. [24] used dynamic shock and ramp compression [28,29] to show that at 330 GPa, iron melts from the hcp phase at 6230±540 K, in agreement with previous experiments in diamond anvil cells by Anzellini et al. [30].
In this work, we present a DFT-MD calculation of the melting curve of iron from 300 to 5000 GPa in order to cover the conditions present in Super-Earth interiors.The melting temperatures are obtained from free energy calculations that are based on a thermodynamic integration (TDI) method.With this technique, we determine which phase, liquid or solid, has lower Gibbs free energy for a given pressure and temperature.For a given pressure, the melting temperature is obtained by interpolating the free energy difference between solid and liquid phases as function of temperature.We also use the TDI method to derive the entropy as a function of temperature and pressure, which we enables to derive isentropes to characterize the temperature profiles in the cores of Super-Earths.
In Sec.II, we describe our computational methods and Gibbs free energy calculations.In Sec.III, we discuss our melting curve and adiabats.We show that our melting temperatures are considerably higher than previous predictions that were extrapolated from data at lower pressure.We also predict that the melting curve is steeper than the adiabats, which implies that the crystallization of iron cores of Super-Earth would always start from the center, provided that core crystallization occurs in a planet's lifetime.In Sec.IV we build models for Super-Earth interiors focusing on planets with up to 1.6 Earth radii and 5.8 Earth masses.For relevant pressure conditions, we find that our melting temperature is much higher than the temperature profiles in available interior models, suggesting that the cores of Super-Earths are completely frozen over their entire lifetime.

A. Ab initio techniques
Using density functional theory molecular dynamics (DFT-MD), we calculate the Gibbs free energy of iron for the solid and the liquid phases at specific pressure-temperature conditions.The melting temperature, T m , is obtained when the Gibbs free energy difference vanishes for a given pressure.For most solid calculations, we assume an hcp crystal structure, but we also performed a few calculations near 300 GPa using the bcc phase and one calculation for the melting temperature of the fcc phase at 5000 GPa.
An alternative approach for calculating the melting temperatures with computer simulations is the Z method, which has been applied to a number of other materials [31][32][33][34][35][36][37].This technique is independent of the Gibbs free energy calculations, and is based on overheating the solid in the microcanonical (N V E) ensemble.The Z method relies on the fact that any solid system that has been sufficiently overheated will spontaneously melt, provided the simulations are long enough.As the latent heat is removed, the temperature of the system drops.If the amount of overheating is carefully calibrated, temperature will drop precisely to the melting temperature, T m .Here we applied the Z method to study the melting of iron at different densities to compare with the melting points we derived from Gibbs free energy calculations.
We perform our DFT-MD simulations with the VASP code [38] using exchange-correlation functional of Perdew, Burke, and Ernzerhof [39].We used a Mg-core pseudopotential of the projectoraugmented wave type [40] with 14 valence electrons per atom (PAW-14) and a core radius of 1.16 Å.Some calculations were done using a Ne-core pseudopotential (PAW-16).Single-particle orbitals have been expanded in plane-waves with a cutoff of 1100 eV in all calculations.The DFT-MD simulations were carried out under the assumption of Born-Oppenheimer approximation.A time step of 0.5 fs was employed and the simulations lasted between 1.0 and 12 ps.We performed preliminary DFT-MD simulation in an hcp cell with 96 atoms using 2 × 2 × 2 and Γ-only k-point grids.We found a that the energy and pressure were underestimated by 150 meV per atom and 14 GPa, respectively by the Γ-only k-point grid for pressures close to 5000 GPa.
All following simulations of both solids and liquids were then performed in a larger cell, containing 144 atoms (except for the bcc phase, where we employed a supercell with 128 atoms) with Γ-only k-point grid in order to converge the thermodynamic properties and prevent dynamic instabilities in the overheated solids.The same parameters were considered for the Z method calculations, with the exception of the simulation time, which was extended up to 8 ps in some cases to ensure the stability of the overheated solids.DFT-MD simulations with 144 atoms and a 2 × 2 × 2 k-point grid showed that the total energy differs by less than 4 meV per atom respect to the Γ-only k-point grid.We also performed calculations with a larger cell of 180 atoms and obtained thermodynamic properties that were consistent with the 144 atom results, which we therefore considered sufficiently well converged for the purpose of this study.The c/a ratio in the hcp supercells was adjusted for every pressure-temperation condition in order to obtain hydrostatic conditions.For liquids, we used cubic cells with 144 atoms and Γ point to sample the Brillouin zone.In addition, DFT-MD simulations were performed with four other pseudopotentials available in VASP, which treat 8 or 16 electrons explicitely with PBE and PW91 functionals.These simulations also used 144 atoms and Γ-point sampling.
We also performed simulations in larger cells of up to 1296 atoms to test the convergence with respect to system size and simulation time.To do this, we trained an on-the-fly machine learning potential, as implementented in VASP 6.We trained the force field on a supercell of 144 iron atoms in the hcp phase using the PAW-16 pseudopotential at 12.938 g/cc and 6000 K. Simulations were carried out at both constant volume (N V T ) and constant pressure (N P T ), obtaining consistent results.After running them for over 25 ps, we observed no significant change in the pressure or free energy of the system with respect to our DFT-MD simulations performed with 144 atoms.We provide more details in the supplementary material.

B. Computation of Gibbs Free energies
Free energy calculations require the knowledge of the entropy, which is not directly accessible from the standard MD simulations.The anharmonic contributions to the free energy of iron are large [41], so a description of the solid phase with quasi-harmonic methods alone would not be appropriate and lead to incorrect melting temperatures.One of the available methods to address this problem is the thermodynamic integration (TDI), which is a general technique to determine the difference in Helmholtz free energy between two systems with potential-energy functions U a (r i ) and U b (r i ).By defining a hybrid potential U λ = U a + λ(U b − U a ), the difference in Helmholtz free energy between the two systems can be computed from where one averages over configurations, r i , generated with forces derived from the hybrid potential.
In order to increase computational efficiency, we adopt a twostep TDI scheme as implemented in previous studies [42][43][44][45][46][47][48].In this scheme, we introduce an intermediate system governed by a classical pair potential, U cl , which we constructed for every density and temperature by matching the forces along the DFT-MD trajectories [49].For every temperature-volume condition, we perform a preliminary DFT-MD simulation and fit a new pair potential to the forces along the computed trajectory.More details on how this method is applied in liquid systems can be found in references [43], [47], and our supplementary material.The TDI procedure for solid systems is discussed in Refs.[42,[45][46][47]50].In the first integration step, one derives the free energy difference, ∆F cl→DFT , between the system interacting with the DFT potential, U DFT , and the system governed by the classical pair potential, U cl .Because the classical forces match the DFT forces well, 5 equally spaced λ points are sufficient to obtain an accurate value of the integral in Eq. ( 2) in this step.
Then a second TDI step is performed in order to obtain ∆F ref→cl , the free energy difference between the classical system and a reference system with known Helmholtz free energy, F ref .For liquids, we chose the ideal gas as the reference system.For solids, we employed an Einstein crystal as the reference system with a combination of two-body and one-body harmonic potentials for the classical system [42,45,47].Recent work has shown that the correction by Frenkel et al. [51] to the free energy of a solid due to a fixed center of mass has been overestimated [52].This correction overestimated of the stability of the solid phase, leading to melting temperatures that were slightly too high as was discussed for materials such as MgO and Be [50,[52][53][54].The alternate correction proposed in Ref. [52] is much smaller than the Frenkel correction and, in practice, it is equivalent to applying no correction because it is smaller than the error bars we obtain.Not applying the Frenkel correction due to a fixed center of mass increases the free energy of the solid which slightly lowers the resulting melting temperature.
Since only classical Monte Carlo simulations are needed, the second integration step is computationally much less expensive (by factor of ∼10 −5 ) than the first step that requires solving the Kohn-Sham equations.DFT calculations.This remains true even though a larger number of simulations is required to accurately compute the integral in Eq. ( 2).The Helmholtz and Gibbs free energies of the DFT system are then obtained from and In order to align the Gibbs free energy of the solid and the liquid at the same pressure, we use the thermodynamic relationship where G 0 = G(P 0 , T ), P 0 = P DFT , P is the target pressure, and V T (P ) is the respective volume of each system along an isotherm of temperature T that we obtain from a separate set of DFT-MD simulations.
For comparison purposes, we also implemented the Weeks-Chandler-Andersen (WCA) potential, given by because it was employed in one earlier study of iron [55].We use it in lieu of our fitted pair potentials when we conduct test calculations for a few cases.This potential adopts the repulsive part of Lennard-Jones potential but removes its attractive part.Mirzaeinia et al. [56] computed the free energy of the WCA liquid by performing a TDI along a path at constant temperature to the low-density limit where the system's free energy is known analytically and corresponds to the ideal gas.This path provides an alternative to the integration path at constant density toward the limit of infinite temperature that we typically employ because our pair potential are finite at the origin.Still, we confirmed the accuracy of the subset of the results in Ref. [56] that are relevant for this study.Following Ref. [55], we set the potential parameters to each temperature an volume such that the reduced temperature T * ≡ k B T /ε = 1.5 and the reduced volume η ≡ πσ 3 /6V = 0.1 to ensure that the system is in the liquid phase.

A. Melting at 330 GPa
Determining the precise melting temperature of iron at the Earth's inner-core boundary conditions, 330 GPa, is crucial for developing models of Earth's interior and understanding its core's crystallization and heat flow.Despite numerous attempts through simulations and experiments, a consensus on the precise melting temperature of pure iron at this pressure has been difficult to reach, as we illustrate in Fig. 1.This figure reveals discrepancies among the predictions for the melting temperature of iron at 330 GPa between different theoretical techniques, such as thermodynamic integration [20,26,27,55,57], two-phase simulations [21,22,58], free-energy based models [59], and the Lindemann law [18], as well as between experiments that have explored similar conditions [23][24][25]30].Conducting experiments and accurate computer simulations has remained a challenge even at lower pressures of the core-mantle boundary [60], and the source of discrepancies among various study is often difficult to identify.[23][24][25]30], while solid symbols represent predictions from ab initio simulations using different methods: TDI [20,26,27,55,57], two-phase simualtions [19,21,22,58], freeenergy based models [59], and the Lindemann law [18].The temperatures from Refs.[27,57] have been interpolated from the values that the authors report at similar pressures.

B e lo n o s
Because of these discrepancies, we have performed careful calculations with our TDI methods for solid and liquid iron at 330 GPa and temperatures ranging from 6000 K to 6800 K.As explained in Sec.II, we fit a new pair potential for every single thermodynamic condition so that our TDI calculations are independent from each other.For the solid, we use the Einstein crystal as the reference system with known free energy.The TDI is conducted in two steps.For the solid phase, we first integrate from the DFT system to a classical system with pair and harmonic Einstein forces.In the second integration step, we gradually turn off the pair forces.
For the liquid, we follow a similar procedure.First we switch from the DFT potential to a classical pair potential and then we turn off the classical forces, leading to a gas of non-interacting particles.To double-check our predictions, we repeated the TDI calculation for the liquid using the WCA potential instead of using our fitted pair potential as the intermediate system, as done in Ref. [55].In Fig. 2, we show an example of the integrand in Eq. ( 2) as a function of λ for both TDI implementations for liquid iron at 330 GPa and 6400 K.While for the pair potential the integrand is linear and varies by ∼0.01 eV/atom (bottom panel), for the DFT→WCA integration the integrand has a large curvature in the same interval and varies ∼18 eV/atom because the WCA potential does not represent the forces between the iron atoms very well.
The DFT→WCA integration can be performed directly as a function of λ if the integrand has been constrained by a sufficient number of points, each requiring computationally expensive DFT-MD simulations.However, fewer integration points may be used if a change of variables λ → λ m is applied with 0 < m 1 in Eq. ( 2).This allows one to work with an integrand that varies more smoothly and obtain the integral via Gaussian quadrature, as done by Sun et al. [55].We sampled the integrand for this particular condition with 12 points, so we are not required to change the integration variable.We were able to obtain a smooth interpolation using cubic spline function for m = 0.5, m = 0.25 and the original integration (m = 1) as we show in Fig. 2. We obtained a consistent value of ∆F WCA→DFT = −3.893± 0.007 eV/atom that varies within the given error bars for all three values of m.
For 330 GPa, Sun et al. [55] reported a melting temperature of 6200 K, which is lower than our value of 6522 K.We primarily attribute this difference to a deviation in the computed Gibbs free energies of the liquid.In the two lower panels of Fig. 3, we show that our liquid Gibbs free energies for PAW-16 pseudopotential are about 35 meV/atom higher than the corresponding values by Sun et al..While we could not extract sufficient details from Sun et al. article to fully explain this difference, we want to point out parts of the free energy calculations that we agree on and others that we do not.We added Fig. 2 for this comparison.
This figure also shows that using the WCA potential as classical system results in a very pronounced curvature of the integrand, which requires both a very careful sampling and a precise integration and interpolation methods to reduce the error in the final integral.Fig. 2 shows that the integrand for the WCA potential varies by 18 eV/atom between λ = 0 to λ = 1, which makes it difficult to control the integration to meV/atom precision level.The lower panel shows when our pair potentials are employed, the integrand is approximately linear and varies by only 0.013 eV/atom.Furthermore, the WCA integral is very sensitive to the quality individual points because one needs to capture the curvature of the function correctly.For example, removing 3 out of 12 points can introduce differences of more than 100 meV/atom in the resulting integral.
We agree on the average energy of the WCA potential [56].Our value of U WCA λ=0 = 33.1 ± 0.6 meV/atom agrees with the 33.9 ± 0.3 meV/atom found by Sun et al. at their smallest value of λ = 10 −6 .They implemented the TDI between the DFT and WCA potentials in a slightly different way, U λ = U WCA + λU DFT , which requires to integrate U DFT λ (open circles and open squares in Fig. 2) instead of our expression, U DFT − U WCA λ (filled symbols).
Sun et al. performed the their TDI calculations using the PAW-8 pseudopotential, which means their values of U DFT λ differs from ours even in the limit of λ = 0, as we can see in the Fig. 2. Sun et al. then corrected their resulting Helmholtz free energies, F 1 , to achieve PAW-16 precision, F 2 , using free energy perturbation theory: On the other hand, as we can see in the bottom panel of Fig. 2, our implementation of TDI, which involves fitting a pair potential that matches the interatomic forces at this specific temperature and pressure, is much less sensitive to the sample quality and interpolation methods because U DFT − U cl λ changes by much less and varies mostly linearly with λ.The choice between using a linear, Gaussian quadrature, or spline interpolation method, and whether we perform the integral directly or using the change of variables in Eq. ( 7), varies the integral by less than 10 −4 eV/atom (∆F cl→DFT = −22.7313± 0.0005 eV/atom).Even if we use two points only, λ = 0 and λ = 1, we can obtain the same value of the integral within 0.3 meV/atom.Therefore, our implementation of TDI is more robust than an integration that uses the WCA potential.
In Fig. 3, we show the resulting Gibbs free energy of liquid and solid iron at 330 GPa that we obtained from TDI results in Fig. 2, using both the WCA and our fitted pair potentials.We repeated these calculations for 6000, 6200, and 6800 K and compare the results between the two TDI implementations for liquids.Our Gibbs free energies of solid iron are in very good agreement with those provided by Sun et al. [55], and a calculation we have done with a larger simulations cell with 240 atoms confirms the convergence values, but the latter approach yields much larger error bars.We compare our results to those from Sun et al. [55].
with respect to the system size.As we can see in the bottom panel, there is remarkable agreement between our two implementations of TDI for liquids.But when we utilize our fitted pair potential, the error bars on the Gibbs free energy are an order of magnitude smaller than for the WCA potential.
However, we observe a systematic offset of ∼35 meV/atom between our calculated Gibbs free energies for liquid iron and those reported by Sun et al [55].Their results rest on a single TDI calculations with the WCA potential for a temperature of T 0 = 6400 K. Values for all other temperatures are calculated along the isobar by utilizing the thermodynamic relationship, G/T = G 0 /T 0 − T T0 H/T 2 dT .So if there is an offset at 6400 K, all other Gibbs free energy values will be shifted, while our calculations for different P -T conditions are independent.
The middle panel of Fig. 3 shows that our solid and liquid free energies become equal at a temperature of 6523 ± 8 K, which is our prediction for the melting temperature at 330 GPa using the PAW-16 pseudopotential.In the same panel we show that when we repeat the calculation with the PAW-14 pseudopotential, ∆G = G liq − G sol increases by a modest amount of 20 meV/atom, which implies a melting temperature increase of approximately 200 K.In Ref. [55], the authors derived a ∆G value of 32 meV from perturbation calculations, which very close to the our result that we with full TDI calculations.The deviation between the two ∆G values is within the error bars of most numerical methods (see Fig. 1).
Using the PAW-14 pseudopotential, we find a melting temperature of 6747 ± 14 K at 330 GPa.Previous free energy calculations of hcp iron performed by Alfè et al. [20,61] used a non-standard PAW-8 pseudopotential that mimicked the inner electrons by introducing a repulsive interaction via a classical pair potential correction in the PW91 functional.Their results point towards a lower melting temperature of ∼ 6350 ± 300 K, with an additional error of 300 K associated to the inherent DFT inaccuracies, including the choice of different pseudopotentials.This is close, but still lower, to our prediction using the PAW-16 pseudopotential.A later work by Sola and Alfè [26] revisited these calculations using diffusion Monte Carlo calculation that incorporate correlation effects more accurately, and obtained a higher temperature of 6900 K ± 400 K for 330 GPa.
We have found, in agreement with other authors [20,55,57,61], that the predicted melting temperature of iron depends on the pseudopotential used, and the discrepancies become particularly pronounced when the 3p electrons are not considered.In order to address this issue in more depth, we performed a detailed study on how much the Gibbs free energies of liquid and solid iron are affected by the choice of the pseudopotential.For four additional VASP pseudopotentials, we redetermined the pressure-density relation along an isotherm, recomputed the free energies and melting temperatures.We performed these calculations for the three pressures of 300, 3000, and 5000 GPa to cover the entire pressure range of interest.We considered two PBE and two PW91 pseudopotentials with 8 and 16 valence electrons per atom, respectively.In the case of the 8-electron pseudopotentials, the 3s and 3p electrons are frozen while they were treated explicitly for the 16-electrons per atom pseudopotentials.We always find consistent melting temperatures when we switched between PBE and PW91 functionals for calculations with pseudopotentials that have the same number of electrons.Furthermore, the predictions with the small-core pseudopotentials, that treat 16 electrons per atom explicitly, were consistent with those from the 14 electrons per atom pseudopotential that we employed for most calculations in this work, and the difference in the resulting melting temperatures are small, about ∼200 K, as we show in Fig. 3.
However, when we switched to the 8-electron PBE and PW91 pseudopotentials, the Gibbs free energies of solid and liquid iron were not affected equally.For these pseudopotentials, the liquidsolid Gibbs free energy difference at 300 GPa was ∼50 meV per atom lower than we had obtained with the 14-electron pseudopotentials.This reduced our predicted melting temperature at 300 GPa by 600 K to approximately 6200 K, which is in good agreement with the PAW-8 calculations by Sun et al. [55] who also used an 8-atom pseudopotential.Shock experiments from Yoo et al. [62] also suggest, as we do, a higher melting temperature.However, at 300 GPa, they report a melting temperature of 6720 K, which is 250 K above our predicted melting temperature using PAW-14 and 470 K below using PAW-16 at this pressure.These results by Yoo et al. [62] were recently reinterpreted by Kraus et al. [24], suggesting a much lower melting temperature at shock conditions that is consistent with their measurements and previous experiments by Anzellini et al. [30].Despite of these differences, we find very good agreement among the predictions from various pseudopotentials for the volume and entropy differences upon melting.

B. Thermodynamic stability of the bcc phase
Other studies have considered that iron melts from the bcc phase at these conditions [27,57,58,63] because it has been argued that this phase may become stable at Earth inner core pressures via a self-diffusion mechanism [64][65][66].However, the bcc phase has not been observed in experiments at these conditions [23][24][25]67] and the size-dependent anomalies of the bcc cell are still controversial [68].In Ref. [27], Belonoshko et al. used thermodynamic integration to compute the free energy differences between the hcp, bcc, and liquid phases of iron at 120 and 360 GPa.While the simulations were done with classical molecular dynamics using the EAM, which does not take the electronic entropy into account explicitly, the contribution from the electrons to the free energy was obtained from ab initio calculations to correct the EAM-based free energies.Their results show that the free energy of the bcc structure is lower than that of hcp by ∆G = G hcp − G bcc = 23 meV/atom at 360 GPa and 6000 K and that the resulting melting temperature of iron in the bcc phase is 7190 K at 360 GPa, higher than the melting temperature of hcp at the same pressure (6800 K).However, our TDI calculations based on DFT-MD at these conditions indicate that the free energy of the bcc structure is considerably higher, hence less stable than hcp, with ∆G = G hcp − G bcc = −65 meV/atom and a melting temperature for hcp of 6922 K at 360 GPa.At 330 GPa, this difference is −57 meV, favoring hcp at Earth inner-core boundary pressures.Therefore, although our calculations show that the hcp is much more stable than the bcc structure, our resulting melting temperatures for both structures are still similar.
In Ref. [57], Y. Sun et al. implemented a different type of thermodynamic integration but used a similar approach to that of Ref. [55], where they performed the TDI calculations with PAW-8 and corrected the free energies to PAW-16 accuracy through free energy perturbation.In contrast to the findings of Ref. [27], and along with our predictions, the authors obtain a higher melting temperature for the hcp phase than for the bcc phase at 360 GPa: 6692 ± 45 K and 6519±80 K, respectively.Near 330 GPa, the authors report melting temperatures of 6357 ± 45 K and 6168 ± 80 K for hcp and bcc, respectively.This study, like ours, also suggest that hcp should be the stable phase of pure iron at core conditions, as this phase has lower free energy and higher melting point than bcc, in agreement with our results.However, the authors of Ref. [57] suggest that PAW-8 underestimates the melting temperatures by more than 600 K compared to PAW-16, while in Ref. [55] the authors suggest that this difference is 400 K. Going beyond free energy perturbation, our TDI calculations using PAW-16 suggest that this difference is actually 466 K, as we obtain a melting temperature of 6060 ± 8 K with PAW-8.This is ∼200 K higher than the melting temperatures reported in Ref. [55] and [57] with the same pseudopotential.Again, the melting temperatures reported in Ref. [57] using PAW-16 are ∼ 200 K lower than our prediction with this pseudopotential.

C. Melting at TPa pressures
Pressures in the cores of Super-Earths can easily reach several TPa [5,8,12,69].Extending the melting curve of iron to such conditions is thus important to determine whether such planets have solid or liquid cores.In Fig. 4  ference, ∆G, between solid and liquid iron at different pressures as a function of temperature.Negative values of ∆G indicate that the liquid phase is more stable than the solid, while ∆G = 0 corresponds to the melting temperature at the given pressure.Since we have determined that the hcp structure has a lower Gibbs free energy at Earth's core conditions, and based on stability arguments of iron phases in Ref. [18], we adopted the hcp structure to perform all other calculations at higher pressures in this work.Fig. 4 shows that ∆G is a sufficiently linear function of temperature.We therefore use a linear interpolation for ∆G(T ) to determine the temperature at which ∆G = 0 for every pressure.The resulting melting points are shown in Fig. 5 along with the isochores we derived using the Z method.Previous ab initio studies of the melting curve of iron are also shown for comparison.First, we notice that the melting temperatures we predict using the Z method are slightly higher, T 1 = 11850 K and T 2 = 25030 K at P 1 = 1008 GPa and P 2 = 4518 GPa, respectively, which agree with our melting curve within an 8% relative error.This also means that iron can withstand significant overheating without melting, as reported in recent shock experiments [70], with critical superheating temperatures 10% to 16% higher than the melting temperatures, but not as high as the critical limit of 23.1% [31].Although the Z method can provide a close upper limit for the melting temperature, it is susceptible to waiting times required for the sample to melt [32], especially in the vicinity of the melting temperature.In addition, it is not clear what the appropriate electronic smearing should be in the microcanonical ensemble simulations of the Z method, which   [18,20,22,26,27,55,58,63] and experiments [23][24][25]30].Isochores of densities 17.44 and 27.23 g cm −3 obtained using the Z method are shown in black filled circles (solid) and open triangles (liquid).Liquids frozen during the simulation are shown in black filled up triangles.Open blue circles correspond to the melting points using the Frenkel correction, which overestimates the melting temperatures.
is particularly important for metals like iron, as the evolution of temperature and pressure of the system depends on this choice.In practice, this means that the ions are not in equilibrium with the electrons, as they follow a Fermi distribution with an associated temperature that is different to that of the ions.If we correct for this effect by considering a smearing of σ = k B T eq , where T eq is the equilibrium temperature that the ions reached in the original Z method calculations, we obtain a melting temperature of 10500 K at 969 GPa for the density of 17.44 g cm −3 , which is in much better agreement with our melting curve that we derived from free energies.Nevertheless, we obtained the specific heat of liquid iron from the slope of the Z curve in energy-temperature space, which yields C liq V = 5.11 k B /atom along the 17.44 g cm −3 isochore that spans through pressures near 1000 GPa, as we see in Fig. 5.This value agrees with the recent experiments by Kraus et al. [24], who estimated a value of 4.2 ± 1.0 k B at similar conditions.At a density of 27.23 g cm −3 , the solid and the liquid have a comparable specific heat of C sol V = 4.94 k B and C liq V = 5.64 k B , respectively.The Z method approach is independent from our TDI method and confirms our melting predictions over a wide pressure interval.
Second, we find that our melting temperatures are higher than the estimations from vibrational frequencies using the Lindemann melting criterion [18].The difference is as large as 5000 K (20%) at 5000 GPa, which shows that iron can withstand significant thermal vibration at high pressure before undergoing melting.At 1500 GPa, our melting line is approximately 750 K higher than the bcc iron melting line of Bouchet et al. [22] who used the two-phase (TP) method to calculate the melting curve up to this pressure.From heat-until-it-melts simulations, the authors conclude that the melting temperature does not significantly depend on whether hcp or bcc crystal structure is used at the pressures of interest and decided to use the bcc structure in all their TP simulations.They acknowledged, however, that calculations with a metastable phase could lead to an underestimation of the melting temperature.While the stability of the bcc structure at Earth's ICB conditions is still a matter of debate [64,66,71,72], the hcp structure has been predicted to be the most stable structure for a wide range of pressures and temperatures [18].
Furthermore, for temperatures as high as 20000 K, we observed that some of our simulation of liquid iron froze spontaneously within a few picoseconds at pressures between 2500 and 5000 GPa.This has also been confirmed by the recent experiments of Kraus et al. [24] that report fast crystallization in the nanosecond time scale.The filled triangles in Fig. 5 mark these conditions.Many of them are at substantially higher temperatures than the extrapolated melting law by Bouchet et al., which underlines that the melting temperature of iron at TPa pressure must be considerably higher than this extrapolation predicts.
At 5 TPa, the Gibbs free energy between the hcp crystal and liquid iron is less than 60 meV/atom at 24000 K and 26000 K, which indicates this temperature is close to the melting line (see Fig. 4).A study of the free energy of iron obtained from phonon calculations in the quasi-harmonic approximation (QHA) [18] predicted an hcp-to-fcc transition at 5.8 TPa, with a steep Clapeyron slope at higher temperatures, implying that the transition should occur at lower pressures for higher temperatures and that both phases should have the same Gibbs free energies at the phase boundary.Therefore, we decided to perform one additional calculation of the Gibbs free energy of iron in the fcc phase at 5 TPa for 27500 K in order to determine whether the fcc phase has a significantly lower Gibbs free energy than the hcp structure, which would have an affect on the predicted melting line.We derived a Gibbs free energy difference of G fcc − G hcp = 56 meV per atom, which means that even at this high temperature the hcp phase is still more stable than the fcc phase, even though predictions from Ref. [18] using QHA indicate that the transition to the fcc phase at 5 TPa should occur at 13300 K. Therefore, we conclude that anharmonic contributions to the free energy play an important role, which is consistent with the underestimated melting temperatures inferred using the Lindemann criterion in Fig. 5 based on QHA.
We also studied finite size effects at this pressure by performing one calculation at T = 30000 K in a larger hcp cell with 180 atoms.The difference in Gibbs free energy between the 144 and 180 atoms cells was found to be 45 meV, which is shown by the yellow triangle in Fig. 4.This means that the effect of size on the melting

H (eV)
FIG. 6. Fractional change in volume, entropy, and enthalpy difference (latent heat) between liquid and solid iron along the melting line as a function of pressure.The black circle and dashed red lines show results from Refs.[73] and [20], respectively.temperature is negligible, as our 144 atoms cell is large enough.
We fitted our melting temperatures from 300-5000 GPa to a Simon equation, with the parameters a = 434.822GPa and c = 1.839 in addition to P 0 = 300 GPa and T 0 = 6469 K.The fractional changes of volumes, entropy, and latent heat of fusion (∆V , ∆S, and ∆H, respectively) along the melting line are shown in Fig. 6 as a function of pressure.At 300 GPa, the volume increase upon melting is 0.116 Å3 per atom, which correspond to a fractional decrease in density of 1.6%.The enthalpy and entropy of melting, respectively, are 1.06 × 10 6 J/kg (0.61 eV/atom) and 1.05 k B /atom at this pressure.These results are an agreement with previous findings at similar conditions [20,73,74].Our entropies of melting, however, decrease with pressure and are slightly higher than the those reported by Kraus et al. [24], who assumed a constant entropy of melting of ∆S = 0.77 k B to derive the melting temperatures between 300 and 1000 GPa.We determined that this value decreased from 1.05 to 0.8 k B /atom over this pressure interval.Novel techniques based on latent heat to detect melting may confirm these values in the near future [75].
As pressure increases, the volume difference between the solid and the liquid decreases and becomes as small as 0.038 Å3 at 5000 GPa, which still represents a fractional difference of 1.21%.Despite this decrease, we always find the solid to be denser than the liquid at the same pressure and temperature, which implies a positive slope for the melting curve.We used these values to calculate the slopes of the melting line from the Clausius-Clapeyron equation dT m /dP = ∆V /∆S.The slopes decrease from 8.6 K/GPa at 300 GPa to 2.6 K/GPa at 5000 GPa, and are consistent with slope of the fitted melting curve in Eq. (9).
In order to study the crystallization of planetary cores, we calculated adiabats of solid and liquid iron.We obtained the Gibbs free The blue line is our fit our melting points according to Eq. ( 9).We also show the core-mantle boundary conditions (CMB) for Uranus (U), Neptune (N), Saturn (blue area), and Jupiter (red area) [9,76] for comparison.The temperature profile of a 10 Earth mass Super-Earth model from Ref. [8] is shown in solid, black line.The 60 GPa isentrope from the rampcompression experiments on iron from Ref. [77] is shown in short dashed line.
energy for both phases along 5 different isotherms, for pressures ranging from 500 to 5000 GPa.The entropies were then obtained from the free energies as S = (E −F )/T , and a spline interpolation was applied to obtain the points of constant entropy.The results are shown in Fig. 7.We observe that the S = 14 k B adiabat intersects the melting line at 1032 GPa and only re-emerges in the solid phase at a much higher pressure 3247 GPa.This implies there exists an extended solid-liquid coexistence region along the adiabat, which has implications for shock and ramp compression experiments that are designed to measure the melting temperature by compressing liquid iron.If ramp waves are employed to avoid shock heating, one can expect such experiments to be nearly adiabatic [24,77,78], such as the recent experiments performed by Kraus et al. [24].Our example illustrate that, up 3247 GPa, such shock experiments would generate a solid-liquid mixtures on the melting line.X-ray diffraction signal would grow with pressure as the solid fraction increases until complete solidification is obtained.This is consistent with the recent measurements by Kraus et al. [24] who observed a solid-liquid mixture when they started the ramp compression from approximately 580 GPa and 8500 K.

IV. RELEVANCE TO THE INTERIOR OF SUPER-EARTHS
Since our melting line is considerably higher than previous predictions, our work will revise some of the assumptions that so far have been made about the temperature distribution in the interiors of Super-Earths.We will now discuss how existing models may have to be adjusted.For planets with 2 Earth masses or more, Gaidos et al. [79] predicted their iron cores to crystallize from the top and suggested iron "snow" near the core-mantle boundary (CMB), which would inhibit convection in their cores.This would shut down the magnetic dynamo and reduce their chances for harboring life [80].Conversely, our ab initio simulations predict the melting line to be steeper than the adiabats up to a pressure of at least 50 Mbar.This means in Super-Earths with up 10.4 Earth masses, the core crystallization will start from the center, like on Earth.We derived this mass estimate by assuming a terrestrial iron fraction of 32.5% and by building interior models following by Seager et al. [5] as well as by Wilson and Militzer [10].
In the following discussion about interior temperature distributions, we will focus on terrestrial planets with up to 1.6 Earth radii and 5.8 Earth masses because observations combined mass-radius models have shown that larger planets must have a low-density outer envelope that is presumably composed of gas or ice [12].To compare with earlier work, we construct models for planets with 2.0 and 5.0 Earth masses.Papuc and Davies [81] reported the average core temperature over an assumed lifespan of 10 billion years.According to their thermal evolution models, average core temperature of a 5.0 (2.0) Earth mass planet would start from 5100 (4300) K and then drop to 4300 (3300) K over a period of 10 Gyr.To compare with these predictions, we constructed a number of interior models based on our computed adiabats for liquid and solid iron and provide some results in Tab.I. Our calculations predict that the core crystallization of a 5.0 (2.0) Earth mass planet would start when average core temperature reaches 12457 (8069) K.The core crystallization would be complete when average temperature reaches 10642 (8069) K.This means our melting calculations imply that the cores of Super-Earths with 2.0-5.0Earth masses would be frozen for their entire 10 Gyr lifespan, if we assume the temperatures predicted by Papuc and Davies are correct.While solid cores are still assumed to cool convectively, there would be no or only a very weak dynamo, as convection in a solid body would produce such small flow speeds that would render the magnetic Reynolds number negligibly small, thus leading to the field decay by magnetic diffusion [82].While we have performed all calculations for pure iron and light elements are known to be present in Earth core, it is reasonable to assume that they would lower the core melting temperature of Super-Earths, but this effect is unlikely to bridge the deviations between our melting predictions and the models by Papuc and Davies that approximately amount to a factor of 2 in temperature.
Noack and Breuer [83] also put together models for the interiors of Super-Earth.For planets with 5 Earth masses, they predict a temperature at the core-mantle boundary (CMB) of 5100 K.According to our calculations, the core crystallization of such planet would set in when the temperature at the CMB reaches 10236 K and already be complete when it reached 8983 K.
Other authors have constructed Super-Earth models with somewhat hotter interiors.Tachinami et al. [84] used mixing length theory to study the thermal evolution of Super-Earth interiors.The temperature profile of 5 Earth mass planet did hardly change over the 10 Gyr lifespan.Respectively, 11300 and 7540 K were predicted for the temperatures in the planet's center and at the CMB.According to our models, core crystallization would start when temperature in the planet's center 15333 K and be already complete when it reaches 12442 K. Wagner et al. [8] also used mixing length theory, and predicted that a 5 Earth masses planet should have a central pressure and temperature of 2000 GPa and 8000 K, respectively.Again these assumed temperatures, imply frozen cores.In Fig. 7, we show the temperature profile of a 10 Earth mass Super-Earth model by the same authors.
Stamenkovic et al. [85] investigated the impact that different viscosity models have on the evolution of Super-Earths.They consider a wide range of initial CMB temperatures from 5100 to 13500 K when they model planets of five Earth masses, which implies some planet models start with frozen cores, some with partially liquid and some with completely liquid cores.If a strongly pressure and temperature dependent viscosity is assumed, it take ∼3 Gyr for the temperature at the CMB to cool from 13500 to 9000 K.At this temperature, our melting calculations again imply that the core of such a planet would be completely frozen.
Recently, Boujibar et al. [16] constructed models for Super-Earth of different masses, core-mantle fractions, and a wide range of interior temperatures in order to determine which planets have partially crystallized cores.Buoyancy effects from core crystallization contribute to convection and the magnetic field generation [86].Results are presented in terms of the temperature at the core-mantle boundary, which is related to the efficiency of retaining the gravitational energy from accretion.We agree with the assumptions and predictions by Boujibar et al. but they based their analysis on the melting line that Morard et al. [63] derived with ab initio simulations and the corresponding fit by Stixrude et al. [58].As we illustrated in Fig. 5, our melting temperatures are similar up to a pressure of 1100 GPa but then are lower for higher pressures.
Currently, the initial temperature profile of a forming planet is not well characterized because no observations have been made.The temperature is controlled by the influx of the planetesimals and radiative energy exchange with the surrounding disk [81,87].Many interior processes, like core-mantle differentiation and radioactive decay matter during the magma ocean phase cooling, are expected to be orders of magnitude faster than after mantle has solidified.Overall, the evolution of terrestrial planets is a complicated process, and as evidenced by the Earth-Venus paradox, it is not a unique process.This is particularly true for rocky exoplanets that have a significant amount of ice, as the presence of rock-ice mixtures can lead to non-layered interiors [88][89][90].With this paper, we aim to improve our understanding of one just element: the state of the iron core.

V. CONCLUSIONS
With DFT-MD simulations, we derived Gibbs free energy of solid and liquid iron and determined its melting curve in the pressure range from 300 to 5000 GPa.We discuss the implications for ramp compression experiments and predict that they generate solid-liquid mixtures over a wide range in pressure.At 5000 GPa, we predict a melting temperature of 25000 K, which is 5000 K higher than predicted from the extrapolation of earlier melting curves.We suggest that the initial temperature profiles in Super-Earth evolution models be investigated in more detail.For some of the published models, our melting results imply that the iron cores start out in solid form and remain frozen for the entire planet's lifetime.In this case, they would not generate a strong magnetic fields but they might still be generated in the mantles of Super-Earths [91,92].On the other hand if the interior temperature profiles are high, the cores of Super-Earths would initially be liquid and would eventually start to crystallize from their centers.Our results also imply, if elemental iron were present in the cores of Jupiter and Saturn, it would occur in solid form (see Fig. 7).However, thermodynamic calculations [45,46] and recent models for Jupiter's interior [93,94] imply that core of giant planets have been eroded and all heavy elements have been mixed with the surrounding metallic hydrogen.6747 K at 330 GPa, which is 688 K higher than the melting temperature we obtain at the same pressure using the PAW-8 pseudopotential.The differences in the predicted melting temperatures between PAW-14 and PAW-16 are small, but very large when compared to predictions from PAW-8 pseudopotential, which shows that the treatment of inner 3p electrons of iron at Earth core pressures and higher is extremely important, as they lead to very different predictions for the melting temperatures of iron at these conditions.We summarize these results in Table III.  is effectively equivalent to applying no correction at all, our melting temperatures became lower than those where the free energy of the solid included the Frenkel correction.We can see the differences in the resulting melting temperatures of iron after including the Frenkel and Navascués correction in Fig. 10.As we can see in Fig. 14, we can already get a well-defined average of these quantities by simulating only 3.5 ps (7000 steps with a timestep of 0.5 fs).Even 2 ps of a DFT-MD simulation are enough to get the energy converged within 1 meV/atom.After generating a force field for solid iron at 290 GPa and 6000 K using this on-the-fly machine learning training process, we were able to simulate 50000 steps in just 12 minutes in one node of 36 cores using the same 144 atoms supercell, which corresponds to a speedup of ×2000.Then, we used the same force field to perform a simulation of iron at the same conditions but in a much larger cell, containing of 1296 atoms (9 × 6 × 6 replicas of an orthorhombic, 4-atoms unit cell of hcp).In both cases, 144 and 1296 atoms, we obtained basically the same mean value of the energy from these simulations, with a difference of less than 16 meV/atom with respect to the smaller 144-atoms simulation, demonstrating convergence with respect to time and system size.The standard deviation decreased considerably when the system size was increased, as expected.FIG.14.Energy and pressure in molecular dynamics simulations of solid iron at 12.938 g cm −3 and 6000 K (∼290 GPa) in the hcp phase.The DFT-MD simulations were performed using the PAW-16 pseudopotential, and the ML potential was trained on these simulations.The inset shows a zoom-in to the fluctuation in a time window of 400 fs.
In Fig. 15, we show how the pressure and energy of iron change as function of the number of iron atoms in the supercell, N .As we can see, a size of 144 atoms is more than enough to converge the energy and pressure with respect to the system size, considering that the associated error in fitting the machine learning potential leads to energy differences about 10 meV/atom.While larger system sizes can be reached by careful supercell design [97], the convergence with respect to system size shown in Fig. 15 demonstrates that 144 is enough to obtain reliable values of pressure and energy.
In addition, we carried simulations at constant pressure in the N P T ensemble to test the stability against fluctuations in the cell volume.Fixing the pressure in our N P T simulations to P 0 = 289.204GPa in the 1296 atoms cell, which was the average pressure in our N V T simulation at 12.938 g cm −3 , resulted in an average density of 12.924 ± 0.0019 g cm −3 , consistent with our N V T simulations.In Fig. 16, we show the fluctuations in the simulation cell dimensions, which show that the simulation cell remained orthorhombic on average.The lattice constants of the hcp lattice
FIG.3.Gibbs free energies per atom of solid (hcp) and liquid iron at P0 = 330 GPa obtained from TDI using the PAW-16 pseudopotential.Top: original energies of each phase.Middle: Gibbs free energy difference ∆G = G liq − G sol .Bottom: We plot the liquid and solid Gibbs free energies with respect to G(P0, T ) ≡ G0 − S0 (T − T0) (with S0 = 13.35kB, T0 = 6000 K, and G0 = 3.939 eV), which is an linear approximation for G liq (P0, T ).TDI calculations of G liq using our fitted pair potential (open red squares) and the WCA potential (open black squares) give consistent values, but the latter approach yields much larger error bars.We compare our results to those from Sun et al.[55].

FIG. 4 .
FIG.4.Gibbs free energy difference between solid and liquid iron at different pressures obtained with the PAW-14 pseudopotential using 144 (except for 108 fcc and 180 hcp).No correction due to a fixed center of mass (Frenkel correction) has been applied to the free energy of the solid.Negative values of ∆G favor melting.Half-filled circle: fcc iron (108 atoms).Yellow-filled triangle: hcp iron (180 atoms).Yellow stars denote the melting temperature at each pressure.

FIG. 7 .
FIG.7.Melting curve and adiabats of solid (dashed black) and liquid (dot-dashed red) iron, where the numbers indicate the respective entropy in kB/atom.The S = 14 kB/atom adiabat is shown in thick, grey line and the open circles indicate the intersection of this adiabat with the melting line.The blue line is our fit our melting points according to Eq. (9).We also show the core-mantle boundary conditions (CMB) for Uranus (U), Neptune (N), Saturn (blue area), and Jupiter (red area)[9,76] for comparison.The temperature profile of a 10 Earth mass Super-Earth model from Ref.[8] is shown in solid, black line.The 60 GPa isentrope from the rampcompression experiments on iron from Ref.[77] is shown in short dashed line.

FIG. 10 .
FIG. 10.Effects of the Frenkel correction of the free energy of an Einstein crystal with a fixed center of mass on the melting temperature.

FIG. 12 .
FIG. 12. Evolution of pressure as a function of time for two typical DFT-MD simulations at low (300 GPa) and high pressures (4000 GPa).Both samples are solid, as depicted by the hcp crystal in the inset.The corresponding average values of pressure are 298.243± 0.155 GPa and 3998.249± 1.071 GPa.

FIG. 15 .
FIG.15.Pressure and internal energy of iron derived from machine learning MD simulations at 6000 K and 12.938 g cm −3 using supercells with different numbers of atoms, N .
[55]4in Ref.[55]), which is 70 meV/atom higher than we obtained with the PAW-16 pseudopotential.Sun et al. reported a smaller value of 53 meV/atom for the Helmholtz free energy difference between PAW-8 and PAW-16 pseudopotentials.
(7)ch does not correct the pressure.Before this correction, Sun et al. reported ∆F Sun WCA→PAW−8 = −3.823eV/atomwithPAW-8pseudopotential(UDFT Ucl (PAW-16, m=1) FIG. 2. Thermodynamic integration (TDI) of liquid iron at 330 GPa and 6400 K.The curves show the integrands of Eqs.(2) and(7).The integrand for our fitted pair potentials is linear function of λ and varies by ∼0.01 eV/atom (lower panel).Conversely, the different integrands for the WCA potential vary by several eV and some show considerable curvature (upper panel).

TABLE I .
Parameters for super-Earths of 2 or 5 Earth masses with terrestrial iron mass fraction of 32.5%.For each mass, we constructed two models: the coldest but still completely liquid core and the hottest but still completely frozen cores.

TABLE III .
Pseudopotentials effects on the melting temperature (no Frenkel correction) at 330 GPa.Our results are compared with findings by Sun et al. if available.