Sequestration of Noble Gases in Giant Planet Interiors

  • Hugh F. Wilson and Burkhard Militzer, Department of Earth and Planetary Science, University of California Berkeley, Berkeley, California 94720, USA and Department of Astronomy, University of California Berkeley, Berkeley, California 94720, USA
Phys. Rev. Lett. 104, 121101
(left) Schematic depiction of the interior of a gas giant (e.g., Jupiter or Saturn) with a layer of H-He immiscibility. Helium-rich droplets form within the immiscibility layer and rain downwards, leading to a slow increase in the helium concentration in the deep interior. Neon is absorbed within the droplets and carried out of the upper atmosphere. (right) P/T curves for Jupiter and Saturn combined with the location of the H-He immiscibility region determined in the work of Morales et al. [14] assuming an overall He atomic molar concentration of 0.0847.(left) Schematic depiction of the interior of a gas giant (e.g., Jupiter or Saturn) with a layer of H-He immiscibility. Helium-rich droplets form within the immiscibility layer and rain downwards, leading to a slow increase in the helium concentratio... Show more

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.

At left, the difference in volume ΔVins between the pure hydrogen/helium cells and the cells with a single Ne/Ar atom added isobarically at 5000 K. At right, the computed difference in free energy ΔGins between the pure H/He cells and the cells with a single Ne/Ar atom added for pressures between 100 and 300 GPa at 5000 K.At left, the difference in volume ΔVins between the pure hydrogen/helium cells and the cells with a single Ne/Ar atom added isobarically at 5000 K. At right, the computed difference in free energy ΔGins between the pure H/He cells and the cells with ... Show more

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<L/2 and zero otherwise. We set L=9.749a.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.

Relationship between the temperature of the immiscibility region and the change in helium concentration which would be required to produce the observed 90% depletion of neon, for ΔGTr values of 2.35±0.45eV. The observed value of 1.2% for ΔXHe is marked with a line.Relationship between the temperature of the immiscibility region and the change in helium concentration which would be required to produce the observed 90% depletion of neon, for ΔGTr values of 2.35±0.45eV. The observed value of 1.2% for ΔXHe is m... Show more

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.

Pair correlation functions g(rRsol) for distances between solvent (H, He) and solute (Ne, Ar) atoms at 200 GPa and 5000 K. The curves are shifted by the effective radius Rsol of the solvent atom in each case, determined from the point where the solvent-solvent g(r) crosses 0.5. Rsol is 0.37 Å for hydrogen and 0.58 Å for helium.Pair correlation functions g(rRsol) for distances between solvent (H, He) and solute (Ne, Ar) atoms at 200 GPa and 5000 K. The curves are shifted by the effective radius Rsol of the solvent atom in each case, determined from the point where the solv... Show more

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.4eV 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.45eV 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.45eV 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).

About the Authors

Image of Hugh F. Wilson
Image of Burkhard Militzer

Related Articles

Synopsis: A Crystal of Light and Atoms
Atomic and Molecular Physics

Synopsis: A Crystal of Light and Atoms

A predicted type of atom-light crystal could host phonon-like excitations, allowing for new ways to simulate the physics of solids.   Read More »

Viewpoint: An Arrested Implosion
Condensed Matter Physics

Viewpoint: An Arrested Implosion

The collapse of a trapped ultracold magnetic gas is arrested by quantum fluctuations, creating quantum droplets of superfluid atoms. Read More »

Focus: Complex Crystals Form from Heterogeneous Particles
Materials Science

Focus: Complex Crystals Form from Heterogeneous Particles

A suspension containing particles with wide-ranging diameters can crystallize into multiple ordered structures. Read More »

More Articles