Partially Diffusive Helium-Silica Compound under High Pressure

Understanding the interiors of the giant planets in our Solar System is a key objective of planetary science and a multidisciplinary challenge combining condensed matter physics, astrophysics, and geophysics. This challenging task demands abundant accurate measurements, accompanied by theoretical models that are used to infer the planetary conditions and fit the available observational constraints. Traditional models^[1–3] describe the outer planets (Jupiter, Saturn, Uranus, and Neptune) with three distinct layers: a gas-rich envelope mainly composed by hydrogen and helium, a denser convective and electrically conductive envelope that yields huge planetary magnetic fields (this layer is mainly composed of metallic fluid hydrogen for Jupiter and Saturn, and of superionic water/ammonia/methane mixture for Uranus and Neptune), together with a compact heavy element central core, with a density discontinuity at the core envelope boundary.

Thanks to updated gravity data from the Juno and Cassini missions and advance in planetary models, more complex models have been developed, in which the planetary interior is inhomogeneous and the core is diluted.^[4–8] For both Jupiter and Saturn, the heavy elements could be gradually distributed or homogeneously mixed with lighter elements (hydrogen and helium) and extended to about half of the planet's total radius.^[6,8] Such extended cores are difficult to explain within standard giant planet formation models,^[9] and therefore it was suggested that Jupiter's fuzzy core could be a result of an energetic giant impact between a large planetary embryo with ten times Earth mass and the proto-Jupiter.^[10]

However, while it is still challenging to construct accurate and unique planetary structure models for the outer giant planets, we can investigate the relevant physical and chemical processes of simple elements at high pressures and high temperatures to guide planetary interiors.^{[5–7,11–13]} Although helium is the most inert element at ambient pressure, it can form stable compounds with other matter under high pressure. For instance, Dong et al.^[14] have experimentally synthesized a novel Na₂He compound at around 113 GPa and 1500 K. In the deep interior of giant planets, the pressures and temperatures may provide sufficient energy for helium to form new compounds^[15–19] together with other ingredients of planet interiors. Among the outer giant planets in our Solar System, another major component of the heavy elements is silica, which undergoes a series of phase transitions with compression^[20–22] and may have a significant influence on the formation and evolution of terrestrial planets.^[23] For the outer planets, the possible mixing or even chemical reaction of helium and silica could be important.

As aforementioned, helium and silica are major components of giant planets, nevertheless whether they can form new stable compounds under pressure and what states they can exist at giant planetary conditions are still open questions. To address these questions, we have systematically studied the helium silica system within the pressure range of planetary core conditions and found four new stable HeSiO₂ phases by crystal structure prediction methods and first-principles calculations. Further molecular dynamics simulations show that the HeSiO₂ compounds may survive inside Saturn with superionic-like helium diffusive phase and in Uranus and Neptune with solid phase, which provides more explication of the formation of diluted core in gaseous planets.

Results. We explore the crystal structures in the helium-silica system in the pressure range of 500–4000 GPa using variable composition structure prediction method. A structure is regarded to be thermodynamically stable if its enthalpy of formation is negative relative to the mixture of the most stable phases of solid helium and silica at the corresponding pressures. As the convex hulls shown in Fig. 1(a), we find that helium and silica will form a new stable HeSiO₂ compound just above 600 GPa and such HeSiO₂ compound inclines to gain its energetic stability with further compression. We have checked the stability of the newly predicted HeSiO₂ compound against different exchange correlation functionals including different van der Waals (vdW) corrections, a full-potential linearized augmented plane wave (FP-LAPW) method, and zero point energy corrections, as shown in Fig. S1 of the Supplemental Material. The results demonstrate that the HeSiO₂ compound can survive under deep planetary pressure.

Zoom In Zoom Out Reset image size

As shown in Fig. 1(b), we find several energetically competitive candidates for the HeSiO₂ compound. A Pnma phase (denoted as Pnma-I), will gain its stability at around 605 GPa, and then, it can transform into a Pmn2₁ phase at ∼ 1100 GPa. We find that the enthalpy of the Pnma-I and Pmn2₁ phases are pretty close, because they share the same point group and have similar crystal structures, as shown in Figs. 1(c) and 1(d). With increasing pressure, the Pmn2₁ phase transforms into another Pnma phase (denoted as Pnma-II), at around 2100 GPa. Finally, a Pnma-III phase, would gain its stability above 2300 GPa. Phonon calculations indicate that the above structures are all dynamically stable and temperature affects Gibbs free energy significantly for Pnma-II phase, as shown in Figs. S2 and S3 of the Supplemental Material. We also calculate the electronic band gaps of these HeSiO₂ compounds and they are completely insulating with wide band gaps between that of He and SiO₂ crystal, as shown in Figs. S4 and S5 of the Supplemental Material. It is clear that the electronic band gaps decrease in He and SiO₂ crystals with increasing pressure. The SiO₂ crystal is predicted to transform into a semiconductor at 4 TPa and could metallize at higher pressure. Noticeably, the HeSiO₂ compounds have an abnormal pressure dependence with an increase in the band gap over the pressure range 500–2000 GPa. Further Bader charge (Table 1 of the Supplemental Material) calculations show that helium atoms obtain some charge (around 0.1 e) from the SiO₂ framework, which indicates that the ionic interaction may play a role in the stability of SiO₂He compounds. Previous work^[22] showed that six-fold silica would transform into a mixed coordination silica with an averaged coordination number of 8 by compression. Interestingly, in our newly found HeSiO₂ compound, all silicon atoms in Pnma-I, Pmn2₁, and Pnma-II phases are seven-fold coordinated with oxygen atoms, while Pnma-III phase is purely eight-fold. These seven-fold and eight-fold configurations are rare in pure silica, which indicates that the inserting of helium atoms can significantly change the packing of silicon and oxygen atoms, as well as the silicon-oxygen bonding behaviors. Meanwhile, we also find that a purely eight-fold HeSiO₄ compound with I422 phase can survive above 400 GPa, as shown in Figs. S6 and S7 of the Supplemental Material. This may shed light on further explorations on chemical coordination in helium compounds.

Since HeSiO₂ compound could exist in high pressure range corresponding to deep interiors of giant planets, we investigated the equation of state of the HeSiO₂ compound, as shown in Fig. 2. In the traditional 3-layer models for of Jupiter and Saturn with distinct layers, their density curves increase sharply at core mantle boundary. Our ab initio calculations show that the density curves of the HeSiO₂ compound have a smooth tendency: increasing by compression, while decreasing by heating. Most importantly, they are located just between that of the core and mantle, which indicates that the newly found HeSiO₂ compound very possibly exists near the core mantle boundary of giant gaseous planets, such as Jupiter and Saturn. Especially when consider diluted core models, helium/hydrogen is expected to erode the core and their density curves would change to a smooth one, which should agree with the density of states of the HeSiO₂ compound even more closely, in the deep interiors of Jupiter and Saturn.

Zoom In Zoom Out Reset image size

In order to gain a better understanding of such HeSiO₂ compound affecting the interior model of giant planets, one should account for both the equation of states calculations and the pressure-temperature phase diagrams, because the internal structure models must be consistent with the phase diagram of the assumed materials and their dynamical behavior. We perform extensive ab initio molecular dynamics (AIMD) simulations at deep planetary conditions to study the dynamical properties of our predicted HeSiO₂ compound. Diffusion coefficients (D) were calculated for the silicon, oxygen, and helium atoms from the averaged mean-square displacements (MSD) to monitor different types of atomic motions (D = ∂ MSD/∂ t). We demonstrated the atomic trajectories from simulations at the initial pressure of 600 GPa for the Pnma-I phase HeSiO₂ compound as an example, as shown in Fig. S8 of the Supplemental Material. At 5000 K, all atoms oscillate around their equilibrium positions, resulting in three horizontal MSD curves with slightly oscillating (D = 0). When heating up to 9000 K, we find that the helium atoms became diffusive freely within the static silica frameworks during the simulations (D_He > D_O = D_Si = 0). These are clearly two different states: the solid phase at 5000 K and the superionic-like partial diffusive phase at 9000 K.

Inspired by this helium diffusive state in the HeSiO₂ compound, we extended the pressure and temperature range of our AIMD simulations to explore the states of the HeSiO₂ compound under the deep interior condition of Saturn and Jupiter, as shown in Fig. 3. Each colored symbol corresponds to an independent simulation to avoid correlation effects, and the pressures and temperatures were obtained from simulations by statistically averaging. Due to the different stable pressure range of HeSiO₂ phases, as well as the core mantle boundary conditions of Saturn (1 TPa and 10000 K) and Jupiter (4 TPa and 20000 K) varying a lot, we separately simulate HeSiO₂ compounds in pressure range of 500–1200 GPa with Pnma-I phase [as shown in Fig. 3(a)] and in pressure range of 2000–5000 GPa with Pnma-III phase [as shown in Fig. 3(b)] up to the melting temperature. It is clear that the superionic-like helium diffusive state is widespread at high pressure and high temperature and exists between the solid phase and fluid phase in the HeSiO₂ phase diagram. With increasing pressure, this helium diffusive state appears at higher temperatures.

Zoom In Zoom Out Reset image size

For comparison, the isentropes of the giant outer planets (Jupiter, Saturn, Uranus, and Neptune) are also plotted in Fig. 3 to represent the pressure temperature profiles at deep interior conditions. We find that the pressure temperature profile of Saturn (magenta) perfectly passes through the superionic-like helium diffusive region of the HeSiO₂ compound, which suggests that such helium diffusive HeSiO₂ compound can exist near the core mantle boundary of Saturn. As for Jupiter, both adiabatic (red) and non-adiabatic (white) pressure temperature profiles are slightly higher than the helium diffusive region. In other words, the HeSiO₂ compound is totally melted in the core region, if we account for a diluted core model of Jupiter. Previous works^[26,27] on typical mantle silicates (MgO, SiO₂, and MgSiO₃) showed that upon melting, the behavior of SiO₂ changes from semi-conducting to semi-metallic, indicating that a magnetic dynamo process would develop in the magma oceans of Super-Earths. We calculate the electronic density of states (DOS) of HeSiO₂ compounds at different temperatures and find that at higher temperature HeSiO₂ compounds may also transform into metallic in fluid state, which may affect the conductivity and therefore the magnetic field generation in Jupiter. The pressure temperature profiles of Uranus and Neptune are much lower than that of Jupiter and Saturn, and in this case the HeSiO₂ compounds may exist in solid form in the deep interior of Uranus and Neptune when helium deposits and erode the core.

We also calculate the elastic and wave velocity properties of the HeSiO₂ compound to explore how it affects the evolution of the interior of Uranus and Neptune. The calculated bulk moduli (K_S), shear moduli (G), compressional velocities (V_P), and shear wave velocities (V_S) of the Pnma phase HeSiO₂ compound, as well as pure helium and silica for comparison, at pressures 500–1000 GPa are shown in Fig. 4. He and HeSiO₂ compound both have almost linear temperature and pressure dependences of bulk moduli and shear moduli. For instance, in the HeSiO₂ compound, the first temperature and pressure derivate of elastic properties almost keep unchanged (∂ K_S/∂ T = –0.013 GPa/K, ∂ G/∂ T = –0.016 GPa/K, ∂ K_S/∂ P = 2.7, and ∂ G/∂ P = 1.1). In contrast, noticeable nonlinear dependences on temperature and pressure are observed for the wave velocities. Comparing with pure silica, the penetration of helium into silica can reduce the compressibility as the previous report in helium-silica mixture.^[29]

Zoom In Zoom Out Reset image size

The chemical composition can also affect the elastic moduli and wave velocities of these planetary matters. Although pure helium has a pretty small bulk moduli and shear moduli, helium has larger compressional velocities and shear wave velocities than silica and HeSiO₂ due to its small density $({V}_{{\rm{P}}}=\sqrt{({K}_{{\rm{S}}}+\displaystyle \frac{4}{3}G)/\rho }$ and ${V}_{{\rm{S}}}=\sqrt{G/\rho }$, where ρ is density). In silica, the Cotunnite-type phase surviving in high temperature region has a small shear moduli and shear wave velocities with respect to other phases, which indicates that an ultralow-velocity zones may exist in the mantle of super-earth. Compared to silica, HeSiO₂ has a slightly smaller bulk moduli but larger shear moduli. Due to the decrease of density by inserting helium, the velocities of HeSiO₂ compound have a larger compressional velocities and shear wave velocities than that of silica. Notably, the temperature does not significantly affect wave velocities. For instance, wave velocities vary 0.16 and 0.23 km/s for compressional velocities and shear wave velocities from 4000 K to 7000 K at 500 GPa and these differences are evidently diminished at high pressure, which decrease to 0.08 and 0.14 km/s at 1000 GPa. Thus, the inserting of helium may evidently increase the wave velocities of silica and our newly predicted HeSiO₂ may affect the model of the deep interior of Uranus and Neptune.

Discussion. The composition of giant planets cannot be measured directly. Instead, their bulk compositions and internal structures can only be inferred indirectly from interior models that fit the available measured physical parameters. Here we used the calculated equations of state (EoS) of our predicted HeSiO₂ compound to derive the density (as shown in Fig. 2.) and the associated pressure temperature phase diagram (as shown in Fig. 3), and then we attempt to sketch internal structure models for Jupiter, Saturn, Uranus, and Neptune.

Our static calculations show that helium and silica can react with each other in the deep interiors of giant planets, and the HeSiO₂ compounds may survive near the core mantle boundary of giant planets. Especially considering that helium may be immiscible with hydrogen at 1–2 MBar, namely "helium rain", it is expected to settle down into deep interior of Jupiter and Saturn. Thus, helium in the deep interior could affect the erosion of the compact heavy-element core, which may result in a gradually expanding core region and forming a diluted one. Here, we only fitted density curves and pressure temperature phase diagram of the HeSiO₂ compound to the representative planetary models, more physical data (such as their masses, radii, gravitational and magnetic fields, 1-bar temperatures, atmospheric composition, and internal rotations) of planets are required to construct a more comprehensive model. Moreover, it should be noted that other elements may influence the stability of silica at high pressure as well, for instance, the hydrogen can also react with silica at high pressure, as we discussed elsewhere.^[12] However, a complete investigation in the H–He–Si–O system or more complicated systems with other elements is beyond the scope of this work.

In summary, using crystal structure prediction and ab initio calculations, we have predicted four phases of HeSiO₂ compound (Pnma-I, Pmn2₁, Pnma-II, and Pnma-III), which can gain their stability in the pressure range of 600–4000 GPa, corresponding to deep interior conditions of giant planets, such as Jupiter, Saturn, Uranus, and Neptune. Due to the spatial occupation of helium atoms, the silicon atoms in these compounds are all bonded to seven/eight oxygen atoms, which are very rare cases in pure silica. Thus, helium and other inert gas atoms can be used as a space filler to design compounds with usual chemical bonding and coordination under high pressure. Equation of states calculations suggests that the density curves of our newly predicated HeSiO₂ compound are close to current models of Jupiter and Saturn, especially when accounting for diluted core models. Furthermore, extensive ab initio molecular dynamics simulations illustrated that the HeSiO₂ compound can survive in the helium diffusive state under the core-mantle boundary conditions of Saturn and metallic fluid state under the core-mantle boundary conditions of Jupiter, which may shed light on the formation or evolution of diluted core of such gaseous giants. Noticeably, under the pressure-temperature conditions in the deep interiors of Uranus and Neptune, the HeSiO₂ compound is found to be in solid form due to lower temperatures. We also carried out elastic and wave velocity properties calculations for the HeSiO₂ compound in the pressure range 500–1000 GPa and found that the inserting of helium increases the compressional and shear wave velocities in the core region of Uranus and Neptune. Our findings may be used to guide future shock wave compression experiments and giant planet interior models, as well as to significantly improve our understanding on giant planets in the solar system and exoplanets.

Methods. We used Magus (machine learning and graph theory assisted universal structure searcher) code to search for the crystal structures, in which we employed the Bayesian optimization^[30] and graph theory^[31] to improve the search efficiency and diversity. We performed extensive crystal structure searches on He_x (SiO₂)_y (x = 1–4, y = 1–4) at 500, 1000, 2000, and 4000 GPa with maximum atom number up to 40. Meanwhile we also performed ternary variable composition crystal structure searches on He_x Si_y O_z . Some compositions are further double checked with extensive fixed composition searches. Each search runs over 25 generations and each generation has a population size of 60 structures. There are 40%–60% of the parents for the evolution of next generation from the lowest enthalpy structures of the last generation, and the left seeds are randomly produced. We also cross checked the searching results with AIRSS^[32,33] combined with CASTEP.^[34]

DFT calculations were performed with the Vienna ab initio simulation package (VASP),^[35] accompanied with the projector augmented-wave (PAW) method.^[36] We chose 3s²3p², 2s²2p⁴, and 1s² as valence electrons for Si, O, and He, and used the generalized gradient approximation (GGA) in the Perdew–Burke–Ernzerhof exchange correlation functional.^[37] All predicted structures were further optimized by the hard version pseudopotentials and employed a plane wave cutoff of 1050 eV and dense Monkhorst–Pack k-point sampling grids with resolutions of 2π × 0.025 Å⁻¹, leading to ionic and cell optimizations with energy and force convergences better than 10⁻⁶ eV and 0.002 eV/Å, respectively. The CASTEP code^[34] was used for the ab initio random structure searching approach searches. Similar convergence parameters as used in the Vienna ab initio simulation package were employed.

Elastic properties at high pressure and temperature are calculated by cij package^[38] based on phonon spectrum from the PHONONPY package.^[39] The static elastic constants are calculated by the stress strain method. We should notice that the calculated results of elastic properties will be influenced by temperature, and calculations based on the quasi-harmonic approximation will be invalid in the superionic/fluid state where some atoms can diffuse freely.

We adopted the canonical NVT ensemble using a Nose–Hoover thermostat^[40] to perform ab initio molecular dynamics (AIMD) simulations in a supercell with 192 atoms for both Pnma-I phase and Pnma-III phase, with Γ-centered k-points sampling, normal version pseudopotentials and a cutoff energy of 720 eV were adopted to ensure energy convergence of better than 10⁻⁵ eV. Each simulation lasts for 12 ps with a time step of 1 fs, and we allowed the first 2 ps for thermalization. We tested the influence from box size and k-mesh choice, as shown in Table 2 of the Supplemental Material. DOS at definite temperatures is estimated by averaging over 10 configurations from molecular dynamics trajectory.

Acknowledgements