Quantum Hydrogen-Bond Symmetrization and High-Temperature Superconductivity in Hydrogen Sulfide
Hydrogen compounds are peculiar as the quantum nature of the proton can crucially affect their structural and physical properties. A remarkable example are the high-pressure phases Goncharov et al. (1996); Loubeyre et al. (1999) of HO, where quantum proton fluctuations favor the symmetrization of the H bond and lower by 30 GPa the boundary between the asymmetric structure and the symmetric one Benoit et al. (1998). Here we show that an analogous quantum symmetrization occurs in the recently discovered Drozdov et al. (2015) sulfur hydride superconductor with the record superconducting critical temperature K at 155 GPa. In this system, according to classical theory Duan et al. (2014, 2015); Errea et al. (2015); Li et al. (); Bernstein et al. (2015), superconductivity occurs via formation of a structure of stoichiometry HS with S atoms arranged on a body-centered-cubic (bcc) lattice. For GPa, the H atoms are predicted to sit midway between two S atoms, in a structure with symmetry. At lower pressures the H atoms move to an off-center position forming a short HS covalent bond and a longer HS hydrogen bond, in a structure with symmetry Duan et al. (2014, 2015); Errea et al. (2015); Li et al. (); Bernstein et al. (2015). X-ray diffraction experiments confirmed the HS stoichiometry and the S lattice sites, but were unable to discriminate between the two phases Einaga et al. (2015). Our present ab initio density-functional theory (DFT) calculations show that the quantum nuclear motion lowers the symmetrization pressure by 72 GPa. Consequently, we predict that the phase is stable over the whole pressure range within which a high was measured. The observed pressure-dependence of is closely reproduced in our calculations for the phase, but not for the phase. Thus, the quantum nature of the proton completely rules the superconducting phase diagram of HS.
The discovery of high-temperature superconductivity in compressed hydrogen sulfide Drozdov et al. (2015) has generated intense interest over the last year, and has led to a number of theoretical studies aimed at understanding the phase diagram of the H-S system as well as the origin of the astonishingly high observed Duan et al. (2014); Bernstein et al. (2015); Duan et al. (2015); Papaconstantopoulos et al. (2015); Errea et al. (2015); Akashi et al. (2015); Nicol and Carbotte (2015); Flores-Livas et al. (2015); Li et al. (). The overall consensus is that HS, the only stable compound formed by hydrogen and sulfur at ambient conditions, is metastable at high pressures and its decomposition gives rise to several H-S compounds. High- superconductivity is believed to occur in a structure with the HS stoichiometry, and is considered to be conventional in nature, i.e., mediated by the electron-phonon interaction Drozdov et al. (2015); Duan et al. (2014); Bernstein et al. (2015); Papaconstantopoulos et al. (2015); Errea et al. (2015); Akashi et al. (2015); Nicol and Carbotte (2015); Flores-Livas et al. (2015). Alternatives to conventional superconductivity have also been discussed Hirsch and Marsiglio (2015). According to structural predictions Duan et al. (2014, 2015); Errea et al. (2015); Li et al. (); Bernstein et al. (2015), HS adopts a rhombohedral form between approximately 112 and 175 GPa, and a cubic at higher pressures. As shown in Fig. 1, the phase is characterized by covalently bonded SH units with a covalent HS bond of length . Each of these H atoms is bonded to the next S atom by a hydrogen HS bond of length . The phase, in contrast, has full cubic symmetry, with so that each H atom resides midway between the two S atoms, as shown in Fig. 1. The structure is nevertheless very close to cubic symmetry, for example, the DFT-relaxed structure, which represents the minimum of the Born-Oppenheimer energy surface (BOES), has a rhombohedral angle of 109.49 at 150 GPa, compared to 109.47 for a perfect bcc lattice. We have verified that imposing a cubic angle on the structure has a negligible effect on the energy difference between the and structures. Consequently, we assume a cubic lattice for both phases in the following.
The bond-symmetrizing second-order transition from to occurs at 175 GPa according to our static lattice calculations. At this pressure, our harmonic phonon calculations show that a -point optical phonon of the high-symmetry phase becomes imaginary, implying that is at a saddle point of the BOES between 112 and 175 GPa, while the phase lies at the minimum. Crystal symmetry guarantees that the transition is of second-order type. As occurs in the high-pressure ice X phase Benoit et al. (1998); Lee et al. (1992, 1993) and other hydrogenated compounds McMahon et al. (1990), the quantum nature of the proton can radically alter the pressure at which the second-order phase transition occurs and, in the present case, can strongly affect the stability of the phase below 175 GPa. Determining the stability ranges of these phases therefore requires the inclusion of vibrational zero-point energy (ZPE) alongside the static BOES energy. However, the presence of imaginary phonon frequencies hinders calculations of the ZPE, since the quasi-harmonic approximation breaks down, and anharmonicity becomes a crucial ingredient.
To elucidate the role of anharmonicity and quantum effects in the pressure range in which the record was observed, we make use of the stochastic self-consistent harmonic approximation (SSCHA) Errea et al. (2013, 2014). The variational SSCHA method was devised for calculating the free energy and phonon spectra while fully incorporating quantum and anharmonic effects, and it is therefore perfectly suited for our purpose. All of the calculations presented here are performed at 0 K. Primitive cells for the and structures contain 4 atoms (1 S atom and 3 H atoms), and therefore a particular nuclear configuration can be described by a 12-dimensional vector containing the atomic coordinates. In the classical limit the ZPE is neglected and the energy of a nuclear configuration is given by the DFT Born-Oppenheimer energy . In the SSCHA, the ZPE is accounted for by approximating the nuclear wave-function by a Gaussian centered on a centroid coordinate , which denotes the average and most probable position of the nuclei. For a given , the width of the Gaussian is obtained by a variational minimization of the expectation value of the sum of the nuclear potential and kinetic energies. In the following analysis it is convenient to split the SSCHA total energy into static and anharmonic-vibrational-ZPE contributions: .
We study the energy landscape along the line defined by , where and are, respectively, the coordinates corresponding to the saddle point and minimum of the BOES, representing the two different symmetries. Here, is a real number describing the reaction coordinate, so that at the centroids are located at the atomic positions of , and at at the atomic positions of . Hence, measures the off-centering of the hydrogen nuclear wave-function and can be associated with the relative coordinate that quantifies the length of the HS hydrogen bond with respect to the symmetric position ( is the lattice parameter). We analyze the curve for a fixed primitive bcc unit cell. As shown in Fig. 2(a) for a cell volume of the curve has a shallow double-well structure favoring the structure by only 5.6 meV/HS. However, after adding the energy calculated with the SSCHA, the full curve shows a clear minimum at , which favors the structure. At a larger volume of 102.11 , which corresponds to a pressure of around 130 GPa, the minimum of is also at , despite the fact that the Born-Oppenheimer well in becomes deeper, as shown in Fig. 2(c). Repeating these calculations for DS, we find that the structure is the most favorable once the ZPE has been included. We therefore conclude that the quantum nature of the nuclei symmetrizes the hydrogen bond and leads to a proton wave-function centered at the atomic positions of for both HS and DS. To eliminate the possibility that the energy minimum occurs beyond the line studied, we performed an unconstrained SSCHA minimization, optimizing both the width of the Gaussians and the centroid positions. The results of this minimization show again that, within stochastic error, the centroid position obtained corresponds to the structure, in which the HS covalent and HS hydrogen bond distances equalize, leading to symmetric hydrogen bonds.
The difference between the vibrational energies of and as a function of the coordinate is weakly dependent on volume. This allows us to interpolate in a wide volume range and estimate the pressure at which the proton wave-function shifts away from the centered position. Our calculations show that this symmetry breaking occurs at 103 GPa in HS and 115 GPa in DS (see Fig. 3). The higher transition pressure in DS is due to weaker quantum effects. This isotope effect is similar to the one observed in the ice VII/ice X transition Song et al. (2003). Considering that below 112 GPa the phase is expected to transform into a very different phase consisting of isolated HS and H molecules with HS stoichiometry Li et al. (), -HS might not be formed. However, DS may adopt the structure at pressures below the transition to the phase.
The quantum proton symmetrization has an enormous impact on the phonon spectra of HS. As mentioned earlier, and shown in Fig. 4, the phonon spectra of -HS have several imaginary modes in the harmonic approximation below 175 GPa. This is analogous to ice X, which has only real positive phonon frequencies once the classical limit predicts symmetrization of the hydrogen bond Caracas (2008); Marqués et al. (2009); Bronstein et al. (2014). On the contrary, the corresponding anharmonic SSCHA phonon spectra show well-behaved phonon dispersion relations with positive frequencies in the pressure range of interest (Fig. 4). The anharmonic renormalization of the phonon energies is huge, especially for the H-S bond-stretching modes in which H atoms move towards the neighboring S atoms, which are precisely those modes which drive the second-order phase transition between the and phases. Therefore, the proximity to the second-order quantum phase transition is the origin of the strong anharmonicity.
While the bond symmetrization in ice X occurs in an insulating system, HS is metallic and the symmetrization strongly affects the superconductivity. Indeed, the calculation of the electron-phonon coupling and the superconducting lend further support to the suggestion that -HS yields the record . We use Wannier interpolated electron-phonon matrix elements in our calculations Calandra et al. (2010) and estimate solving the isotropic Migdal-Eliashberg equations. The phonon frequencies and polarizations that enter the electron-phonon calculations are calculated using the SSCHA. The results obtained for the structure using anharmonic phonon frequencies agree well with experimental measurements of for HS and DS and correctly capture the observed increase in with decreasing pressure. We also find an isotope coefficient for HD substitution of at 210 GPa and at 155 GPa in good agreement with experiment (see Fig. 4(c)). The electron-phonon coupling constant , which scales with the phonon frequencies as , is enhanced with decreasing pressure due to the overall softening of the phonon modes. This explains the smooth decrease of with increasing pressure. Between approximately 130 and 150 GPa the increase in is compensated by the decrease in the average phonon frequency and saturates. For the structure we also present SSCHA calculations keeping the centroids at the position. We find a rapid drop of with decreasing pressure as in previous harmonic calculations Akashi et al. (2015), in stark contrast to the experimental data. Therefore, the observed high- superconductivity cannot be explained by HS in the phase. The sudden drop in measured for DS below 150 GPa Drozdov et al. (2015), which is not present in HS, could arise from the need for higher temperatures to anneal the DS sample or could indicate the symmetry breaking that we predicted at 115 GPa. Indeed, the predicted transition pressure depends on the choice of the exchange correlation functional. Even if our choice of the PBE exchange-correlation functional Perdew et al. (1996) appears appropriate, based on agreement between the experimentally observed equation of state Einaga et al. (2015) and DFT calculations, we cannot exclude a small error on the transition pressure.
Supercell calculations for the SSCHA Errea et al. (2013, 2014) and linear response calculations Baroni et al. (2001) were performed within DFT and the generalized gradient approximation functional Perdew et al. (1996) as implemented in the Quantum ESPRESSO Giannozzi et al. (2009) code. We used ultrasoft pseudopotentials Vanderbilt (1990), a plane-wave cutoff energy of 60 Ry for the kinetic energy and 600 Ry for the charge density. The charge density and dynamical matrices were calculated using a 32 Monkhorst-Pack shifted electron-momentum grid for the unit cell calculations. This mesh was adjusted accordingly in the supercell calculations. The electron-phonon coupling was calculated by using electron and phonon momentum grids composed of up to randomly displaced points in the Brillouin zone. The isotropic Migdal-Eliashberg equations were solved using 512 Matsubara frequencies and .
The SSCHA calculations were performed using a 333 supercell for both HS and DS in the phase, yielding dynamical matrices on a commensurate 333 -point grid. The difference between the harmonic and anharmonic dynamical matrices in the 333 phonon momentum grid was interpolated to a 666 grid. Adding the harmonic matrices to the result, the anharmonic dynamical matrices were obtained on a grid. These dynamical matrices were used for the anharmonic electron-phonon coupling calculation. The SSCHA calculations for were performed with a 222 supercell. For consistency, the vibrational energies presented in Fig. 2 were also calculated using a 222 supercell. The electron-phonon calculations for were, however, performed with the SSCHA dynamical matrices interpolated to a 666 grid from the 222 mesh.
The curves in Fig. 2 were obtained as follows. was calculated for and with the SSCHA. With the SSCHA calculation at , we extracted with no further computational effort. Considering that the derivative of the curve at vanishes by symmetry, we can get straightforwardly a potential of the form . The curves presented in Fig. 2 were obtained in this way. The extra point obtained at for HS at (see Fig. 2(a)) confirmed the validity of the fitting procedure. The BOES energies were calculated for many points yielding an accurate fitting curve. Fig. 3 was obtained using a polynomial interpolation of the BOES in the volume range shown and adding the curves that are practically independent of volume.
We acknowledge financial support from the Spanish Ministry of Economy and Competitiveness (FIS2013- 48286-C2-2-P), French Agence Nationale de la Recherche (Grant No. ANR-13-IS10-0003-01), EPSRC (UK) (Grant No. EP/J017639/1), Cambridge Commonwealth Trust, National Natural Science Foundation of China (Grants No. 11204111, 11404148, and 11274136), and 2012 Changjiang Scholars Program of China. Computer facilities were provided by the PRACE project AESFT and the Donostia Internatinal Physics Center (DIPC).
- A. F. Goncharov, V. V. Struzhkin, M. S. Somayazulu, R. J. Hemley, and H. K. Mao, Science 273, 218 (1996).
- P. Loubeyre, R. LeToullec, E. Wolanin, M. Hanfland, and D. Hausermann, Nature 397, 503 (1999).
- M. Benoit, D. Marx, and M. Parrinello, Nature 392, 258 (1998).
- A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
- D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian, and T. Cui, Sci. Rep. 4 (2014).
- D. Duan, X. Huang, F. Tian, D. Li, H. Yu, Y. Liu, Y. Ma, B. Liu, and T. Cui, Phys. Rev. B 91, 180502 (2015).
- I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114, 157004 (2015).
- Y. Li, H. Liu, Y. Zhang, Y. Ma, C. J. Pickard, J. Nelson, R. J. Needs, I. Errea, M. Calandra, and F. Mauri, Submitted, arXiv:1508.03900 .
- N. Bernstein, C. S. Hellberg, M. D. Johannes, I. I. Mazin, and M. J. Mehl, Phys. Rev. B 91, 060511 (2015).
- M. Einaga, M. Sakata, T. Ishikawa, K. Shimizu, M. Eremets, A. Drozdov, I. Troyan, N. Hirao, and Y. Ohishi, ArXiv e-prints (2015), arXiv:1509.03156 [cond-mat.supr-con] .
- L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics, 3rd ed., Vol. 5 Statistical Physics (Butterworth-Heinemann, 1980).
- D. A. Papaconstantopoulos, B. M. Klein, M. J. Mehl, and W. E. Pickett, Phys. Rev. B 91, 184511 (2015).
- R. Akashi, M. Kawamura, S. Tsuneyuki, Y. Nomura, and R. Arita, Phys. Rev. B 91, 224513 (2015).
- E. J. Nicol and J. P. Carbotte, Phys. Rev. B 91, 220507 (2015).
- J. A. Flores-Livas, A. Sanna, and E. K. U. Gross, ArXiv e-prints (2015), arXiv:1501.06336 [cond-mat.supr-con] .
- J. Hirsch and F. Marsiglio, Physica C: Superconductivity and its Applications 511, 45 (2015).
- C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev. Lett. 69, 462 (1992).
- C. Lee, D. Vanderbilt, K. Laasonen, R. Car, and M. Parrinello, Phys. Rev. B 47, 4863 (1993).
- M. I. McMahon, R. J. Nelmes, W. F. Kuhst, R. Dorwarth, R. O. Piltz, and Z. Tun, Nature 348, 317 (1990).
- I. Errea, M. Calandra, and F. Mauri, Phys. Rev. Lett. 111, 177002 (2013).
- I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B 89, 064302 (2014).
- M. Song, H. Yamawaki, H. Fujihisa, M. Sakashita, and K. Aoki, Phys. Rev. B 68, 014106 (2003).
- R. Caracas, Phys. Rev. Lett. 101, 085502 (2008).
- M. Marqués, G. J. Ackland, and J. S. Loveday, High Pressure Research 29, 208 (2009).
- Y. Bronstein, P. Depondt, F. Finocchi, and A. M. Saitta, Phys. Rev. B 89, 214101 (2014).
- M. Calandra, G. Profeta, and F. Mauri, Phys. Rev. B 82, 165111 (2010).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
- P. Giannozzi et al., J. Phys. Condens. Matter 21, 395502 (2009).
- D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).