MISCIBILITY CALCULATIONS FOR WATER AND HYDROGEN IN GIANT PLANETS

The past decade has been a period of extraordinary discoveries in the field of exoplanets.³ For the first time, we now have quantitative estimates for occurrence rates of different types of planets in our galaxy. A recent study by Petigura et al. (2013) that focused on planets with periods up to 100 days in the Kepler sample demonstrates the prevalence of Neptune and sub-Neptune exoplanets. In this sample, the planets with a radius larger than 1.5 Earth radii have a mean density (Weiss & Marcy 2014) that implies that they are composed of both heavy elements (rocks and ices) and gas (hydrogen and helium; for the mass–radius relation, see Seager et al. 2007). Despite the absence of sub-Neptune planets in our solar system, we may be able to place constraints on the formation process of ice giants by studying the interior and the evolution of Uranus and Neptune. Both planets have sizable gaseous envelopes surrounding cores composed of heavier elements.

The core accretion model (Pollack et al. 1996; Bodenheimer & Lissauer 2014; Helled & Bodenheimer 2014) suggests that the giant planets form in two distinct phases. First, rocky and icy planetesimals accumulate to form a dense core. Once a critical core mass has been reached, a runaway accretion of hydrogen–helium gas sets in and lasts until all gas from the planet's neighborhood has been depleted. Most often it has been assumed that very little gas is present when the initial core forms (Pollack et al. 1996). Because of this assumption, planetary interior models typically assume a dense core of rock and ice with a sharp transition to a gaseous outer envelope that is composed of hydrogen, helium, and a small fraction of heavier elements (Stevenson 1985; Hubbard et al. 1995; Guillot 1999, 2005; Hubbard 1999; Militzer et al. 2008). While the heavier-element fraction in Jupiter's envelope has been measured in situ by the Galileo entry probe (Wong et al. 2004), there is no direct measurement that characterizes the state of a giant planet core. It is thus plausible that some mixing of gas and icy planetesimals has occurred when a giant planet formed and that it may have persisted until today.

The mean density of Uranus and Neptune, on the other hand, suggests that they have a much higher fraction of heavy elements than Jupiter and Saturn (Hubbard 1984). Based on cosmological abundances, it is assumed that water is among the dominant species, followed by methane and amonia, even though only trace amounts of it have been detected spectroscopically in the atmospheres of Uranus and Neptune (Encrenaz 2003). This paper focuses on water–hydrogen mixtures because of the prevalence of water in our solar system and the possibility that gases and ices mix when giant planets form.

A recent study (Fortney et al. 2013) focused on the late accretion of planetesimals when a large, gaseous envelope is present. It was predicted that even large planetesimals of 100 km in diameter may disintegrate upon entry and the material would be distributed throughout the gas envelope. Depending on the mixing properties of hydrogen and heavier elements, this material may either remain in the envelope or gradually settle onto the existing core. We will study hydrogen–water mixture in this article because little information is available about how well hydrogen mixes with other elements at high pressure. The phase separation of hydrogen and helium has been predicted to occur in the interiors of Saturn (Stevenson & Salpeter 1977; Fortney & Hubbard 2004; Morales et al. 2009; Soubiran et al. 2013) and also Jupiter (Wilson & Militzer 2010). The resulting release of gravitational energy has been named as a reason for the observed excess in Saturn's luminosity. It would be of interest to know whether a similar process involving the separation of water–hydrogen mixtures could operate in the interiors of ice giant planets. If a phase separation occurs, the envelope would be progressively depleted of heavy elements, which affects the interior structure and the luminosity of a giant planet.

On the contrary, it is possible that the core of an initially differentiated planet would dissolve into the surrounding layer of hot, dense hydrogen. Results from recent ab initio simulations predict that the rocky and icy components of the cores in Jupiter and Saturn are miscible in metallic hydrogen (Wilson & Militzer 2012a; Wahl et al. 2013; González-Cataldo et al. 2014). While a mixture state is thermodynamically preferred, it is not known how fast the cores of giant planets erode because gravitational forces counteract the advection of heavy elements. This may give rise to a semi-convective regime that has already been suggested to occur in the gas giant interiors (Leconte & Chabrier 2012, 2013).

Using recent updates of the gravitational moments of Uranus and Neptune (Jacobson 2007, 2009), revised interior models have been constructed. Helled et al. (2011) fitted a density profile to these data and showed that they could be satisfied by a single layer of hydrogen and helium with a compositional gradient of heavier elements. On the other hand, Nettelmann et al. (2013) matched the available constraints by constructing a typical three-layer model assuming a rocky core, surrounded by an intermediate layer of water in a liquid or superionic state (Cavazzoni et al. 1999; French et al. 2009; Wilson et al. 2013), and a gaseous outer envelope. While Helled's model assumes that water and hydrogen are completely miscible throughout the planet interior, Nettelmann et al. conversely assume that both fluids would not mix at the interface of the two layers at approximately 10 GPa and 2000 K. This underlines why the mixing properties of water and hydrogen are important for giant planet interiors.

Both aforementioned models assume an adiabatic, fully convective behavior for each layer. This assumption is not consistent with the observed luminosity of Uranus (Pearl et al. 1990; Pearl & Conrath 1991). Podolak et al. (1990) suggested a stably stratified interior model instead. Nevertheless, a convective layer of a conducting material is needed to sustain a magnetic dynamo and to produce the quadrupolar fields observed for Uranus and Neptune (Ness et al. 1986, 1989; Connerney et al. 1987, 1991). Water is assumed to be the dominant species in the layer where magnetic fields of ice giant planets are generated. If water and hydrogen mix, a gradient of composition may be introduced into the planet's interior during formation. It would thus be possible for semi-convection to operate in ice giant planets, which would reduce the overall heat flux and allow for a vigorous convection in some layers.

Some experimental data for water–hydrogen mixtures are available. Seward & Franck (1981) studied mixtures up to 0.25 GPa and observed a phase separation below 650 K. Based on this result, the presence of water clouds in the deep atmosphere of the ice giants has been inferred (Fegley & Prinn 1986). More recently, Bali et al. (2013) investigated the properties of water–hydrogen in the presence of minerals under deep-Earth conditions up to 2.5 GPa. Hydrogen was released due to chemical reactions between water and the minerals. The analysis of micro-inclusions of fluid within the minerals showed that phase separation of hydrogen and water was possible at temperatures below 1200 K. These findings favor differentiated interior models for Uranus and Neptune.

At the present time, no experimental data are available for higher pressure and temperature, but materials under such conditions can be studied efficiently with the ab initio simulations that we used throughout this paper. Results from such simulations have been shown to agree very well with shock wave measurements for hydrogen (Lenosky et al. 1997; Knudson et al. 2004; Loubeyre et al. 2012). Similarly, results of shock wave experiments of water have been matched with ab initio simulations (Knudson et al. 2012). Recently, improvements have been made to compute the dielectric constant of water with such methods (Pan et al. 2013). Here we use the same technique to study water–hydrogen mixtures from 2 to 70 GPa and from 1000 to 6000 K (see Figure 1). For various mixing ratios, we compute the equation of state (EOS), in particular the pressure and internal energy as a function of density and temperature. In addition, we derive the Gibbs free energy of mixing by performing a thermodynamic integration. We show that the mixtures behave closely to an ideal mixture. The entropy of mixing is the dominant term in the Gibbs free energy of mixing, indicating that no phase separation occurs at all conditions under consideration.

We also analyze changes in ionic species that are present in the mixture as a function of the pressure and the temperature. Finally, we comment on some implications of computed miscibility properties. We show, for instance, that a sharp transition from a hydrogen-rich envelope to a water-rich phase is thermodynamically unstable and that a three-layer picture for the icy giants is a simplification. Furthermore, we compare the mass–radius relationship for Neptune-like and sub-Neptune planets depending on their differentiation. We found that for a given radius and mass, a fully differentiated planet has a much lower hydrogen content than a homogeneously mixed planet.

Our investigations of the hydrogen–water mixture rely on ab initio molecular dynamics (MD) simulations in which the nuclei are treated as classical particles while the electrons are considered quantum mechanically using density functional theory (DFT; Hohenberg & Kohn 1964; Kohn & Sham 1965). We used the Vienna Ab initio Simulation Package (Kresse & Furthmüller 1996). The time step of the MD simulations was set to 0.2 fs, which is short enough to describe the molecular vibrations accurately. Each simulation was performed for a minimum of 0.5 ps and up to 7 ps for the lowest temperatures. To keep the temperature constant, we employed a Nosé thermostat (Nosé 1984, 1991). For the electrons, we used a Fermi–Dirac distribution within a finite temperature scheme (Mermin 1965). We used projector augmented wave pseudopotentials (Blöchl 1994) with a cutoff radius of ${r}_{\mathrm{cut}}=0.8\;{a}_{0}$ for hydrogen and ${r}_{\mathrm{cut}}=1.1\;{a}_{0}$ for oxygen. We used a plane-wave basis energy cutoff at 1100 eV. To sample the Brillouin zone, we used the Baldereschi point (Baldereschi 1973), except for pure water, for which a 2 × 2 × 2 Monkhorst–Pack K-point grid (Monkhorst & Pack 1976) was used. We performed a few simulations with finer K-point grids and found consistent results within the statistical error bars. For the exchange-correlation functional, we chose the generalized gradient approximation of the Perdew, Burke, and Ernzerhof (PBE) type (Perdew et al. 1996) because it gave reasonable results for pure hydrogen (Caillabet et al. 2011; Loubeyre et al. 2012) and pure water (French & Redmer 2009; Knudson et al. 2012). We also verified our predictions by performing additional simulations with the van der Waals density functional by Dion et al. (2004) and Klimeš et al. (2011). The resulting mixing properties were in agreement with our PBE predictions at 1000 K. A more detailed comparison of the PBE and van der Waals functionals is given in the recent work by Santra et al. (2013).

To explore the phase separation at a given pressure and temperature, it is necessary to determine the Gibbs free energy. Standard MD only provides the internal energy and the pressure, but not the entropy of the system. Therefore, we performed a thermodynamic integration using an auxiliary classical pair potential (de Wijs et al. 1998; Morales et al. 2009; Wilson & Militzer 2010, 2012a, 2012b; McMahon et al. 2012; Wahl et al. 2013) because it provides the Helmholtz free energy difference ${F}_{1}-{F}_{0}$ between two systems characterized by two different potentials ${U}_{0}(\{{\boldsymbol{r}}\})$ and ${U}_{1}(\{{\boldsymbol{r}}\})$:

$$\begin{eqnarray}&&{F}_{1}-{F}_{0}={\displaystyle \int }_{0}^{1}{\langle {U}_{1}-{U}_{0}\rangle }_{\lambda }d\lambda ,\end{eqnarray} \tag{ 1 }$$

where the parameter λ defines a hybrid potential ${U}_{\lambda }={U}_{0}+\lambda ({U}_{1}+{U}_{0})$. The $\langle .{\rangle }_{\lambda }$ means that we take the average over a trajectory computed with the potential ${U}_{\lambda }$. We performed the integration in two distinct steps. First, we integrated between the potential as given by the DFT, ${U}_{\mathrm{DFT}}$, and a set of classical pair potentials, ${U}_{\mathrm{cl}}$, that we constructed by matching the forces on configurations taken from a DFT trajectory that was computed beforehand for a particular temperature, density, and composition .

Zoom In Zoom Out Reset image size

While pair potentials provide a sufficiently good description of water at high pressure for our thermodynamic integration technique to work well (Wilson & Militzer 2012a; Wilson et al. 2013), special care needs to be taken to describe molecular hydrogen. Pair potentials have a deep minimum that represents the intramolecular binding. If such attractive pair potentials are used without any repulsive many-body term, unphysical chains and clusters of hydrogen nuclei form. In Militzer (2013), a repulsive many-body potential was constructed and a stable thermodynamic integration technique was obtained for molecular and partially dissociated hydrogen. Since here we are dealing with water–hydrogen mixtures that lead to a more diverse set of short-lived chemical species, we pursued a different and simpler approach. We constructed sets of non-bonding pair potentials that we have developed and tested recently (Wahl et al. 2015). By removing all attractive parts from the pair potentials, we prevented the formation of unphysical clusters. Hydrogen molecules still form gradually as we switch from the non-bonding potentials to the DFT forces. We determined that using between 5 and 7 λ points was still sufficient to accurately calculate the integral in Equation (1). Examples of the evolution of this average as a function of λ are shown in Figure 2.

Zoom In Zoom Out Reset image size

To derive the Helmholtz free energy of the classical system, ${F}_{\mathrm{cl}}$, we performed another thermodynamic integration to a reference system with known Helmholtz free energy, F₀. For this paper, we used the ideal gas as reference. For this integration, we used many more λ-steps as it only requires classical simulations, which are 10⁵ times faster than DFT computations.

With this two-step integration approach, we are able to derive the Helmholtz free energy for the DFT system at different densities, temperatures, and concentrations. Adding the PV term, we obtain the Gibbs free energy, ${G}_{\mathrm{DFT}}={F}_{\mathrm{DFT}}+{P}_{\mathrm{DFT}}V$, where ${P}_{\mathrm{DFT}}$ is the pressure given by the DFT calculation and V the volume of the simulation cell.

We considered different concentrations of water molecules, x,

$$\begin{eqnarray}&&x=\displaystyle \frac{{N}_{{{{\rm H}}}_{2}{{\rm O}}}}{{N}_{{{{\rm H}}}_{2}{{\rm O}}}+{N}_{{{{\rm H}}}_{2}}},\end{eqnarray} \tag{ 2 }$$

where N_i is the number of molecules of type i in the simulation cell. We performed simulations for the following ${N}_{{{{\rm H}}}_{2}{{\rm O}}}$:${N}_{{{{\rm H}}}_{2}}$ ratios: 48:0, 40:16, 32:32, 24:48, 16:64, and 0:92. Except when noted otherwise, the concentrations refers to the total contents in the simulation cell regardless of what chemical species form at various pressures and temperatures. In Figure 3, we show a snapshot from a simulation for x = 0.5 and ${N}_{{{{\rm H}}}_{2}{{\rm O}}}$:${N}_{{{{\rm H}}}_{2}}$ = 32:32, which implies that 32 oxygen and 128 hydrogen atoms were present in the simulation cell. When we varied the water concentration, we always replaced one water molecule with two hydrogen molecules so that the volume of the simulation cell would not change too much. We performed simulations in the pressure range from 2 to 70 GPa along four isotherms at 1000, 1500, 2000, and 6000 K.

Zoom In Zoom Out Reset image size

From the DFT MD and the thermodynamic integration, we extracted pressure, internal energy, and Gibbs free energy as a function of temperature, density, and concentration. The results are given in Table 1.

In Figures 4–6, we plot various thermodynamic quantities for different pressures and concentrations at 6000 K. We used a spline interpolation for the pure water and pure hydrogen system. We compared our simulation results for the mixtures with an ideal mixing approximation using our simulation results for H₂ and H₂O and an additive volume law at constant pressure and temperature. The ideal mixing approximation was found to reproduce our simulation results for the mixtures well. Some small deviations of the order of a few percent (up to 10% locally), in particular for the density and for the internal energy, can be identified, however. Nevertheless, we conclude that the ideal mixing approximation is robust and will be sufficiently accurate for the construction of most planetary interior models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We computed the Gibbs free energy of mixing:

$$\begin{eqnarray}{{\rm \Delta }}G(x,P,T) & = & G(x,P,T)-x\;{G}_{{{{\rm H}}}_{2}{{\rm O}}}(P,T)\\ & & -(1-x){G}_{{{{\rm H}}}_{2}}(P,T),\end{eqnarray} \tag{ 3 }$$

where G is a Gibbs free energy per molecule and G_i is the Gibbs free energy of the pure system of molecule i. A homogeneous mixture is the thermodynamically preferred state if

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\partial }^{2}{{\rm \Delta }}G}{\partial {x}^{2}}\right|}_{P,T}\geqslant 0,\end{eqnarray} \tag{ 4 }$$

which implies that ${{\rm \Delta }}G$ is a convex function of x at a given pressure and temperature. If the Gibbs free energy difference shows a partial or a full concavity, the homogeneous system is unstable and undergoes a partial or complete phase separation.

In Figure 7, we plot the Gibbs free energy of mixing for different temperature and pressure conditions. Most of the curves are convex, and when they are not, a convex curve can still be drawn within the error bars. We therefore conclude that, for the whole set of conditions we explored, a homogeneous mixture of hydrogen and water in any proportion is the thermodynamically stable state. We were not able to identify an indication for a phase separation at any condition that we explored.

Zoom In Zoom Out Reset image size

For a given pressure and temperature, we can split the Gibbs free energy of mixing into three terms, ${{\rm \Delta }}G={{\rm \Delta }}E+P{{\rm \Delta }}V-T{{\rm \Delta }}S$. The internal energy term, ${{\rm \Delta }}E$, represents the interaction between the different species in the mixture. The $P{{\rm \Delta }}V$ term measures deviations in density from an ideal mixture. Finally, $-T{{\rm \Delta }}S$ is the full entropy of mixing including ideal and non-ideal contributions.

From the example at 6000 K and 70 GPa shown in Figure 8, we can infer that the contributions to the Gibbs free energy of mixing from the internal energy and the $P{{\rm \Delta }}V$ term are small, while the entropy is the dominant term by far. The computed entropy can be very well approximated by the ideal entropy of mixing of water and hydrogen molecules ${S}_{\mathrm{id}}=-{k}_{{{\rm B}}}[x\mathrm{ln}x+(1-x)\mathrm{ln}(1-x)]$. We find this approximation to work very well even if the system is almost fully dissociated at, e.g., 6000 K and 70 GPa (see Figures 9 and 10). If one tries to use the entropy of mixing for two atomic systems instead, the agreement with the simulation results is inferior (Figure 8). We conclude that at high temperature when many short-lived species are present, the attraction between ionic species is still sufficiently large so that the molecular entropy of mixing is a much better approximation than relying on a mixture of atoms.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We performed a systematic study of the different contributions in the Gibbs free energy of mixing and always found that the entropy is the dominant term that can be matched well with an ideal mixing approximation for molecules. The largest deviations from this approximation arise for high water concentration. The presence of hydrogen appears to slightly alter the dissociation fraction of water molecules in the mixtures. This is the strongest non-ideal mixing effect that we identified. Similarly, the presence of helium atoms appears to increase the stability of hydrogen molecules when various hydrogen–helium mixtures are compared for a given pressure and temperature (Vorberger et al. 2007).

In addition to studying various thermodynamic functions, we also determined the different chemical species present in the mixture for the different pressure, temperature, and concentration. We used a similar species analysis method to that employed by Vorberger et al. (2007). At each time step, we computed the distances between the nuclei. If a pair of nuclei remained closer than a given distance for a minimum duration (10 times the vibration period of molecular hydrogen ${\tau }_{{{{\rm H}}}_{2}}$ = 7.6 fs), we considered the two nuclei bound. From pairs of bound nuclei, we built larger chemical species if the bonds remained contiguous. Using the first minimum in the radial distribution functions (Soubiran & Militzer 2014), the following distance limits were derived: ${l}_{{{\rm H}}-{{\rm H}}}=1.0$ Å, ${l}_{{{\rm O-H}}}=1.3$ Å, and ${l}_{{{\rm O}}-{{\rm O}}}=1.9$ Å.

Figures 9 and 10 show examples of the chemical compositions of the mixture. We observed that H, H₂, O, HO, and H₂O make up more than 99% of the observed species in the system. A few transient bigger molecules were identified also. The temperature has a strong impact on the degree of dissociation, as expected. While at 2000 K and below the mixture is almost fully molecular, at 6000 K we observe a system that is nearly fully dissociated. We thus expect a much higher electric conductivity at 6000 K because the dissociation is generally associated with ionization, especially if hydrogen is present (Collins et al. 2001).

For both temperatures, the fraction of dissociated species appears to increase with pressure. Water molecules tend to split into hydroxide and hydrogen ions or even to fully dissociate, entering into the regime where no stable chemical bonds exist.

The computed EOS of the homogeneous mixtures shows only small deviations from the ideal mixing approximation. We thus conclude that an ideal mixing law may provide quite reasonable estimates for the purpose of planetary modeling in general as long as no phase separation is expected and that very accurate EOSs for the pure systems are used.

Moreover, we do not find a phase separation for the hydrogen–water mixtures we studied. As was observed for other mixtures at higher pressure (Wilson & Militzer 2012a; González-Cataldo et al. 2014), the entropy is nearly ideal and contributes the most to the free energy of mixing, which stabilizes the homogeneous mixture.

We have to stress here that our miscibility predictions differ slightly from what has been inferred from high-pressure experiments by Bali et al. (2013; see Figure 1 for comparison). In these experiments, samples of different minerals were compressed up to a 2.5 GPa in the presence of water. Due to chemical reactions, hydrogen was released. In some conditions a phase separation of water and hydrogen occurred, even for temperatures above 1000 K. While these experiments may be relevant for the Earth's interior, the conditions are quite different from the envelope of the giant planets, where no minerals are present. We therefore suggest that diamond anvil cell experiments of water–hydrogen mixtures in a mineral-free environment should be performed.

With knowledge of the experimental work by Bali et al. (2013), Aranovich (2013) constructed a van Laar mixing model using the molar volume of the pure species and a van Laar parameter fitted on earlier experiments at the kbar regime. This model predicts the hydrogen–water mixture to phase separate at even higher temperatures than suggested by Bali et al. (2013). When we constructed a van Laar mixing model, we were not able to obtain a good fit to our ab initio Gibbs free energies, even though the molar volumes in our simulations of the pure species agree quite well with the values that Aranovich used. To match our ab initio data, the van Laar model would need to be extended by including non-ideal mixing effects for the volume, and the van Laar mixing parameter W would need to be made temperature and pressure dependent.

The fact that our results indicate the water–hydrogen mix at 2–70 GPa has implications for our understanding of ice giant planets. The generic three-layer model (Nettelmann et al. 2013) with a hydrogen-rich outer envelope, a water-rich intermediate layer, and a rocky core may be a simplification. The boundary between the hydrogen and water layers is at about 10 GPa and 2000 K in Uranus and Neptune, which falls into the parameter regime that we explored with our simulations. Yet our results show that such a sharp boundary between the two layers is thermodynamically unstable. Convection will efficiently mix the two layers unless the system approaches a semi-convective state.

Therefore, we cannot rule out the three-layer picture completely because the timescale of the mixing has to be taken into account. In the core accretion model (Pollack et al. 1996), solid planetesimals are accreted first, followed by a runaway gas accretion. The core is assumed to differentiate quickly, letting the lighter components rise to the surface of the core. Based on our results, water would start to mix with the gas, but the timescale is then a key parameter. If the mixing is only diffusive, the planet remains in a nearly differentiated state because the diffusive timescale of water in hydrogen can be as long as 10¹² yr based on the diffusion coefficient estimates in Soubiran & Militzer (2014) and the size of the envelope given by Nettelmann et al. (2013). On the other hand, if there is a vigorous and sustained convective activity, then in a few convective timescales (on the order of a 100 yr timescale; Hubbard 1968b) the planet should be homogenized, which would be incompatible with the multiple-layers assumption. Nevertheless, Leconte & Chabrier (2012, 2013) showed that, depending on the conditions, semi-convection may set in. If a gradient of composition exists in the interior of a planet, gravity may prevent an efficient convection of heavy elements. In this case, a series of alternating diffusive and convective layers are predicted to occur. Although the long-term evolution of such a semi-convective state is not fully understood, it is a possible mechanism for maintaining fairly steep compositional gradients. The planet would then have a water-rich deep inner envelope and a gradual, semi-convective transition to a hydrogen-rich envelope.

Finally, we want to discuss the importance of the water–hydrogen miscibility for exoplanet interior modeling. Given a mass and radius observation for a sub-Neptune planet, we find that the inferred hydrogen fraction depends significantly on whether water and hydrogen occur in mixed form. In Figure 11 we compare the interior properties of hypothetical planets composed of water and hydrogen only. While Valencia et al. (2010) found very similar mass–radius relationships for homogeneously mixed and differentiated, two-layer planets composed of iron and silicates, here we find much larger deviations because hydrogen and water have very different compressibilities. For a planet of 10 Earth masses, we find that the radius varies between 2.6 and 10.5 Earth radii for pure water and pure hydrogen planets. Respectively, the central pressure varies between 5.8 and 0.07 Mbar.

Zoom In Zoom Out Reset image size

In Figure 11, we also compare the hydrogen fraction for fully mixed and differentiated water–hydrogen planets for a given radius. If a planet with 10 Earth masses and 4 Earth radii were detected, one would infer a hydrogen fraction of only 8% (0.8 Earth masses) if one assumed a differentiated interior. The central pressure would be 5.6 Mbar, and the pressure at the water–hydrogen boundary would be 0.12 Mbar, which is far below the molecular-to-metallic transition pressure in pure hydrogen (Vorberger et al. 2007). On the other hand, if one assumed a homogeneously mixed interior structure, the hydrogen mass fraction would increase to 25% (2.5 Earth masses), and the central pressure would decrease to 1.4 Mbar, which is sufficiently high for hydrogen molecules to dissociate. This reaction would also increase the electrical conductivity of the mixture and further the generation of magnetic fields.

The inferred hydrogen fraction increases significantly if a mixed interior is assumed because hydrogen gas is very compressible and then some hydrogen fluid is exposed to much higher pressure. For differentiated planets of 10 Earth masses, we find that hydrogen is never exposed to more than 0.31 Mbar regardless of the planet's radius. The central pressure was found to decrease with increasing hydrogen content because water is diluted with a material of smaller density.

A difference of 0.8 and 2.5 Earth masses in hydrogen contents also has implications for our understanding of the environment where the planet formed.

Using DFT MD simulations, we showed that in the range of pressure from 2 to 70 GPa and temperature from 1000 to 6000 K, the water–hydrogen mixtures behave closely to an ideal mixture. We found that no phase separation is expected for this parameter range. A homogeneous mixture is always thermodynamically preferred for all concentrations. Our simulation results are in disagreement with water–hydrogen mixture experiments in mineral cells (Bali et al. 2013) in the lowest pressure–temperature regime explored by our calculations. We would suggest that high-pressure experiments of hydrogen–water mixture should be performed in a mineral-free environment in order to better constrain the phase diagram in this regime.

Since we predict no phase separation for water and hydrogen, this has consequences for the ice giant planets. If a planet already has a mixed water–hydrogen layer, the thermodynamic properties, in particular the entropy of mixing, prevent any differentiation from occurring.

If a planet has two separate water and hydrogen layers, this implies two things. First, the icy materials must have been delivered early on when very little gas was present, or the icy planetesimals must have been sufficiently large so that they penetrated through the existing atmosphere. Both possibilities would be consistent with the core accretion model.

Furthermore, such a planet could not be fully convective; otherwise, water and hydrogen would mix assuming that the pressure at the interface is in the 2–70 GPa range. It is possible, however, for the planet to remain predominantly differentiated because the interior may assume a semi-convective state, in which water and hydrogen would not mix efficiently. Because the long-term dynamics of this semi-convection is not sufficiently well understood, it is difficult to model the evolution of planetary interiors on the gigayear timescale when compositional gradients are present.

We also showed that, unlike in the case of rock–iron mixtures, the mixing of water and hydrogen can drastically increase the estimates of the hydrogen content of sub-Neptune exoplanets. This effect has to be taken into account when their interior structure and evolution are modeled.

This work was supported by NASA and NSF. Computers at NAS and NCCS were used.