First-Principles Studies of the Metallization and the Equation of State of Solid Helium
The insulator-to-metal transition in solid helium at high pressure is studied with first-principles simulations. Diffusion quantum Monte Carlo (DMC) calculations predict that the band gap closes at a density of g/cm and a pressure of terapascals, which is % higher in density and higher in pressure than predicted by density functional calculations based on the generalized gradient approximation (GGA). The metallization density derived from GW calculations is found to be in very close agreement with DMC predictions. The zero point motion of the nuclei had no effect on the metallization density within the accuracy of the calculation. Finally, fit functions for the equation of state are presented and the magnitude of the electronic correlation effects left out of the GGA approximation are discussed.
At low pressure, helium is an inert gas that exhibits a very large electronic excitation gap of eV and has the highest ionization energy of all atoms, eV. This is because helium has no core electrons, so its valence electrons are bound more strongly than in heavier atoms where screening effects play a role. Given such a strong binding, extreme pressures are needed to reach metallization. In fact, after neon boettger (); needs06 (), solid helium is expected to have the highest metallization pressure among all elemental solids.
Metallic solid helium is expected to be present in the outer layers of white dwarfs (WD) chab05 (). After the initial star has exhausted all its nuclear fuel, it sheds its outer layer and leaves behind a dense carbon-oxygen core of the size of the earth that is surrounded by an envelope of pure helium, hydrogen, or a mixture. The fossil star then spends the remaining of its lifetime cooling until vanishing luminosity. Measuring the luminosity of the oldest WD would therefore constrain the age of the galaxy Oswalt (), which qualifies WDs as stellar chronometers.
Extracting the correct physics from WDs depends on how consistent the cooling models are with the observed luminosity isern (). Characterizing helium at high pressure is important because its properties regulate the heat transport across the outer layers. The metallization transition is important because it marks the point where the heat transfer switches from electronic conduction in interior WD layers to photon diffusion in the exterior.
Most WD models rely on semi-analytical descriptions in the chemical picture chab05 () where one treats helium as a collection of stable atoms, ions, and free electrons interacting via approximate pair potentials. While such approaches work well at low density, they cannot describe the complex interactions in a very dense system, and a more fundamental description is required instead. First-principle methods, such as density functional theory and diffusion quantum Monte Carlo, are necessary as they provide a full accounting of quantum and statistical laws that govern the electrons and nuclei.
Recent calculations by Kietzmann et al. redmer () relied on DFT and the generalized gradient approximation (GGA) gga () to calculate the electrical conductivity and locate the insulator-to-metal transition in dense fluid helium. Kowalski et al. Kowalski () went beyond GGA to study the EOS, the electrical and optical properties of fluid helium up to g/cm. Stixrude and Jeanloz loz08 () studied the band gap closure in fluid helium over a wide range of densities including conditions of giant planet interiors.
In this Letter, we study the metallization in solid helium at densities above g/cm using several first-principle simulation techniques. Our intention is to give an assessment of the accuracy of the widely used GGA for calculating total energies and band gaps at extreme conditions. For this reason, we use the accurate, but expensive diffusion quantum Monte Carlo (DMC) method and compare with the GW approximation to correct the GGA band gaps. We also derive the equation of state (EOS) and study the effect of the zero point motion of the nuclei on the band gap closure. We expect our first-principle study to serve as a guide for future laser experiments that are planned to extend the EOS measurements to very high pressures nellisloz ().
Our DFT calculations were performed with the ABINIT plane-wave basis code abinit (). The electron-nuclei interactions were treated by a local Troullier-Martin norm-conserving pseudopotential fhi98pp () with a core radius of 0.4 a.u. We use Perdew-Burke-Ernzerhof GGA gga () for the exchange-correlation functional. We worked with an 8x8x8 Monkhorst-Pack k-point grid and a plane wave energy cutoff of Ha.
Mao et al. mao93 () have demonstrated experimentally that solid helium still remains in the hcp structure up to high pressures ( GPa), apart from a limited fcc loop along the melting line at low temperatures. So we adopt the hcp solid phase and optimize the cell geometry iteratively until the pressure has converged to within .
The band structure in the inset in Fig. 1 shows that solid helium in the hcp structure exhibits an indirect band gap. The lowest unoccupied state occurs at the wave vector. The highest occupied state is at , which is very close to the M wave vector. Subsequently, we approximate the excitation gap by the difference in energy between the valence M point and the conduction point. As Fig. 1 shows, the band gap decreases almost linearly with density. GGA predicts the band gap closure at a density of g/cm, which corresponds to a pressure of TPa.
Kohn-Sham GGA is known to systematically underestimate the band gap. DMC is expected to predict the band gap width more accurately, because it explicitly includes electronic correlation effects. In fact, DMC has been used successfully to describe the electronic ground state in weakly as well as in strongly correlated systems foulkes (). Excited states were also calculated reliably with DMC Martin94 (); needs06 (). One needs a large supercell to describe all correlation effects, which comes at a high computational cost since DMC scales as high as , where is the number of electrons.
Our DMC simulations were performed with the CASINO code casino (). In DMC, the Schrödinger equation is solved stochastically by simulating branching and diffusion in imaginary time. A trial wave function enters in the propagation of electronic configurations and we use the fixed-node approximation to avoid the fermion sign problem. Our trial wave function is of the Slater-Jastrow form. The parameters in the Jastrow factor were optimized by variance minimization. They comprise electron-electron, electron-nucleus and electron-electron-nucleus terms. that the We use GGA orbitals for the Slater part. We also keep the same pseudopotential as in DFT calculations. We picked a conservative high energy cutoff of Ha. For efficiency, the orbitals are represented numerically using blip functions alf04 (). Our results are well converged with a time step of a.u. To calculate the band gap, we promote one electron, either spin-up or -down, from the valence band to the conduction band (). The calculations for the excited state used the same Jastrow parameters as for the ground state calculation because the DMC energy depends only on the nodes of the trial wavefunction that are determined by the Slater determinant. To include the and M wave vectors into one DMC band gap calculation, we have chosen rectangular 2x1x1 or 4x2x2 supercells with, respectively 16 or 128 electrons. We also performed simulations with up to 3x3x3 triangular supercells with 108 electrons to determine the ground state energies and to study their finite-size dependence.
Figure 1 shows the DMC band gaps to be larger than the GGA results, as expected. The DMC curves for the 2x1x1 and 4x2x2 rectangular supercells lie parallel to the GGA curve, hence showing a gap correction that is independent of density. DMC simulations in a 4x2x2 rectangular cell predict a metallization density of g/cm. We consider the 4x2x2 gap results to be converged with respect to system size because the ground state energy agrees well with triangular 3x3x3 supercells and the remaining finite size corrections are small.
The 4x2x2 DMC band gaps agree very well with the GW results, which also show a linear behavior with density similar to all previous calculations. The GW approximation has proven to be a reliable method for correcting the GGA band gaps in a variety of materials johnson (). Our GW band gap corrections are calculated within an accuracy of eV after converging the number of bands () and plane waves ( Ha). In comparison, our metallization density is significantly higher than the linear-muffin-tin-orbitals prediction of g/cm for helium in the fcc phase marvin ().
We tried improving the nodes of the trial wave function in DMC by adding a backflow correction to the Slater determinant. This method introduces further correlation to the trial wave function by replacing the electron coordinates in the determinant by a set of collective coordinates. The DMC total energy decreased slightly but the band gap did not change within error bars ( eV).
We also studied whether the zero point motion of the nuclei has any effect on the metallization density because one might expect the disorder introduced by the zero point motion to reduce the metallization density, as already noted in the fluid Kowalski (); loz08 (). We generated series of configurations with path integral Monte Carlo (PIMC) pimc () simulations for different temperatures between 500 and 5000 K. The helium atoms interact with an effective pair potential that we constructed by matching the forces forcematching () of a density functional molecular dynamics (DFT-MD) simulation at 5000 K. We verified the accuracy of the potential by comparing the original DFT-MD pair correlation function with that of classical Monte Carlo simulations. We then computed the GGA band gaps for the PIMC configurations and compare with perfect hcp lattice results. Within error bars, the zero point motion had no effect on the band gap.
DMC is computationally intensive which limits the size of the supercell used to simulate the infinite solid. Our biggest supercells consist of 3x3x3 triangular ( electrons) and 4x2x2 rectangular ( electrons) primitive unit cells. We correct the energies by considering the following finite size effects. To first order, the independent-particle finite size effects (IPFSE) dominate needs99 (). The error arises from an incomplete k-point sampling of the Brillouin zone in the DMC supercell. In DFT, the computational cost is directly proportional to the number of k-points because the electrons are treated as independent particles. In DMC, however, bigger supercells are needed to include more k-points. So the limitation is more severe since the computational cost in DMC scales with the number of particles in the supercell as O(aN + bN). Our IPFSEs range from eV/el up to eV/el with increasing density in our 3x3x3 triangular cell. We also include kinetic energy corrections of the order of eV/el, which are due to long range correlations chiesa ().
We correct for Coulomb finite size effects (CFSE) that are due to the long range interaction between charged particles and their periodic images. These effects decay with system size as and tend to lower the energy slightly. We reduce these effects by using the model interaction potential (MPC) instead of the Ewald interaction foulkes (). The difference between the total ground state energies computed with these two interactions is eV/el for the largest volume. The error increases with density because the periodic image charges are closer in smaller supercells. We obtain very similar band gaps with the Ewald and the MPC interactions. This is due to the cancellation of CFSE errors when taking energy differences to calculate the band gap.
After correcting the DMC ground state energies for finite size effects, we compare with DFT calculations in Fig. 2. We use two ways to estimate the correction energy that is missing in GGA. First we report the DFT-DMC difference directly for our larger 3x3x3 supercell. Secondly we extrapolate the DMC energies to infinite size as function of . We derive an uncertainty of the resulting DMC energies by comparing linear and quadratic extrapolation. The resulting error bars in Fig. 2 are small except for the highest density. In the density range of consideration, the correlation energy error in GGA is approximately constant, 0.36 eV/el.
The correlation error in GGA becomes less important with increasing density because it stays constant while the energy increases with compression. This trend is illustrated in Fig. 2 where we relate the DFT-DMC energy difference to the energy needed to compress helium to density, . The plotted ratio, , approaches unity in the low density limit. In the opposite high density limit, the graph shows how it tends to zero as helium approaches the state of a one-component plasma with a neutralizing background. The kinetic energy of homogeneous electron gas dominates the correlation and eventually the Coulombic energy terms. GGA is expected to describe this limit well since the exchange-correlation functional was derived from DMC simulations of the homogeneous electron gas exc ().
Since the corrections to the GGA energy appear to be independent of density, there are no corrections to pressures derived from GGA. We were able to represent our zero-temperature static lattice GGA pressure data in the insulating regime by a Vinet EOS curve Vinet () with the parameters cm/mol as zero-pressure volume, GPa as bulk modulus, and as its derivative. The fit reproduced our GGA data from 3 to 1200 GPa with an accuracy of 3%. The comparison with experiments loubeyre () was studied in great detail earlier lowPsim ().
It was not possible to extend the Vinet fit into the metallic regime. Instead, we adopt a fit based on the parametrization of the homogeneous electron gas energy. This fit includes the kinetic, Coulombic, exchange as well as correlation terms and, in this case also ionic contributions, . In units of GPa and cm/mol, the coefficients are , , , and where the leading coefficient is taken from the free Fermi gas. The fit reproduces our DFT data points from 1.2 to 1600 TPa within . In Fig. 3, we show the pressure over a large density range. In the high density limit, the correlation effects decrease and the DFT pressure approaches the free homogeneous electron gas behavior.
In conclusion, we have demonstrated that solid helium reaches a metallic state at an extreme pressure of TPa, which is significantly larger than predicted by standard GGA method. For WD interiors, this implies that the inner layer of metallic helium is thinner and the outer region where photon diffusion dominates the heat transport is larger than previously predicted.
With quantum Monte Carlo we have shown that standard GGA methods underestimate the band gap in solid helium by eV, which translates into an underestimation of the metallization pressure by . The GW band gap corrections are in good agreement with DMC calculations, which offers the possibility of using GW for correcting the band gaps derived from GGA simulation of fluid helium at high pressure and to make more realistic comparisons with shock wave measurements of conductivity and reflectivity.
Finally, we determined the equation of state and presented a fit. We analyzed the correlation effects that are missing in GGA and demonstrated that their importance decreases relative to the total energy with increasing density as helium approaches the state of a one-component plasma with a rigid neutralizing background.
This work was supported by NSF and NASA. We thank the Cambridge group for the CASINO code. NERSC provided computational resources. We acknowlegde useful discussions with E. Quataert, P. Chang, and R. Jeanloz.
- (1) J. C. Boettger, Phys. Rev. B 33 6788 (1986).
- (2) N. D. Drummond and R. J. Needs, Phys. Rev. B. 73, 024107 (2006).
- (3) C. Winisdoerffer and G. Chabrier, Phys. Rev. E, 71, 026402, (2005); V. E. Fortov et al, JETP, 97, 259, (2003).
- (4) T. D. Oswalt, J. A. Smith, M. A. Wood, P. Hintzen, Nature 382, 692, (1996).
- (5) J. Isern, E. Garcìa-Berro, M. Hernanz, G. Chabrier, ApJ, 528, 397, (2000); Brad M. S. Hansen, ApJ, 520, 680, (1999); S. K. Leggett, M. T. Ruis and P. Bergeron, ApJ, 497, 294, (1998).
- (6) A. Kietzmann et al., Phys. Rev. Lett. 98, 190602 (2007).
- (7) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (8) P. M. Kowalski et al., Phys. Rev. B 76, 075112 (2007).
- (9) L. Stixrude and R. Jeanloz. Submitted to Proc. Nat. Ac. Sci. (2008).
- (10) W. J. Nellis et al., Phys. Rev. Lett. 53, 1248 (1984); J. Eggert et al, Phys. Rev. Lett. 100, 124503 (2008).
- (11) X. Gonze, et al, Zeit. Kristallogr. 220, 558-562 (2005).
- (12) M. Fuchs and M. Scheffler, Comput. Phys. Commun. 119, 67 (1999).
- (13) P. Loubeyre et al,Phys. Rev. Lett 71, 2272 (1993).
- (14) W. M. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- (15) L. Mitàš and R. M. Martin, Phys. Rev. Lett. 72, 2438 (1994).
- (16) B. R. Oppenheimer, N. C. Hambly, A. P. Digby, S. T. Hodgkin, D. Saumon, Science, 292 698 (2001).
- (17) R. J. Needs, M. D. Towler, N. D. Drummond and P. Lopez Rios, CASINO 2.1 (Univ. of Cambridge, 2005).
- (18) D. Alfè, M. J. Gillan, Phys. Rev. B 70, 161101(R) (2004).
- (19) K. A. Johnson and N. W. Ashcroft, Phys. Rev. B 58 15548 (1998).
- (20) D. A. Young, A. K. McMahan, and M. Ross, Phys. Rev. B 24 5119 (1981).
- (21) D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- (22) P. R. C. Kent et al., Phys. Rev. B 59, 1917 (1999).
- (23) S. Chiesa, D. M. Ceperley, R. M. Martin, M. Holzmann,Phys. Rev. Lett. 97, 076404 (2006).
- (24) D. M. Ceperley and B. J. Alder, Phys. Rev. lett. 45, 566 (1980).
- (25) P. Vinet et al., Phys. Rev. B 35, 1945 (1987).
- (26) P. Loubeyre et al., Phys. Rev. Lett. 71, 2272 (1993).
- (27) C. P. Herrero, J. Phys. Condens. Matter 18 3469 (2006); C. Cazorla and J. Boronat ibid. 20 015223 (2008)
- (28) S. Izvekov et al., J. Chem. Phys. 120, 10896 (2004).