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

TL;DR

This study uses ab initio simulations to analyze how hydrogen, helium, and silicon dioxide mixtures behave under giant planet interior conditions, informing models of core erosion and element mixing.

Contribution

It provides new computational data on the structural, diffusion, and viscosity properties of these mixtures at relevant pressures and temperatures.

Findings

Identified changes in chemical behavior with pressure and temperature.

Provided diffusion coefficients for species in the mixture.

Estimated time scales for core erosion and mixing processes.

Abstract

Recent observations of Jupiter and Saturn provided by spacecraft missions, such as Juno and Cassini, compel us to revise and improve our models of giant planet interiors. Even though hydrogen and helium are by far the dominant species in these planets, heavy elements can play a significant role in the structure and evolution of the planet. For instance, giant-planet cores may be eroded by their surrounding fluid, which would result in a significantly increased concentration of heavy elements in the hydrogen-helium envelope. Furthermore, the heavy elements could inhibit convection by creating a stabilizing gradient of composition. In order to explore the effects of core erosion, we performed ab initio simulations to study structural, diffusion and viscosity properties of dense multi-component mixtures of hydrogen, helium, and silicon dioxide at relevant pressure-temperature conditions. We computed radial distribution functions to identify changes in the chemical behavior of the mixture and to reveal dissociation trends with pressure and temperature. The computed diffusion coefficients of the different species as well as the viscosity provide constraints for the time scale of the dynamics of the core erosion and the mixing of its constituents into the envelope, which will help improve planetary models.

Figures

Equations

$$D_{\alpha }=\frac{1}{3N_{\alpha }}\sum _i^{N_{\alpha }}{\int }_0^{+\infty }⟨{\text{v}}_i(\tau )\cdot {\text{v}}_i(0)⟩\text{d}\tau ,$$

$$\eta =\frac{V}{3k_{\text{B}}T}\sum _{\{ij\}}{\int }_0^{+\infty }⟨{\sigma }_{ij}(\tau ){\sigma }_{ij}(0)⟩\text{d}\tau ,$$

$${\gamma }_{\text{v}}(t)=a_0{e}^{-t/{\tau }_0}+\sum _{l=1}^{l_{\text{max}}}{e}^{-t/{\tau }_l}\left[{\mu }_l\cos ({\omega }_{l}t)+{\nu }_l\sin ({\omega }_{l}t)\right]$$

$$D_{\alpha }=a_0{\tau }_0+\sum _l{\tau }_l\frac{{\mu }_l+{\nu }_l{\omega }_l{\tau }_l}{1+{\omega }_l^2{\tau }_l^2}$$

$${\gamma }_{\sigma }(t)=A{e}^{-t/{\tau }_1}+B{e}^{-t/{\tau }_2}\sin (\omega t),$$

$$\frac{k_{\text{B}}T\eta }{V}=A{\tau }_1+\frac{B{\tau }_2}{1+{\omega }^2{\tau }_2^2}.$$

$$g_{\alpha \beta }(r)=\frac{V}{4\pi r^2{N}_{\alpha }{N}_{\beta }}⟨\sum _{i_{\alpha },j_{\beta }}\delta (r-|{\text{r}}_{i_{\alpha }}-{\text{r}}_{j_{\beta }}|)⟩,$$

$$\delta \sim 0.34R{e}^{-1/4}L,$$

Full text

††preprint: APS/123-QED

I Introduction

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 require the accretion of heavy rocky or icy core materials before the planets can further accrete gas Pollack et al. (1996). 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 planetFortney et al. (2013). 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 – H2O, SiO2, 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 Wilson and Militzer (2012a, b); Wahl et al. (2013); González-Cataldo et al. (2014). 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 into 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 Wong et al. (2004), 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 Nettelmann et al. (2012); Helled and Guillot (2013); Militzer and Hubbard (2013); Hubbard and Militzer (2016); Militzer et al. (2016). 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 Militzer (2013) 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 Soubiran and Militzer (2016).

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 Stevenson (1982). 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 Leconte and Chabrier (2012, 2013). 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 French et al. (2012), 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 SiO2. 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.

II Methods

II.1 Simulation methods

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 Kresse and Furthmüller (1996). 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 Militzer (2013) and added from 2 to 4 SiO2 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 Nosé (1984, 1991).

At each time step a DFT calculation was performed to determine the electronic density. We used the Kohn-Sham scheme Kohn and Sham (1965) at finite temperature Mermin (1965) and populated the eigenstates using a Fermi-Dirac distribution. Following previous work on H-He mixtures Militzer (2013), we used the Perdew, Burke and Ernzerhof (PBE) exchange correlation functional Perdew et al. (1996) which uses the generalized gradient approximation (GGA). To improve the efficiency of the calculation, we used projector augmented wave (PAW) pseudo-potentials Blöchl (1994) 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 Baldereschi (1973) to sample the Brillouin zone, which was found to be sufficient in the case of H-He mixtures Militzer (2013). We investigated the properties of the mixtures from 5000 K to 15000 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.

II.2 Calculation of the ionic transport properties

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 $\alpha$ by a Green-Kubo formula Frenkel and Smit (2002); Danel et al. (2012),

[TABLE]

where $N_{\alpha}$ is the number of particles of type $\alpha$ and $\textbf{v}_{i}(\tau)$ is the velocity vector of the $i$th particle at time $\tau$. The brackets represent an ensemble average over multiple origins along the full length of the simulation Allen and Tildesley (1987).

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

[TABLE]

where $V$ and $T$ are the volume and the temperature respectively, $k_{\textrm{B}}$ the Boltzmann constant and $\sigma_{ij}(\tau)$ 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 Meyer et al. (2014) and fitted the autocorrelation function with a multi-time scale function. For the velocity autocorrelation function $\gamma_{\textbf{v}}(t)=\frac{1}{3N_{\alpha}}\sum_{i}^{N_{\alpha}}\langle\textbf{v}_{i}(t)\cdot\textbf{v}_{i}(0)\rangle$, we used the following expression:

[TABLE]

The parameters $a_{0}$, $\tau_{0}$, $\{\mu_{l},\nu_{l},\tau_{l},\omega_{l}\}_{l}$ were determined with a least-square fit. We used parameters up to $l_{\textrm{max}}=1$ for H and He. For Si and O we needed to go up to $l_{\textrm{max}}=3$, but, for some cases, the multi-time scale fit did not converged 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:

[TABLE]

For the autocorrelation function of the stress tensor $\gamma_{\sigma}(t)=\frac{1}{3}\sum_{\{ij\}}\langle\sigma_{ij}(t)\sigma_{ij}(0)\rangle$, we used a simpler functional form:

[TABLE]

where $A$, $B$, $\tau_{1}$, $\tau_{2}$ and $\omega$ are fit parameters. Inserting Eq. 5 into Eq. 2, we obtained the following expression for the viscosity,

[TABLE]

As an example, we show in Fig. 1 the fit of the stress-tensor autocorrelation function of a 220 H, 18 He, 3 SiO2 mixture at 15000 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 allow us a much more closer fit of the autocorrelation function.

III Results and discussion

III.1 Structure of the mixture

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 Collins et al. (2001); Morales et al. (2010). It has been shown that adding helium stabilized the molecular phase near the molecular-to-metallic transition Vorberger et al. (2007). In order to investigate the influence of adding SiO2 to the system, we computed the pair-correlation function between particles of type $\alpha$ and $\beta$:

[TABLE]

where $r$ is a variable that typically varies from 0 to half the size of the simulation box, $\textbf{r}_{i_{\alpha}}$ is the position vector of the ${i_{\alpha}}$th particle of type $\alpha$ 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 Soubiran and Militzer (2015a, b).

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 H2 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 Caillabet et al. (2011) and in H-He mixtures Vorberger et al. (2007). It is very interesting to note that the peak height slightly decreases as more SiO2 entities are inserted. Unlike helium, SiO2 does not stabilize the H2 molecules. At 10000 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 as 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 Soubiran and Militzer (2015b). 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 H2 molecules, as was previously predicted for water-hydrogen mixtures Soubiran and Militzer (2015a, b). The existence of H-O bonds also explains why the peak in the H-H pair-correlation function decreases with the number of inserted SiO2. As the number of free oxygen atoms increases, the number of H-O bonds increases breaking H-H bondsbbecause 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 H2 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.

III.2 Diffusion properties

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.2. We only report calculations for the highest concentration in SiO2 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 decreases as the pressure increases because the particles interact more strongly at higher density. It is also interesting to note that at 10000 and 15000 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. French et al. (2012), 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 of waterSoubiran and Militzer (2015b). It is safe to assume that the dissociation of H2 is also at play in the present case.

The calculations 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 15000 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.

III.3 Viscosity

It is usually assumed that the deep interior of giant planets is fully convecting Militzer et al. (2016) but the possible erosion of the core could inhibit the convection Leconte and Chabrier (2012, 2013) 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.2.

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*-4* to 10*-3* 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\sim 10^{12}$ for the convection in Jupiter Guillot et al. (2004). The values that we obtained for the viscosity are in agreement with H-He mixture calculations in similar conditions French et al. (2012). 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 SiO2. In the dilute limit, the properties of the multi-component mixture are mostly the same as that of the H-He mixture. The heavy elements barely influence the bulk viscosity of the system.

III.4 Erosion time-scale

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 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 King et al. (2013) provide scaling laws for the viscous boundary layer thickness for non-rotating systems:

[TABLE]

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=10^{12}$ for the deep interior of Jupiter Guillot et al. (2004)) and higher Prandtl numbers (above 1 for the simulations but closer to $10^{-2}$ for Jupiter French et al. (2012)). But comparison with experiments in other range of parameters give confidence for an extrapolation of this scaling law towards more extreme regimes such as the ones in giant planets. If we take $L\sim 10^{7}$ m for the eddy size of the convective cell Guillot et al. (2004), we find a boundary layer thickness of the order of $\delta\sim 3$ km. Using the estimate from Fig. 7 we can estimate the diffusion coefficient to be close to 0.1 mm2/s, which results in a typical diffusion time scale of $3\times 10^{6}$ years. Such a time scale would suggest a nearly complete erosion of the core for a planete 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.

IV Conclusion

Our ab initio simulations of H-He-SiO2 mixtures under conditions relevant for giant planet interiors showed a smooth dissociation of hydrogen as the pressure and the temperature increase. Unlike helium, SiO2 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 Stevenson and Salpeter (1977); Morales et al. (2009); Lorenzen et al. (2011); Soubiran et al. (2013). 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 mm2/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.

V Acknowledgments

We are very thankful to Bruce Buffett for the instructive discussions on the convection in planetary interiors. The authors achnowlegde 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.