Stability and metallization of solid oxygen under high pressure

TL;DR

This study uses first principles theory to map oxygen's phase diagram under high pressure, revealing a stable metallic phase and resolving longstanding discrepancies between theory and experiment.

Contribution

It introduces a comprehensive theoretical analysis including anharmonic effects and exchange interactions, clarifying oxygen's phase stability and metallization under high pressure.

Findings

Identification of a stable metallic P6_3/mmc phase at high pressures and temperatures.

Resolution of previous disagreements on epsilon-O8 stability and transition pressures.

Accurate modeling of anharmonic and thermal effects crucial for phase predictions.

Abstract

The phase diagram of oxygen is investigated for pressures from 50 to 130~GPa and temperatures up 1200 K using first principles theory. A metallic molecular structure with the $P6_3/mmc$ symmetry ($\eta^{'}$ phase) is determined to be thermodynamically stable in this pressure range at elevated temperatures above the $\epsilon$(${O_8}$) phase. Long-standing disagreements between theory and experiment for the stability of $\epsilon$(${O_8}$), its metallic character, and the transition pressure to the $\zeta$ oxygen phase are resolved. Crucial for obtaining these results are the inclusion of anharmonic lattice dynamics effects and accurate calculations of exchange interactions in the presence of thermal disorder.

Figures

Full text

Oxygen has received a great deal of attention because it is a fundamental element, one of the most abundant on earth and the only one known with a diatomic molecule that carries a magnetic moment  Lundegaard et al. (2006). It has a rich high pressure ($P$) phase diagram with multiple solid phases exhibiting diverse physical properties  Freiman and Jodl (2004); H rl (1962); LeSar and Etters (1988); Lundegaard et al. (2009); Meier and Helmholdt (1984); Schiferl et al. (1983); Gorelli et al. (2002); Goncharenko et al. (2004); Desgreniers et al. (1990); Shimizu et al. (1998); Zhu et al. (2012); Nomura et al. (2014). One of the most interesting among them is the $\epsilon$-phase, which is stable over a large pressure range. It has been studied extensively both theoretically  Serra et al. (1998); Gebauer et al. (2000); Neaton and Ashcroft (2002) and experimentally by X-ray diffraction Johnson et al. (1993); Weck et al. (2002) and spectroscopic measurements  Akahama and Kawamura (1996); Gorelli et al. (1999); Akahama and Kawamura (2000); Agnew et al. (1987). Despite previous theoretical works suggesting that it has $C2/m$ symmetry  Johnson et al. (1993); Weck et al. (2002); Akahama et al. (1995) and its strong infrared absorption Freiman and Jodl (2004) alluding to $O_{2}$ molecules forming larger units Gorelli et al. (1999); Gebauer et al. (2000), only recent experiments have defined the exact structure as $O_{8}$ clusters ($\epsilon$(${O_{8}}$))  Lundegaard et al. (2006); Fujihisa et al. (2006); Meng et al. (2008). Upon compression to 96 GPa, $\epsilon$(${O_{8}}$) tranforms to the metallic $\zeta$ phase Desgreniers et al. (1990); Akahama et al. (1995), which is even superconducting at 0.6 K  Shimizu et al. (1998). Experimentally, Goncharov et al.  Goncharov et al. (2011) proposed an $\eta^{{}^{\prime}}$ phase in the pressure range of 44 to 90 GPa and at temperatures ($T$) near 1000 K. They suggested it to be an isostructure of the $\eta$-O phase previously proposed at low pressures Goncharenko (2005); Lundegaard et al. (2009).

Most of the previous theoretical studies have been limited to 0 K and show a significant disagreement between the calculated and measured $\epsilon$-$\zeta$ transition pressures, as well as the structures of the $\epsilon$ and $\zeta$ phases Freiman et al. (2018). Ma et al.  Ma et al. (2007), suggested $C2/m$ as the best candidate for $\zeta$-$O_{2}$ but it remains a matter of debate. The $\epsilon$-$\zeta$ transition occurs at only 35 GPa in their study. Theoretical confirmation for the $\eta^{{}^{\prime}}$ phase has not been reported yet. Moreover, the experimental evidence for its stability is not conclusive, and its nature and phase boundaries are unknown.

Here we report results on the phase diagram of solid oxygen for pressure up to 130 GPa, at both 0 K and finite $T$. First, we resolve the existing inconsistencies between theory and experiment regarding $\epsilon$ ($O_{8}$) and the $\epsilon$-$\zeta$ transition. Next, we focus on elevated temperatures where we show that inclusion of anharmonic effects and accurate exchange energy calculations in the presence of thermal disorder are both crucial for determining the stability of $O_{2}$ phases.

We start by examining oxygen at 0 K in the $P$ range of 10-130 GPa and consider the relevant structures $\eta^{{}^{\prime}}$ Goncharov et al. (2011), $\epsilon$ ($O_{8}$) Lundegaard et al. (2006), and $\zeta$ ($C2/m$) Ma et al. (2007). Density-functional theory calculation (DFT)  Kohn and Sham (1965) are performed with ABINIT Gonze et al. (2009), using Troullier-Martins Pseudopotentials Troullier and Martins (1991) a plane-wave expansion with a 80-Hartree cut off, and k-point grids of $16^{3}$, $4^{3}$, and $12\times 12\times 10$ for $\eta^{{}^{\prime}}$, $\epsilon$ ($O_{8}$), and $\zeta$ ($C2/m$) phases, respectively, ensuring enthalpies convergence to better than 1 meV/atom. In order to test the effect of the exchange correlation functional approximation of DFT on the relative stability of these structures, calculations with both the generalized gradients approximation (GGA) and the local density approximation (LDA) were performed (Fig. 1). The GGA transition between $\epsilon$ ($O_{8}$) and $\zeta$($C2/m$) is at 35 GPa, which is in agreement with the previous theoretical results Ma et al. (2007). Note that during structural optimization at $P>50$ GPa, we have restricted the occupation of the electronic states of $\epsilon$(${O_{8}}$) in order to prevent its spontaneous transformation to $\zeta$($C2/m$) as reported by Ma et al.  Ma et al. (2007). Within LDA, the $\epsilon$(${O_{8}}$) – $\zeta$($C2/m$) transition is 25 GPa and it is clear that $\epsilon$(${O_{8}}$) is the most sensitive to the choice of exchange correlation functional. To examine the significance of spin polarization, we have also performed spin-polarized calculations within GGA and LDA (Fig. 1(c)). The results show negligible effects on the enthalpies, which is in agreement with previous studies reporting that the $O_{2}$ spin is suppressed at pressures above 10 GPa Goncharenko (2005).

It is well known that the GGA and LDA introduce errors in the ground state energy calculations that depend on the electronic properties of the system; the overall tendency is to favor better metals. Therefore, the fact that O2 undergoes metallization in the pressure range of interest and the notable differences among the electronic properties of the competing structures ($\epsilon$(${O_{8}}$) is an insulator while $\zeta$($C2/m$) and $\eta^{{}^{\prime}}$ are metallic; see Supplementary material), raise the question of whether there is a significant non-cancellation of LDA and GGA errors. To examine this, we have carried out hybrid exchange calculations within the Heyd-Scuseria-Ernzerhof approximation (HSE06)  Heyd et al. (2006) as implemented in VASP  Kresse and Hafner (1993). Full structural optimizations within HSE06 were performed using a 1000 eV plane-wave cut-off and slightly reduced k-point grids: ($12\times 12\times 4$) for $\eta^{{}^{\prime}}$ and ($6\times 6\times 6$) for $\zeta$($C2/m$), ensuring convergence of relative enthalpies to better than 2 meV/atom. Additionally, HSE06 electronic band structures were calculated on both HSE06 and GGA relaxed structures.

The results with HSE06 corrections are shown in Fig. 2. The relative enthalpy of $\epsilon$($O_{8}$) is most strongly affected. At around 90 GPa, it is lowered by as much as 53 and 44 meV/atom relative to $\eta^{{}^{\prime}}$ and $\zeta$($C2/m$), respectively. These corrections are sufficient to make $\epsilon$($O_{8}$) the preferred structure in the $P$-$T$ region where it has been observed experimentally. The large effect on $\epsilon$($O_{8}$) is understood by the fact that (at 0 K) this phase is insulating at $P<107$ GPa whereupon further compression it metalizes by way of band overlap. The computed closure of the bandgap at 107 GPa and 0 K is in a good agreement with the experimental observation of metallization at 96 GPa Desgreniers et al. (1990); Akahama et al. (1995); indeed, thermal effects are expected to close the gap at lower pressure. Our analysis (see Supplementary material) indicates that in order to reproduce the experimentally observed metallization of $\epsilon$($O_{8}$) at $P>96$ GPa, it is necessary to perform both structural optimizations and electronic band structure calculations using hybrid exchange. To summarize so far, we are able to obtain the low-$T$ thermodynamic stability and non-metallic character of oxygen in the $\sim$50–100 GPa range in agreement with the experimental observations Lundegaard et al. (2006). However, note that our results do not support the $\epsilon$($O_{8}$)–$\zeta$($C2/m$) transition at around 107 GPa. Even though the two structures are close in energy, this suggests that $C2/m$ may not be the correct structure of the $\zeta$-$O_{2}$ phase.

Having established the phase properties of $O_{2}$ at low $T$, we turn our attention to its thermodynamic stability at elevated $T$. For this, we have first computed the phonon dispersions of the $\eta^{{}^{\prime}}$ and $\epsilon$-$O_{8}$ phases at 50 GPa using Density-Functional Perturbation Theory (DFPT) Gonze (1997) as employed in ABINIT Gonze et al. (2009), with the same convergence parameters used to calculate their enthalpies. For the $\epsilon$(${O_{8}}$) phase, a $4^{3}$ q-point phonon grid for computing the dynamical matrices was sufficient to achieve convergence for Helmholtz free energies and entropy better than 1 meV/atom. For the $\eta^{{}^{\prime}}$ phase, the dynamical matrices were computed on $3^{2}\times 5$ , $3^{2}\times 9$, $5^{3}$, and $7^{2}\times 5$ q-point grids, in all cases yielding imaginary phonon frequencies (see Fig. 3(a)). This could mean that the structure is mechanically unstable and is probably the reason why it was not found in previous structure search studies. However, depending on the nature of the soft phonon modes, it can be stabilized at finite $T$ or even at 0 K due to quantum zero point motion. Hence, we have examined more closely the phonon modes where the instability is most pronounced using the frozen phonon method. The results in Fig. 3(b) show that the unstable mode is indeed a shallow double well potential. The appearance of a double well potential associated with a structural instability is similar to what has been observed in other elements such as $Ca$  Teweldeberhan et al. (2010). However, compared to $Ca$, the potential barrier in $\eta^{{}^{\prime}}$-$O_{2}$ is much higher. Nevertheless, this mode can be thermally stabilized at $T$ of around several hundred K and is likely to contribute significantly to the entropy of the $\eta^{{}^{\prime}}$ phase at elevated $T$. For determining the thermodynamic stability of $\eta^{{}^{\prime}}$-$O_{2}$, it is therefore necessary to go beyond the quasi-harmonic approximation.

Gibbs free energies ($G$) at finite $T$ were computed using first principle molecular dynamic (FPMD) simulations in the $\sim$49–70 GPa $P$ range and $T=500$, 800 and 1200 K, using finite-$T$ DFT Kohn and Sham (1965) within PBE-GGA Perdew et al. (1997) and VASP Kresse and Hafner (1993). The simulations were carried out with 300-atom supercells, the $\Gamma$ k-point, a 6-electron projector augmented wave pseudopotential, and a 900 eV plane-wave cut-off in the canonical $NVT$ (constant number of particles, $N$, volume $V$, and $T$) ensemble using Born-Oppenheimer dynamics with a Nos$\acute{e}$-Hoover thermostat. For each $V$ and $T$, the system was initially equilibrated within $2$ ps and then ran for additional $6$ ps or more using a 0.75 fs ionic time-step. $G$ was calculated as $G=E_{0}+P_{0}V+P_{\mathrm{ph}}V+U_{\mathrm{ph}}-TS$, where $E_{0}$ and $P_{0}$ are the 0 K DFT energy and pressure, $U_{\mathrm{ph}}$ and $P_{\mathrm{ph}}=-\frac{\partial F_{\mathrm{ph}}}{\partial V}|_{N,T}$ are the phonon internal energy and pressure, $S$ is the entropy, and $F_{\mathrm{ph}}=U_{\mathrm{ph}}-TS$. Here $U_{\mathrm{ph}}$ and $S$ are obtained by integrating vibration density of states (VDOS), which for $\eta^{{}^{\prime}}$ are calculated by taking a Fourier transforms of velocity autocorrelation functions (VACF). Although $S$ is calculated using a harmonic partition function, the VDOS from PFMD represent thermally renormalizes phonons and capture most of the anharmonic free energy Teweldeberhan et al. (2010). For $\epsilon$($O_{8}$), where the harmonic approximation is sufficient, VDOS are calculated from DFPT.

A comparison of the VDOS of the $\eta^{{}^{\prime}}$ and $\epsilon$(${O_{8}}$) phases is shown in Fig. 4(a), from which the anharmonicity of the former is evident. The enthalpies of the two phases as a function of $T$ are shown in Fig. 4(b). Within GGA, the enthalpy of $\eta^{{}^{\prime}}$ remains lower in the entire temperature range of interest. The lower frequency modes of $\eta^{{}^{\prime}}$-$O_{2}$ slightly lower its $H$ relative to that of $\epsilon$(${O_{8}}$), but the effect is small and after adding the HSE06 correction computed earlier, the enthalpy of $\epsilon$(${O_{8}}$) becomes lower at all temperatures. The soft modes have much more pronounced effect on the entropy. The $TS$ terms are potted as a function of temperature in Fig. 4(c). As expected, the $\eta^{{}^{\prime}}$ phase has higher entropy. At 750 K, its $TS$ is 37 meV/atom higher than that of the $\epsilon$(${O_{8}}$) phase. In the context of relative stabilities of molecular crystals, this is a relatively large value. However, it is still not sufficient to overcome the enthalpy differences between the two phases and to make $\eta^{{}^{\prime}}$ preferred at elevated temperatures.

In order to find a clue that may solve the problem for the finite temperature stability, we examine the local structural oder of the $\epsilon$(${O_{8}}$) and $\eta^{{}^{\prime}}$ phases. The atomic arrangements in the two structures are shown in the Insets of Fig. 5. In both cases, the molecules are arranged in layers, with their bonds perpendicular to the the layers. Fig. 5(a) shows a distribution of distances between the molecular center of masses (CM) of the 0 K crystals. The CM-CM distributions look quite different and it is clear that $\eta^{{}^{\prime}}$ is the more symmetric structure. However, if we examine the atomic arrangements shown in the Insets of Fig. 5, we see that locally the $\epsilon$(${O_{8}}$) phase can be viewed just as a distortion of the hexagonal molecular arrangements found in the $\eta^{{}^{\prime}}$ phase. This observation suggest that the introduction of thermal disorder may bring the average local order of the two structures closer to each other. CM-CM distance distributions computed from finite-temperature FPMD trajectories (Fig. 5(b)) show that this is indeed the case. Namely, the short range orders of the two structures become similar, on average, in the presence of thermal disorder. As seen in the Insets in Fig. 5, the difference of the electronic density of states between the two structures near the Fermi energy also diminish with temperature.

The implication of the above observations is that the hybrid exchange corrections to the energies computed on ideal 0 K crystals may not be adequate at fine $T$. Indeed, exchange interactions are in principle short-range. In the HSE06 implementation here, the effective range of the Hartree-Fock exchange is around 4 A (hence the CM-CM range shown in Fig. 5). Therefore, if the local orders of the two structures become similar at elevated $T$, then the differences between their HSE06 corrections are also expected to diminish with temperature. This, in turn, will make $\eta^{{}^{\prime}}$-$O_{2}$ more competitive at finite $T$ compared to the 0 K case. We have therefore performed HSE06 calculations on atomic configurations taken from the FPMD trajectories at 1200 K. For each structure, 5-10 equally spaced (in time) configurations were taken randomly from the trajectories and their energies computed within GGA-PBE and HSE06 with exactly the same simulation parameters. We verified that the fluctuations in the energy differences between HSE06 and GGA-PBE are negligible, which indicates that the GGA-PBE ensemble is sufficient for this analysis.

The calculations reveal that at 50 GPa and 1200 K the hybrid exchange correction of $\eta^{{}^{\prime}}$-$O_{2}$ relative to $\epsilon$(${O_{8}}$) is only 9 meV/atom – much smaller than the 0 K value of 60.8 meV/atom. The larger entropy of $\eta^{{}^{\prime}}$ is now sufficient to compensate for this smaller correction and $\eta^{{}^{\prime}}$ is thus the preferred phase at high $T$ (Fig. 6). We have recomputed the Gibbs free energies of the two structures as a function of $T$ by simply interpolating the HSE06 correction between the 0 and 1200 K values. The result shows that $\eta^{{}^{\prime}}$-$O_{2}$ becomes thermodynamically stable at temperatures above about 780 K at 50 GPa, consistent with the experimental observations by Goncharov et al. Goncharov et al. (2011).

In conclusion, we have shown that $\epsilon$(${O_{8}}$) is stable and non-metallic in the $\sim$ 50–100 GPa range, in agreement with measurements. Computing the structural and electronic properties of the competing oxygen structures using a beyond GGA exchange functional is essential for this result. Furthermore, we have determined that the experimentally proposed $\eta^{{}^{\prime}}$ structure is mechanically unstable at low $T$ within a classical ion dynamics treatment due to soft phonon modes. However, it is stabilized at finite $T$ where the thermally renormalized phonons also contribute to it having a relatively large entropy. In the presence of thermal disorder, the differences between the local structural order and electronic properties of $\epsilon$(${O_{8}}$) and $\eta^{{}^{\prime}}$ diminish. The interplay of all these factors - anharmonicity, exchange effects, and thermal disorder, results in $\eta^{{}^{\prime}}$ becoming the thermodynamically stable phase at elevated $T$. At 50 GPa, we predict the transition to take place at around 780 K.

We thank Prof. R. Hoffmann and Dr. V. Askarpour for helpful discussions. This work was supported by NSERC, Acenet and LLNL. S.A.B. performed work at LLNL under the auspices of the US Department of Energy under contract No. DE-AC52-07NA27344.