AB INITIO EQUATIONS OF STATE FOR HYDROGEN (H-REOS.3) AND HELIUM (He-REOS.3) AND THEIR IMPLICATIONS FOR THE INTERIOR OF BROWN DWARFS

2.1.1. DFT-MD Data

2.1.2. EOS Table Composition

The several components of our hydrogen EOS table (H-REOS.3) are shown in Figure 1. The orange region covers states of matter that represent the greatest challenge for many-particle theory since strong correlations and quantum effects occur that lead to pressure-induced dissociation and ionization processes. An appropriate method to describe such states is the DFT-MD technique; see the previous subsection.

To cover the remaining parts of the ρ–T plane, we connected the DFT-MD data to EOS data derived from chemical models such as the FVT for the molecular and dissociated fluid (see Juranek et al. 2002), the SCvH-EOS (Saumon et al. 1995) for the partially ionized plasma at lower densities, and the Chabrier–Potekhin (CP) model for fully ionized matter (Chabrier & Potekhin 1998). Note that interpolations between two EOS models had only to be performed at the interface between the DFT-MD data and the region labeled with "path integral Monte Carlo (PIMC) based interpolation" in Figure 1; see also below. At all remaining interfaces, the points of the different EOS models merge smoothly into each other. In particular, the CP model connects almost perfectly to the DFT-MD data at high densities; see Figure 2. For low temperatures (black: 100 K isotherm) as well as for high temperatures (red: 100 kK isotherm) the transition at 70 g cm⁻³ is performed with a deviation ⩽0.15% in the thermal EOS P(ρ, T). We find the same accuracy for intermediate temperatures in P(ρ, T) as well as for the transition within the caloric EOS u(ρ, T).

Zoom In Zoom Out Reset image size

However, the low-density connection to the DFT-MD data turned out to be more difficult. Up to 10 kK we used the FVT data which lead to similar results as the SCvH-EOS but can be connected more smoothly to the DFT-MD data. The subtle details appear in regions where significant ionization occurs, e.g., in the transition region left of the DFT-MD data above 10 kK; see Figure 1. The choice between the FVT+ model (FVT including ionization; see Holst et al. 2007) and the SCvH model came out in favor of the latter for the following reason: We compared the PIMC results from Hu et al. (2011) with the 30 kK, 50 kK, and 100 kK isotherms of the FVT+ and the SCvH-EOS; see Figure 3. We plot the ratio of pressure and density versus the density. In this P/ρ-representation, non-ideal contributions to the EOS can nicely be seen as deviations from horizontal lines that represent the ideal gas law (P/ρ = RT/M with the molar mass M).

Zoom In Zoom Out Reset image size

It turns out that the accurate PIMC data (triangles) clearly favor a connection from the SCvH data (dashed lines) to the DFT-MD data (filled circles); the FVT+ data (dotted lines) underestimate the pressure substantially. For 30 kK and 50 kK, the PIMC data are very close to the SCvH results in the low-density range. When the SCvH data start to underestimate the pressures at densities above 0.1 g cm⁻³, the PIMC results turn into the direction of the DFT-MD data—both ab initio methods agree there nicely. Thus, we have performed a smooth interpolation from the SCvH data via the PIMC data to the DFT-MD results for the final H-REOS.3 isotherms (solid lines) from 20 kK to 150 kK; see interpolation regime in Figure 1. In detail, the following data points of the H-REOS.3 table are generated using natural cubic spline interpolations: 0.05–0.1 g cm⁻³ at 20 kK, 0.03–0.1 g cm⁻³ at 30 kK, 0.06–0.2 g cm⁻³ at 50 kK, 0.08–0.2 g cm⁻³ at 75 kK, 0.3–0.5 g cm⁻³ at 100 kK, and 0.1–0.5 g cm⁻³ at 150 kK. For temperatures above 150 kK, hydrogen is fully ionized and well described by the CP model. This can be connected to the SCvH data at lower densities without violating thermodynamic consistency; see Section 2.4.

The helium EOS (He-REOS.3) consists of three major parts; see Figure 4: DFT-MD EOS data cover the intermediate and high-density range of the density-temperature plane, where correlation and degeneracy effects are important and have to be treated accurately, while for temperatures up to 10 kK and low densities we use a virial EOS based on an ideal gas model and very accurate virial coefficients; see Section 2.2.2. It turned out that the helium EOS model of Saumon et al. (1995; SCvH-He) connects smoothly to our DFT-MD data for most T–ρ points and is thus the most appropriate one for our purpose. Interpolations were only performed between 60 kK and 300 kK since ionization of the He atoms occurs in the transition region between the SCvH and DFT-MD data. In the following, we describe the methods used for generating the EOS data, discuss the composition of the final EOS table, and compare with other ab initio data for helium.

Zoom In Zoom Out Reset image size

2.2.1. DFT-MD Data

2.2.2. The Virial and High-temperature EOS

For the intermediate region between the strongly correlated quantum regime treated with DFT-MD and the ideal gas limit, we introduced correction terms to the ideal gas equation. They describe weak coupling between the particles and result in a smooth transition between the DFT-MD data and the ideal gas EOS. Such correction terms can be expressed by virial coefficients B_i(T). The respective virial EOS is composed of an ideal part P^id and a correlation part P^cor as follows:

$$\begin{equation} P(\rho,T) = P^{\mathrm{id}} + P^{\mathrm{cor}} = \frac{\rho RT}{M} \left[ 1 + \sum \limits _{i=2}^\infty B_i(T)\left(\frac{\rho }{M}\right)^{i-1}\right] \;. \end{equation} \tag{ 1 }$$

R is the universal gas constant and M the molar mass of helium. Note that in this notation the ideal part ρRT/M is without ionization, while we finally use a more elaborated ideal gas model with ionization; see below. The contribution of the virial coefficients to the specific internal energy u = U/m is calculated via the fundamental connection between the thermal and caloric EOS,

$$\begin{equation} -\rho ^2\left(\frac{\partial u}{\partial \rho }\right)_T = T\left(\frac{\partial P}{\partial T}\right)_{\rho } -P \quad. \end{equation} \tag{ 2 }$$

Using Equations (1) and (2), we obtain

$$\begin{equation} u(\rho,T) = u^{\mathrm{id}}-\frac{RT^2}{M}\left[ \sum \limits _{i=2}^\infty \frac{\mathrm{d}B_i(T)}{\mathrm{d}T}\frac{1}{i-1}\left(\frac{\rho }{M}\right)^{i-1} \right] \,. \end{equation} \tag{ 3 }$$

We investigated several sets of virial coefficients from B₂(T) up to B₅(T) (Slaman & Aziz 1991; Bich et al. 2007; Cencek et al. 2012; Shaul et al. 2012a, 2012b) to find the best transition from the ideal gas to the DFT-MD data. The virial coefficients finally used are given in Table 4.

Table 4. Virial Coefficients B_i taken from Shaul et al. (2012b) by Interpolation of their given Data (Labeled with a) and from Bich et al. (2007) where no Interpolation was needed (Labeled with b) which we used for the Virial Part of our Helium EOS

T	B2	B3	B4	B5
(K)	(cm3 mol−1)	(cm6 mol−2)	(cm9 mol−3)	(cm12 mol−4)
60a	7.46483	161.676	1553.11	15478.4
100a	10.4982	147.007	1322.62	9857.49
200a	11.6929	121.309	882.001	4458.25
300a	11.5403	104.351	650.567	2547.36
600b	10.651	78.73	⋅⋅⋅	⋅⋅⋅
1000b	9.5497	59.44	⋅⋅⋅	⋅⋅⋅
2000b	7.9556	38.60	⋅⋅⋅	⋅⋅⋅
3000b	7.0330	29.16	⋅⋅⋅	⋅⋅⋅
6000b	5.5459	17.14	⋅⋅⋅	⋅⋅⋅
10000b	4.5542	11.05	⋅⋅⋅	⋅⋅⋅

Download table as:  ASCIITypeset image

For temperatures below 600 K, we used the classical results for B₂(T) − B₅(T) from Shaul et al. (2012b) derived from state-of-the-art Mayer-sampling Monte Carlo calculations. For temperatures between 600 K ⩽ T ⩽ 10 kK, we found smoothest connections using the ab initio virial coefficients B₂(T) and B₃(T) of Bich et al. (2007). They derived their results from a high-quality 18 parameter ab initio interaction potential for helium; see Hellmann et al. (2007) for details.

We calculated the ideal parts P^id and u^id of the virial EOS from the free energy of an ideal plasma model as described in Förster et al. (1992):

$$\begin{equation} F = Nk_BT\left[\sum \limits _{z=0}^2 \alpha _z\mathrm{ln}\left(\frac{n_z\lambda _z^3}{\sigma _z}\right) -1\right]. \end{equation} \tag{ 4 }$$

The index z represents the z-fold charged atoms with a relative fraction of α_z = n_z/n with respect to the total atom density n = N/V = n₀ + n₁ + n₂. The thermal wavelength of the atoms is given by λ_z = h(2πm_zk_BT)^−1/2 and the internal partition functions σ_z were calculated using the Planck–Larkin convention (Planck 1924; Larkin 1960):

$$\begin{eqnarray} \sigma _z &=& \exp \left(-\sum \limits _{z^\prime =0}^z \frac{I_{z^\prime }}{k_BT} \right)\\ && \times \sum _m g_{z,m}\left[ \exp \left(\frac{-E_{z,m}}{k_BT}\right)-1+\frac{E_{z,m}}{k_BT} \right],\nonumber \\ z&=&0,1,\quad \sigma _2 = \exp \left(-\sum \limits _{z^\prime =0}^2 \frac{I_{z^\prime }}{k_BT} \right). \end{eqnarray} \tag{ 5 }$$

The ionization energies I_z, the energy levels E_{z, m}, and the corresponding statistical weights g_{z, m} are taken from the NIST database (Kramida et al. 2012). Note that here is one difference to the SCvH-He EOS (see below), where only the ionization energies I_z are considered in the internal contribution to the free energy while the E_{z, m} terms are neglected.

The relative fractions of the atoms α_z are determined by coupled Saha equations

$$\begin{equation} \mu _{0} = \mu _{1} + \mu _{e^-} ,\quad \mu _{1} = \mu _{2} + \mu _{e^-} \;, \end{equation} \tag{ 7 }$$

with the boundary condition of charge neutrality of the system: n_e = n₁ + 2n₂. P^id and u^id are then obtained by differentiation of the specific free energy (f = F/m):

$$\begin{equation} P = \rho ^2\left(\frac{\partial f}{\partial \rho } \right)_T ,\quad u = f-T\left(\frac{\partial f}{\partial T} \right)_\rho \;. \end{equation} \tag{ 8 }$$

For temperatures above 10 kK, helium starts to ionize. The additional ideal contributions to the pressure and the internal energy from the unbound electrons lead to a non-horizontal behavior of the ideal gas isotherms; see Figures 7 and 8.

It turns out that the helium EOS of Saumon et al. (1995; SCvH-He) is the most appropriate one above 10 kK for our purposes: describing the partial and fully ionized helium with interaction and correlation effects and a smooth transition into the DFT-MD data. The full description of their EOS can be found in their paper, while we only focus on the main points for a better understanding of the composition of our EOS table. The free energy of their model can be written as

$$\begin{equation} F(n_0,n_1,n_2,n_e,V,T) = F_{\mathrm{id}} + F_{\mathrm{conf}} + F_{\mathrm{DH}} + F_{\mathrm{int}} \;. \end{equation} \tag{ 9 }$$

F_id + F_int, the ideal and internal part, is in principle the same expression as Equation (4), but with a simpler approximation for the internal partition function—it only takes into account the ionization energies I_z and treats electrons on the level of the Fermi integral in order to include degeneracy effects. The contribution F_conf is due to the interaction of the helium atoms and calculated via fluid perturbation theory. The interaction of charged particles is treated using the Debye–Hückel approximation F_DH. The fully ionized plasma for densities above 3 g cm⁻³ is described by a screened one-component plasma model (Chabrier 1990).

We extended the SCvH-He isotherms down to 10⁻¹⁰ g cm⁻³ (smallest density of the virial EOS) using the ideal plasma model that connects almost perfectly to the SCvH-He data. This is why the respective red area in Figure 4 is labeled as "SCvH-EOS+ideal plasma limit." As for hydrogen, the connection of the partially ionized plasma to the DFT-MD data becomes difficult when significant ionization takes place in the transition region. For helium this is the case between 60 kK and 300 kK where we applied spline interpolation for generating a smooth transition between the SCvH-He and the DFT-MD data; see Figure 4.

2.2.3. EOS Table Composition

Having presented the constitutive parts of our EOS tables in the previous section, we now explain how we constructed the final EOS tables. The requirements of an accurate description of the underlying physics, of preserving thermodynamic consistency, and obtaining smooth transitions between the various models could always be met within a reasonable compromise. The isotherms of each model are shown in Figures 5–8. The pressure isotherms are again displayed in the P(ρ, T)/ρ representation that reduces the ideal contribution to a constant offset for each isotherm and stresses the non-ideal contributions to the EOS.

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

As mentioned above, for temperatures below 20 kK, the transition from the ideal gas to the DFT-MD data has been performed via the virial EOS; see Figures 5 and 6 where the contribution of the virial coefficients are in principle those parts of the isotherms that differ from the ideal gas (dotted) horizontal lines below 0.1 g cm⁻³ above 300 K. There, the virial expansion connects right into the lowest DFT-MD data. Taking into account pressure and energy isotherms below 300 K, it is reasonable to use the virial EOS up to 0.3 g cm⁻³ because DFT-MD simulations are challenging there. For example, at 60 K the virial expansion shows the right negative slope in internal energy, indicating the transition to the solid, which is covered by the DFT-MD data above 0.3 g cm⁻³. However, the respective DFT-MD pressures are partly below the ideal gas that is not indicated by the virial expansion. That is why the virial expansion is favored in this regime.

As shown in Figures 5 and 6, we find a very good agreement with the ab initio data from Militzer (2009; open triangles) except small deviations in pressure for the two lowest densities below 10 kK. EOS data from the SCvH-He model (dashed lines) overestimate the pressure significantly below 3000 K and densities around 0.1 g cm⁻³ and overestimate the internal energy for all isotherms up to 10 kK at the same densities. Note that the SCvH isotherms merge into the corresponding DFT-MD data for temperatures above 10 kK at the transition densities. This is why we prefer this model to extend our EOS table into the partial and fully ionized regime of the density-temperature plane; see Figures 7 and 8, on which we focus our discussion now.

The ideal gas model includes, by definition, no interaction effects, and hence no pressure ionization at higher densities. This leads to a false prediction of an atomic system for increasing density caused by the higher recombination probabilities within the Saha Equation (7), leading to a drop of the pressure and energy isotherms (dotted curves). However, the ideal plasma model is exact in the low-density limit and therefore serves as the extension of the ideal SCvH-He data down to 10⁻¹⁰ g cm⁻³.

The most difficult connection between the SCvH-He and DFT-MD data occurred for temperatures between 60 kK and 300 kK; see Figure 4. This is caused by the ionization of helium and the effects of recombination and pressure ionization for higher densities. While we used the SCvH-He data at the 20 kK and 30 kK energy isotherms up to 0.5 g cm⁻³, we had to apply interpolations, in particular for the pressure isotherms from 60 kK on. In detail, the following data points of the He-REOS.3 table are generated using natural cubic spline interpolations: 0.01–0.1 g cm⁻³ at 60 kK, 0.05–0.3 g cm⁻³ at 100 kK, 0.1–3 g cm⁻³ at 200 kK, and 0.1–3 g cm⁻³ at 300 kK. This interpolated data can be identified in Figures 7 and 8, where the solid He-REOS.3 curves deviate from the SCvH-He data (dashed lines). These data underestimate the pressures systematically in the relevant region predicted from our simulations (circles) and those of Militzer (2009; open triangles), which coincide.

As shown in Figures 7 and 8, the transition from the SCvH-He model to the DFT-MD data for temperatures above 300 kK proceeds very smoothly, while all ab initio data (circles and open triangles) are located on the final He-REOS.3 isotherms (solid lines). Thus, no interpolations were performed here. The 6 MK and 10 MK isotherms are pure SCvH-He data which proceed almost horizontally in the representations and thus are nearly perfect ideal gas isotherms.

The zero point u₀ of the specific internal energy u(ρ, T) can be chosen arbitrarily. We decided to fix our EOS tables to the limiting case of the perfect ionized ideal plasma at very high temperatures and very low densities. One can simply show that this limiting case is material dependent and obeys the following equation:

$$\begin{equation} u=\frac{3}{2}\frac{Nk_BT}{m}=\frac{3}{2}\frac{k_BT}{\rho }(n_e+n_i) =\frac{3}{2}\frac{(Z+1)RT}{M_{{\rm at}}}. \end{equation} \tag{ 10 }$$

N is the total particle number of the system, k_B the Boltzmann constant, n_e and n_i the electron- and ion-densities of the system, Z the atomic number, and M_at the molar mass of the atomic system. For hydrogen with Z = 1 and M_at = 1 g  mol⁻¹, we fixed the zero point to u₀(10 MK) = 249.433 MJ g⁻¹. For helium with Z = 2 and M_at = 4 g    mol⁻¹, we obtained u₀(10 MK) = 93.538 MJ  g⁻¹. Therefore, at the lowest densities at 10 MK in the H-REOS.3 and He-REOS.3 tables, the specific internal energy has the respective u₀ values; see Tables 2 and 3.

A sensible test for the accuracy of our REOS.3 tables is checking thermodynamic consistency that is given by Equation (2). To get a relative deviation with respect to the pressure of the system, we divide Equation (2) by P and subtract unity which leads to

$$\begin{equation} \frac{T}{P}\left(\frac{\partial P}{\partial T} \right)_\rho + \frac{\rho ^2}{P}\left(\frac{\partial u}{\partial \rho } \right)_T - 1= \Delta \:. \end{equation} \tag{ 11 }$$

For instance, Δ = 0 denotes perfect thermodynamic consistency while Δ = 0.1 implies a violation of 10%. The derivatives that are needed for the evaluation of Equation (11) are obtained from a more resolved EOS grid generated via cubic spline interpolation. In detail we used natural cubic splines for interpolations along isotherms because the spacing between the original density grid is small, unlike the original temperature grid that has certain jumps, e.g., for helium with 6 kK, 10 kK, 20 kK. Here we applied cubic Akima splines that do not tend to overshoot interpolating this rougher grid.

The check of thermodynamic consistency for the H-REOS.3 and He-REOS.3 is shown in Figures 9 and 10, respectively. The color code represents areas where thermodynamic consistency is ensured within a certain accuracy. For instance, a performance with errors ⩽1% is given in black, ∼10% in red, and ∼100% in yellow. As expected, most inconsistencies occur at the boundary between two EOS models at 10 kK where ionization occurs and in the interpolation regimes. Note that the SCvH model contains thermodynamic inconsistencies within their interpolation areas as well. Interestingly, we obtain red areas in particular in the transition region from 10 kK to 20 kK at low densities for both EOS tables. This can be explained directly by our criterion. This transition region is within or near the ideal gas regime but at the points where ionizations becomes increasingly important, leading to a slope in the pressure and energy isotherms; see Figures 5–8. This slope has a strong influence on the isochores with respect to the interpolation scheme that takes into account the two neighboring points for each data point. For instance, for helium we have the 6 kK, 10 kK, 20 kK, 30 kK, and 60 kK isotherms, which is a rather large spacing for interpolation and obtaining the quite sensitive derivative (∂P/∂T)_ρ; see Equation (11). This term has to compensate the (∂u/∂ρ)_T term, which is zero within the ideal gas model because its caloric EOS is independent of the density. Hence, the red areas in the ideal gas regime should be considered as intrinsic effects of our consistency criterion, producing ostensible inconsistencies. However, they are not important for the calculation of isentropes for GPs and BDs. Remarkably, these thermodynamic inconsistencies are absent at the transition region from the SCvH model to the EOS of CP in the Hydrogen EOS; see Figures 1 and 9. Since the hydrogen is fully ionized there, no significant changes in the slope of the neighboring isotherms occur.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we conclude that in those parts of the H-REOS.3 and He-REOS.3 which are relevant for calculating isentropes for GPs and BDs (see Figures 1 and 4), thermodynamic consistency is ensured to be better than 5%, in most cases better than 1%.

High-pressure experiments using static diamond anvil cells or dynamic shock compression techniques are a crucial test for each EOS. Since an extensive comparison of our DFT-MD data for hydrogen with a variety of experiments is given in Becker et al. (2013), we will focus here on the helium EOS. In contrast to hydrogen, only a few high-pressure experiments for helium are published. The cold curve at 300 K has been investigated by several static anvil compression experiments (Lallemand & Vidal 1977; Mills et al. 1980; Loubeyre et al. 1993; Mao et al. 1988; Eggert et al. 2008). We find a very good agreement with our He-REOS.3 and these experimental results. In particular, our data coincide with the experiments below 0.3 g cm⁻³. At higher densities, our pressures are systematically higher than in the experiments with a maximum deviation from the Loubeyre data in the solid of 10% at 1.5 g ccm⁻¹.

To test our EOS at higher temperatures more relevant to the BD regime, we compare to the principal Hugoniot experiments of Nellis et al. (1984), Eggert et al. (2008) and Celliers et al. (2010). Of particular interest are the gas gun shots from Nellis et al. (1984) that yielded a pressure of 16 GPa at 0.41 g cm⁻³ and ∼12, 000 K. This is an EOS point very close to the adiabat of Gliese-229b; see Section 4. The experimental data for the principal Hugoniot together with theoretical predictions from our He-REOS.3 (solid) and the DFT-MD and PIMC results of Militzer (2009; dashed) are shown in Figure 11. Note that the data of Eggert et al. (2008; circles) rely on the EOS of the reference material (quartz) used in the experiment. This EOS has been recalibrated by Knudson & Desjarlais (2009), leading to a density correction of the original helium results. Celliers et al. (2010) estimated this correction to be 10%, leading to the shift from the filled to the shaded circles in Figure 11.

Zoom In Zoom Out Reset image size

The ab initio predictions agree very nicely for lower and high pressures, but our He-REOS.3 is stiffer at intermediate pressures which yields a smaller maximum compression; see inset of Figure 11. Both theoretical curves coincide with the experimental results of Nellis et al. (1984) but only with the shifted highest compression of Eggert et al. (2008). This issue might be solved with a detailed reanalysis of their data based on the new quartz standard. We point out that the relevant experimental data for shallow depths of BDs (Nellis et al. 1984) are matched by our He-REOS.3.

The major constituents of GPs and BDs are hydrogen and helium and a small fraction of heavier elements. Hence, interior models are based on a mixture of their EOS. As we provide two separate EOS tables for hydrogen and helium, we apply the additive volume rule to obtain a linear mixture H-He EOS (LM-EOS) with the mass fraction of helium (Y) and hydrogen (X = 1 − Y):

$$\begin{eqnarray} \frac{1}{\rho _{\rm mix} (P,T)} & = & \frac{X}{\rho _{H} (P,T)} + \frac{Y}{\rho _{{\rm He}} (P,T)} \;, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} u_{\rm mix}(P,T) & = & Xu_H(P,T) + Yu_{{\rm He}}(P,T) \;. \end{eqnarray} \tag{ 13 }$$

Taking into account a representative EOS for the heavier elements (Z), such as water for GPs or a scaled helium EOS as done within the BD calculations in Section 4, Z/ρ_Z(P, T) has to be added to the right-hand side of Equation (12) and Zu_Z(P, T) to Equation (13). The LM-EOS can then be applied for arbitrary compositions X/Y/Z.

The alternative approach is to perform ab initio simulations for a mixture of hydrogen and helium (and perhaps heavier elements) with the desired fractions to obtain the EOS for given values of X/Y/Z; see Lorenzen et al. (2009). This leads to the RM-EOS which includes non-linear mixing effect in a genuine way. In a recent paper, Militzer & Hubbard (2013) provide EOS data for a real mixture with a helium fraction of Y = 0.2466. Pressure and energy isotherms of their original DFT-MD data (open triangles) and their interpolated data (circles) are shown in Figures 12 and 13 together with our LM isotherms for this Y value. Note that the 100 kK RM isotherm and all RM data below 0.2 g cm⁻³ are extrapolated. Furthermore, we do not see a first-order phase transition in the original RM data for hydrogen around 1 g cm⁻³ below 2000 K as has been predicted recently (Lorenzen et al. 2010; Morales et al. 2010b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find a remarkably good agreement of our LM-REOS.3 with the original RM data, especially in the region above 10 kK that is relevant for the interior of BDs. As an example, we show the adiabat of Gliese-229b (dashed violet curve in Figure 12). We conclude that non-linear mixing effects are not relevant for ρ–T conditions that occur in BDs but may be important for GPs as stated in Militzer & Hubbard (2013); see the next section. Therefore, we apply our LM-REOS.3 for modeling the interior of BDs.

Another important phenomenon that occurs in real H–He mixtures is the demixing of both components at high pressures. This issue has been discussed for decades (Stevenson & Salpeter 1977b; Lorenzen et al. 2011; Morales et al. 2013b). It affects the interiors of rather cold objects such as Saturn and presumably Jupiter (Stevenson & Salpeter 1977a; Fortney 2004). Caused by demixing, helium droplets can form and sink into the planet, thereby releasing gravitational energy as an additional energy source. This process can explain the luminosity problem of Saturn (Fortney & Hubbard 2003), where homogeneous evolution models predict an age of ∼2.5 Gyr in contradiction to the age of the solar system of 4.56 Gyr. This problem could be solved alternatively assuming a layered convection for the interior structure of Saturn as proposed by Leconte & Chabrier (2013).

However, a demixing effect is not included in a linear mixture of our EOS data. Therefore, we predict a thermodynamic and mass range for homogeneous objects consisting of hydrogen and helium with a solar helium fraction where demixing will occur in their interiors. The results are based on the demixing regions calculated from ab initio simulations by Lorenzen et al. (2011) and Morales et al. (2013b), who predict the phase separation above 1 Mbar for temperatures below 8000 K. These conditions are fulfilled in the interiors of objects with a mass >0.125 M_Jup and a temperature at the 100 bar level of T(100 bar) ⩽ 700 K. For comparison, the temperatures at 100 bar for the massive GP KOI-889b and the BDs Corot-3b, Gliese-229b, and Corot-15b, which we will investigate in Section 4, are 1100 K and 2500 K, well above the demixing region.

The REOS.3 tables were built with the main purpose of calculating interior models and MRR for BDs. Due to the lack of observational constraints on the internal structure other than mass and radius, we assume a BD to be a homogeneous sphere. Since there is an overlapping mass regime of BDs and GPs (Hennebelle & Chabrier 2008; Padoan & Nordlund 2004), we also assume that the GP KOI-889b consists of a homogeneous layer, in contrast to the three-layer model applied for Jupiter; see the previous section.

The pressure-density profile of BDs and GPs is given by an adiabat that is fixed by an atmospheric boundary condition P_at; see Section 4.2. The adiabat is derived from a linear mixture EOS with certain amounts of hydrogen (X), helium (Y), and heavier elements (Z). The interior profiles of an object with given mass M and radius R are derived by integrating the two following structure equations using a fourth-order Runge–Kutta scheme:

$$\begin{eqnarray} \frac{dr}{dm} & = & \frac{1}{4\pi r^2\rho } \;, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \frac{dP}{dm} & = & -\frac{Gm}{4\pi r^4} \;. \end{eqnarray} \tag{ 15 }$$

G is the gravitational constant and the Lagrangian coordinate m is the mass within a shell of thickness r. The boundary conditions are given by r(0) = 0, r(M) = R, and P(M) = P_at. The aim of these calculations is to study the influence of the EOS on the interior profiles by comparing results using the new REOS.3 tables and the SCvH EOS tables. Therefore, we do not consider the rather complex evolution of these objects with their characteristic early fusion processes; see Burrows et al. (2001). For simplicity, our atmospheric profiles assume the same age of 5 Gyr for all objects studied here.

There are several BDs and many massive GPs with measured mass and radius (see Figure 16), from which we selected the GPs KOI-889b (Hébrard et al. 2013) and WASP-18b (Maxted et al. 2013) and the following BDs: Corot-3b (Deleuil et al. 2008), Corot-15b (Bouchy et al. 2011b), WASP-30b (Anderson et al. 2011), KOI-205b (Díaz et al. 2013), Kepler-39b (Bouchy et al. 2011a), and LHS-6343c (Johnson et al. 2011). We also plotted two predictions for the mass and radius of the first bona fide BD Gliese-229b (Nakajima et al. 1995) made by Marley et al. (1996) and Allard et al. (1996).

In the following, we will study four of these objects with increasing mass in detail: KOI-889b, Corot-3b, Gliese-229b, and Corot-15b. We use the input parameter mass M in Jovian masses M_J, radius R in Jovian radii R_J, and metallicity Fe/H of the host star from the literature as listed in Table 5.

Table 5. Literature Data of our Studied Objects

Object	Mass	Radius	Fe/H
	(MJ)	(RJ)
KOI-889b	9.98 ± 0.5	1.03 ± 0.06	−0.07 ± 0.15
(Hébrard et al. 2013)
Corot-3b	21.66 ± 1	1.01 ± 0.07	−0.02 ± 0.06
(Deleuil et al. 2008)
Corot-15b	63.3 ± 4.1	1.$12_{-0.15}^{+0.3}$	0.1 ± 0.2
(Bouchy et al. 2011b)
Gliese-229ba	46.$2_{-14.8}^{+11.8}$	0.$87_{-0.07}^{+0.11}$	−0.2 ± 0.4

Note. ^aNote that mass and radius for Gliese-229b are derived from the fitting formulae given in Marley et al. (1996) while Fe/H = −0.2 ± 0.4 is taken from Schiavon et al. (1997).

Download table as:  ASCIITypeset image

BDs and GPs host complex and actively circulating molecule-dominated atmospheres and, depending on the temperature, grains that lead to dusty atmospheres; see Burrows et al. (2011), Showman & Kaspi (2013), Zhang & Showman (2014), Basri (2000), and Baraffe (2014) for reviews.

Measured emitted spectra of the objects are generally necessary to obtain realistic atmospheric profiles, but evolutionary models can also be used to assess atmospheric structure. For Gliese-229b, where accurate spectra have been obtained, but the mass and radius are poorly known, (Allard et al. 1996; Marley et al. 1996; Oppenheimer et al. 1998; Saumon et al. 2000), we use the atmospheric profile from Marley et al. (1996); see Figure 15. Since we focus on the effect of the EOS on the interior profile, we make some simplifying assumptions for the other objects, KOI-889b, Corot-3b, and Corot-15b. We neglect irradiation from the host star and set their age to 5 Gyr. Given the measured masses, we extracted the respective effective temperatures from the models of Baraffe et al. (2003). Based on the measured surface gravities and estimated effective temperatures for these objects, we used the radiative-convective atmosphere model of Marley et al. (1996) and Fortney et al. (2008) to calculate the temperature structure, and the depth at which the radiative atmosphere becomes convective and stays convective. This is the tie point for the deep atmospheric structure.

Zoom In Zoom Out Reset image size

The resulting atmospheric profiles are shown as solid lines in Figure 15. We show the transition points, where the temperature profiles become adiabatic. These points define the atmospheric boundary condition P_at as well. We obtain for Corot-15b P(2000 K) = 40 bar, for Gliese-229b P(1800 K) = 52 bar, for Corot-3b P(1500 K) = 74 bar, and for KOI-889b P(1000 K) = 58 bar. The dashed lines in Figure 15 represent the respective adiabats, which determine the interiors. While this is a simple prescription, since we are only interested in interior differences between EOS given reasonable upper boundary conditions, this treatment certainly suffices.

We have calculated MRR for homogeneous and adiabatic objects with an atmospheric boundary like Gliese-229b of P(1800K) = 52 bar, a helium content of Y = 0.27, and a metallicity of Z = 2%. The results for the different EOS are shown as green curves (solid: REOS.3, dashed: SCvH) in Figure 16. We obtain systematically higher radii using the REOS.3 with a maximum deviation from the SCvH-EOS of ∼6% at 1 M_J. For masses above 20 M_J, the deviation remains nearly constant at ∼2.5%. This is in agreement with Militzer & Hubbard (2013), who obtained slight radius enhancement using their ab initio EOS data compared to the SCvH EOS for planets with masses between 0.5 and 2 Jovian masses. This finding can now be extended up to 70 Jovian masses and interpreted directly in terms of the EOS.

Zoom In Zoom Out Reset image size

The main deviations between REOS.3 and SCvH arise within the WDM regime, where strong correlations and quantum effects are important and dissociation and ionization processes occur. There, the DFT-MD EOS data are not as compressible as the SCvH ones; see the lower maximum compression at the principal Hugoniot curve in Section 1. Since the objects with lower masses around 1 M_J are dominated by WDM, the more compressible SCvH data lead to significantly smaller radii than the REOS.3 data. At higher densities and/or temperatures the system is fully ionized and becomes increasingly degenerate, which is accurately described by both EOS leading to similar MRR.

The next aim was to investigate the variation of the radius with respect to the metallicity. For simplicity, the heavier elements (Z) are represented by a scaled helium EOS (ρ_Z(P, T) = 4 × ρ_He(P, T), u_Z(P, T) = u_He(P, T)/4). We increased the fraction of Z within the mixture and calculated the respective radii for the mean value of the mass of our four selected objects and the maximum errors of these masses; see Table 5. Each object has a different adiabat, fixed by the different atmospheric boundary; see Figure 15.

The results are shown in Figure 17, where the metallicities are given in terms of the metallicity of the host star. The horizontal solid lines represent the mean values and the color shaded areas indicate the error bars of the objects observed radii.

Zoom In Zoom Out Reset image size

The solid black lines represent our calculated radii with respect to the mean value of the mass using the REOS.3 data. Dashed black lines indicate the radii for the maximum observed uncertainties in total mass. The red solid and dashed lines illustrate the corresponding results using the SCvH-EOS.

Recall that for every object we try just one atmospheric boundary, meaning one tie point in P–T for the interior adiabat. This allows us to readily compare how different EOS affect the radius, while the remaining parameters are not changed.

In all cases, we find a larger radius for the REOS.3 compared to SCvH. For Gliese-229b, the REOS.3 curve for the mean mass intersects the mean radius at a metallicity of Z = 1.56%. The SCvH radii have no value for the mean radius but are well located within its error bars. Due to the large error bars in mass, REOS.3 and SCvH results overlap and there are many possible metallicities up to 9% for the mean radius. For the other three objects, the accuracy of the metallicity derivation is, of course, tied to the simple choice of the atmospheric boundary condition, but the trends in radius with metallicity are clear. In the case of KOI-889b, we find disjoint solutions with respect to the underlying EOS. The REOS.3 result of Z ∼ 3% for the mean values of mass and radius is again slightly larger than the SCvH result. The curves for Corot-3b are disjoint as well. Here the temperature dependence is nicely illustrated. For our chosen atmospheric conditions of T(74 bar) = 1500 K we find no solution for the mean radius. However, with a slightly higher initial temperature of T(52 bar) = 1800 K, as for Gliese-229b, we obtain a metallicity of ∼2% with REOS.3 (dotted line in the Corot-3b frame). In the case Corot-15b, we fail to reach even the shaded error bar region of its radius with both EOS. This BD is very large for its mass compared to WASP-30b or LHS-6343c; see Figure 16. Such a size can be explained by assuming a younger object age (since objects contract with time) or within advanced evolutionary calculations that assume a more opaque atmosphere slowing the cooling and contraction of the object; see Burrows et al. (2011).

Finally, we present interior profiles of Gliese-229b and Corot-3b calculated with Z = 2% and of KOI-889b with a metallicity of Z = 4% for the given masses and appropriate atmospheric conditions; see above. The results are shown in Figures 18 and 19. All these figures contain curves for the temperature in units of 10 kK (red), the pressure in Gbar (black), and the density using REOS.3 data (solid) and the SCvH EOS (dashed). The more massive an object, the higher the central values; see the right panel of Figure 18 for objects with a composition and an atmospheric boundary like Gliese-229b. The respective MRR of these objects is shown in Figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While for the selected GP and BDs the central temperatures are similar but slightly higher using REOS.3, the central pressures of the SCvH results are higher by ∼10% due to the higher densities. Here one can see that the pressure is more influenced by the density than by the temperature as typical for degenerate matter. The higher pressures and densities of the SCvH results are a consequence of the predicted smaller radii, hence the objects are more compact.