text:
# High Resolution X-Ray Spectroscopy of Seyfert 2 Galaxies

## 1. Introduction

In standard models of active galactic nuclei, Seyfert 2 galaxies are interpreted as systems viewed at sufficiently high inclination angles that the central continuum source is obscured by the surrounding dusty torus (e. g., Krolik 1999). However, residual soft X-ray emission is observed from many of these systems. This must originate in the unobscured nuclear vicinity, or perhaps in other regions of the galaxy. While one expects some soft X-ray emission to be produced by reprocessing of the strong continuum flux in the circumnuclear environment, collisionally-heated gas associated with starburst regions might also provide an important contribution (e. g., Wilson et al. 1992).

High resolution X-ray spectroscopy provides the least ambiguous means of distinguishing between these two interpretations. Discrete spectral features formed in collisional plasmas are qualitatively different from those formed in photoionized gas (Liedahl 1999). In collisional plasmas, emission lines are primarily produced by electron impact excitation from the ground state, which favors strong resonance transitions. Fe L-shell transitions are especially prominent. In photoionized plasmas, the electron temperature and density are too low for excitation by electron impact, and the lines are produced mostly via radiative cascades following recombination, which favors the lowest-lying multiplet levels, thereby producing forbidden transitions. Also, Fe L lines are much weaker than hydrogenic and heliumlike K-shell transitions of lower Z elements. Another characteristic difference involves the radiative recombination continua (RRC). In collisional plasmas, these are broad weak features, barely discernible above the bremsstrahlung continuum. In photoionized plasmas, the RRC are much narrower and far more prominent (Liedahl and Paerels 1996).

The grating experiments on Chandra and XMM-Newton have provided the first X-ray spectra of Seyfert 2’s with sufficient resolution and sensitivity to investigate these issues in detail. Here we present Chandra HETG spectra of Markarian 3 and XMM-Newton RGS spectra of NGC 1068. The data are consistent with all emission arising in photoionized gas. We can use these spectra to derive sensitive constraints on the structure of the circumsource environment.

## 2. Markarian 3

The Chandra HETG spectrum of Markarian 3 is shown in Figure 1 (Sako et al. 2000). A variety of narrow features are readily discernible and are explicitly labeled in the figure. These include Ly series transitions of hydrogenic O, Ne, Mg, Si, S, and Fe, heliumlike complexes of O, Ne, Mg, Si, S, Ar, and Fe, near-neutral fluorescence K$`\alpha `$ lines of Mg, Si, S, and Fe, and a few weak Fe L lines for a range of charge states (Fe XVII through XXIV).

The He$`\alpha `$ complexes are especially interesting. For O VII, the forbidden line is clearly the strongest, and the relative line ratios agree with the values expected for pure recombination (Porquet & Dubau 2000). For both Ne IX and Si XIII, the resonance lines are roughly comparable to the forbidden lines, but the ratios are still not as high as expected for electron impact excitation. This suggests possible contributions from both recombination and collisional excitation emission from the same plasma or spatially-distinct plasmas. The presence of the Fe L features might also bolster that interpretation.

However, we are firmly convinced that this is not what is happening in this spectrum. Both the heliumlike resonance lines and the Fe L lines are not produced collisionally, but instead are photoexcited by the central continuum radiation. Two major clues lead us to this interpretation: (1) The Fe L lines are too weak relative to the heliumlike resonance lines to be consistent with electron impact excitation, unless we significantly reduce the iron abundance in the putative collisional component. (2) The specific Fe L lines observe are all 3d–2p and 3p–2s. These have high oscillator strengths, and are the most prominent features expected for photoexcitation. There are no accompanying 3s– 2p transitions, which have much lower oscillator strengths. In collisional plasmas, the two sets of transitions have near comparable line intensities.

At low to intermediate line optical depths, photoexcitation and recombination provide comparable contributions to the X-ray line spectrum of photoionized gas. This interpretation suffices for Markarian 3. We see no conclusive evidence for a contribution from collisional gas.

## 3. NGC 1068

The XMM-Newton RGS spectrum of NGC 1068 is shown in Figure 2 (Kinkhabwala et al. 2001a). We see a very line-rich spectrum dominated by hydrogenic and heliumlike lines (C through Si), accompanied by weak Fe L-shell transitions. In addition, we see very prominent RRC. All RRC for the hydrogenic and heliumlike charge states of C, N, O, and Ne are detected.

The narrow RRC provide unambiguous evidence that most of the X-ray emission we are observing is produced in photoionized gas. This conclusion is supported by the He$`\alpha `$ complexes, in which the forbidden line is brighter than the resonance. The ratios of the forbidden line intensities to the RRC intensities agree with expectations for recombination cascades.

Nevertheless, as for Markarian 3, recombination alone does not provide a good quantitative fit to the entire NGC 1068 spectrum. The heliumlike resonance lines are still too bright relative to the forbidden lines, and the Fe L lines are brighter than would be expected for purely recombining gas. Here again, we can explain all of these features by allowing for additional contributions from photoexcitation. Conclusive evidence for this comes from the hydrogenic line series: The higher series lines are observed to be significantly brighter relative to Ly$`\alpha `$ than is expected for either recombining gas or collisionally-excited gas. This is a natural consequence of photoexcitation at intermediate optical depths due to saturation of the highest oscillator strength lines.

We have developed a fully self-consistent model for the NGC 1068 spectrum in terms of a photoionized medium which correctly accounts for both the photoexcitation and radiative recombination contributions (Kinkhabwala et al. 2001a, Behar et al. 2001). Individual column densities and velocity distributions for each ion are left as free parameters. The model provides an amazingly good fit to these high statistics data. Fits to the earlier Chandra HETG spectra of Markarian 3, NGC 4151 (Ogle et al. 2000), and the Circinus galaxy (Sambruna et al. 2000) as well as to a recent observation of NGC 4507 yield similar results (Kinkhabwala et al. 2001b). In particular, we find that the Ogle et al. (2000) conclusion that the NGC 4151 spectrum indicates the presence of hot collisional gas is incorrect. All of the features in the spectra of all Seyfert 2 galaxies so far observed with Chandra and XMM are consistent with fluorescence of cold gas and photoexcitation and recombination cascades in photoionized gas.

## 4. Conclusion

The high resolution X-ray spectra of Seyfert 2 galaxies display a wealth of emission features produced by fluorescence, photoexcitation, and recombination cascades in radiatively-driven gas. There is no evidence for measureable contributions from hot collisionally-ionized gas in any of the sources which have been studied to date. Further observations of this kind should be very helpful in elucidating the nature of the circumsource medium in active galactic nuclei.

## References

Behar, E., et al. 2001, in ASP Conf. Ser. Vol. TBD, Mass Outflow in Active Galactic Nuclei: New Perspectives, ed. Crenshaw, D.M., Kraemer, S.B. & George, I.M. (San Francisco: ASP), TBD

Kinkhabwala, A., et al. 2001(a,b), in preparation

Krolik, J. 1999, Active Galactic Nuclei, Princeton University Press

Liedahl, D. A. 1999, in X-ray Spectroscopy in Astrophysics, ed. van Paradijs, J., & Bleeker, J. A. M. (Springer), 189

Liedahl, D. A., & Paerels, F. 1996, ApJ, 468, L33

Ogle, P. M., et al. 2000, ApJ, 545(2), L81

Porquet, D., & Dubau, J. 2000, A&AS, 143, 495

Sako, M., Kahn, S. M., Paerels, F., & Liedahl, D. A. 2000, ApJ, 543, L115

Sambruna, R. M., et al. 2001, ApJ, 546(1), L13

Wilson, A. S., Elvis, M., Lawrence, A., & Bland-Hawthorn, J. 1992, ApJ, 391, L75


format: text

text:
# Direct experimental verification of applicability of single-site model for angle integrated photoemission of small 𝑇_K concentrated Ce compounds

## Abstract

Bulk-sensitive high-resolution Ce 4f spectra have been obtained from 3d $``$ 4f resonance photoemission measurements on La1-xCexAl2 and La1-xCexRu2 for $`x=0.0,0.04,1.0`$. The 4f spectra of low-Kondo-temperature ($`T_\mathrm{K}`$) (La,Ce)Al2 are essentially identical except for a slight increase of the Kondo peak with $`x`$, which is consistent with a known increase of $`T_\mathrm{K}`$ with $`x`$. In contrast, the 4f spectra of high-$`T_\mathrm{K}`$ (La,Ce)Ru2 show a Kondo-like peak and also a 0.5 eV structure which increases strongly with $`x`$. The resonance photon-energy dependences of the two contributions are different and the origin of the 0.5 eV structure is still uncertain.

The single-impurity Anderson model (SIAM) has been used with great success to describe angle integrated 4f photoemission spectra of dense-Kondo Ce compounds, both the low energy Kondo features and the high energy ionization excitations, after proper subtraction of a surface spectrum , in spite of claims to the contrary . Theory support for such claims might come from Nozières ‘exhaustion’ idea, that in the simple Kondo lattice there are not enough conduction electrons to screen out the magnetic moments as in the SIAM, essentially because only electrons located within the Kondo temperature ($`T_\mathrm{K}`$) of the Fermi level ($`E_\mathrm{F}`$) can take part in the Kondo screening . The exhaustion problem is thus most serious in the lowest-$`T_\mathrm{K}`$ systems, and would imply a decreased $`T_\mathrm{K}`$ in the lattice system relative to the impurity system. In this work, we have made a direct experimental comparison of 4f spectra of concentrated and dilute Ce compounds by performing bulk-sensitive high-resolution Ce 3d $``$ 4f resonant photoemission spectroscopy (RPES) of La1-xCexAl2 and La1-xCexRu2, for which the $`T_\mathrm{K}`$’s are 5 K and larger than 1000 K, respectively, in the concentrated limit. They have the same cubic Laves structure in which the number of nearest Ce neighbors of a Ce ion is only four. For $`x=0.04`$, 85% of the Ce ions are completely isolated from other Ce ions if next-nearest Ce-Ce couplings are neglected, so that the 4f spectrum is dominated by the single-site contribution.

Polycrystalline samples of La1-xCexAl2 and La1-xCexRu2 ($`x=0.0,0.04,1.0`$) were prepared by arc melting and characterized by x-ray diffraction. Some elemental Ru impurity was found in the dilute (La,Ce)Ru2 sample. Although the Ru impurity contribution should be removed by the RPES extraction of the Ce 4f spectrum, nonetheless the impurity renders our results for the (La,Ce)Ru2 system to be tentative. Ce 3d $``$ 4f RPES was performed at beam line BL25SU of SPring-8 . All samples were fractured to reveal clean surfaces in a vacuum of $`2\times 10^{10}`$ Torr. The overall energy resolution was about 100 meV and the sample temperature was kept at 20 K.

Figure 1 shows a comparison of the Ce 4f spectra of concentrated and dilute (La,Ce)Al2 obtained from the RPES spectra by a usual method . The sharp two-peak structure near $`E_\mathrm{F}`$ is a generic feature of low-$`T_\mathrm{K}`$ Ce compounds, consisting of the tail of a Kondo resonance (KR) and its spin-orbit replica (SOR) near 0.25 eV. Other structures at higher binding energy are the material-specific ionization features of the hybridization matrix element $`[V(ϵ)]^2`$ between Ce 4f and conduction electrons, within the SIAM description. There is perfect agreement of the two spectra except for an overall reduction of the KR and its SOR in the dilute system. The known reduc-

tion of $`T_\mathrm{K}`$ with decreasing Ce concentration , due to volume expansion, is opposite to expectations from Nozières ‘exhaustion’ idea but easily explains the small spectral change within a SIAM description, including the stronger intensity reduction of the KR relative to its SOR. The results are a direct experimental confirmation of the applicability of a single-site description of the angle integrated Ce 4f spectrum in this low-$`T_\mathrm{K}`$ system, contrary to objections  that have been made.

The Ce 4f spectra of concentrated and dilute (La,Ce)Ru2 are presented in Fig. 2. In great contrast to the (La,Ce)Al2 system, the two spectra differ greatly, especially by a 0.5 eV structure in the concentrated system. At present the origin of the 0.5 eV structure is unclear. The difference might well imply the need to go beyond a single-site model in the high-$`T_\mathrm{K}`$ system due to strong intersite couplings arising from large Ce 4f-Ru 4d hybridization . But this conclusion may be premature. As shown in the inset, the off-resonance spectrum has a strong $`x`$-dependence, that a peaky 0.5 eV structure in the concentrated system is shifted to $`E_\mathrm{F}`$ for the dilute system. Our calculation (not shown) of the 4f spectrum using the Gunnarsson-Schönhammer method  and taking $`[V(ϵ)]^2`$ to be the off-resonance spectrum, does show the 0.5 eV structure in the concentrated system, though with an intensity that is rather smaller than in the experimental spectrum. The fact that the 0.5 eV structure in CeRu2 is found to become much stronger when the spectrum is obtained at a slightly higher photon energy motivates yet another explanation, that it arises from an incoherent Auger transition, as clearly observed in CeFe2 . Further experimental investigation is now in progress.

This work was supported by the U.S. NSF at the University of Michigan (No. DMR-99-71611), by the U.S. DoE at the University of Michigan (No. DE-FG-02-90ER45416), and by the U.S. NSF at the University of California at San Diego (No. DMR-97-05454).


text:
# The transient X–ray pulsar 4U 0115+63 from quiescence to outburst through the centrifugal transition

## 1. Introduction

Accreting collapsed stars in transient X–ray binaries are subject to very large variations of mass inflow rate and provide a laboratory to test the physics of accretion over a range of different regimes that is unaccessible to persistent sources (see e.g. Parmar et al. 1989; Campana et al. 1998). The highly magnetic ($`B10^{12}`$ G), spinning neutron stars that are hosted in Hard X–ray Transients (HXRTs) have long been suspected to be in the so-called propeller regime when quiescent. In this regime accretion onto the neutron star surface is inhibited by the centrifugal action of the rotating magnetosphere (Illarionov & Sunyaev 1975). On the other hand there is little doubt that accretion onto the neutron star surface takes place in these systems while in outburst, as they display very similar properties to those of persistent X–ray pulsars.

Relatively little is known on the low luminosity states of HXRTs. Observations of their quiescent state are still sparse (Campana 1996; Negueruela et al. 2000) and there are only hints that the transition from the lowest outburst luminosities to quiescence occurs in a sudden fashion (see e.g. the case of V 0332+53, Stella et al. 1986). According to models of accretion onto rotating magnetic neutron stars (e.g. Illarionov & Sunyaev 1975), the transition from the accretion to the propeller regime (or vice versa) should take place over a very limited range of mass inflow rates and give rise to luminosity variation in the order $`1000`$ for neutron star spinning at a few second period. The observation of a paroxysmal luminosity variation when a HXRT emerges out of quiescence or fades away at the end of an outburst therefore holds a great potential as a diagnostic of the different regimes experienced by a rotating magnetic neutron star.

In this paper we report on BeppoSAX observations of the transient X–ray pulsar 4U 0115+63, which reveal the quiescent state of the source and, more crucially, provide the first convincing evidence for the centrifugal transition regime of any X–ray transient.

## 2. 4U 0115+63

Like most HXRTs, 4U 0115+63 hosts a magnetic neutron star orbiting a Be star in a moderately eccentric orbit ($`e0.3`$). The spin and orbital periods are $`P_0=3.62`$ s and $`P_{\mathrm{orb}}24.3`$ d (Cominsky et al. 1978; Rappaport et al. 1979). The system has been extensively studied during its frequent outbursts ($`20`$ have been detected so far), which occasionally reach a peak luminosity of $`10^{38}\mathrm{erg}\mathrm{s}^1`$ and usually last about a month (Bildsten et al. 1997; Campana 1996). The neutron star magnetic field is measured with good accuracy, $`B_0=1.3\times 10^{12}[(1+z)/1.2]`$ G, through the cyclotron features detected in its X–ray spectrum (Santangelo et al. 1999 and references therein). Here $`(1+z)=1+GM/(c^2R)1.2`$ (we scale the neutron star mass and radius as $`M=1.4M_{1.4}M_{}`$ and $`R=10R_6`$ km, respectively) is the gravitational redshift at the neutron star surface. The source distance was determined to be $`(8\pm 1)d_8`$ kpc, such that the X–ray flux to luminosity conversion is correspondingly accurate (Negueruela & Okazaki 2001).

Before the present study, 4U 0115+63 was not detected in its quiescent state to a limit of $`6\times 10^{33}\mathrm{erg}\mathrm{s}^1`$ (0.4–6 keV; Campana 1996). At the end of the 1990 outburst a quick decrease of the luminosity was observed below a level of $`5\times 10^{36}\mathrm{erg}\mathrm{s}^1`$, possibly indicating the transition to the onset of the centrifugal barrier (Tamura et al. 1992).

## 3. BeppoSAX Observations

A 35 ks observation of 4U 0115+63 was carried out with the Italian/Dutch satellite BeppoSAX (Boella et al. 1998) on 1999 August 3 06:37:04 UT – August 4 01:41:32 UT, covering an orbital phase range (0.94–0.97) close to periastron. The source count rate varied by a very large amount (from $`\stackrel{>}{}0.02`$ to $`5`$ c s-1) showing a steep increase by a factor of $`\stackrel{>}{}250`$ in $`15`$ hr (see Fig. 1; note that the Y-axis is logarithmic). Superposed to this trend, variations of up to a factor of $`20`$ in less than one hour were clearly present.

Inspection of the RossiXTE ASM light curve of 4U 0115+63 shows that the dramatic flux increase revealed by the BeppoSAX observation was not followed by a source outburst over the following days; the source remained instead below a daily average of $`2\times 10^{36}\mathrm{erg}\mathrm{s}^1`$.

In consideration of the extreme variability, we divided the entire observation in different time intervals (for spectral analysis) as well as intensity intervals (for temporal analysis111Pulsations were detected also at the lowest intensity interval, even though the source spectrum could not be meaningfully characterized; this is why time intervals were used instead in the spectral analysis.), in order to characterize the source properties at different luminosity levels. Energy spectra from the BeppoSAX LECS (0.1–4 keV), MECS (1.8–10 keV) and PDS (15–200 keV) were accumulated in each of these intervals. In all cases the source spectrum from the three instruments was well fit by a model consisting of an absorbed power law with high-energy cut-off plus absorption, with photon index of $`\mathrm{\Gamma }1.0`$, cut-off energy $`E_{\mathrm{cut}}12`$ keV, $`e`$folding energy $`E_{\mathrm{fold}}8`$ keV. The column density amounts to $`N_H1.7\times 10^{22}\mathrm{cm}^2`$ and it is consistent with the galactic value. This model fits also nicely the spectrum of the entire observation and is compatible also with the source spectrum measured during outburst (see e.g. Nagase 1989). The best fit parameters are reported in Table 1 for each interval (see also Gastaldello et al. 2001, in preparation). We note that all spectral parameters are nearly constant throughout the observation (see Table 1). Variations are observed only in the column density ($`4\%`$ probability of getting a higher $`\chi ^2`$ by chance from a constant distribution) and in the power law index (probability of $`1\%`$). The main contribution to the variations of the power law index comes from the first interval. If we exclude this value, we obtain a probability of $`10\%`$. Marginal evidence was found for the second cyclotron harmonic at a centroid energy of $`22.5\pm 2.5`$ keV in the spectrum from the entire observation ($`95\%`$ significance; Santangelo et al. 1999 measured $`24.16\pm 0.07`$ keV during outburst). Note that the first harmonic falls just in the gap between the MECS and the PDS spectra.

In this letter we estimate the source bolometric luminosities by using the 0.1–200 keV unabsorbed luminosities as derived from the best fit parameters of the spectrum from the entire observation and the MECS count rates in the 1.8–10 keV band. The corresponding conversion factor was $`1\mathrm{c}\mathrm{s}^1=8.5\times 10^{35}\mathrm{erg}\mathrm{s}^1`$. The luminosity inferred in this way varied between $`2\times 10^{34}`$ and $`5\times 10^{36}d_8^2\mathrm{erg}\mathrm{s}^1`$ (see Fig. 2).

Pulsations at the $`3.62`$ s neutron star spin period were detected in all intensity intervals. The pulsed fraction (semi-amplitude of modulation divided by the mean source count rate) in the 1.8–10 keV energy band was determined for each interval after folding the data at the best period. Values ranged between $`29\%`$ and $`54\%`$ with the lowest value corresponding to the lowest intensity interval (see Fig. 2).

A further 86 ks BeppoSAX observation was carried out on 2000 August 13 21:33:06 UT – August 16 00:22:36 UT around apastron (orbital phase of 0.42–0.51). The source was very weak and remained at a virtually constant level of $`2\times 10^3`$ c s-1 in the MECS; it was not detected in the LECS ($`<2\times 10^3`$ c s-1). The spectrum was softer than in the August 1999 observation. Fitting a power-law model to the MECS data and fixing the column density to the best fit value of the periastron observation yielded a photon index of $`\mathrm{\Gamma }=2.6_{0.8}^{+1.0}`$ (at $`68\%`$ c.l., $`\chi _{\mathrm{red}}^2=0.2`$ for 2 d.o.f.). The uncertainties in the spectral parameters translate also into a fairly large uncertainty in the inferred 0.1–200 keV unabsorbed luminosity, which is $`(0.63)\times 10^{33}\mathrm{erg}\mathrm{s}^1`$. A pure blackbody model gave also a reasonable fit for a temperature of $`kT=0.7_{0.2}^{+0.3}`$ keV ($`68\%`$ c.l., $`\chi _{\mathrm{red}}^2=0.5`$).

The $`3.62`$ s pulsations were not detected; a $`3\sigma `$ upper limit of $`30\%`$ on the pulsed fraction was derived (see also Campana et al. 2001 in preparation).

In order to facilitate the comparison with the source properties while in outburst we also analyzed the data from the 48 ks BeppoSAX observation that took place on 1999 March 26, during the outburst decay. This observation covered an orbital phase interval of 0.61–0.65. The source light curve is also shown in Fig. 1 for comparison; the average 0.1–200 keV luminosity was $`8\times 10^{37}\mathrm{erg}\mathrm{s}^1`$. The source spectrum could be described by a power law with a high energy cut-off with $`\mathrm{\Gamma }0.8\pm 0.1`$, $`E_{\mathrm{cut}}9.4\pm 0.5`$ keV, $`E_{\mathrm{fold}}16\pm 2`$ keV and $`N_H(1.5\pm 0.1)\times 10^{22}`$ cm-2. To obtain an acceptable fit also an iron line and a cyclotron line feature had to be included.

## 4. Evidence for the transition from the propeller to the accretion regime

The factor of $`\stackrel{>}{}250`$ luminosity variation of 4U 0115+63 during the August 1999 close to periastron observation, when the source luminosity ranged between $`2\times 10^{34}`$ and $`5\times 10^{36}\mathrm{erg}\mathrm{s}^1`$ in $`15`$ hr is the most extreme ever seen in a HXRT. Previously the steepest flux increases detected from this source (and other HXRTs) were those associated to the rise to an outburst peak, involving variations of up to a factor of $`34`$ on a comparable time scale. On the contrary during the August 2000 observation close to apastron the source luminosity was low and nearly constant around a level of $`2\times 10^{33}\mathrm{erg}\mathrm{s}^1`$.

One possibility to explain the extreme variability in the August 1999 observation is that the source was subject to factor of $`\stackrel{>}{}250`$ variations in the mass inflow rate, which in turn gave rise to a comparable variation in accretion luminosity (under the hypothesis that accretion onto the neutron star surface took place unimpeded also for luminosities as low as $`10^{34}\mathrm{erg}\mathrm{s}^1`$). This possibility faces serious problems. Firstly, models of Be star disks/winds predict a neutron star mass capture rate variation of a factor of $`5`$ at the most for a binary system such as 4U 0115+63 over an orbital phase interval of 0.94–0.97 (Raguzova & Lipunov 1998). Moreover, within the direct accretion scenario (required by interpretation above), the issue as to why extreme variability manifests itself over a range of intermediate luminosities (between $`2\times 10^{34}`$ and $`5\times 10^{36}\mathrm{erg}\mathrm{s}^1`$), while the source variability is much less pronounced for both higher and lower luminosities, remains unaddressed.

The behaviour of 4U 0115+63 has a much more natural explanation in terms of the different regimes that are expected for a rotating magnetic neutron star subject to a variable mass capture rate. The inflow of matter towards an accreting magnetic star is dominated by gravitational forces and possibly mediated by an accretion disk down to the magnetospheric boundary at $`r_\mathrm{m}`$. From this radius matter can proceed inwards to the star surface only if the centrifugal force due to magnetospheric drag is weaker than gravity. This translates into the condition that the radius at which a test particle in Keplerian circular orbit corotates with the central object \[the so-called corotation radius, $`r_{\mathrm{cor}}=(GMP^2/4\pi ^2)^{1/3}`$, with $`P`$ the neutron star spin period and $`G`$ the gravitational constant\] is larger than $`r_\mathrm{m}`$ (Illarionov & Sunyaev 1975). In this direct accretion regime, the matter to radiation conversion efficiency is high $`L(R)=GM\dot{M}/R`$ and gives rise to luminous X–ray sources (here $`\dot{M}`$ is the mass inflow rate). If $`r_\mathrm{m}>r_{\mathrm{cor}}`$ the drag by the rotating magnetosphere is super-Keplerian, such that the centrifugal force exceeds gravity and accretion is inhibited. However, as the inflowing matter reaches $`r_\mathrm{m}`$, an accretion luminosity of $`L(r_\mathrm{m})GM\dot{M}/(2r_\mathrm{m})`$ must be released (the factor of 1/2 comes from assumption that the flow down to $`r_\mathrm{m}`$ is mediated by a disk). $`L(r_\mathrm{m})`$ is expected to scale approximately as $`\dot{M}^{9/7}`$. Therefore, above and below the transition to the propeller regime a smooth, close to linear, relationship between the accretion luminosity and $`\dot{M}`$ is expected.

An important prediction is that the accretion luminosity across the transition from direct accretion to the propeller regime (corresponding to $`r_\mathrm{m}r_{\mathrm{cor}}`$) or vice versa should be characterized by a sudden luminosity jump

$`L_{\mathrm{lim}}(R)/L_{\mathrm{lim}}(r_\mathrm{m})`$ $``$ $`2r_{\mathrm{cor}}/R`$
$`2\left({\displaystyle \frac{GMP^2}{4\pi ^2R^3}}\right)^{1/3}`$
This is a factor of $`800`$ in the case of 4U 0115+63 (Corbet 1996; Campana & Stella 2000). In practice the transition separating the two regimes is expected to take place over a finite, though small, range of mass inflow rates around $`\dot{M}_{\mathrm{lim}}`$, such that a very steep dependence of the accretion luminosity on $`\dot{M}`$ ensues temporarily.

Since the magnetospheric radius changes in response to variations of the mass inflow rate, the absolute luminosity at which the centrifugal barrier is expected to close can be determined based on models of the interaction between the inflowing matter and the rotating magnetosphere. In a simple spherical accretion approximation $`r_\mathrm{m}`$ is expected to scale as $`\dot{M}^{2/7}`$ (Davidson & Ostriker 1973). A similar dependence is obtained in the detailed models of the disk/magnetosphere that were developed by (e.g.) Ghosh & Lamb (1979) and Wang (1996) in the regime in which the disk is dominated by gas pressure (the one relevant to the case of 4U 0115+63). In order to account for the predictions of different models, we adopt the scaling above and introduce a correction factor $`\xi `$, the ratio of the magnetospheric radius determined on the basis of a given model to that of simple spherical accretion. We obtain

$$L_{\mathrm{lim}}(R)1\times 10^{37}\xi ^{7/2}B_0^2P_0^{7/3}M_{1.4}^{2/3}R_6^5\mathrm{erg}\mathrm{s}^1$$
(2)
where $`\xi =1`$ for spherical accretion (by definition), $`\xi =0.5`$ in the model by Ghosh & Lamb (1979) and $`\xi 1`$ in the model by Wang (1996).We take this range as representative of the accuracy with which $`r_\mathrm{m}`$ can be predicted by current models ($`M_{1.4}`$ and $`R_6`$ are the neutron star mass and radius in units of 1.4 solar masses and $`10^6`$ cm, respectively. See e.g. Campana et al. 1998). Since the source distance (and, therefore, luminosity) and neutron star magnetic field are fairly accurately measured in the case of 4U 0115+63, we obtain $`L_{\mathrm{lim}}(R)\xi ^{7/2}\mathrm{\hspace{0.17em}10}^{37}\mathrm{erg}\mathrm{s}^1`$; correspondingly $`L_{\mathrm{lim}}(r_\mathrm{m})\xi ^{7/2}\mathrm{\hspace{0.17em}2}\times 10^{34}\mathrm{erg}\mathrm{s}^1`$. The corresponding lines for $`\xi =1`$ (continuous) and $`\xi =0.5`$ (dashed) are shown in Fig. 1. It is immediately apparent that the variations of the August 1999 observation fall just in the luminosity range of the transition between the propeller and the direct accretion regime, where very large luminosity variations are expected in response to modest changes of the mass inflow rate. On the contrary the quiescent state luminosity of the August 2000 observation, during which the source flux was approximately constant, lies in the propeller regime, while the luminosity during the outburst observation on March 1999 is well in the range of the direct accretion regime.

## 5. A Simple Model for the Centrifugal Transition Regime

In this section we develop a simple model for the centrifugal transition regime and compare its predictions to the results from the August 1999 observation of 4U 0115+63. In general, we express the source accretion luminosity in the transition regime as the sum of two contributions: (a) the luminosity of the disk extending down to the magnetospheric boundary, $`L_{\mathrm{disk}}`$; we assume that mass flows through the disk at rate $`\dot{M}`$, which is equal to the rate at which mass is captured at the neutron star accretion radius; (b) the luminosity released within the magnetosphere $`L_{\mathrm{mag}}`$ by the fraction $`f`$ of the mass inflow rate that effectively accretes onto the neutron star surface. We have

$`L`$ $`=`$ $`L_{\mathrm{disk}}+L_{\mathrm{mag}}=L(r_\mathrm{m})+f(L(R)L(r_\mathrm{m}))`$
$`=`$ $`GM\dot{M}[1/2r_\mathrm{m}+f(1/R1/2r_\mathrm{m})]`$
This approximation the direct accretion regime corresponds to $`f=1`$ and the propeller regime to $`f=0`$.

In the case of 4U 0115+63, being $`2r_{\mathrm{cor}}800R`$, $`fL(R)>(1f)L(r_\mathrm{m})`$ for $`f\stackrel{>}{}10^3`$, i.e. the luminosity produced by matter accreting onto the neutron star surface is dominant over most of the centrifugal transition. This has two important implications. Firstly, being dominated by the release of energy within the magnetosphere, the emitted X–ray spectrum should be similar to that observed in the direct accretion regime (i.e. source outbursts) for luminosities $`10^{37}\mathrm{erg}\mathrm{s}^1`$ and remain nearly unchanged across most of the centrifugal transition. We note also that according to models of accretion columns onto magnetic neutron stars, the spectrum is virtually insensitive to accretion rate variations as long as the optical depths remains $`<1`$ (Nagel 1981). Only for luminosities of $`\stackrel{<}{}10^{34}\mathrm{erg}\mathrm{s}^1`$ (corresponding to $`f\stackrel{<}{}10^3`$), when $`L_{\mathrm{disk}}`$ becomes non-negligible, substantial spectral changes are to be expected. The apparent (relative) stability of the BeppoSAX spectra during the August 1999 observation are consistent with this expectation. Secondly, as the transition regime involves only a relatively small variation of $`r_\mathrm{m}`$ the geometry of accretion close to the polar caps should show only minor changes (Wang & Welter 1981; Parmar et al. 1989). Therefore the pulse fraction is expected to remain essentially unaltered as long as $`L_{\mathrm{mag}}`$ dominates. On the other hand, $`L_{\mathrm{disk}}`$ is expected to be unpulsed and the source pulsed fraction should decrease close to bottom of the transition regime. This is also consistent with the results from the August 1999 observation of 4U 0115+63 (see Table 1 and Fig. 2).

We explored a simple model to determine $`f`$ and the source behaviour across the centrifugal transition: we considered a magnetic dipole field the axis of which is tilted relative to the neutron star rotation axis by an angle $`\chi `$. We determined the azimuthal ($`\varphi `$-)dependence of the magnetospheric boundary in the disk plane by equating the ram pressure of radially free-falling matter with the local magnetic pressure. The magnetospheric boundary takes an elongated shape given by

$`r_\mathrm{m}(\dot{M},\varphi )`$ $`=`$ $`\xi B^{4/7}R^{12/7}`$
$`[1+3\mathrm{sin}^2\chi \mathrm{sin}^2\varphi ]^{2/7}(2GM\dot{M}^2)^{1/7}`$
(Jetzer, Strässle & Straumann 1998). The minimum radius $`r_\mathrm{m}(\dot{M},0)`$ corresponds also to $`\chi =0`$, the approximation usually adopted in models of the disk-magnetosphere interaction. The maximum radius $`r_\mathrm{m}(\dot{M},\pi /2)`$ is only a factor of $`[1+3\mathrm{sin}^2\chi ]^{2/7}1.49`$ larger. If $`r_{\mathrm{cor}}<r_\mathrm{m}(\dot{M},0)`$ the magnetospheric boundary is larger than the corotation radius for any $`\varphi `$ and the propeller regime applies (i.e. $`f=0`$). If on the contrary $`r_{\mathrm{cor}}>r_\mathrm{m}(\dot{M},\pi /2)`$ every point on the boundary is within the corotation radius and the standard accretion regime applies ($`f=1`$). In the intermediate regime, $`r_\mathrm{m}(\dot{M},0)<r_{\mathrm{cor}}<r_\mathrm{m}(\dot{M},\pi /2)`$, only the fraction $`f=\mathrm{\Delta }\varphi /2\pi `$ of the magnetospheric boundary for which $`r_{\mathrm{cor}}>r_\mathrm{m}(\dot{M},\varphi )`$ leads to accretion onto the neutron star surface and the transition regime applies. The accretion luminosity versus mass inflow rate curves have been calculated by using this model with $`\xi =1`$. Different curves correspond to different values of the angle $`\chi `$, the sharpest transitions occurring for the lowest values of $`\chi `$, as expected (see Fig. 3).

We adopted the model above also to describe the evolution of the source pulse fraction versus luminosity (see Fig. 2). If $`L_{\mathrm{disk}}`$ is unpulsed, the observed pulsed fraction $`k`$ can be expressed as $`kk_{\mathrm{mag}}f/(R/2r_\mathrm{m}+f)`$, where $`k_{\mathrm{mag}}`$, the pulsed fraction of $`L_{\mathrm{mag}}`$ alone, can effectively be regarded as a free parameter. $`k`$ depends on $`\xi `$ (which sets the onset luminosity of the transition regime) and $`\chi `$ (which determines the range of the mass inflow rates over which the transition regime applies). Here too the dependence on $`\chi `$ is such that for low values the transition is sharp, whereas for large values the transition occurs over a mass inflow rate variation by $`4`$ (see Fig. 3).

Modelling the time evolution of the source luminosity during the August 1999 observation would require knowing the time evolution of $`\dot{M}`$. We approximated $`\dot{M}(t)`$ with the mass capture rate by the neutron star in its orbital motion within the equatorial disk of the Be star companion over the phase interval 0.94–0.97. We assumed a 10 km s-1 nearly Keplerian disk-outflow velocity and a density profile $`r^3`$ (as in Raguzova & Lipunov 1998). This model predicts a $`\dot{M}`$ variation of $`\stackrel{<}{}3`$ over the phase interval covered by the August 1999 observation. We remark that this is among the largest mass capture rate variations that wind models predict. The centrifugal transition model described above with $`\chi =35^\mathrm{o}`$ and $`\xi =1`$, was then used to calculate the line shown in Fig. 2; this reproduces reasonably well the overall behaviour of the light curve. Accordingly a luminosity variation of a factor of $`\stackrel{>}{}250`$ is produced in response to a factor of $`2`$ variation in $`\dot{M}`$: clearly the centrifugal transition works like a very efficient amplifier. Moreover, these values of $`\chi `$ and $`\xi `$ are among those that match nicely the pulsed fraction versus luminosity variation (see Fig. 2). By adopting the same model for the equatorial wind of the Be star over the phase interval of the August 2000 observation (0.42–0.51), the predicted $`\dot{M}`$ is a factor of $`10`$ smaller (a value for which the neutron star is well in the propeller regime) giving rise to a luminosity in good agreement with the observed value.

## 6. Conclusions

The BeppoSAX observations of the transient X–ray pulsar 4U 0115+63 revealed for the first time the presence of the extreme variations as the source approached periastron. The luminosities encompassed by these variations are within the range predicted by modelling of the centrifugal transition, which separates the propeller from the direct accretion regime over a small interval of mass inflow rates. On the contrary during a quiescent state observation the source was most probably in the propeller regime, as its luminosity was lower still and nearly constant (however see below).

Before the present work the evidence for the onset of the centrifugal barrier was based on the sudden steepening of the outburst decay below a luminosity of $`10^{36}\mathrm{erg}\mathrm{s}^1`$ in V0332+53 (Stella et al. 1986) and 4U 0115+63 (Tamura et al. 1992). In both cases, the source became quickly undetected in the relevant (collimator) instruments, such that the transition towards the propeller regime could not be seen. It was also proposed that several other sources in their quiescent (or low) state host a neutron star in the propeller regime: among these are the HXRTs A 0538–66 (Campana 1997; Corbet et al. 1997) and A 0535+26 (Negueruela et al. 2000). The evidence reported here is far more convincing in that dramatic source flux variations of the kind expected in the centrifugal transition regime were observed for the first time. These variations occur within the luminosity interval predicted by models of disk-magnetospheric interaction for the case of 4U 0115+63, where the neutron star spin and magnetic field, as well as distance are fairly accurately measured. The source flux level and absence of sizeable variability in the faint state away from periastron is also in agreement with basic expectations for the regime in which the centrifugal barrier is fully closed. Yet, owing to poor statistics a thermal like spectrum cannot be ruled out. The analogy with neutron star soft X–ray transients suggests that such a spectral component, if present, might be due to reemission of heat in the inner crust caused by pycnonuclear reactions (Brown, Bilsten & Rutledge 1998; Campana et al. 1998; Colpi et al. 2001). The inferred blackbody luminosity222Note that the spectrum emerging from a hydrogen atmosphere would be different from a pure blackbody; yet bolometric corrections are relatively small (factor of 2–3 depending on the spectrum) and are ignored here. is in $`10^{33}\mathrm{erg}\mathrm{s}^1`$ to be compared with an expected deep crustal heating luminosity of $`5\times 10^{33}\mathrm{erg}\mathrm{s}^1`$ under the hypothesis that 4U 0115+63 has accreted at a time average rate of $`5\times 10^{15}\mathrm{g}\mathrm{s}^1`$ (as suggested by the RossiXTE ASM light curves, http://xte.mit.edu/lcextrct/) for some $`10^4`$ years. Independent of the origin of the quiescent emission from 4U 0115+63, the data presented in this paper provide substantial new evidence in favor of the centrifugal transition regime.

This research has made use of SAXDAS linearized and cleaned event files (Rev.2.0) produced at the BeppoSAX Science Data Center. This work was partially supported by ASI, Co-fin and CNAA grants.


text:
# Iron under Earth’s core conditions: Liquid-state thermodynamics and high-pressure melting curve from ab initio calculations

## I Introduction

The last few years have seen important progress in calculating the thermodynamic properties of condensed matter using ab initio techniques based on density-functional theory (DFT) . There has been particular attention to the thermodynamics of crystals, whose harmonic free energy can be obtained from phonon frequencies computed by standard DFT methods . The ab initio treatment of liquid-state thermodynamics is also important, and thermodynamic integration has been shown to be an effective way of calculating the DFT free energy of liquids . These developments have made it possible to treat phase equilibria, including melting properties, by completely ab initio methods. We report here DFT free-energy calculations on high-pressure/high-temperature liquid iron, which we combine with earlier results on the solid  to obtain the complete melting curve and the variation of the volume and entropy of melting along the curve. We also present results for some key thermodynamic properties of the liquid, which we compare with data from shock experiments and other sources. The general methods developed here may be useful for other problems involving phase equilibria under extreme conditions. A brief report of this work was presented earlier .

The properties of high-pressure/high-temperature Fe are of great scientific importance because the Earth’s core consists mainly of Fe, with a minor fraction of light impurities . The melting curve is particularly important, since it provides one of the very few ways of estimating the temperature at the boundary between the liquid outer core and the solid inner core . Because of this, strenuous efforts have been made to measure the melting curve , but the extreme pressures and temperatures required ($`p330`$ GPa, $`T6000`$ K) make the experiments very demanding. Ab initio calculations therefore have a major role to play, and several independent attempts to obtain the melting curve using different ab initio strategies have been reported recently . The rather unsatisfactory agreement between the predictions makes a full presentation of the technical methods all the more important.

The calculation of melting properties using ab initio free energies was pioneered by Sugino and Car  in their work on the melting of Si at ambient pressure. Related methods were subsequently used by de Wijs et al.  to study the melting of Al. In both cases, thermodynamic integration (see e.g. Ref. ) was used to obtain the ab initio free energy from the free energy of a simple reference system, and we follow the same strategy here. The other recent calculations  on the high-pressure melting of Fe employed ab initio methods in a different way. Free energies were not calculated, but instead an empirical parameterised form of the total-energy function was fitted to DFT total energies calculated for representative configurations of the solid and liquid. The empirical energy function was then used in molecular dynamics simulations of very large systems containing coexisting solid and liquid.

The detailed DFT techniques used in this work are identical to those used in our work on h.c.p. Fe . In particular, we use the generalised gradient approximation (GGA) for exchange-correlation energy, in the form known as Perdew-Wang 1991 , which reproduces very accurately a wide range of experimental properties of solid iron, as noted in more detail elsewhere . We also use the projector-augmented-wave (PAW) implementation of DFT , which is an all-electron technique similar to other standard implementations such as full-potential augmented plane waves (FLAPW) , as well as being closely related to the ultrasoft pseudopotential method . We have used the VASP code , which is exceptionally stable an efficient for metals, with the implementation of an extrapolation of the charge density which increases the efficiency of molecular dynamics simulations by almost a factor of two .

The calculation of melting properties demands very high precision for the free energies of the two phases, as emphasised elsewhere . The required precision is set by the value of the entropy of melting, and one finds that in order to calculate the melting temperature to within 100 K the non-cancelling error in the free energies must be reduced to $`10`$ meV/atom. The use of identical electronic-structure methods in the two phases is clearly necessary; but it is certainly not sufficient, since the detailed free-energy techniques differ in the two phases. In the solid, we relied heavily on harmonic calculations, whereas the liquid-state calculations rely on relating the free energy to that of a reference liquid. It is therefore essential to reduce the statistical-mechanical errors below the tolerance, and we aim to demonstrate that this has been achieved.

In the next Section, we summarise the technical methods, and Sec. 3 then reports our results for the DFT free energy of liquid Fe over a wide range of thermodynamic states. Sec. 4 presents our calculated melting properties, which we compare with experimental results and the predictions of other ab initio calculations. Our free-energy results have been used to compute a variety of other thermodynamic quantities for the liquid, and we compare these in Secs. V and VI with direct shock measurements as well as published extrapolations of other experimental data. In the final Section, we give further discussion and a summary of our conclusions.

## II Techniques

The key thermodynamic quantity calculated in this work is the ab initio Helmholtz free energy $`F`$, which, with the statistical mechanics of the nuclei treated in the classical limit, is:

$$F=k_\mathrm{B}T\mathrm{ln}\left\{\frac{1}{N!\mathrm{\Lambda }^{3N}}𝑑𝐑_1\mathrm{}𝑑𝐑_N\mathrm{exp}\left[\beta U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_N;T_{\mathrm{el}})\right]\right\},$$
(1)
where $`𝐑_i`$ ($`i=1,\mathrm{}N`$) are the positions of the $`N`$ nuclei, $`\mathrm{\Lambda }=h/(2\pi Mk_\mathrm{B}T)^{1/2}`$ is the thermal wavelength, with $`M`$ the nuclear mass and $`\beta =1/k_\mathrm{B}T`$. The quantity $`U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_N;T_{\mathrm{el}})`$ is the DFT electronic free energy calculated with the $`N`$ nuclei fixed at positions $`𝐑_1,\mathrm{}𝐑_N`$. This is given, following the Mermin formulation of finite-temperature DFT , by $`U_{\mathrm{AI}}=ETS`$, where the DFT energy $`E`$ is the usual sum of kinetic, electron-nucleus, Hartree and exchange-correlation terms, and $`S`$ is the electronic entropy, given by the independent-electron formula: $`S=k_\mathrm{B}T_{\mathrm{el}}_i[f_i\mathrm{ln}f_i+(1f_i)\mathrm{ln}(1f_i)]`$, with $`f_i`$ the thermal (Fermi-Dirac) occupation number of orbital $`i`$. In exact DFT, the exchange correlation (free) energy $`E_{\mathrm{xc}}`$ has an explicit dependence on $`T_{\mathrm{el}}`$ but we assume here that $`E_{\mathrm{exc}}`$ has its zero-temperature form. It was shown in Ref.  that quantum corrections to the classical approximation to the free energy are negligible in the high-temperature solid, and the same will be true a fortiori in the liquid.

Our earlier work  should be consulted for technical details of the PAW implementation of DFT that we use. We note here that under Earth’s core conditions it is not accurate enough to neglect the response of the $`3p`$ electrons, and in principle these should be explicitly included in the valence set along with the $`3d`$ and the $`4s`$ electrons. However, as shown earlier , the computational effort of including 3$`p`$ electrons explicitly can be avoided with almost no loss of precision by mimicking the effect of the 3$`p`$ electrons by an effective pair potential. The procedure for constructing this pair potential was described in Ref. . The pair potential used here is exactly the same as we used in our thermodynamic calculations on the h.c.p. solid, so that a good cancellation of any residual errors is expected. To avoid any possible doubt on this score, we have also done spot checks on the effect of including both 3$`s`$ and 3$`p`$ electrons explicitly in the valence set, as reported in Sec. IV. The outermost core radius in our PAW calculations is 1.16 Å. At Earth’s core pressures, the atoms in both liquid and solid come closer than the diameter of the ionic cores, which therefore overlap. We will show in Sec. IV that this too has only a very small effect on the free energies.

The present calculations, like those reported earlier on the h.c.p. solid, make use of ‘thermodynamic integration’ , which is a completely general technique for determining the difference of free energies $`F_1F_0`$ of two systems whose total-energy functions are $`U_1`$ and $`U_0`$. The basic idea is that $`F_1F_0`$ represents the reversible work done on continuously and isothermally switching the energy function from $`U_0`$ to $`U_1`$. To do this switching, a continuously variable energy function $`U_\lambda `$ is defined as:

$$U_\lambda =(1\lambda )U_0+\lambda U_1,$$
(2)
so that the energy goes from $`U_0`$ to $`U_1`$ as $`\lambda `$ goes from 0 to 1. In classical statistical mechanics, the work done in an infinitesimal change $`d\lambda `$ is:

$$dF=dU_\lambda /d\lambda _\lambda d\lambda =U_1U_0_\lambda d\lambda ,$$
(3)
where $`_\lambda `$ represents the thermal average evaluated for the system governed by $`U_\lambda `$. It follows that:

$$F_1F_0=_0^1𝑑\lambda U_1U_0_\lambda .$$
(4)
We use this technique to calculate the ab initio free energy $`F_{\mathrm{AI}}`$ of liquid Fe by identifying $`U_1`$ as the ab initio total energy function $`U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_N)`$ and $`U_0`$ as the total energy $`U_{\mathrm{ref}}(𝐑_1,\mathrm{}𝐑_N)`$ of a simple model reference system, whose free energy $`F_{\mathrm{ref}}`$ can be calculated. Then $`F_{\mathrm{AI}}`$ is given by:

$$F_{\mathrm{AI}}=F_{\mathrm{ref}}+_0^1𝑑\lambda U_{\mathrm{AI}}U_{\mathrm{ref}}_\lambda .$$
(5)
In practice, we calculate $`U_{\mathrm{AI}}U_{\mathrm{ref}}_\lambda `$ for a suitable set of $`\lambda `$ values, and perform the integration numerically. The average $`U_{\mathrm{AI}}U_{\mathrm{ref}}_\lambda `$ is evaluated at each $`\lambda `$ using constant-temperature ab initio molecular dynamics, with the time evolution generated by the total-energy function $`U_\lambda =(1\lambda )U_{\mathrm{ref}}+\lambda U_{\mathrm{AI}}`$.

As explained in more detail elsewhere , the computational effort needed to perform the thermodynamic integration is greatly reduced if the fluctuations of the energy difference $`\mathrm{\Delta }UU_{\mathrm{AI}}U_{\mathrm{ref}}`$ are small, for two reasons. First, the amount of sampling needed to calculated $`\mathrm{\Delta }U_\lambda `$ to a given precision is reduced; second, the variation of $`\mathrm{\Delta }U_\lambda `$ as $`\lambda `$ goes from 0 to 1 is reduced. In fact, if the fluctuations $`\delta \mathrm{\Delta }U\mathrm{\Delta }U\mathrm{\Delta }U_{\mathrm{AI}}`$ are small enough, one can neglect this variation and approximate $`F_{\mathrm{AI}}F_{\mathrm{ref}}+\mathrm{\Delta }U_{\mathrm{AI}}`$, with the average taken in the ab-initio ensemble. If this is not good enough, the next approximation is readily shown to be:

$$F_{\mathrm{AI}}F_{\mathrm{ref}}+\mathrm{\Delta }U_{\mathrm{AI}}+\frac{1}{2k_\mathrm{B}T}(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}.$$
(6)
Our task is therefore to search for a model $`U_{\mathrm{ref}}`$ for which the fluctuations $`U_{\mathrm{AI}}U_{\mathrm{ref}}`$ are as small as possible.

The problem of mimicking the fluctuations of ab initio energy $`U_{\mathrm{AI}}`$ in high-$`p`$/high-$`T`$ liquid Fe using a reference system was studied in detail in a recent paper . We showed there that a $`U_{\mathrm{ref}}`$ consisting of a sum of pair potentials:

$$U_{\mathrm{ref}}=U_{\mathrm{th}}+U_{\mathrm{pair}},$$
(7)
in which

$$U_{\mathrm{pair}}=\frac{1}{2}\underset{ij}{}\varphi (𝐑_i𝐑_j),$$
(8)
can be arranged to mimic the fluctuations of $`U_{\mathrm{AI}}`$ very precisely, if we choose $`\varphi (r)`$ to be a repulsive inverse-power potential $`\varphi (r)=B/r^\alpha `$, with suitable values of $`B`$ and $`\alpha `$. In this expression for $`U_{\mathrm{ref}}`$, we have included a term $`U_{\mathrm{th}}`$ which depends on thermodynamic state, but does not depend on the positions $`𝐑_i`$. We define this $`U_{\mathrm{th}}`$ so as to minimize the mean square value of $`\mathrm{\Delta }U`$, which requires the following condition:

$$\mathrm{\Delta }U_{\mathrm{AI}}U_{\mathrm{AI}}U_{\mathrm{pair}}U_{\mathrm{th}}_{\mathrm{AI}}=0.$$
(9)
It might at first seem surprising that such a simple $`U_{\mathrm{ref}}`$ can reproduce $`U_{\mathrm{AI}}`$ accurately, since it does not explicitly represent the metallic bonding due to partial filling of the $`3d`$-band, which is absorbed into the term $`U_{\mathrm{th}}`$. However, the only requirement on $`U_{\mathrm{ref}}`$ is that the fluctuations $`\delta \mathrm{\Delta }U`$ should be small for the given liquid state, and we shall demonstrate below that this is indeed achieved. We comment further on this question in Sec. VII.

The free energy of the reference system can be expressed as:

$$F_{\mathrm{ref}}=U_{\mathrm{th}}+F_{\mathrm{pair}},$$
(10)
where $`F_{\mathrm{pair}}`$ is the free energy associated with $`U_{\mathrm{pair}}`$. An important advantage of our chosen form of $`U_{\mathrm{ref}}`$ is that $`F_{\mathrm{pair}}`$ depends non-trivially only on a single thermodynamic variable, rather than depending separately on temperature $`T`$ and atomic number density $`n`$. Let $`F_{\mathrm{pair}}^xF_{\mathrm{pair}}F_{\mathrm{pg}}`$ be the excess free energy of the reference system, i.e. the difference between $`F_{\mathrm{pair}}`$ and the free energy $`F_{\mathrm{pg}}`$ of the perfect gas at the given $`T`$ and $`n`$. Then the quantity $`f_{\mathrm{pair}}^xF_{\mathrm{pair}}^x/Nk_\mathrm{B}T`$ depends only on the dimensionless thermodynamic parameter $`\zeta `$ defined as:

$$\zeta =Bn^{\alpha /3}/k_\mathrm{B}T.$$
(11)
This simplifies the representation of $`F_{\mathrm{ref}}`$, since we can write:

$$F_{\mathrm{ref}}=U_{\mathrm{th}}+F_{\mathrm{pg}}+Nk_\mathrm{B}Tf_{\mathrm{pair}}^x(\zeta ).$$
(12)
We also expect the representation of $`U_{\mathrm{th}}`$ to be simple. Note first that since $`U_{\mathrm{th}}`$ is defined so that $`\mathrm{\Delta }U_{\mathrm{AI}}=0`$, then from Eqn (9) we must have:

$$U_{\mathrm{th}}=U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}.$$
(13)
Now if the fluctuations of $`U_{\mathrm{AI}}U_{\mathrm{pair}}`$ are indeed small, then the value of $`U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}`$ should be very close to the value of $`U_{\mathrm{AI}}U_{\mathrm{pair}}`$ evaluated at zero temperature with all atoms on the sites of the perfect h.c.p. lattice having the same density as the liquid. Denoting by $`U_{\mathrm{AI}}^0`$ and $`U_{\mathrm{pair}}^0`$ the values of the zero-temperature ab initio and pair-potential energies for the perfect lattice, and defining $`U_{\mathrm{th}}^0U_{\mathrm{AI}}^0U_{\mathrm{pair}}^0`$, we can then write:

$$U_{\mathrm{th}}=U_{\mathrm{th}}^0+\delta U_{\mathrm{th}},$$
(14)
where $`\delta U_{\mathrm{th}}`$ will be a small quantity depending weakly on volume and temperature. The accurate computation and representation of $`U_{\mathrm{AI}}^0`$ were discussed in Ref. , and the accurate computation of $`U_{\mathrm{pair}}^0`$ is clearly trivial, so that the treatment of $`U_{\mathrm{th}}^0`$ is straightforward. The small difference $`\delta U_{\mathrm{th}}U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}U_{\mathrm{th}}^0`$ is evaluated from the AIMD simulations described in Sec. III C.

We conclude this Section by summarising our route to the calculation of the ab initio free energy $`F_{\mathrm{AI}}`$ of the liquid. We employ Eqn (6), ignoring the higher-order fluctuation terms. Recalling that $`\mathrm{\Delta }U_{\mathrm{AI}}=0`$, and combining Eqns (12) and (14), we have:

$`F_{\mathrm{AI}}`$ $``$ $`F_{\mathrm{ref}}+(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}/2k_\mathrm{B}T`$ (15)
$`=`$ $`U_{\mathrm{th}}^0+\delta U_{\mathrm{th}}+F_{\mathrm{pg}}+Nk_\mathrm{B}Tf_{\mathrm{pair}}^x(\zeta )+(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}/2k_\mathrm{B}T.`$ (16)
We now turn to the calculation of the reduced free energy $`f_{\mathrm{pair}}^x(\zeta )`$ of the reference system and the small quantities $`\delta U_{\mathrm{th}}`$ and $`(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}`$. We shall also give evidence that with our chosen reference model the higher-order fluctuation terms omitted from Eqn (16) are indeed negligible.

## III Free energy of the liquid

### A Inverse-power reference system

The PAW calculations used to validate the inverse-power reference system are those reported in Ref. . They consist of a set of AIMD simulations performed at 16 thermodynamic states covering the temperature range $`30008000`$ K and the pressure range $`60390`$ GPa. All the simulations were performed on a 67-atom system using $`\mathrm{\Gamma }`$-point sampling, with a time step of 1 fs. We stress that such a small system with such limited sampling cannot be expected to yield very precise results for thermodynamic quantities, and our only purpose here is to demonstrate the adequacy of the reference system. At each thermodynamic state, the system was equilibrated using the reference system itself, and AIMD data were then accumulated for a time span of 5 ps.

We showed in Ref.  that the inverse-power model, with parameters $`\alpha =5.86`$ and $`B`$ chosen so that for $`r=2`$ Å the potential $`\varphi (r)`$ is 1.95 eV, reproduces very closely the ab initio liquid for the state $`T=4300`$ K, $`\rho =10700`$ kg m-3. We have studied the strength of the $`\delta \mathrm{\Delta }U`$ fluctuations for all 16 thermodynamic states, using exactly the same reference model for all states, and we report in Table I the normalized strength of these fluctuations, which we characterize by the quantity $`\sigma \left[(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}/N\right]^{1/2}`$. Two points should be noted: First, $`\sigma `$ is small, since its typical value of 100 meV is markedly smaller than the typical thermal energies $`k_\mathrm{B}T`$ (258 meV at the lowest temperature of 3000 K). Once $`\sigma `$ is as small as this, little is gained by further improvement of the reference system. Second, $`\sigma `$ does not vary strongly with thermodynamic state, so that the reference system specified by the values of $`\alpha `$ and $`B`$ given above can be used for all the thermodynamic states of interest here.

### B Free energy of reference system

In order to cover the range of thermodynamic states of liquid Fe that interests us, we need accurate values of the excess free energy of the reference system for $`\zeta `$ values going from 2.5 to 5.0. There have been many studies of the thermodynamic properties of inverse-power systems, including one on the free energy of the liquid $`1/r^6`$ system , but since these do not provide what we need we have made our own calculations of $`f_{\mathrm{ref}}^x(\zeta )`$ for the $`1/r^{5.86}`$ case. Our strategy is to start with standard literature values for the excess free energy of the Lennard-Jones (LJ) liquid (we use the results reported in Ref. ), and to use thermodynamic integration to go from the LJ system to the inverse-power system, so that $`U_0`$ and $`U_1`$ in Eqns (23) represent the LJ and inverse-power total energies respectively. In doing this, our target was to keep technical errors small enough so that the final free energy $`F_{\mathrm{ref}}`$ is correct to better than 5 meV/atom.

We note the following technical points. The calculations were done at a standard volume per atom, usually taken to be 8.67 Å3, with the temperature chosen to give the required value of $`\zeta `$. Ewald techniques were used to avoid cutting off the inverse-power potential at any distance – we regarded this as essential, since a cut-off would compromise the scaling properties of the reference system. The classical molecular dynamics simulations used to compute $`U_1U_0_\lambda `$ were done using the constant-temperature technique, with each atom taken to have a mass of 55.86 a.u., and the time-step set equal to 1 fs. For each thermodynamic state, we are free to choose any convenient values for the LJ parameters $`ϵ`$ and $`\sigma `$. Our criterion for choosing these is that the fluctuations of $`U_1U_0`$ should be kept reasonably small, but with the proviso that the initial LJ system must be in the liquid state. In many cases, we have checked for consistency by using different $`ϵ`$ and $`\sigma `$ values. Since we require $`f_{\mathrm{ref}}^x(\zeta )`$ in the thermodynamic limit of infinite system size, we have made careful checks on size effects. Tests on systems containing up to 499 atoms show that size errors in $`f_{\mathrm{ref}}^x(\zeta )`$ are less than 1 meV/atom, and this is small enough to ensure that $`F_{\mathrm{ref}}`$ has a precision of better than 5 meV. Most of this error arises from the error in the literature values of the LJ free energy. As a further check on our techniques, we have done calculations on the $`1/r^6`$ system at selected thermodynamic states, and compared with the free energy results of Laird and Haymet .

We have checked our procedures by repeating most of the calculations using the perfect gas as reference system, so as to be free from possible errors in the free energy of the LJ system. For these calculations we used a different form for $`U_\lambda `$, namely

$$U_\lambda =(1\lambda ^2)U_0+\lambda ^2U_1,$$
(17)
where $`U_1`$ is the total energy of the inverse power system and $`U_0=0`$ is that of the perfect gas. Using this functional form Eq. (4) becomes

$$F_1F_0=_0^1𝑑\lambda \mathrm{\hspace{0.17em}2}\lambda U_1U_0_\lambda .$$
(18)
The advantage of using this different functional form for $`U_\lambda `$ is that the value of the integrand does not need to be computed for $`\lambda =0`$, where the dynamics of the system is that of the perfect gas. In this case, since there are no forces in the system, there is nothing to prevent the atoms from overlapping, and the potential energy $`U_1`$ diverges. Not computing the integrand at $`\lambda =0`$ only partially solves this problem, since for small values of $`\lambda `$ the forces on the atoms are small, the atoms can come close together, and the potential energy $`U_1`$ fluctuates violently. However, we found that by performing long enough simulations, typically 1 ns, we could calculate the integral with an accuracy of ca. 1 meV/atom. These calculations with the perfect-gas reference system give excess free energies of the inverse power system that are systematically 5 meV/atom lower than those obtained using the LJ reference system. Our belief is that the discrepancy arises from a small systematic error in the free energies given in Ref. .

After all these tests, calculations of $`f_{\mathrm{ref}}^x(\zeta )`$ were done at a regularly spaced set of $`\zeta `$ values at intervals of 0.25, and we found that the results could be fitted to the required precision by the following 3rd-degree polynomial:

$$f_{\mathrm{ref}}^x(\zeta )=\underset{i=0}{\overset{3}{}}c_i\zeta ^i.$$
(19)
The values of the coefficients are: $`c_0=1.981,c_1=5.097,c_2=0.1626,c_3=0.009733`$.

### C From reference to full ab initio

To achieve our target precision of 10 meV/atom in the ab initio free energy $`F_{\mathrm{AI}}`$ of the liquid, two sources of error must be studied: system size effects and electronic $`k`$-point sampling. An important point to note is that these errors only affect the small terms $`\delta U_{\mathrm{th}}`$ and $`(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}/2k_\mathrm{B}T`$ in Eqn (16), since $`f_{\mathrm{ref}}^x(\zeta )`$ refers already to the infinite system, and $`k`$-point errors in $`U_{\mathrm{th}}^0`$ are negligible. We also study the validity of neglecting the higher-order fluctuation terms in Eqn (16).

We focus first on the quantity $`\delta U_{\mathrm{th}}`$ in Eqn (16). To study size errors in this quantity, we calculated the thermal average $`U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}`$ for a range of system sizes. These test calculations were done on systems of up to 241 atoms at $`V=8.67`$ Å3/atom and $`T=4300`$ K using $`\mathrm{\Gamma }`$-point sampling. The preparation and equilibration of these systems were done using the inverse-power reference system. Since the latter so closely mimics the ab initio system for the 67-atom cell, it should provide a well equilibrated starting point for ab initio simulation of larger systems. The duration of all the ab initio simulations after equilibration was 1 ps. The results of these tests are summarised in Table II, where we report the value of $`\delta U_{\mathrm{th}}`$ per atom, i.e. the quantity $`\delta U_{\mathrm{th}}/N[U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}U_{\mathrm{th}}^0]/N`$ (see Sec. II). Since $`U_{\mathrm{th}}^0/N`$ is independent of the system size, the variation of the reported quantity arises solely from size dependence of $`U_{\mathrm{AI}}U_{\mathrm{pair}}_{\mathrm{AI}}/N`$. We see that with $`125`$ atoms $`\delta U_{\mathrm{th}}/N`$ is converged to better than 5 meV/atom, and that already with 67 atoms the size error is of the order of 10 meV/atom.

We tested for $`k`$-point errors in $`\delta U_{\mathrm{th}}`$ by performing calculations using both four and 32 Monkhorst-Pack  sampling points. Since explicit AIMD calculations with so many $`k`$-points would be extremely expensive, we use the following procedure. From an existing $`\mathrm{\Gamma }`$-point simulation we take a set of typically 10 atomic configurations separated by 0.1 ps. The ab initio total energies of these configurations calculated with the different $`k`$-point samplings are then compared. For sampling with four $`k`$-points, we did calculations on systems of up to 241 atoms, but the heavier calculations with 32 $`k`$-points were done only on the 67-atom system. The results of these tests for the thermodynamic state $`V=8.67`$ Å3/atom and $`T=4300`$ K are also reported in table II, where we see that for the smallest system containing 67 atoms the difference with respect to a calculation with the $`\mathrm{\Gamma }`$-point only is $`9`$ meV/atom, but as the number of atoms in the cell is increased above $`125`$ the difference becomes negligible. The result for the calculation with 32 $`k`$-points is identical to the one with four $`k`$-points and is not reported in the table. We also found that fluctuations of the energy differences between the calculations done with the $`\mathrm{\Gamma }`$-point only and those with 4 $`k`$-points are extremely small.

Similar, but less extensive, tests of system-size and $`k`$-point errors have also been performed at the state $`V=6.97`$ Å3/atom, $`T=6000`$ K, and we find that the variation of these errors with system size is numerically almost the same as before. The indication is therefore that $`\delta U_{\mathrm{th}}`$ can be obtained to a precision of ca. 5 meV/atom from simulations on systems of 125 atoms or more. Unfortunately, it is not practicable yet to do all our AIMD simulations with this system size, and in practice we have computed $`\delta U_{\mathrm{th}}`$ from $`\mathrm{\Gamma }`$-point simulations on the 67-atom system, and corrected the results by adding 10 meV/atom, which from the present evidence appears to the almost constant error in the $`\mathrm{\Gamma }`$-point 67-atom results.

As expected, the numerical values of $`\delta U_{\mathrm{th}}`$ are small, and depend weakly on temperature and pressure across the range of thermodynamic states of interest. We find that they can be represented to within $`3`$ meV/atom by a sum of third-degree polynomials in $`V`$ and $`T`$:

$$\delta U_{\mathrm{th}}/N=\underset{i=0}{\overset{3}{}}\left(a_iV^i+b_iT^i\right),$$
(20)
with the following fitting parameters (units of eV, Å and K): $`a_0=0.649`$; $`a_1=4.33\times 10^2`$; $`a_2=4.19\times 10^3`$; $`a_3=6.48\times 10^5`$; $`b_0=0.296`$; $`b_1=6.51\times 10^5`$; $`b_2=7.46\times 10^9`$; $`b_3=2.07\times 10^{13}`$.

To test the validity of neglecting the higher-order fluctuation terms omitted from Eqn (16), we have performed full thermodynamic integration for four different thermodynamic states, the first three with $`V=8.67`$ Å3/atom and $`T=4300`$, 6000 and 8000 K and the fourth with $`V=6.97`$ Å3/atom and $`T=8000`$ K, using the five equally spaced $`\lambda `$ values 0, 0.25, 0.5, 0.75 and 1.0. These calculations were done using $`\mathrm{\Gamma }`$-point sampling on the system of 67 atoms. We have seen that this system size is not big enough to yield the required precision for $`F_{\mathrm{AI}}`$, but it should certainly be enough to test the adequacy of the second-order formula. In Table III we report a comparison between the results obtained from the integral using the five $`\lambda `$ values and those from the second order formula, and we see that they are practically indistinguishable. The Table also indicates that the term $`(\delta \mathrm{\Delta }U)^2_{\mathrm{AI}}/2k_\mathrm{B}T`$ is rather insensitive to thermodynamic state and can be approximated to the required precision by setting it equal to 10 meV/atom. We have used this constant value in evaluating the ab initio free energy by Eqn (16).

## IV Melting properties

From our parameterized formulas for the ab initio Helmholtz free energies $`F(V,T)`$ of the h.c.p. solid (Ref. ) and the liquid (present work), we immediately obtain the Gibbs free energies $`G(p,T)F(V,T)+pV`$, and for each pressure the melting temperature $`T_m`$ is determined as the $`T`$ at which the latter free energies are equal for the solid and liquid. The resulting melting curve is reported in Fig. 1 for pressures from 50 to 350 GPa. On the same plot, we show the ab initio melting curve reported very recently by Laio et al. . We also compare with experimental melting curves or points obtained by shock experiments or by static-compression using the diamond anvil cell (DAC). DAC determinations of the melting curve of Fe and other transition metals have been performed by several research groups . The early results of Williams et al.  lie considerably above those of other groups, and are now generally discounted. This still leaves a range of ca. 400 K in the experimental $`T_\mathrm{m}`$ at 100 GPa. Even allowing for this uncertainty, we acknowledge that our melting curve lies appreciably above the surviving DAC curves, with our $`T_\mathrm{m}`$ being above that of Shen et al.  by ca. 400 K at 100 GPa. We return to this discrepancy below.

Shock measurements should in principle be able to fix a point on the high-pressure melting curve at the thermodynamic state where melting first occurs on the Hugoniot. However, temperature is notoriously difficult to measure in shock experiments. The temperatures obtained by Yoo et al.  using pyrometric techniques are generally regarded as being too high by at least 1000 K. This has been confirmed by our recent ab initio calculations  of Hugoniot temperature for h.c.p. Fe. We therefore disregard their data point on the melting curve. In the shock measurements of Brown and McQueen  and Nguyen and Holmes , no attempt was made to measure temperature, which was estimated using models for the specific heat and Grüneisen parameter; the approximate validity of these models is supported by our ab initio calculations  on h.c.p. Fe. However, the identification of the Hugoniot melting point has been hampered by the possible existence of a solid-solid transition. In their measurements of sound velocity on the Hugoniot, Brown and McQueen  believed that they had observed a solid-solid transition as well as a separate melting transition. The new shock results of Nguyen and Holmes  using improved techniques indicate that there is no solid-solid transition, and we place greater weight on their Hugoniot melting point. We plot in Fig. 1 the point reported by Brown and McQueen  as lying on the melting curve, though for the reasons just explained, we are cautious about accepting it. We also plot the point obtained from the measurements of Nguyen and Holmes . The pressure of 221 GPa is taken directly from their measurement of the onset of melting, while the temperature at this point is taken from our calculation of the Hugoniot temperature of the h.c.p. solid at this pressure, as reported in Ref.  (see also following section).

We now consider possible sources of error in our DFT calculations. First, we recall that even with the best available GGA for exchange-correlation energy the low-temperature $`p(V)`$ relation for h.c.p. Fe is not in perfect agreement with experiment. This has been shown by a number of independent calculations using all-electron techniques  as well as pseudopotential  and PAW  techniques, all of which agree closely with each other. Roughly speaking, the pressure is underpredicted by ca. 10 GPa at near-ambient pressures and by ca. 8 GPa in the region of 300 GPa. The pressure error can be thought of as arising from an error in the Helmholtz free energy, so that the true free energy $`F_{\mathrm{true}}`$ can be written as $`F_{\mathrm{true}}=F_{\mathrm{GGA}}+\delta F`$, where $`F_{\mathrm{GGA}}`$ is our calculated free energy and $`\delta F`$ is the correction. If we take the pressure error $`\delta p(\delta F/V)_T`$ to be linear in the volume, then $`\delta F`$ can be represented as $`\delta F=b_1V+b_2V^2`$, where $`b_1`$ and $`b_2`$ are adjustable parameters determined by least-squares fitting to the experimental pressure. If we now neglect the temperature dependence of $`\delta F`$, and simply add $`\delta F(V)`$ to the calculated free energies of solid and liquid, this gives a way of gauging our likely errors. We find that this free-energy correction leads to a lowering of the melting curve by ca. 350 K in the region of 50 GPa and by ca. 70 K in the region of 300 GPa.

The second error source we consider is the PAW implementation, and specifically our choice of the division into core and valence states, and the PAW core radii. As mentioned earlier, at Earth’s core pressures the $`3p`$ electrons, and to a lesser extent the $`3s`$ electrons, must be treated as valence states. Moreover, the choice of the maximum PAW core radius may also affect the calculations, because under such high pressures and temperatures the atoms come so close that the cores overlap. These errors may affect the melting curve if they fail to cancel between the liquid and the solid. To check both these possible problems, we have performed trial PAW calculations with the much smaller core radius of 0.85 Å and with both $`3s`$ and $`3p`$ states in the valence set; with this choice of core radius the overlap of the cores in the liquid and the high temperature solid is almost negligible. We have then used Eq. (6) to calculate the free energy difference between the systems described with the two PAW approximations, repeating the calculations for both the liquid and the solid. To do that we have drawn two sets of 30 statistically independent configurations from two long simulations performed with the original PAW approximation on the solid and the liquid at $`V=7.18`$ Å3/atom and $`T=6700`$ K. As expected, we find a significant shift in the total electronic (free) energies. This shift is almost constant, thus validating the use of Eq. (6), but the important result is that it is almost the same for the liquid and the solid, the two numbers being $`F_{\mathrm{hard}}^lF_{\mathrm{soft}}^l=0.210`$ eV/atom $`F_{\mathrm{hard}}^sF_{\mathrm{soft}}^s=0.204`$ eV/atom. Here, $`F_{\mathrm{hard}}^l`$ is the free energy calculated with small core and $`3s`$ and $`3p`$ states in valence, and $`F_{\mathrm{soft}}^l`$ the free energy with large core and the $`3s`$ and $`3p`$ frozen in the core, plus the effective pair-potential; the superscripts $`s`$ and $`l`$ indicate the solid and the liquid respectively. The effect is small, and stabilises the liquid by 6 meV/atom, which has the effect of shifting the melting curve down by $`60`$ K.

As we show in Fig. 1, if we include both these corrections they bring our low-temperature melting curve into quite respectable agreement with the DAC measurements of Shen et al., while leaving the agreement with the shock points of Nguyen and Holmes essentially unaffected. There is still a considerable discrepancy with the DAC curve of Boehler  and the ab initio results of Laio et al.

We now turn to the changes of volume and entropy on melting. Our calculated volume of melting (volume of liquid minus volume of coexisting h.c.p. solid at each pressure expressed as a percentage of the volume of the solid at that point) is plotted as a function of pressure in Fig. 2. We also show the melting volume predicted by the ab initio calculations of Laio et al.  at the pressure 330 GPa, and it is encouraging to note that their value of 1.6 % is quite close to ours. The free-energy correction discussed above makes only a small difference to the calculated volume of melting: at 50 GPa the correction makes the volume of melting increase from 5.0 to 5.8 %, while at 300 GPa it is affected by less than 0.1 %. The most striking feature of our results is the steep decrease of $`\mathrm{\Delta }V`$ by a factor of about three in the range from 50 to 200 GPa, and its approximate constancy after that.

Our predicted entropy of melting $`\mathrm{\Delta }S_m`$ (entropy per atom of liquid minus entropy per atom of coexisting solid) is plotted as a function of pressure in Fig. 3, where we also show the ab initio value of Laio et al.  at 330 GPa. The agreement of our value (1.05 $`k_\mathrm{B}`$) with theirs (0.86 $`k_\mathrm{B}`$) is reasonably close. The entropy of melting also decreases with increasing $`p`$, but more moderately than $`\mathrm{\Delta }V/V`$, the decrease between 50 and 200 GPa being only 30 %. We note the relevance to the slope of the melting curve, given by the Clausius-Clapeyron relation: $`dT_m/dp=\mathrm{\Delta }V/\mathrm{\Delta }S`$. (This relation is satisfied identically by our results, since they are all derived from free energies.) The strong decrease of $`dT_m/dp`$ between 50 and 200 GPa and its approximate constancy thereafter is mainly due to the variation of $`\mathrm{\Delta }V/V`$.

## V Hugoniot properties

Since shock experiments are the only direct way of obtaining thermodynamic information for high-$`p`$/high-$`T`$ liquid Fe, it is important to test our predictions against the available shock data. The data that emerge most directly from shock experiments consist of a relation between the pressure $`p_\mathrm{H}`$ and the molar volume $`V_\mathrm{H}`$ on the so-called Hugoniot line, which is the set of thermodynamic states given by the Rankine-Hugoniot formula :

$$\frac{1}{2}p_\mathrm{H}(V_0V_\mathrm{H})=E_\mathrm{H}E_0,$$
(21)
where $`E_\mathrm{H}`$ is the molar internal energy behind the shock front, and $`E_0`$ and $`V_0`$ are the molar internal energy and volume in the zero-pressure state ahead of the front. The pressure-volume and temperature-pressure relations on the Hugoniot are straightforwardly obtained from our ab initio calculations: for a given $`V_\mathrm{H}`$, one seeks the temperature $`T_\mathrm{H}`$ at which the Rankine-Hugoniot relation is satisfied, and from this one obtains $`p_\mathrm{H}`$ (and, if required, $`E_\mathrm{H}`$). In experiments on Fe, $`V_0`$ and $`E_0`$ refer to the zero-pressure b.c.c. crystal. We obtain $`E_0`$ directly from GGA calculations that we performed on the ferromagnetic b.c.c. crystal, as described earlier , but we use the experimental value of $`V_0`$. The slight shift produced by using instead the theoretical value of $`V_0`$ was noted earlier . Melting in shock experiments is usually detected by monitoring the sound velocity , which shows marked discontinuities of slope along the Hugoniot. In a simple melting transition, there are discontinuities at two characteristic pressures $`p_s`$ and $`p_l`$, which are the points where the solid and liquid Hugoniots meet the melting curve. Below $`p_s`$, the material behind the shock front is entirely solid, while above $`p_l`$ it is entirely liquid; between $`p_s`$ and $`p_l`$, the material is a two-phase mixture.

We present in Fig. 4 our calculated $`T_\mathrm{H}(p_\mathrm{H})`$ Hugoniot curve for the liquid, together with our curve for the solid reported earlier  and our ab initio melting curve. Without the free energy correction $`\delta F`$ of Sec. IV, we find $`p_s=`$ 229 and $`p_l=`$ 285 GPa. The very recent shock data of Nguyen and Holmes  give values of 221 and 260 GPa respectively, so that our $`p_s`$ value is very close to theirs, and our $`p_l`$ value is also not very different.

If the correction $`\delta F`$ is included in calculating the melting curve, then for consistency it must be included also in the solid and liquid Hugoniots. It is straightforward to obtain the corrected $`p_\mathrm{H}`$ and $`E_\mathrm{H}`$ as a function of $`V_\mathrm{H}`$ for the two phases. But in the Rankine-Hugoniot equation we also need $`E_0`$ for the b.c.c. crystal, and this will be subject to a correction similar to $`\delta F`$, but of unknown size. To supply the missing information, we add to the ab initio energy of b.c.c. Fe a correction term $`\delta F_{\mathrm{bcc}}`$, which we represent as $`c_1+c_2V`$. The constants $`c_1`$ and $`c_2`$ are fixed by requiring that the equilibrium volume of the b.c.c. crystal and the low-temperature transition pressure between the b.c.c. and h.c.p. phases be correctly given. The resulting ‘corrected’ $`T_\mathrm{H}(p_\mathrm{H})`$ Hugoniots of the solid and liquid are reported in Fig. 4. The shifts in the curves are of about the same size as those discussed earlier  for the solid when we replaced the calculated b.c.c. volume $`V_0`$ in the Rankine-Hugoniot equation by its experimental value, and are an indication of the inherent uncertainty due to DFT errors. The corrected values $`p_s=`$ 243, $`p_l=`$ 298 GPa are now in somewhat poorer agreement with the experimental values. This is a rather sensitive test of DFT errors, since the shallow angle at which the Hugoniot curves cross the melting line amplifies the effect of the errors.

We now turn to our liquid-state results for $`p_\mathrm{H}(V_\mathrm{H})`$ compared with the shock data of Brown and McQueen  (Fig. 5), including for completeness our results for the solid reported earlier . We report results both with and without the free energy correction $`\delta F`$, using the experimental b.c.c. volume $`V_0`$ in the Rankine-Hugoniot equation in both cases. We mark on the Figure the volumes above which the shocked material is entirely solid and below which it is entirely liquid. Above the upper volume, we report our calculated h.c.p. Hugoniots, and below the lower volume the liquid Hugoniots. In the interval between them, we linearly mix the two. We note that the $`\delta F`$ correction makes little difference to the liquid Hugoniot, which lies above the experimental values by ca. 3 %.

Shock experiments on Fe have given values for the adiabatic sound speed $`v_S=(K_S/\rho )^{1/2}`$ of the liquid, with $`K_S`$ the adiabatic bulk modulus and $`\rho `$ the mass density. Fig. 6 shows our ab initio values for $`v_S`$ of the liquid as a function of pressure on the Hugoniot, both with and without the $`\delta F`$ correction, compared with the shock data of Refs. . Up to the pressure of $`260`$ GPa, the experimental points refer to the solid or the two-phase region, so it is in the liquid region above this pressure that the comparison is significant. In that region, our agreement with the experimental data is close, the discrepancies being $`2`$ and $`<1`$ % for our uncorrected and corrected $`v_S`$ values respectively.

We conclude this Section by reporting results for the Grüneisen parameter $`\gamma `$ on the liquid Hugoniot. This parameter is defined as $`\gamma V(p/E)_V=\alpha K_TV_\mathrm{m}/C_v`$, with $`\alpha `$ the volume expansion coefficient, $`K_T`$ the isothermal bulk modulus, $`C_v`$ the constant-volume molar specific heat, and $`V_\mathrm{m}`$ the molar volume. Assumptions or estimates of its values have played a key role in constructing parameterised equations of state for Fe. Our calculated $`\gamma `$ on the liquid Hugoniot is almost exactly constant, varying in the narrow range from 1.51 to 1.52 as $`p`$ goes from 280 to 340 GPa.

## VI Thermodynamics of the liquid

Although directly measured data on high-$`p`$/high-$`T`$ liquid Fe all come from shock experiments, attempts have been made to combine these data with measurements at lower $`p`$ and $`T`$ using parameterised models for quantities such as $`K_S`$, $`\gamma `$ and $`C_v`$ to estimate thermodynamic properties away from the Hugoniot curve . These attempts have been crucial in trying to understand how the properties of the Earth’s liquid core deviate from those of pure liquid Fe. We present here a brief comparison with these experimentally based extrapolations for the two quantities that determine the seismic properties of the outer core: the density $`\rho `$ and the adiabatic bulk modulus $`K_S`$.

Since the outer core is in a state of turbulent convection, the variation of its thermodynamic properties with depth is expected to follow an adiabat. We therefore present our comparisons on adiabats specified by their temperature $`T_{\mathrm{ICB}}`$ at $`p=330`$ GPa, which is the pressure at the inner-core/outer-core boundary (ICB) . We choose the two temperatures $`T_{\mathrm{ICB}}=5000`$ and 7000 K, because the results of Sec. IV indicate that the melting temperature at the ICB pressure lies between these limits. Our comparisons (Tables IV and V) show that the uncorrected ab initio density is very close (within a few tenths of a percent) to the extrapolated experimental data at $`p=150`$ GPa, and in slightly poorer agreement (within $`1.5`$ %) at $`p=350`$ GPa. As expected, the free-energy correction lowers the predicted density, resulting in larger discrepancies with experiment of 1.5 % and 2.5 % at $`p=150`$ and 350 GPa respectively. Our uncorrected ab initio $`K_S`$ values also agree more closely with the experimental data, being typically within 2 %, while the corrected predictions disagree with the data by up to 8 %. However, given the closer agreement between ab initio and experiment on the Hugoniot (Sec. V), it is possible that some of the disagreements may be due to deficiencies in the experimental extrapolation.

## VII Discussion and conclusions

In assessing the reliability of our results, we consider three sources of error: first, the uncontrolled DFT errors inherent in the GGA for exchange and correlation energy; second, the controllable errors in the detailed electronic-structure implementation of GGA, and specifically in the use of PAW to calculate the total ab initio (free) energy $`U_{\mathrm{AI}}(𝐑_1,\mathrm{}𝐑_N;T_{\mathrm{el}})`$ for each set of atomic positions $`𝐑_1,\mathrm{}𝐑_N`$; third, the statistical-mechanical errors, including system-size effects. We have endeavoured to reduce errors of the third kind below 10 meV/atom for the liquid. In our earlier free-energy calculations on the h.c.p. solid , the corresponding error was estimated as $`15`$ meV/atom. Taking these errors together, and recalling that the resulting error in melting temperature $`T_m`$ is roughly the combined free-energy error divided by Boltzmann’s constant, we find an expected $`T_m`$ error of ca. $`\pm 300`$ K. We have also attempted to control errors of the second kind by changing the division between core and valence states and by reducing the core radius. These tests suggest that the associated error in $`T_m`$ is probably no more than ca. $`\pm 100`$ K. The inherent DFT errors are more difficult to quantify, but we have demonstrated that the known discrepancies in the low-temperature $`p(V)`$ relation for h.c.p. Fe almost certainly lead to an overestimate of $`T_m`$ by ca. 350 K at 50 GPa and ca. 70 K at 300 GPa, and we have corrected for this. We have also seen the significant shifts in the Hugoniot curves resulting from DFT errors. We believe the remaining uncertainty in $`T_m`$ from this source could be as much as 300 K.

Our attempts to correct for DFT errors give a melting curve which is in quite good agreement with the recent measurements of Shen et al. , and with estimates based on shock data ; the methods used to estimate temperature in the shock experiments are also supported by our ab initio results for Grüneisen parameter $`\gamma `$ and specific heat $`C_v`$ . Our melting curve is still above the experimental data of Boehler  by $`800`$ K in the pressure region up to ca.  100 GPa. We cannot rule out the possibility that some of this discrepancy is due to our DFT errors. Our substantial disagreement with the ab initio melting curve of Laio et al.  must be due to other reasons. We are currently working with authors of Ref.  to discover the cause of the disagreement, and we hope to report on this in the future.

A key part of our strategy for eliminating system-size errors in the calculated free energies is the use of an empirical reference model which accurately reproduces the fluctuations of total energy. At first sight, the use of a reference model based on a purely repulsive pair potential might seem surprising, since it does not explicitly include a description of metallic bonding. An empirical reference model (there called an ‘optimised potential model’) is also used in the work of Laio et al. , though they use it in a different way from us. Their optimised potential model is a form of the ‘embedded atom model’ (EAM) , which explicitly includes metallic bonding. As described in our earlier work , we have investigated the consequences of using the EAM as a reference model. We showed there that for present purposes fluctuations of the bonding energy are negligible, and that under these circumstances the EAM is almost exactly equivalent to a model based on repulsive pair potentials. We also showed that there is no numerical advantage in using the EAM as reference model for the calculation of free energies. The use of different reference models per se therefore appears to have nothing to do with the current disagreement between ab initio melting curves.

The agreement of our ab initio results with the limited data from shock experiments on the liquid is reasonably satisfactory. In particular, our predicted Hugoniot relation $`p_\mathrm{H}(V_\mathrm{H})`$ is almost as good as we found earlier for the solid. The adiabatic sound velocity of the liquid is also predicted to within $`13`$ %, the discrepancy depending on whether or not we attempt to correct for DFT errors. The good agreement for the Grüneisen parameter $`\gamma `$ is also encouraging. Our results for the h.c.p. solid  indicated that $`\gamma `$ varies little with pressure or temperature for $`100<p<300`$ GPa and $`4000<T<6000`$ K, and has a value of ca. 1.5. Our present results indicate that the same is true of the liquid.

The ab initio free-energy techniques outlined here could clearly be adapted to a wide range of other problems, so that melting curves could be calculated for many materials, including those of geological interest, like silicates. We have recently completed ab initio calculations of the melting curve of aluminium up to pressures of 150 GPa, which are in excellent agreement with static-compression and shock data, as will be reported elsewhere .

In conclusion, we have shown how ab initio free-energy calculations based on thermodynamic integration can be used to obtain the melting curve and the volume and entropy of melting of a material over a wide pressure range. We have emphasised that the key requirement on the reference system used in thermodynamic integration is that it faithfully mimics the fluctuations of ab initio energy in thermal equilibrium. Our ab initio melting curve of Fe over the pressure range $`50350`$ GPa agrees fairly well with experimental data obtained from both static-compression and shock techniques, but significant discrepancies remain to be resolved. Our ab initio predictions for quantities obtained directly from shock experiments, including the Grüneisen parameter of the liquid, agree closely with the measured data in most cases.

## Acknowledgments

The work of DA was supported by NERC Grant GST/02/1454 to G. D. Price and M. J. Gillan, and by a Royal Society University Research Fellowship. We thank NERC and EPSRC for allocations of time on the Cray T3E machines at Edinburgh Parallel Computer Centre and Manchester CSAR service, these allocations being provided by the Minerals Physics Consortium (GST/02/1002) and U.K. Car-Parrinello Consortium (GR/M01753). Calculations were also performed at the UCL HiPerSPACE Centre funded by the Joint Research Equipment Initiative. We thank Dr. G. Kresse for technical assistance with the PAW calculations.