# Sequestration of Noble Gases in Giant Planet Interiors

Phys. Rev. Lett. 104, 121101

Jupiter is the most extensively probed and best understood of the giant planets, but many questions regarding its detailed composition, formation, and interior structure remain unanswered. One issue of major importance to structural models is the question of whether hydrogen and helium mix homogeneously throughout the planet or whether a layer of hydrogen-helium immiscibility exists deep within the interior [1–3]. In the immiscibility layer, helium would form dense droplets which would rain down into the deeper interior, resulting in a gradual and ongoing transfer of helium from regions above the immiscibility layer to regions below (Fig. 1). Such a layer almost certainly exists in Saturn, as evident from the observed depletion of helium from its upper atmosphere (compared to protosolar values) and the apparent excess luminosity of the planet [4]. For the hotter interior of Jupiter the case is less clear since there is no measurable excess luminosity and the observed helium depletion from the upper atmosphere is quite small (0.234 by mass compared to 0.274 in the protosolar nebula [5–7]). Theoretical attempts to determine the pressure or temperature range in which H and He are immiscible using successively more sophisticated levels of theory [3,8–14] have produced quite different results; however, recent work [13,14] provides a hydrogen-helium immiscibility line which is very close to the Jupiter isentrope in the 100–300 GPa region.

The strongest evidence for H-He immiscibility may in fact come from the depletion of neon. Jupiter’s upper atmosphere was found by the Galileo entry probe to be extremely deficient in neon: although neon makes up $1/600$ by mass in the solar system it comprises only $1/6000$ by mass in Jupiter’s upper atmosphere [5]. Prior to these measurements, Roulston and Stevenson [15] proposed that hydrogen-helium immiscibility could lead to neon depletion on the assumption that neon would preferentially dissolve in the helium-rich phase in the immiscibility layer. This would lead to Jupiter’s neon content being gradually carried down within the helium droplets and concentrating in the deep interior. There is, however, a lack of direct experimental evidence for whether Ne will indeed preferentially dissolve in the helium phase as proposed. The pressure-temperature conditions corresponding to phase separation can currently only be attained in shock-wave experiments lasting only tens of nanoseconds [16]. Is it also not known why the chemically similar noble gas argon is not seen to be depleted but instead is present at slightly above-solar concentrations ( $≈1.6$ times solar [5]) comparable to most other detectable trace heavy elements in Jupiter, and whether this indicates that the depletion mechanism acting upon neon does not act upon argon or whether it indicates a very high initial argon concentration. In order to resolve these issues, here we present ab initio free energy calculations with a view to determining the solubility behavior of neon and argon in H and He phases at pressures corresponding to the postulated immiscibility region.

The distribution of a trace species between coexisting phases is dependent upon the Gibbs free energy of transfer, $ΔGTr$, being the change in $G$ when an atom of the trace species is moved from one phase to the other at constant $P$ and $T$, in this case

Here we compute $ΔGTr$ for Ne and Ar in pure H and He within the density functional theory molecular dynamics (DFT-MD) framework, within a temperature and pressure range of 100–300 GPa and 3000 to 7000 K. Determination of free energies from MD is a nontrivial problem for which a number of methods have been developed. We use a two-step coupling constant integration (CCI) approach [17] similar to that recently applied by Morales et al. [14], in which the Gibbs free energy of the DFT system is determined by adiabatically transforming the system in two steps: (a) from the DFT system to a system of atoms governed by a classical potential and (b) from the classical system to a noninteracting gas.

The CCI method provides a general scheme for computing the Helmholtz free energy difference between two systems governed by potential energy functions $U1$ and $U2$. Constructing an artificial system with potential $Uλ=(1-λ)U1+λU2$,

where at each integration point, the average is taken over a sample of configurations obtained in the $Uλ$ system. Since the difference in Gibbs free energy between a system at two different pressures can be found by the thermodynamic integration $G(P2)-G(P1)=∫P1P2VdP$, we performed all CCI calculations at pressures close to 200 GPa and then integrated the equation of state for each system outwards to obtain $G$ values at pressures from 100 to 300 GPa.

The first part of the CCI was the integration from the noninteracting system to the classical system. The classical potential used was a pair potential of a modified Yukawa form [14]:

for $r and zero otherwise. We set $L=9.749 a.u.$, then fitted $a$ and $b$ to the $g(r)$ functions of the DFT systems at high pressure [18]. The integration from the noninteracting to classical system used 16 $λ$ values. We checked potentials obtained by a force-matching method and found that the numerical results achieved for the free energy of the final system were not altered by the different potential.

The second part was the integration from the classical forces to the DFT forces. We found that the variation in $⟨UDFT-Uclassical⟩$ was sufficiently close to linear in $λ$ to allow a fit with only three $λ$ points to be used—checks against calculations with five $λ$ points resulted in discrepancies smaller than 0.1 eV. All DFT simulations were performed using the VASP code [19]. We used 128 H atoms or 64 He atoms and a single $Ne/Ar$ atom per cell. We used pseudopotentials of the projector-augmented wavefunction type [20], the exchange-correlation functional of Perdew, Burke, and Ernzerhof [21], a cutoff of 1000 eV and eight $k$ points in the Monkhorst-Pack grid. All MD simulations used a time step of 0.4 fs and were run for 5000 time steps, with the first 1000 steps discarded for equilibration.

The total Gibbs free energy computed for each system is a sum of five terms:

where $P0$, $V0$ are the pressure and volume at which the $F$ values were computed, and $P1$ is the pressure of interest. The CCI procedure was undertaken at pressures within 1% of 200 GPa and $VdP$ integration was used to correct the values back to 200 GPa exactly. Pressure-volume curves were obtained from a series of five MD simulations on each system at pressures spaced from 100 to 300 GPa, and by fitting the resulting data points with a piecewise power law fit.

In order to validate this method, we also performed simulations via an alternative free energy calculation method based on the particle insertion formalism of Widom [22]. Using only the $Γ$ point for Brillouin Zone sampling, we computed the free energies associated with the insertion of Ne and Ar into He and H cells at volumes corresponding to a Wigner-Seitz radius for the electrons of 2.4 bohr radii, then integrated along isotherms to obtain values of $ΔGTr$ which were then compared with CCI-computed values. The results were found to agree within the relevant error bars. Since the CCI method is more computationally efficient we applied it for the computation of the final, eight $k$-point results. We also estimated the quantum correction to the classical free energy resulting from the fluctuations around the classical trajectories of the nuclei. Using the first term of a Wigner-Kirkwood expansion in $ℏ$ [23], we estimate the free energy correction at 5000 K and 200 GPa to be of the order of 0.01 eV per atom or less, and can consequently be neglected.

The computed values of $ΔGTr$ for neon and argon are given in Table I. For neon at 200 GPa we find $ΔGTr$ values of approximately $-2.4 eV$ with temperature variation from 3000 to 7000 K producing only a small variation in $ΔGTr$. The negative sign here indicates a preference for solubility in helium. Thermodynamic integration of the pressure-volume curves as shown in Fig. 2(a) from 200 GPa shows that the helium preference becomes more pronounced with increasing pressure. For argon, we see contrasting behavior: at 200 GPa and 5000 K we have a $ΔGTr$ of $+5.1±0.45 eV$ with the positive sign now indicating a preference for the hydrogen phase. The magnitude of the preference for hydrogen solubility increases somewhat with both temperature and pressure.

Following the work of Roulston and Stevenson [15], we expect the rate at which neon is removed from the upper envelope to be related to the loss rate of helium by

This implies that the observed depletion of neon will be approximately given by

where $XQ,0$ and $XQ,1$ refer to the original (protosolar) and present-day molar atomic concentrations of species $Q$ in the upper Jovian atmosphere, respectively. Based on the measurements of Von Zahn et al. [6] for the current helium concentration and the estimate of Lodders [7] for the protosolar concentration, we find a helium depletion $XHe0-XHe1$ value of approximately 1.2%. Combining this with values of $ΔGTr$ of $-2.35±0.45 eV$ for neon partitioning, we obtain the relationship between $T$ and neon depletion shown in Fig. 3. The observed value of approximately 0.1 for the neon depletion ratio corresponds to $T$ values of between approximately 4000 and 6000 K. This is consistent with the expected temperature of the immsicibility region. The computed value of $ΔGTr$ is thus consistent with the assumption that the observed depletions of both helium and neon are due entirely to helium rain within the hydrogen-helium immiscibility layer.

For argon, the positive value of $ΔGTr$ implies that Ar will be almost entirely excluded from the He phase. Since the helium phase remains only a very small portion of the planet this will lead only to a miniscule enhancement in the argon content of the upper atmosphere. This implies that the measured concentration, approximately 1.6 times the solar value [5], should be close to the true argon concentration of the planet as a whole.

The difference in solubility behavior between neon and argon invites further examination. The difference in the free energies of insertion is governed primarily by the volume change $ΔVins$ associated with the insertion at constant pressure of the noble gas atom into the pure-solvent cell. As shown in Fig. 2(a) the effective volume of neon is larger by $0.86 Å3$ at 200 GPa and 5000 K in hydrogen than in helium, while argon shows the opposite trend, being larger by $0.73 Å3$ in helium than hydrogen. In Fig. 4 we plot the pair correlation function $g(r)$ for the solvent atoms surrounding each species of noble gas atom. The $g(r)$ curves have been shifted by $Rsol$, the effective radius of the solute atom derived from the point at which the $g(r)$ for H-H or He-He crosses 0.5. There is a clear difference in the exclusion behavior of neon and argon, with the small-distance tail of the H-Ar curve allowing a closer effective approach than in helium, in contrast to of the H-Ne and He-Ne curves where helium approaches more closely.

As a possible interpretation, we note that in the $P$, $T$ range of interest the H atom is essentially ionized [24,25] while the He atom retains its electrons. A helium atom thus is repelled from the $Ne/Ar$ atom by the electron-electron interactions dominated by Pauli exclusion, whereas hydrogen atoms may more easily penetrate the outer shells and are repelled primarily by core-core repulsion. The Ar atom’s additional electron shell thus gives it a larger effective volume to exclude the helium atom, but much less so the hydrogen. If this model is correct then we would expect the noble gases krypton and xenon to likewise exhibit a preference for hydrogen solubility, a behavior consistent with their observed nondepletion the upper atmosphere [5].

In this work we have considered only pure helium and pure hydrogen phases. In practice the helium phase will have very little hydrogen, but the hydrogen-dominant phase still contains some helium [14]; however, we do not expect this to qualitatively change the results. Another limitation not considered in this study is whether the partitioning coefficient changes as the neon concentration in helium increases; it should be noted that the required molar concentration of neon in the pure-helium phase will be quite large. We also cannot exclude the possibility that neon forms its own phase; however, we consider this unlikely due to the small initial Ne concentration.

These results strongly support the existence of hydrogen-helium phase separation in Jupiter as an explanation for the observed Ne depletion. We have also shown that argon will be preferentially excluded from a helium droplets explaining the observed lack of depletion of this element. Further work to more accurately determine the location of the hydrogen-helium immiscibility line at low pressure and low temperature would allow us to make a quantitative estimate of the neon concentration in Saturn to be tested by future missions. Furthermore, neon may be added as a tracer in laboratory experiments to detect the phase separation of H-He mixtures because neon scatters x rays more strongly.

## References

1. E. E. Salpeter, Astrophys. J. 181, L83 (1973).
2. D. J. Stevenson, Phys. Rev. B 12, 3999 (1975).
3. D. J. Stevenson and E. E. Salpeter, Astrophys. J. Suppl. Ser. 35, 221 (1977).
4. D. J. Stevenson, Science 208, 746 (1980).
5. H. B. Niemann et al., Science 272, 846 (1996).
6. U. von Zahn, D. M. Hunten, and G. Lehmacher, J. Geophys. Res. 103, 22815 (1998).
7. K. Lodders, Astrophys. J. 591, 1220 (2003).
8. D. J. Stevenson and E. E. Salpeter, Astrophys. J. Suppl. Ser. 35, 239 (1977).
9. D. M. Straus, N. W. Ashcroft, and H. Beck, Phys. Rev. B 15, 1914 (1977).
10. W. B. Hubbard and H. E. DeWitt, Astrophys. J. 290, 388 (1985).
11. J. E. Klepeis, K. J. Schafer, T. W. Barbee, and M. Ross, Science 254, 986 (1991).
12. O. Pfaffenzeller, D. Hohl, and P. Ballone, Phys. Rev. Lett. 74, 2599 (1995).
13. W. Lorenzen, B. Holst, and R. Redmer, Phys. Rev. Lett. 102, 115701 (2009).
14. M. A. Morales et al., Proc. Natl. Acad. Sci. U.S.A. 106, 1324 (2009).
15. M. S. Roulston and D. J. Stevenson, Proceedings of the EOS Abstr. Fall Meeting, December 1995, San Francisco (American Geophysical Union, Washington, DC, 1995), Vol. 76, p. 343.
16. M. D. Knudson, D. L. Hanson, J. E. Bailey, C. A. Hall, and J. R. Asay, Phys. Rev. Lett. 90, 035505 (2003).
17. D. Alfé, M. J. Gilland, and G. D. Price, Nature (London) 405, 172 (2000).
18. The parameters ($a$, $b$) obtained for the pair potential between each pair of species were as follows, using atomic units: H-H:(1.0,2.5), He-He:(5.5,1.6), H-Ne:(8.0,1.86), He-Ne:(8.76,1.34904), H-Ar:(13.0,2.119), He-Ar:(16.5,1.36125).

19. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
20. P. E. Blochl, Phys. Rev. B 50, 17953 (1994).
21. J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
22. B. Widom, J. Chem. Phys. 39, 2808 (1963).
23. E. L. Pollock (personal communication).
24. B. Militzer, M. Magro, and D. Ceperley, Contrib. Plasma Physics 39, 151 (1999).
25. S. T. Weir, A. C. Mitchell, and W. J. Nellis, Phys. Rev. Lett. 76, 1860 (1996).

## Related Articles

Atomic and Molecular Physics

### Viewpoint: Bose Polarons that Strongly Interact

Researchers have used impurities within a Bose-Einstein condensate to simulate polarons—electron-phonon combinations in solid-state systems. Read More »

Particles and Fields

### Synopsis: Deuterons Spin Together for a Thousand Seconds

Researchers set a new record for the in-plane spin-alignment lifetime of deuterons circulating in a magnetic storage ring. Read More »

Quantum Physics

### Synopsis: Testing Quantum Physics with Neutrinos

An experiment similar to the Bell inequality test confirms that neutrino oscillation is a quantum physics effect that is incompatible with alternative classical models. Read More »