On the possibility of metastable metallic hydrogen

TL;DR

This paper investigates whether metallic hydrogen can exist in a metastable state at lower pressures using first-principles simulations, finding it likely not to be stable at ambient conditions, which impacts its practical applications.

Contribution

It provides the first computational evidence on the metastability of metallic hydrogen at lower pressures, addressing a key question for its practical use.

Findings

Metallic hydrogen is metastable at high pressures.

It is unlikely to be metastable at ambient conditions.

Implications for fundamental physics and applications are discussed.

Abstract

Metallic hydrogen is expected to exhibit remarkable physics. Examples include high-temperature superconductivity and possible novel types of quantum fluids. These could have revolutionary practical applications. The pressures required to obtain metallic hydrogen, however, are expected to be significant. For practical, and terrestrial applications, a key question is therefore whether this phase is metastable at lower pressures. In this work, this possibility is investigated, using first-principles simulations. The results show that metallic hydrogen is metastable, but strongly suggest that it is not so to ambient conditions. Implications of these results for fundamental physics, and also practical applications of metastable metallic hydrogen are discussed.

Figures

Equations

$$E_{\text{lattice}}=E_{\text{static}}+E_{\text{vib}},$$

$${\omega }_j(\mathbf{q}{)}^2\ge 0$$

Full text

In 1935, Wigner and Huntington predicted Wigner and Huntington (1935) that sufficient pressure would dissociate hydrogen molecules; and any Bravais lattice of such atoms would be metallic. Initial interest in metallic hydrogen was primarily related to astrophysical problems Baraffe et al. (2010). Eventually, and to more recently though, there has been significant interest in this phase at relatively low temperatures McMahon et al. (2012). This can be attributed to the remarkable properties that are expected. This includes, for example, high-temperature superconductivity Ashcroft (1968); McMahon and Ceperley (2011a, 2012). The possibility of a quantum liquid ground-state has also been suggested Bonev et al. (2004); in which case, hydrogen may have quantum ordered states that represent novel types of quantum fluids Babaev et al. (2004). Applications of the (expected) remarkable physics could revolutionize several fields; such as energy, technology, and rocketry Silvera and Cole (2010).

The pressures involved to dissociate hydrogen molecules are expected to be significant though [$447(3)$ GPa McMinis et al. (2015)]. Particularly important for practical, and terrestrial applications is therefore whether metallic hydrogen is metastable. This possibility was suggested in 1972 Brovman et al. (1972); Kagan and Brovman (1972). However, it has received little Geng et al. (2012) to no attention since then; both at low (to zero) pressures and higher (but below molecular dissociation).

The theoretical capabilities available in 1972, however, limited the calculations Brovman et al. (1972) to an approximate, perturbative treatment; and consideration of only one- (and few two-) atom lattices. Theoretical capabilities have significantly progressed McMahon et al. (2012) though; properties can now be accurately calculated from first principles. Considering this problem with modern theoretical methods is therefore critical for addressing this gap in knowledge; and also providing important information for experiments. This is particularly timely, with recent claims Dias and Silvera (2017) (see also Refs. Goncharov and Struzhkin, 2017; Eremets and Drozdov, 2017; Loubeyre et al., 2017; Silvera and Dias, 2017) of the production of metallic hydrogen at low temperatures in static, diamond anvil cell experiments.

The answer to the question of metastability presupposes solution to a series of interrelated problems:

1. (1)

Determination of the minimum-energy crystal structure(s). This includes proof that they lie at stationary points (with respect to all parameters). 2. (2)

Proof of (dynamic) stabilities. 3. (3)

Analysis of the relation between the ground- and metastable-state structure(s) (assuming that such exist), and the molecular phase. This includes that of the processes that are undergone during transitions between them, and also when the pressure is removed. 4. (4)

Determination of the lifetime of this phase.

In this Article, we present results from first-principles simulations in support of this answer. Particular focus is given to problems ($1$) and ($2$) above; and indirectly to ($3$). A discussion of these results is given, and conclusions are drawn; including implications for fundamental physics, and also practical applications. A discussion of the methods used follows the main text. A Supplementary Information (SI) accompanies this work.

I Results

I.1 Body-Centered Tetragonal Structures

The most promising candidate structures for metastable metallic hydrogen (at least at higher pressures; and likely what one would find experimentally) are the most stable ones just above molecular dissociation. While these are not known experimentally, first-principles simulations McMahon and Ceperley (2011b) suggest that they can be represented by a body-centered tetragonal (BCT) structure (possibly including lower-symmetry distortions Geng et al. (2012)); as shown in Fig. 1.

In this representation, there are five structures of note (in atomic hydrogen), defined by the axial ratio $c/a$; three of which have been considered before: those of $\beta$-Sn ($c/a<1$), diamond ($c/a=\sqrt{2}$), and Cs-IV ($c/a>1$), and the other two with $c/a\ll 1$ and $\gg 1$. Their structural parameters are reported in the SI. An analysis of (only) these structures reveals several general, qualitative features about metastable metallic hydrogen.

Consider the phase diagram (at [math] K), below molecular dissociation; shown in Fig. 2.

The pressure range can be (qualitatively) divided into three regions:

1. (a)

$>300$ GPa: the relative enthalpies of all of the structures remain relatively flat; the only exception being the diamond structure (discussed in more detail below). 2. (b)

$200$–$300$ GPa: the enthalpies of the $\beta$-Sn, diamond, and Cs-IV structures become nearly degenerate; while those with $c/a\ll 1$ and $\gg 1$ begin to decrease and increase, respectively. 3. (c)

$<200$ GPa: the $\beta$-Sn, diamond, and Cs-IV structures become degenerate (precisely at $100$ GPa); while those with $c/a\ll 1$ and $\gg 1$ sharply decrease and increase, respectively.

This division (as meaningful) is supported by the results of additional calculations; presented and discussed below, and also in Section I​ B​ $1$ of the SI.

Additional insight (including into Fig. 2) is obtained by considering the energy of the BCT structure over the entire range of $c/a$; shown in Fig. 3.

As will be supported by the results of additional calculations (below), there are two especially important points that are suggested by these results:

First: the collapse of energy barriers in the potential-energy surface (PES) of (atomic) metallic hydrogen, between $300$ and $200$ GPa. Consider the diamond structure. At relatively high pressures, this structure forms an energy barrier (a maximum) between those of $\beta$-Sn and Cs-IV. While Fig. 2 shows that its relative enthalpy decreases continuously with decreasing pressure, it is between $300$ and $200$ GPa that it (as an energy barrier) does so appreciably. And, by $100$ GPa, it has collapsed completely, leaving only a single minimum in this region of $c/a$ ($\approx\sqrt{2}$). Notice in Fig. 2 that the $\beta$-Sn, diamond, and Cs-IV structures here become degenerate in enthalpy.

Second: the tendency towards the formation of molecules, below $200$ GPa. Consider the structures with $c/a\ll 1$ and $\gg 1$. The sharp increases in energy as $c/a\rightarrow 0$ and $\rightarrow\infty$ are due to the Coulomb repulsion resulting from Brillouin planes with like charges that approach each other. This (approach) though is also related to the formation of molecules; resulting in a minimum of energy, prior to its sharp increase. This explains the sharp decrease in energy of the $c/a\ll 1$ structure. Note the nearest-neighbor distance (at [math] GPa) of $0.9900$ Å.

The ability of atoms to approach each other is the result of the collapse of (atomic) energy barriers (discussed above). While this occurs, there is also a hardening of molecular ones. Notice the relative heights (and their trends with pressure) of the maxima separating the $c/a\ll 1$ and $\gg 1$ structures from the atomic ($c/a\approx\sqrt{2}$) region.

I.2 Structure Prediction

While the BCT structures (discussed above) are the most promising candidates at pressures near molecular dissociation, there is no guarantee that they are representative of those (potentially) so at lower ones. Consider the lattice energy $E_{\text{lattice}}$ of metallic hydrogen; in the adiabatic approximation:

[TABLE]

where $E_{\text{static}}$ is the static (electronic ground-state energy, with a fixed lattice) contribution, and $E_{\text{vib}}$ is that resulting from atomic vibrations. In metallic hydrogen (and unlike ordinary metals), $E_{\text{static}}$ favors anisotropic structures; see Ref. Brovman et al., 1972 for a discussion. $E_{\text{vib}}$, on the other hand, favors symmetric ones. The subtle (and important) details of these terms, including their relative magnitudes, depends on the pressure. At low pressures, $E_{\text{lattice}}$ is determined largely by $E_{\text{static}}$. It is with increasing pressure that $E_{\text{vib}}$ becomes (especially) important.

Searches were performed for metallic (atomic and/or mixed atomic/molecular) structure(s) with the lowest $E_{\text{static}}$ at [math] GPa. These revealed several candidate structures of (potentially) metastable metallic hydrogen. An analysis of the most stable ones over the entire pressure range further supports the results above; such is provided in Section I​ B​ $1$ of the SI. At present, however, it is relatively low pressures that are of main interest.

Below $200$ GPa, all structures found show a tendency towards the formation of molecules. This is consistent with the (second) point suggested by Fig. 2 (above). It is supported though by a detailed analysis (below) of the structures themselves.

In searches with more than a few atoms, only molecular or mixed atomic/molecular structures were found. In addition, for the latter, as the number of atoms increases, there is a greater proportion of molecules (to atoms). This trend is also reflected in their relative enthalpies. Results and discussion for the most stable of these structures are provided in Section I​ B of the SI. That this trend occurs means that such structures may not be at stable stationary points.

The results above suggest that any metastable structure(s) of metallic hydrogen are (completely) atomic. Such were found in searches with a few or less atoms. The most stable (and nearly degenerate in energy) ones are shown in Fig. 4.

Notice that they are all very similar in their symmetries.

The tendency towards the formation of molecules though can be seen, even in these structures. Notice the close proximity of atoms (into the page); this is approximately $0.9914$ Å. This is analogous to the results for the $c/a\ll 1$ BCT structure (discussed above). This means that these structures too may not be at stable stationary points; confirmed, in results presented below. Further results and discussion about these structures themselves is provided in Section I​ B​ $3$ of the SI.

I.3 Dynamic Stabilities

A crystal is dynamically (mechanically) stable, if it executes a stable vibrational motion about its equilibrium configuration Born and Huang (1954). By solving the equations of motion for the normal modes of vibration, this is equivalent to the condition that their frequencies $\omega$ satisfy:

[TABLE]

for all phonon branches $j$ and wavevectors $\mathbf{q}$. (In other words, they are real; for an imaginary frequency means that the system, subject to a small displacement, will disrupt exponentially with time.) Dynamic stabilities were determined by calculating the phonon density of states $F(\omega)$ (along both the real and imaginary axes), for each structure.

Consider first the BCT structures; Fig. 5 shows $F(\omega)$ of Cs-IV.

Note that all other BCT structures are dynamically unstable; see Section I​ C​ $1$ of the SI, for results and discussion about the $\beta$-Sn structure.

Figure 5LABEL:sub@fig:PHDOS_tetragonal:stable shows pressures from $500$ (down) to $250$ GPa. The absence of imaginary frequencies proves dynamic stability (at least to these pressures). Consider also the (qualitative) changes in the spectra. To $350$ GPa, they (mostly) retain their shape, with frequencies simply shifting to lower values. Such is expected for a structure that remains qualitatively the same (in terms of its underlying PES). By $250$ GPa, though, while there continues to be a shift to lower frequencies, differences become apparent.

Figure 5LABEL:sub@fig:PHDOS_tetragonal:unstable shows pressures below $250$ GPa. By $200$ GPa, a small (dynamic) instability appears. This then increases, with decreasing pressure. The spectra remain (qualitatively) similar though (at least until [math] GPa). This behavior, along with that seen in Fig. 5LABEL:sub@fig:PHDOS_tetragonal:stable (discussed above), is consistent with a transition region between $300$ and $200$ GPa, where the underlying PES of hydrogen changes significantly.

Consider now the candidate structures; Fig. 6 shows $F(\omega)$ for those shown in Fig. 4.

The structures exhibit very similar spectra (as expected, given their similarities). Notice their tails that extend to very high imaginary frequencies. While these have about an order of magnitude smaller magnitude, an expanded view reveals peaks. Therefore, while the most stable candidate (atomic) structures (Fig. 4) are at stationary points, they are unstable ones. This is further supported by additional results; see Section I​ C​ $2$ of the SI. Note that also presented and discussed therein are dynamic stabilities of these structures at higher pressures.

II Discussion

The possibility of metastable metallic hydrogen has been investigated in detail. This was considered from two approaches, by an analysis of: the BCT structures (stable above molecular dissociation); and the prediction of new candidate ones at [math] GPa. The results obtained are self-consistent, and provide strong evidence regarding this possibility.

The results in Section I are concisely summarized by the processes undergone on the PES of hydrogen, as the pressure is removed (following molecular dissociation):

1. (a)

At (relatively) high pressures ($>300$ GPa): the PES of (atomic) metallic hydrogen is well defined; and it remains qualitatively similar below molecular dissociation. 2. (b)

At intermediate pressures ($200$–$300$ GPa): the energy barriers of metallic hydrogen on the PES collapse; while those of the (completely) molecular phase harden. 3. (c)

At lower pressures ($<200$ GPa): the PES of molecular hydrogen is well defined.

From these results, we conclude that (atomic) metallic hydrogen is metastable; but only to $300$–$200$ GPa. Below $200$ GPa, metallic hydrogen has no region of stability.

Metastable metallic hydrogen, even at these pressures, may have important practical applications. Such pressures are significantly less than those expected for molecular dissociation McMinis et al. (2015), and are possible to study under static conditions Eremets (1996). The importance though depends on its properties.

There are therefore several important, open questions. One mentioned, are the properties of metastable metallic hydrogen; for example, superconductivity. Whether this phase is even realizable though depends on its lifetime; and such properties on its finite-temperature stability. Answering these, and related questions, will be the subjects of future work.

III Methods

III.1 Electronic-Structure Calculations

All calculations were performed using the Quantum ESPRESSO (QE) density-functional theory code Giannozzi et al. (2009). The pseudopotential method based on the projector augmented-wave method Kresse and Joubert (1999) was used 111We used the pseudopotential H.pbe-kjpaw_psl.0.1.UPF from http://www.quantum-espresso.org/., to replace the bare Coulomb potential of the proton. The Perdew–Burke–Ernzerhof generalized gradient approximation exchange and correlation functional Perdew et al. (1996) was used. QE is based on a plane-wave basis set, and we used kinetic-energy cutoffs of $57.5$ Ry for the wavefunction and $345.5$ Ry for charge density and potential. For Brillouin-zone sampling, at least $32{\times}32{\times}32$ (unshifted) $\mathbf{k}$ points were used. The smearing scheme of Methfessel–Paxton Methfessel and Paxton (1989) was used for Brillouin-zone integrations, with a smearing width of $0.02$ Ry. These choices give a total convergence, in energy, for example, to better than $0.5$ meV$/$proton.

III.2 Geometry Optimizations

Stationary points of structures were found by performing constant-pressure geometry optimizations. These were done using the Broyden–Fletcher–Goldfarb–Shanno algorithm Fletcher (1987), as implemented within QE. Energies, forces, and pressures were converged to $10^{-5}$ Ry, $10^{-4}$ Ry$/$a.u., and $0.5$ kbar, respectively.

III.3 Structure Prediction

Calculations were based on the ab initio random structure searching approach Pickard and Needs (2011). These were performed using the electronic-structure calculations and geometry optimizations, as described above. For the former, the basis-set cutoffs were reduced, to $46$ and $221$ Ry, respectively; as was the Brillouin-zone sampling, to $8{\times}8{\times}8$ (shifted) $\mathbf{k}$ points. For the latter, the energy and force convergence-criteria were relaxed by an order of magnitude. These are (essentially) default values, capable of generating candidate structures for further analysis.

Random structures were constructed by generating random primitive (lattice) vectors, renormalizing the volume (to the expected average), and proton configurations. Searches were performed over unit cells containing $1$ to $10$ atoms. In particular, geometry optimizations (at [math] GPa) were performed; for each unit-cell size, $1000$ trials. Symmetries of the (relaxed) structures were then determined with the FINDSYM program Stokes and Hatch (2005).

III.4 Dynamic Stabilities

Phonons were calculated using density-functional perturbation theory, as implemented within QE. (Phonon) density of states $F(\omega)$ were calculated over a grid of $4{\times}4{\times}4$ phonon wavevectors $\mathbf{q}$. Such a $\mathbf{q}$-point grid is estimated to be twice as dense (in each direction) as sufficient to accurately calculate $F(\omega)$ of metallic hydrogen McMahon and Ceperley (2011b, a, 2012).