AB INITIO EQUATION OF STATE FOR HYDROGEN–HELIUM MIXTURES WITH RECALIBRATION OF THE GIANT-PLANET MASS–RADIUS RELATION

The semi-analytical model by Saumon & Chabrier (1992, hereafter SC) for the equation of state (EOS) of hydrogen and its extension to hydrogen–helium mixtures (Saumon et al. 1995) were very successful and have been used in numerous calculations for the interiors of giant planets. However, with the development of ab initio computer simulation techniques, many uncontrolled approximations can now be avoided, simplifications inherent to analytical EOS models and severely limiting their predictive capabilities in the regime of high density and low temperature where interactions between particles are strong. Relying solely on analytical methods, it is difficult to determine the ionization state of the different chemical species that are present in the dense fluid.

Ab initio simulations allow one to study a fully interacting system of particles and to determine its properties by deriving the electronic states explicitly for every configuration of nuclei. No parameters are adjusted to match experimental data, but ab initio simulations still rely on approximations to solve the Schrödinger equation. However, they are not specific to the particular material or the pressure–temperature conditions under consideration.

In this paper, we rely on density functional molecular dynamics (DFT-MD) simulations that have been employed before to study hydrogen (Lenosky et al. 2000; Militzer et al. 2001; Desjarlais 2003; Bonev et al. 2004; Nettelmann et al. 2008, 2012; Morales et al. 2010; Caillabet et al. 2011; Collins et al. 2012), helium (Militzer 2006, 2009; Stixrude & Jeanloz 2008), and hydrogen–helium mixtures (Vorberger et al. 2007b, 2007a; Militzer et al. 2008; Militzer 2009; Hamel et al. 2011). While the computation of the pressure and the internal energy is straightforward from DFT-MD simulations, the entropy is not directly accessible. However, an accurate knowledge of the adiabats of hydrogen–helium mixtures at high pressure is of crucial importance for the determination of the temperature profile, the density, and the thermal energy budget in the interior of a giant planet. In 2008, two groups constructed Jupiter interior models from DFT-MD simulations (Militzer et al. 2008; Nettelmann et al. 2008). While the derived pressures and internal energies can be considered to be more reliable than those predicted by the SC model, both papers predicted very different interior temperature profiles for Jupiter (Militzer & Hubbard 2009). Using ab initio thermodynamic integration techniques (TDI), Militzer (2013) recently showed that the work by Nettelmann et al. (2008) overestimated the temperature at Jupiter's core–mantle boundary (CMB) by 2450 K (15%) while we underestimated it by 2870 K (18%) in Militzer et al. (2008). The revised temperature for Jupiter's CMB is 16,150 K and the corrections to the SC EOS model are in fact only −350 K.

At the conditions of Jupiter's CMB, hydrogen is metallic and characterized by a high degree of electronic degeneracy. Such a degenerate state is described rather well by the SC model. However, when we applied the TDI technique to explicitly determine the entropy over a wide range of pressure–temperature conditions, we identified a number of discrepancies between the DFT-MD results and the SC model. Near the molecular-to-metallic transition, our simulations predict a significant shift of the adiabat toward higher densities. At high temperature where electronic excitations matter, our computed entropies are higher than those of the SC model. We also do not perfectly reproduce the SC entropies in the molecular regime at low density.

Rather than providing a separate hydrogen and helium EOS and relying on the linear mixing approximation (Saumon et al. 1995; Nettelmann et al. 2008), we computed the EOS over a wide range of density–temperature conditions for a representative mixing ratio of N_He = 18 helium atoms in N_H = 220 hydrogen atoms, corresponding to a helium mass fraction of Y = 0.245, which is close to the solar value. This means that the nonideal mixing effects are fully incorporated. In Vorberger et al. (2007b), we showed, for example, that the presence of helium makes the hydrogen molecules more stable and reduces the dissociation fraction at given pressure and temperature. Even if other mixing ratios become of interest, as the result of helium rain (Stevenson & Salpeter 1977; Morales et al. 2009; Lorenzen et al. 2009; Wilson & Militzer 2010; McMahon et al. 2012), one is still better off by starting from an EOS for a typical hydrogen–helium mixture and then perturbing the mixing ratio by a comparatively small amount. Increasing or decreasing the helium fraction requires knowledge of a helium or hydrogen EOS, respectively. For the helium EOS, we recommend our first-principles computation (Militzer 2009) because it provides simulation data points for the pressure, P; internal energy, E; Helmholtz free energy, F; and entropy, S, over a wide parameter range and a thermodynamically consistent free energy fit for interpolation. For available hydrogen EOS work, we refer to the recent review by McMahon et al. (2012), but there has also been a considerable theoretical effort to compute the hydrogen EOS with semi-analytical techniques (Dharma-wardana & Perrot 2002; Kraeft et al. 2002; Rogers & Nayfonov 2002; Safa & Pfenniger 2008; Ebeling et al. 2012; Alastuey & Ballenegger 2012). If the perturbation in the helium fraction is sufficiently small, then one may use the SC EOS for simplicity.

We base our ab initio entropy calculations on our recent article (Militzer 2013), where we showed how the TDI technique can be extended to study molecular hydrogen and how it can be applied efficiently to determine the entropy at high temperature where electronic excitations matter. The TDI technique allows one to determine the difference in the Helmholtz free energy between two interacting many-body systems at fixed density and temperature (Morales et al. 2009; Wilson & Militzer 2010, 2012a, 2012b; McMahon et al. 2012; Wahl et al. 2013). We apply this method to determine the free energy difference between the DFT simulations and a system of classical forces that we construct:

$$\begin{equation} F_{\rm DFT} - F_{\rm cl} = \int _0^1 d\lambda \; \left< V_{\rm KS} - V_{\rm cl} \right>_\lambda. \end{equation} \tag{ 1 }$$

The angle brackets represent an average over trajectories governed by forces that are derived from a hybrid potential energy function, V_λ = λV_KS + (1 − λ)V_cl. V_cl is the potential energy of the classical system and V_KS is the Kohn–Sham energy (Kohn & Sham 1965). The presence of electronic excitations leads to an intrinsic contribution to the entropy and affects the forces on the nuclei (de Wijs et al. 1998) that need to be derived from the Mermin free energy (Mermin 1965), Ω = V_KS − TS_el. We combined both contributions into the following expression for the ab initio entropy (Militzer 2013):

$$\begin{equation} TS = \left< V_{\rm KS} \right> + \left< K_{\rm ion} \right> - \int _0^1 d\lambda \; \left< \Omega - V_{\rm cl} \right>_\lambda - F_{\rm cl}. \end{equation} \tag{ 2 }$$

<V_KS > includes contributions from partially occupied excited states. The λ integration was performed using five independent MD simulations with λ equally spaced between 0 and 1. To make this integration process efficient, we construct the pair potentials of the classical system to match the DFT forces as closely as possible (Izvekov et al. 2003). The computation of classical free energy is performed with Monte Carlo methods by thermodynamic integration to a system of noninteracting particles.

All simulations were performed with the VASP code (Kresse & Furthmüller 1996) with pseudopotentials of the projector-augmented wave type (Blöchl 1994) and a plane wave basis set cutoff of at least 1000 eV. The Perdew–Burke–Ernzerhof exchange-correlation functional (Perdew et al. 1996) was used throughout, but it has recently been shown that simulations based on the local density approximation yielded very similar results for Jupiter's deep interior (Militzer 2013). In the same article, we also performed a combined finite-size and k point analysis that demonstrated that simulations with 256 electrons and the zone-average point k = (1/4, 1/4, 1/4) are sufficiently accurate. All results that we report in this article were thus obtained with 220 hydrogen and 18 helium atoms in periodic boundary conditions.

We used an MD time step 0.2 fs, except for temperature of 50,000 K and above where we used a time step of 0.1 fs to accurately capture the more rapid collisions between particles at elevated temperatures. All standard DFT-MD simulations that we performed to determine P and E were 2.0 ps long, except at the highest temperatures, where 1.0 ps were found to be sufficient because the auto-correlation times are short and the error bars are small. All simulations were initialized with positions and velocity vectors from converged MD simulations at nearby densities and temperatures. This allowed us to run the TDI simulations for only 0.5 ps at each λ point.

We also adjusted the number of orbitals in the calculations to accommodate the partial occupation of excited electronic states according to the Mermin functional (Mermin 1965). The number of orbitals was increased until the error in the integral of the Fermi function was reduced to less than 10⁻⁵. This required many orbitals at high temperature and low density. As many as 816 orbitals were used, a significant increase in the computational cost over the 128 needed for ground-state calculations. This is the primary reason why we omitted simulations that would lead to entropy values above approximately 12.5 k_b el.⁻¹ The regime of higher temperatures can be studied much more efficiently with path integral Monte Carlo (PIMC) simulations because the computational cost of this alternative first-principles simulation technique scales like 1/T. PIMC simulations have been applied to hydrogen (Pierleoni et al. 1994; Magro et al. 1996; Militzer et al. 1999; Militzer & Ceperley 2000; Militzer & Graham 2006; Hu et al. 2010, 2011), helium (Militzer 2006, 2009), and hydrogen–helium mixtures (Militzer 2005) at high pressure and temperature and most recently also to study the EOS of carbon and water (Driver & Militzer 2012).

As Figure 1 indicates, we did not perform any TDI calculations at very high densities (r_s ⩽ 1.75) and temperatures below 5000 K because hydrogen–helium mixtures will exhibit signs of spontaneous phase separation in the course of a typical DFT-MD simulation (Militzer 2013) under such conditions. The immiscibility of hydrogen and helium at high pressure leads to helium rain on giant planets. This process is the accepted explanation for the observed luminosity excess of Saturn (Stevenson & Salpeter 1977) and has recently been invoked to explain the depletion of neon in Jupiter's atmosphere (Wilson & Militzer 2010). Because of constraints on the size and time of DFT-MD simulations, one can only observe the spontaneous separation of the mixture into a helium-rich and hydrogen-rich phase for conditions that are sufficiently deep in the region of immiscibility. The determination of pressure–temperature boundaries of this region can more accurately be derived with Gibbs free energy calculations. The need to accurately determine the non-ideal entropy of mixing was the primary motivation for Morales et al. (2009) to adopt the TDI technique to study hydrogen–helium mixtures. While our understanding of this process is still incomplete, its characterization goes beyond the scope of this study, which reports results for the mixed phase only, even though the separate state may have a lower Gibbs free energy for some high density points at 5000 K. We also revisit the question of a first order phase transition from molecular to metallic hydrogen because all ab initio simulations report that a discontinuous change in thermodynamic variables only occurs at temperatures well below those in the interiors of giant planets (Militzer et al. 2008; Morales et al. 2013).

Zoom In Zoom Out Reset image size

Finally, we need to discuss how we add the quantum effects in the motion of nuclei to our DFT-MD results that treat all nuclei classically (Militzer 2013). For the hydrogen–helium mixtures in giant planets, it is the vibrational state of the hydrogen molecules that matters for temperatures up to ∼5000 K while the rotational quantization is negligible for a planet's deep interior. A rigorous treatment of the nuclear motion with path integral molecular dynamics is very expensive, but while this paper was under review, the first results for the pressure and conductivity of dense hydrogen have been reported (Morales et al. 2013). An accurate method for the determination of the entropy is still not available. In this study, we use an approximate method to derive the quantum correction due to molecular vibrations from the intramolecular pair potentials, V_mol, that we have already derived. We derive the difference between the quantum and the classical entropy by solving a 1D problem through exact diagonalization. The resulting eigenstates are used to compute the canonical partition function that is compared with the classical analogue, an integral over the Boltzmann factor. The quantum corrections for all thermodynamic variables are derived from the comparison of both partition functions. Since our intermolecular potentials, V_mol, are derived from the DFT-MD forces at a particular density and temperature, our approach captures some of the interaction effects among hydrogen molecules and with helium atoms. The magnitude of this correction is discussed in Militzer (2013). We have added the zero-point correction to all EOS results that we discuss in the following section.

We report the computed EOS in the form of a table, a series of figures, and in analytical form as a two-dimensional fit of the free energy. In Table 1, we provide the thermodynamic functions that directly follow from analysis of the DFT-MD trajectories. The pressure and internal energy were computed for 391 different density–temperature points (see Figure 1). The 1σ error bars correspond to the statistical uncertainty that arises from the finite length of the MD simulations. For 131 points in Table 1, the thermodynamic integration was performed with five λ points and the free energy and entropy are reported in addition. Only counting the production runs that led to results in Table 1, the total CPU time consumed for this project amounted to 850,000 core-hours on Intel Nehalem processors. This is equivalent to using 100 cores for an entire year, which is a considerable amount of computer time by today's standards.

In Figures 2–9, we plot the internal energy, pressure, Helmholtz free energy, and entropy as a function of temperature and density. Every circle corresponds to a particular DFT-MD simulation listed in Table 1, without any interpolation being performed. The dashed lines are the results of the most common version of the analytical SC EOS model where the different thermodynamic functions have been smoothly interpolated across the molecular-to-metallic transition in hydrogen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To accommodate the wide parameter range of our simulations, we plot the different thermodynamic functions on logarithmic scale. Since these functions depend strongly on density and temperature, we added a second panel where we removed most of this dependence by introducing a scale factor equal to r_s or T raised to some power. Here r_s is the Wigner–Seitz radius that specifies the density of the system according to $({4 \pi }/{3}) r_s^3=V/N_e = n^{-1}$. The number density, n, is derived from the number of electrons, N_e = (N_H + 2N_He), per unit volume V. The mass density is given by ρ = n(N_Hm_H + N_Hem_He)/(N_H + 2N_He), where m_H and m_He are the masses of the hydrogen and helium atoms. The rescaling of the ordinate makes it easier to identify the deviations from the SC model while our simulation results can still be reproduced easily. The ordinates are plotted in atomic units. Lengths including r_s are given in Bohr radii (a₀ = 5.29177209 × 10⁻¹¹ m), energies in Hartrees (4.35974380 × 10⁻¹⁸ J) per electron (el.), and entropies are specified in units of k_b el.⁻¹, where k_b is Boltzmann's constant.

Figure 2 shows a comparison between the internal energies from DFT-MD simulations with the predictions of the SC EOS model. For a low density of r_s = 2.4, excellent agreement is found for a temperature range from 1000 to 20,000 K. Hydrogen gradually changes from a molecular state to an atomic state in this temperature interval and, from the good agreement, one may conclude that the thermally activated dissociation of molecules is well described in the SC model. However, above 20,000 K, the SC model predicts a strong and artificial increase in the internal energy that is the result of an inaccurate description of electronic excitations. This deviation was first identified by Militzer & Ceperley (2001) when predictions from the SC model were compared with PIMC simulations. Figure 3 shows that this deviation is present at 20,000 K for the whole density range under consideration and extends to much higher temperatures also.

Figure 2 shows that the favorable agreement between DFT-MD results and SC predictions below 20,000 K continues to hold up to a density of r_s = 1.6. When the internal energy is compared for a higher density of r_s = 1.0 or 1.2 where hydrogen is metallic, one finds that DFT-MD results and SC predictions are offset by a nearly constant amount.

The internal energy curves of r_s = 1.6 and 2.4 appear to cross over in Figure 2 at a temperature of 27,000 K, which is consistently predicted by DFT-MD results and the SC model. Figure 3 shows that this is a consequence of internal energy exhibiting a minimum when plotted at constant temperature as a function of density. At high density, the internal energy sharply rises because of Pauli exclusion effects between the electrons. In the low density limit, the internal energy rises also because the ionization fraction increases as a result of the increased gain in entropy that is associated with electrons becoming free particles.

In Figure 4, we compare the pressure predicted from the DFT-MD simulation with the SC model. At a high density of r_s = 1.0 where the hydrogen–helium mixture is metallic, we find fairly good agreement over the entire temperature range. This implies that the deviation that we identified for the internal energy in this regime varies slowly with density and does not significantly affect the pressure in the SC model.

At a low density of r_s = 2.2 and 2.4, we found good agreement up to a temperature of 5000 K. At this temperature, we see a small decrease in slope in the DFT-MD data that is missing in the predictions of the SC model. We attribute this slope change to the dissociation of molecules in the DFT-MD simulations. At 20,000 K, the SC model predicts a significant decrease in slope that is not present in the DFT-MD data. This slope change in the SC predictions can again be attributed to an inaccurate description of ionization, which leads to deviations over the whole density range under consideration (Figure 5). At an intermediate density of r_s = 1.6 close to the molecular-to-metallic transition, we find that the SC model overestimates the pressure for temperatures up to about 20,000 K and underestimates for higher temperatures. The deviations around 100 GPa, 5000 K, and r_s = 1.6 (0.75 g cm⁻³) are of particular significance. The DFT-MD simulations predict pressures that are much lower than those of the SC model. This leads to a significant departure in the resulting adiabats. Its implication for the interiors of giant planets will later be further analyzed.

In Figures 6 and 7, the Helmholtz free energy from DFT-MD simulations and the SC model are compared. In general, the agreement appears to be much better than for other thermodynamic functions that are derivatives of it. Still, one finds that the SC model overestimates the free energy in the metallic regime, mirroring the deviations that we have discussed for the internal energy.

In Figure 8, the entropies at different densities are compared as a function of temperature. At a very high density of r_s = 1.0, very good agreement between DFT-MD results and the SC model is found up to 50,000 K. For lower densities, the SC model predicts a sharp entropy increase at 20,000 K, which is again a result of the treatment of ionization effects. This trend is not confirmed by the DFT-MD simulations. One also finds significant deviations at lower temperatures, in particular around r_s = 1.6. Even at a relatively low density of r_s = 2.2, the agreement is not perfect. From 4000 to 20,000 K, the SC model underestimates the entropy and it overestimates the entropy for lower temperatures. Figure 9 shows that such deviations persist over a wider density range. In principle, one expects a non- or weakly interacting gas of hydrogen molecules and helium atoms to be perfectly described by the SC model. However, the density that we can efficiently study with DFT-MD simulations does not yet appear to be low enough for the deviations to decay to zero.

We fitted our ab initio results for P, E, F, and S in Table 1 with a two-dimensional spline function that represents the Helmholtz free energy in terms of temperature, T, and electron density, n. By construction, this fit is thermodynamically consistent. We employ the same functional form that we used to represent the free energy of hot, dense helium in Militzer (2009), except the splines here are functions of n rather than log (n). Table 2 provides the free energy as well as the required derivatives on a number of (n, T) knot points. Atomic units are used throughout.

Table 2. Coefficients of Free Energy Fit for the Equation of State

	rs	T	Coefficient
	(a.u.)	(a.u.)	(a.u.)
$\frac{\partial F}{\partial T}$	1.75	0.001013	−8.67942 × 10−1
F	1.75	0.001013	−9.23942 × 10−2
F	1.75	0.003131	−9.76575 × 10−2
F	1.75	0.006512	−1.12223 × 10−1
F	1.75	0.011911	−1.44120 × 10−1
F	1.75	0.020533	−2.07646 × 10−1
F	1.75	0.034302	−3.23222 × 10−1
F	1.75	0.056290	−5.29287 × 10−1
F	1.75	0.091403	−8.92177 × 10−1
F	1.75	0.147476	−1.53056
F	1.75	0.237021	−2.64345
F	1.75	0.380018	−4.58672
$\frac{\partial F}{\partial T}$	1.75	0.380018	−1.41575 × 10+1

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

To evaluate the fit for (n*, T*), we first construct a separate one-dimensional cubic spline function, F_n(T), for every density on a grid ranging from r_s = 3.581 to 0.536 (0.0670–20.0 g cm⁻³). At every density, the free energy is given on a number of temperature knots and its first derivative, (∂F/∂T)|_n, is specified for the highest and lowest temperatures. We construct a similar one-dimensional spline function that represents (∂F/∂n)|_T(T) at the smallest and largest density. We then evaluate all these spline functions at T* and construct a one-dimensional spline function $F_{T^*}(n)$ from the free energy values and its first derivatives at the boundaries. This provides us not only with a straightforward way to obtain the free energy at every (n, T) point but we can also derive the pressure and entropy by taking analytical derivatives

$$\begin{equation} P=n^2\left.\frac{\partial F}{\partial n}\right|_T \;\;{\rm and}\;\; S={-}\left.\frac{\partial F}{\partial T}\right|_n \;\;. \end{equation} \tag{ 3 }$$

The internal energy and Gibbs free energy then follow from E = F + TS and G = F + PV.

When we constructed this fit, we made sure every EOS point in Table 1 was well reproduced. We extended the domain of the fit a bit beyond the range of the DFT-MD data. This leads to a smoother representation of the data in the interior of the domain and also allows us to gradually approach the SC EOS in the limit of low density. As Figure 10 shows, we were able to smoothly match the SC adiabats for entropy values from 6 to 10 and again for 13 and 14 k_b el.⁻¹ A disagreement remains for S = 11 and 12 k_b el.⁻¹, but the SC EOS is not thermodynamically consistent in the regime of 10,000 to 20,000 K and no attempt was made to reproduce those adiabats. So far, only the one exoplanet HAT-P-32b (Hartman et al. 2011) appears to have an internal entropy in excess of 11 k_b el.⁻¹; see Figure 12.

Zoom In Zoom Out Reset image size

We also find a significant disagreement in the high density limit between our DFT-MD adiabats and the predictions of the SC model. Starting with entropy values of S = 10 k_b el.⁻¹, the DFT-MD results predict the adiabats reach states of higher temperatures and higher pressure for a given density.

The most significant result of Figure 10 is the deviations along the S = 7 k_b el.⁻¹ adiabat. The DFT-MD simulations predict a decrease in the slope of the adiabat exactly where the hydrogen molecules dissociate (Militzer et al. 2008). Since the SC model interpolates between separate atomic/metallic and molecular thermodynamic descriptions, it has no predictive power in the regime of pressure dissociation where all the different species interact strongly.

In Figure 11, we compare different predictions for the interior adiabat of a Jupiter-like planet that were all constructed to pass through Jupiter's actual temperature at 1 bar, 166.1 K. A constant hydrogen–helium ratio of 220:18 without heavier elements was consistently assumed in all calculations to facilitate a direct comparison. The measured hydrogen–helium ratio for Jupiter's atmosphere is 220:15 (von Zahn et al. 1998), while the protosolar hydrogen–helium ratio (presumed equal to Jupiter's bulk ratio) corresponds to about 220:21 (Lodders 2003); thus, our adopted ratio is a compromise between these values. This implies that Figure 11 only approximately represents Jupiter's interior because typical models include both a small fraction of heavier elements and a non-constant helium fraction for the molecular and metallic layers in order to match various observational constraints such as mass, radius, and the gravitational moments J₂, J₄, and J₆. Such modifications change Jupiter's interior pressure–temperature profile by a small amount, on the order of 3%. For Saturn, a larger change is expected but nevertheless the interior pressure–temperature profiles for all objects discussed in this article were constructed with the simplifying assumption of a constant hydrogen–helium ratio of 220:18 without the consideration of heavier elements.

Zoom In Zoom Out Reset image size

The comparison in Figure 11 reveals sizable discrepancies (Militzer & Hubbard 2009) in temperature between the different curves, although one might expect all DFT-MD-based curves to agree in the displayed T–P region because the underlying EOS data there were generated by the same method. See, however, the discussion below. At a low pressure of 10 GPa, the Nettelmann et al. (2008) work predicted temperatures that were 5.7% (210 K) higher than our new TDI results. At a pressure of 4000 GPa near Jupiter's CMB, the deviation amounts to 19% (3300 K). In Nettelmann et al. (2012), the temperature was overestimated by 7.5% (300 K) and 15% (2450 K) at 10 and 4000 GPa, respectively. The associated revision of the adiabats introduced a more pronounced reduction in slope into the region of molecular dissociation from 20–70 GPa. In Militzer et al. (2008), the temperature was underestimated at 4000 GPa by 18% (2450 K).

There are three main reasons for the deviations that all relate to approximations that are avoided in the TDI method: (1) inaccuracies in the analytical EOS, (2) approximations when the adiabats are derived from the ab initio pressures and internal energies, and (3) linear mixing approximation when the properties of a mixture are derived from simulations of pure hydrogen and helium.

(1) Only the TDI method presented here allows one to directly determine absolute entropies. In Militzer et al. (2008) and Nettelmann et al. (2008, 2012), the authors instead relied on the pressures and internal energies from standard DFT-MD simulations to determine the shape of the adiabats in the region from 10–4000 GPa where such simulations work efficiently. A separate analytical technique was needed to determine the starting point of the adiabat at ∼10 GPa. The analytical fluid variational theory used in Nettelmann et al. (2008, 2012) to anchor the adiabats is not compatible with our TDI results, which explains the deviations of the adiabats at 10 GPa. Thus, the discrepancies at low pressure are not a consequence of the ab initio simulations but rather the result of inaccuracies in the anchor point, which then also affect the adiabats at higher pressure that are constructed from ab initio results. In Militzer et al. (2008), the SC model was used as a starting point, which has relatively good agreement with the TDI results. (2) When the adiabats are constructed from the ab initio pressures and internal energy (Militzer 2009), (∂T/∂V)_S = −T(∂P/∂T)_V / (∂E/∂T)_V , one needs P and E on a fine grid of density–temperature points for interpolation. In practice, one can only perform a finite number of simulations and the results have finite size and statistical uncertainties. An insufficient number of simulations is the reason why the temperature on Jupiter's adiabat at 4000 GPa was underestimated in our 2008 work and why the shape of the adiabat was not accurately determined in Nettelmann et al. (2008). Finally (3), the linear mixing approximation in Nettelmann et al. (2008, 2012) has some impact on the adiabats. However, this is a small effect compared to points (1) and (2). With the knowledge of our TDI results, N. Nettelmann & R. Redmer (2013, private communication) changed the anchor point and constructed the selected adiabat as shown in Figure 11. The agreement improved substantially but a small deviation of −3% at 10 GPa and +2.5% at 4000 GPa remained, which can be attributed to points (2) and (3).

The vexing problem of radius anomalies of transiting giant planets (Burrows et al. 2007) has continued with the addition of more objects (Laughlin et al. 2011). Figure 12 shows measurements posted in the online Exoplanet Encyclopaedia³ (Schneider et al. 2011) as of late 2012 and different theoretical curves that we will discuss below. Briefly, the problem arises from a significant population of exoplanets that have radii too large to be explained by thermal distension from retained primordial heat, and there is a further population with radii well below those expected for primarily H–He composition even at zero temperature. Any point that falls below the S = 6 k_b el.⁻¹ curve can be explained by invoking the presence of a rocky core and/or admixture of heavy elements, which reduces the radius for given mass (Miller & Fortney 2011). But it is not so simple to classify the population of anomalously distended giant exoplanets, for the degree of distension depends on such factors as the planet's age and degree of irradiation from the host star (Fortney & Nettelmann 2010) and possible additional heating mechanisms such as ohmic dissipation (Batygin et al. 2011).

Zoom In Zoom Out Reset image size

In order to most clearly exhibit the differences in predicted radius between the DFT-MD simulations and the analytical SC model, we model a planet of mass M as an H–He object of constant entropy S and fixed helium fraction of Y = 0.245 with neither a rocky core nor heavy element component in the gas envelope. Since we do not have DFT-MD simulation data at very low densities, we switch back to the SC model below 0.0670 g cm⁻³, the lowest density of the free energy fit to our DFT-MD data.

As is well known (Chandrasekhar 1957), formally, such an object has a precisely defined radius where the temperature T, mass density ρ, and pressure P simultaneously go to zero. In a real object, the ideal-gas outer layers cannot persist in such an isentropic state and instead the temperature reaches a finite limit set by the effective temperature for the radiation balance in the outer layers. The radius as measured by transit observations also depends on sources of slant opacity in these outer layers. As explored by Burrows et al. (2007), the value of such a radius can vary by several 0.1 R_J, depending on the atmospheric model. Adjusting the parameters controlling the atmosphere can move a model closer to agreement with objects of inflated radii, but sometimes a mismatch remains. A major result of the present paper is that differences in the EOS alone can also lead to radius changes of several 0.1 R_J. Because we consider only strict-adiabatic models, the deviations that we point out are entirely due to differences in the EOS at high pressure. In Figure 13, we exhibit these differences for various entropies. For entropy values up to 9 k_b el.⁻¹, the DFT-MD calculations consistently predict smaller planet radii than the SC model, which is a direct consequence of the density enhancement on the adiabats around 100 GPa as illustrated in Figures 10 and 11. Figure 10 also shows that the DFT-MD and SC adiabats for 10, 11, and 12 k_b el.⁻¹ cross over in the density range from 0.15 to 0.7 g cm⁻³. This is the reason why the DFT-MD calculations predict larger planet radii than the SC model for massive planets with M > 0.5M_J but significantly smaller radii for light planets. The deviations between the DFT-MD and SC predictions in Figure 13 reach values up to approximately 0.4 Jupiter radii.

Zoom In Zoom Out Reset image size

The exoplanet HD 209458b (middle open data point in Figure 12) fortuitously falls near an entropy S ≈ 9.5 k_b el.⁻¹ and mass M ≈ 0.7M_J where ΔR ≈ 0. Nevertheless, the interior T–ρ and P–ρ profiles in Figures 14 and 15 differ significantly between the DFT-MD and SC EOSs. These figures also compare interior profiles for simplified (pure H–He mixtures on an adiabat) models of Jupiter and Saturn.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 16 shows the differences in evolutionary behavior of our simplified planetary models. This figure plots the value of the central temperature, T_central, versus the central density, ρ_central, for a range of adiabats and masses. During the evolution of a planet of constant mass, the central density increases monotonically while the central temperature exhibits a maximum. During the initial contraction, the temperature in the center increases at first as the material is subjected to increasing pressure. When a degenerate interior state is reached, the contraction ceases and the whole planet starts to cool. According to DFT-MD simulations, the maximum temperature reached is up to 10,000 K lower than predicted by the SC model. This deviation may have consequences for the evolution of cores in giant planets that remain to be explored.

Zoom In Zoom Out Reset image size

This paper provides an equation state table for hydrogen–helium mixtures in giant-planet interiors that was derived from ab initio computer simulations. The combination with an efficient thermodynamic integration technique enabled us to calculate the absolute entropy and free energy directly, in addition to pressure and internal energy that follow from standard simulations.

Our complete EOS table with 391 density–temperature points as well as a thermodynamically consistent free energy fit is included in this publication so that our EOS can be easily incorporated in future models for giant-planet interiors.

We have identified significant deviations for the SC EOS models. The new DFT-MD EOS causes low-entropy giant-planet models (S ⩽ 8 k_b el.⁻¹) to shrink in comparison to SC models by up to 0.08 Jupiter radii. However, for hot giant planets with mass exceeding 0.5 M_J and interior entropy values ranging from 10–12 k_b el.⁻¹, the DFT-MD simulations predict significantly larger radii. The correction to the SC model reaches 0.4 Jupiter radii for the hottest planets. Thus, the revision suggests that some of the most inflated giant exoplanets are at lower entropies than was previously inferred. Our revision could ameliorate the "inflated giant exoplanet" discrepancy to some extent, but perhaps not for HD209458b. The matter is to be revisited with detailed evolutionary calculations based on our revised EOS.

This work has been supported by NASA and NSF. Computational resources at NCCS were used. We thank N. Nettelmann and R. Redmer for sharing their adibats for a 220:18 hydrogen–helium ratio.