Properties of hydrogen, helium, and silicon dioxide mixtures in giant planet interiors

Gas giant planets, such as Jupiter and Saturn, are predominantly comprised of hydrogen and helium, existing in various phases. The core accretion model for giant planet formation requires the accretion of heavy rocky or icy core materials before the planets can further accrete gas.¹ Even after gas have been accreted, the planet also continues to capture planetesimals, but their fate remains unclear as they may dissolve within the gaseous envelope or reach the core of the planet.² However, the present-day distribution of these heavy elements in gas giants is not completely understood. Furthermore, the planet may not be completely differentiated into a core and an envelope. Ab initio Gibbs free energy calculations using thermodynamic integration (TDI) showed that all the main species—H₂O, SiO₂, MgO, and Fe—that could be found in the core are miscible in the metallic hydrogen that composes the deep region of the envelope and engulfs the core.^3–6 Thus, the core could be progressively eroded by the surrounding H-He envelope, but the time scale of this erosion process remains uncertain. If the time scale is short, then the cores of Jupiter and Saturn would most likely be completely dissolved today, resulting in a heavy-element-enriched H-He fluid envelope. In contrast, if the erosion time scale is long compared to the age of the planet, the core remains mostly intact, resulting in a H-He envelope of near solar composition.

Except for oxygen, the Galileo entry probe measured a heavy element concentration (metallicity) in the outer envelope of Jupiter that is three times higher than solar metallicity,⁷ which is an indication that heavy elements are mixed in with the H-He fluid. The Cassini and Juno missions will measure the gravitational fields of Saturn and Jupiter with high precision, which will help to constrain the mass distribution in the planet. Combining insights from such data with improved planetary models, it is possible to infer the distribution of the heavy elements and to constrain the size of a core.^8–12 An important step forward has already been achieved by computing more accurate equations of state (EOS) for the relevant constituents. We have performed ab initio TDI calculations to determine the EOS of H-He mixtures¹³ and we investigated the influence of heavy materials on the thermodynamic properties of the envelope by performing simulations of mixtures of H, He, and heavy elements.¹⁴ 

The redistribution of heavy elements of the core into the envelope may also strongly depend on the dynamics of the fluid right on top of the core. Even if the erosion were fast, the mixing of the erosion products into layers directly above the core may be a very slow and difficult process due to effects of gravity. Convective forces may be too weak to dredge up heavy materials into the outer layers of the planet.¹⁵ It is thus possible that a semi-convective or even stratified layer could occur above the core, significantly slowing down the mixing of the heavy elements with the envelope.^16,17 The behavior of the heavy-element-bearing fluid can be readily determined by characterizing its transport properties. The transport properties of hydrogen and helium mixtures without heavy elements have already been partially investigated numerically,¹⁸ but the effect of the inclusion of heavy elements needs to be addressed.

In this article, we report results from ab initio simulations of multi-component mixtures of hydrogen, helium, and silicone dioxide. We focus on the dilute limit of a few percent in number of SiO₂. We explore the microscopic structure of the system using pair-correlation functions in order to identify chemical bonds. We determine the diffusion coefficients of the different species and the viscosity of the fluid. These quantities are crucial for determining the dynamics near the core-envelope boundary.

The results presented in this article rely on a series of molecular dynamics (MD) simulations based on density functional theory (DFT), using the Vienna ab initio simulation package.¹⁹ We set up the simulation in a cubic cell with periodic boundary conditions encompassing 220 hydrogen and 18 helium atoms following Militzer's set up¹³ and added from 2 to 4 SiO₂ entities. We used a 0.2 fs time step for a total duration of a minimum of 1 ps. The simulations were performed at constant volume and constant temperature using a Nosé thermostat.^20,21

At each time step, a DFT calculation was performed to determine the electronic density. We used the Kohn-Sham scheme²² at finite temperature²³ and populated the eigenstates using a Fermi-Dirac distribution. Following the previous work on H-He mixtures,¹³ we used the Perdew, Burke, and Ernzerhof (PBE) exchange correlation functional²⁴ which uses the generalized gradient approximation (GGA). To improve the efficiency of the calculation, we used projector augmented wave (PAW) pseudo-potentials²⁵ including a frozen core for oxygen and silicon atoms only. We used a 1200 eV energy cutoff for the plane-wave basis set. The number of bands was adjusted to the concentration, density, and temperature so that the spectrum of fully and partially occupied eigenstates was fully covered. We used the Baldereschi point²⁶ to sample the Brillouin zone, which was found to be sufficient in the case of H-He mixtures.¹³ We investigated the properties of the mixtures from 5000 K to 15 000 K and from 25 GPa to 2000 GPa covering the metallic region of the interior of gas giant planets nearly up to the core boundary.

From the DFT-MD trajectory, it is possible to extract information on the dynamics of the system. For instance, we can use the fluctuation-dissipation theorem and its applications on the autocorrelation functions to determine the ionic transport properties of the mixtures. The autocorrelation function of the velocity is related to the diffusion coefficient of species α by a Green-Kubo formula^27,28

where $Nα$ is the number of particles of type α and $vi(τ)$ is the velocity vector of the ith particle at time τ. The brackets represent an ensemble average over multiple origins along the full length of the simulation.²⁹ 

Similarly, it is possible to determine the viscosity by computing the autocorrelation function of the off-diagonal components of the stress-tensor

where V and T are the volume and the temperature, respectively, $kB$ is the Boltzmann constant, and $σij(τ)$ is the ij off-diagonal component of the stress-tensor, with ij in {xy, yz, xz} in Cartesian coordinates.

Because of the small number of heavy element atoms and the finite length of the simulations, the autocorrelation functions and thus the integrals can become very noisy. In order to minimize the effect of the noise, we followed Meyer et al.'s methodology³⁰ and fitted the autocorrelation function with a multi-time scale function. For the velocity autocorrelation function $γv(t)=13Nα∑iNα⟨vi(t)·vi(0)⟩$⁠, we used the following expression:

The parameters a₀, τ₀, ${μl,νl,τl,ωl}l$ were determined with a least-squares fit. We used parameters up to $lmax=1$ for H and He. For Si and O, we needed to go up to $lmax=3$⁠, but, for some cases, the multi-time scale fit did not converge and we had to employ the exponential decay only (first term in Eq. (3)). Once the parameters have been determined, the integration of Eq. (1) yields

For the autocorrelation function of the stress tensor $γσ(t)=13∑{ij}⟨σij(t)σij(0)⟩$⁠, we used a simpler functional form

where A, B, τ₁, τ₂, and ω are fit parameters. Inserting Eq. (5) into Eq. (2), we obtained the following expression for the viscosity:

As an example, we show in Fig. 1 the fit of the stress-tensor autocorrelation function of a 220 H, 18 He, 3 SiO₂ mixture at 15 000 K and 1600 GPa. The agreement between the fit and the actual autocorrelation function is excellent giving confidence on the value of the viscosity despite the high level of noise on the integral of the autocorrelation function also plotted in Fig. 1. The use of the multi-time scale allows us a much more closer fit of the autocorrelation function.

The simulations we performed cover a large range of conditions and different behaviors are expected for the mixture. For instance, at low pressure and temperature, hydrogen is molecular, but it dissociates and becomes metallic when the temperature or the pressure increases.^31,32 It has been shown that adding helium stabilized the molecular phase near the molecular-to-metallic transition.³³ In order to investigate the influence of adding SiO₂ to the system, we computed the pair-correlation function between particles of types α and β

where r is a variable that typically varies from 0 to half the size of the simulation box, $riα$ is the position vector of the $iα$-th particle of type α, and the brackets denote a time average over the duration of the MD simulation. The pair-correlation function is a standard tool to characterize the microscopic arrangement around each type of particle and can help determine the existence of molecular bonds.^34,35

Fig. 2 shows that the pair-correlation function between hydrogen atoms changes significantly as the pressure increases for a fixed temperature of 5000 K. Over the entire pressure range, there is an exclusion radius of about 1 bohr. At 25 GPa, the sharp peak at 1.4 bohr represents H₂ molecules meaning that the system is not entirely dissociated. But this peak disappears as the pressure increases indicating a complete dissociation of the molecules. This is consistent with observations in pure hydrogen³⁶ and in H-He mixtures.³³ It is very interesting to note that the peak height slightly decreases as more SiO₂ entities are inserted. Unlike helium, SiO₂ does not stabilize the H₂ molecules. At 10 000 K, as plotted in Fig. 3, the 1.4 bohr peak is completely smoothed out, which means that hydrogen is fully dissociated. At both temperatures, we see that the second peak shifts towards smaller radii as pressure increases, which is simply an effect of compression and the increase in density.

The O-H pair-correlation function, plotted in Fig. 4 for a temperature of 5000 K, shows similar features to the H-H pair-correlation function. We observe a sharp peak, around 2 bohr at 25 GPa, that progressively disappears as pressure increases. But this peak is visible longer than for H-H. This peak indicates that hydrogen bonds with oxygen to form hydroxide or water molecules. The existence of such molecules was confirmed by a cluster analysis of the simulations similar to the one we performed for hydrogen-water mixtures.³⁵ At 25 GPa, we find that 45% of the oxygen atoms are chemically bonded forming hydroxide molecules and 10% are in water molecules. These molecules progressively dissociate but are slightly more stable than H₂ molecules, as was previously predicted for water-hydrogen mixtures.^34,35 The existence of H-O bonds also explains why the peak in the H-H pair-correlation function decreases with the number of inserted SiO₂. As the number of free oxygen atoms increases, the number of H-O bonds increases breaking H-H bonds because the O-H bond is more stable than the H-H bond. This results in a dissociation of molecular hydrogen. The oxygen does not stabilize H₂ molecules, like helium. Instead, it chemically reacts with hydrogen.

The H-Si pair-correlation functions at 5000 K, plotted in Fig. 5, do not exhibit the same molecular features. There is an exclusion radius of nearly 2 bohr but the first peak does not exhibit the features of a molecular system. We can assume that hydrogen does not form stable bonds with silicon under these conditions.

The dynamics of the core erosion in giant planets can be slightly impacted by the diffusion properties of the heavy elements. If the diffusion is fast, eroded materials can mix efficiently with the surrounding H-He fluid and the erosion is likely to be rapid. But if the diffusion is slow, it may render the erosion process extremely inefficient and thus there could still be a core in giant planets despite the thermodynamic predictions favoring the complete dissolution of the core into the envelope.

To compute the diffusion properties, we calculated the autocorrelation function of the velocities as indicated in Section II B. We only report calculations for the highest concentration in SiO₂ because of the small number of heavy elements, which prevents us to have a good statistics. We do not expect, however, the diffusion coefficient to be dependent on the concentration of heavy elements at least in the diluted limit.

The diffusion coefficients of H and He are plotted in Fig. 6. They show a very smooth evolution as a function of pressure and temperature. As the temperature increases, the diffusion is easier because the thermal velocity is higher. But the diffusion coefficients decrease as the pressure increases because the particles interact more strongly at higher density. It is also interesting to note that at 10 000 and 15 000 K, we have roughly a factor of 2 difference between the diffusion coefficient of helium and that of hydrogen. But it increases to a factor of almost 4 at 5000 K. The calculated coefficients are in reasonable agreement with the diffusion coefficient for H-He mixtures as computed by French et al.,¹⁸ indicating a marginal influence of the heavy elements on the diffusion of H and He. There is also a sort of plateau in the diffusion coefficient of hydrogen below 100 GPa at 5000 K. We observed a very similar feature in the diffusion of hydrogen in hydrogen water mixtures and attributed it to the dissociation of hydrogen and water.³⁵ It is safe to assume that the dissociation of H₂ is also at play in the present case.

The calculation of the diffusion coefficients for the heavy elements is much more challenging because of the poor statistics due to the small number of heavy atoms but also because the decorrelation time is much longer as these particles move much more slowly, and we would need longer simulations to have a complete description of the autocorrelation function. Because of that, the fit is of lower quality and we were not always able to obtain a converged fit using the multi-time scale decomposition as in Eq. (3). For some simulations, we had to use a simple exponential decay. Therefore, the diffusion coefficients presented in Fig. 7 can only be considered as an order of magnitude.

The trends observed on H and He are also visible for the diffusion coefficients of Si and O despite the scatter in the curves in Fig. 7. The pressure has a negative effect on the diffusion, while temperature makes the diffusive process easier. The oxygen and the silicon seem to have similar diffusion coefficients at low temperature, but at 15 000 K, oxygen diffuses faster than silicon as expected by the difference of mass. Overall, we observe a slower diffusion of the heavy elements compared to H and He as expected. However, the diffusion is only slower by one order of magnitude at most for the range of parameters we explored.

It is usually assumed that the deep interior of giant planets is fully convecting,¹² but the possible erosion of the core could inhibit the convection^16,17 and a semi-convection pattern could set in instead. One key parameter in the stability analysis is the viscosity of the fluid. By computing the autocorrelation function of the off-diagonal components of the stress-tensor, we can determine the viscosity of the mixture as explained in Section II B.

As shown in Fig. 8, the viscosity evolves smoothly as a function of the temperature and the pressure. The range of values goes from 10⁻⁴ to 10⁻³Pa s, which indicates that the convection of the metallic phase in giant planets is most likely to be turbulent and not laminar. The expected Reynolds number is of the order of $Re∼1012$ for the convection in Jupiter.³⁷ The values that we obtained for the viscosity are in agreement with H-He mixture calculations in similar conditions.¹⁸ While the dependence in temperature is unclear at low pressure, at higher pressure, we observe that the viscosity slightly decreases as the temperature increases but the difference is not significant. Because of the size of the errorbars and some uncertainty on the goodness of the fit of the stress-tensor autocorrelation function for every conditions investigated, we are unable to identify any influence from the addition of SiO₂. In the dilute limit, the properties of the multi-component mixture are mostly the same as those of the H-He mixture. The heavy elements barely influence the bulk viscosity of the system.

Based on the results we obtained, it is possible to do an estimate of a possible time scale for the erosion of the core in a Jupiter-like planet. If we neglect the time it takes for an atom to overcome the energy barrier when it separates from the solid phase to go to the surrounding fluid, we can estimate the erosion time-scale as the time necessary for the atom to diffuse through the boundary layer of the convective cell on top of the core. Indeed, a typical convective cell is homogeneously mixed except in its viscous boundary layers. In these layers, most of the mixing is performed through diffusion. By estimating the diffusion time of the heavy elements through the boundary layer, we can determine the time it takes for the atom to reach the region of the fluids where it is actually well-mixed and thus how fast the core can dissolve into the envelope. Numerical simulations of Rayleigh-Bénard convective cells³⁸ provide scaling laws for the viscous boundary layer thickness for non-rotating systems

where L is the size of the convective cell and Re is the Reynolds number. The simulations were performed for systems that are quite different from gas giants, with lower Reynolds number than for Jupiter convection for instance (⁠$Re=1012$ for the deep interior of Jupiter³⁷) and higher Prandtl numbers (above 1 for the simulations but closer to $10−2$ for Jupiter¹⁸). But comparison with experiments in other range of parameters gives confidence for an extrapolation of this scaling law towards more extreme regimes such as the ones in giant planets. If we take $L∼107$ m for the eddy size of the convective cell,³⁷ we find a boundary layer thickness of the order of $δ∼3$ km. Based on the values displayed in Fig. 7, we can estimate the diffusion coefficient to be close to 0.1 mm²/s, which results in a typical diffusion time scale of $3×106$ years. Such a time scale would suggest a nearly complete erosion of the core for a planet of 4.5 billion years. But we assumed here a full and efficient convection. Yet, if a semi-convective pattern sets in, the efficiency of the mixing could be significantly decreased and the erosion process also inhibited. More work in this area is thus necessary.

Our ab initio simulations of H-He-SiO₂ mixtures under conditions relevant for giant planet interiors showed a smooth dissociation of hydrogen as the pressure and the temperature increase. Unlike helium, SiO₂ does not stabilize molecular hydrogen because oxygen reacts with hydrogen to form water and hydroxide molecules. The location of the dissociation and the metallization of hydrogen is of first importance because it is related to the demixing of hydrogen and helium.^39–42 If hydrogen is stabilized or is less metallic, the demixing occurs at higher pressures and, thus, deeper in the planet. A more detailed analysis of the electronic properties is, therefore, necessary to identify the effect of the heavy elements on the metallization.

The calculation of the diffusion coefficients of the different species shows that hydrogen and helium are at most only marginally affected by the presence of the heavy elements. We were able to estimate the diffusion coefficients for oxygen and silicon in the H-He metallic phase and obtained values between 0.01 and 1 mm²/s. We also obtained values for the viscosity ranging between 0.1 and 1 mPa s indicating a turbulent behavior of the convection in giant planets. These ionic transport properties are typically extremely challenging to determine in high-pressure experiments, yet crucial for the determination of the dynamics of the deep interior of a planet.

Based on our calculated transport properties, we estimate the typical time scale for heavy elements to migrate from the core to the envelope, by diffusion through the boundary layer, to be of the order of a few million years. This is short on a geological time-scale and suggests that the core of giant planets could be entirely eroded if the convection is efficient. However, if semi-convection sets in, the mixing of the heavy elements into the envelope could be significantly slowed down, which would stabilize the core. Additional work is required in this area to better characterize the dynamics of the envelope near the core boundary.

We are very thankful to Bruce Buffett for the instructive discussions on the convection in planetary interiors. The authors acknowledge support from the U.S. National Science Foundation (AST 1412646) and from the U.S. Department of Energy (DE-SC0010517 and DE-SC0016248). Calculations were performed on Pleiades from NAS as well as on Comet and Stampede from the XSEDE program.