High temperature superconductivity in sulfur and selenium hydrides at high pressure

Due to its low atomic mass hydrogen is the most promising element to search for high-temperature phononic superconductors. However, metallic phases of hydrogen are only expected at extreme pressures (400 GPa or higher). The measurement of a record superconducting critical temperature of 190 K in a hydrogen-sulfur compound at 200 GPa of pressure[1], shows that metallization of hydrogen can be reached at significantly lower pressure by inserting it in the matrix of other elements. In this work we re-investigate the phase diagram and the superconducting properties of the H-S system by means of minima hopping method for structure prediction and Density Functional theory for superconductors. We also show that Se-H has a similar phase diagram as its sulfur counterpart as well as high superconducting critical temperature. We predict SeH3 to exceed 120 K superconductivity at 100 GPa. We show that both SeH3 and SH3, due to the critical temperature and peculiar electronic structure, present rather unusual superconducting properties.


PACS numbers:
Under high pressure conditions, insulating and semiconducting materials tend to become metallic, because, with increasing electronic density, the kinetic energy grows faster than the potential energy. As metallicity is a necessary condition for superconductivity, superconductivity becomes more likely under pressure [2,3]. Wigner and Huntington [4], already in 1935 suggested the possibility of a metallic modification of hydrogen under very high pressures. Ashcroft and Richardson predicted [5,6] hydrogen to become metallic under pressure and also the possibility to be a high temperature superconductor. The high critical temperature (T C ) of hydrogen [7][8][9] is a consequence of its low atomic mass leading to high energy vibrational modes and in turn to a large phase space available for electron-phonon scattering to induce superconductivity [10]. However, the estimated pressure of metallization [11,12] is beyond the current experimental capabilities and it has been a challenge to confirm this hypothesis [13][14][15][16].
It was only recently that hydrogen-rich compounds started to be explored as a way to decrease the tremendous pressure of metallization in pure hydrogen [17], essentially performing a chemical precompression. The first system explored experimentally was silane (SiH 4 ) [18]. Soon after, many other precompressed hydrogen rich materials have been explored experimentally [19][20][21][22] and theoretically [23][24][25][26][27][28][29][30][31][32][33][34][35]. The importance of a systematic search for a crystalline ground state has been put in evidence for disilane (SiH 6 ), where structures enthalpically higher lead to transition temperatures of the order of 130 K. Interesting structures have been proved not to be the global minimum and for the correct ground state was found a rather moderate electron-phonon coupling and T C of 25 K [36]. In agreement with experimental evidence.
Recently it was reported that sulfur hydride (SH 2 ), when pressurized, becomes metallic and superconduct-ing. For pressures above 180 GPa an extremely high transition temperature of about 190 K was measured [1]. This T C is higher than in other superconductors known so far, including cuprates and pnictides. The experimental evidence is supported theoretically [32,35,37], and crystal prediction methods suggest that the system becomes superconducting with a SH 3 stoichiometry. In this work we re-investigate extensively the S-H phases with state of the art ab-initio material search minima hopping methods [38][39][40] (MHM) and compute the superconducting properties with the completely parameter-free Density Functional Theory for Superconductors (SCDFT). We also extend the analysis to the Se-H system, predicting a fairly similar phase diagram and comparable superconducting properties.

METHODS
Electronic and phononic structure calculations are performed within density-functional theory as implemented in the two plane-wave based codes abinit [41], and espresso [42] within the local density approximation LDA exchange correlation functional. The core states were accounted for by norm-conserving Troullier-Martins pseudopotentials [43]. The pseudopotential accuracy has been checked against all-electron (LAPW+lo) method as implemented in the elk code (http://elk.sourceforge.net/). In order to predict the ground state structure of sulfur/selenium hydride compounds we use the minima hopping method [38][39][40] for the prediction of low-enthalpy structures. This method has been successfully used for global geometry optimization in a large variety of applications [44][45][46]. Given only the chemical composition of a system, the MHM aims at finding the global minimum on the enthalpy surface while gradually exploring low-lying structures. Moves on the enthalpy surface are performed by using variable cell shape molecular dynamics with initial velocities approximately chosen along soft mode directions. We have used 1,2,3 formula units of H 3 S and H 3 Se at selected pressures of 100, 150, 200, 250 and 300 GPa. The relaxations to local minima are performed by the fast inertia relaxation engine [47] and both atomic and cell degrees of freedom are taking into account. Final structural relaxations and enthalpy calculations were performed with the vasp code [48]. The plane-wave cutoff energy was set to 800 eV, and Monkhorst-Pack k-point meshes with grid spacing denser than 2π × 0.01 Å −1 , resulting in total energy convergence to better than 1 meV/atom. Superconducting properties have been computed within density-functional theory for superconductors (SCDFT) [49][50][51]. This theory of superconductivity is completely ab-initio, fully parameter-free and proved to be rather accurate and successful in describing phononic superconductors [52][53][54][55]. It allows to compute all superconducting properties including the critical temperature and the excitation spectrum of the system [56].

CRYSTAL STRUCTURE PREDICTION AND ENTHALPIES
Experimentally little is known on the high pressure stability and composition of the S-H system and, to the best of our knowledge, nothing is known about the Se-H. Therefore we investigate their low temperature phase diagram by means of the MHM for the prediction of lowenthalpy structures. Computed enthalpies as a function of pressure are reported in Fig. 1. We consider the H 3 S stoichiometry as well as its elemental decomposition (sul-fur + hydrogen), its decomposition into H 2 S + hydrogen and H 4 S -hydrogen [57]. At low pressure we find the Cccm structure up to 95 GPa and the R3m (β-Po-type) rhombohedral structure between 95 and 150 GPa. Above 150 GPa, we confirm [32,37] the cubic Im3m (bcc) as the most stable lattice.
In a similar way we have studied the Se-H phase diagram. Chemically, selenium is known to have very similar physical properties to sulfur and this system is not an exception. The enthalpies of the phases found in our MHM runs are shown in Fig. 1b. Once again we use the Cccm structure as reference since, as in the case of sulfur, it is the most stable at low pressures and up to 70 GPa. Between 70 GPa to 100 GPa, we find that the H 2 Se + hydrogen decomposition is more stable than the H 3 Se stoichiometry. H 3 Se returns to be the most stable composition above 100 GPa and at least up to 250 GPa.
Therefore from our analysis both systems in the range 50 GPa to 250 GPa show, with increasing pressure, two phase transitions. The S-H system, always stable in the H 3 S stoichiometry, has a first order phase transition from Cccm to R3m at ∼100 GPa, then the R3m rhombohedral distortion decreases continuously up to 150 GPa where it transforms into the Im3m cubic structure. The Se-H system at low pressure is also stable in the H 3 Se stoichiometry but becomes unstable to a phase separation into H 2 S + hydrogen in the range from 70 GPa to 100 GPa. Above 100 GPa another discontinuous phase transition occurs, directly into the Im3m cubic structure. Note that 100 GPa is also the pressure below which the Im3m structure would distort into the β-Po R3m, therefore depending on experimental conditions this rhombohedral phase may occur as a metastable one. The sequence of transformation is highlighted in Fig. 1 by means of shaded areas.

ELECTRONIC AND PHONONIC PROPERTIES OF H3SE AND H3S AT HIGH PRESSURE
We focus now on the properties of H 3 S and H 3 Se in the pressure range of stability of the Im3m structure. The two materials present very similar properties. At 200 GPa electronic band dispersions and Fermi surfaces are barely distinguishable, as seen in Fig. 2. And in the range of pressure between 100 to 200 GPa there are no significant changes in the electronic properties apart from the overall bandwidth that increases with pressure. An important aspect of the electronic structure is the presence of several Fermi surface sheets, with no marked nesting features and with Fermi states both at low and high momentum vector. At small momentum (close to the Γpoint, center of the Brillouin zone in Fig. 2) there are three small Fermi surfaces (only the green larger one can be seen in the figure, smaller ones being inside it). However, these provide a rather small contribution to the total density of states (DOS) at the Fermi level which mostly comes from the two larger Fermi surface sheets. These are of hybrid character, meaning that their Kohn-Sham (KS) states overlap both with H and S/Se states (the overlap is expressed in the figure by the color-scale of the band lines), suggesting that they will be coupled with both hydrogen and S/Se lattice vibrations (more details on this point will be given in the next section). Overall the DOS shows a square root behavior of the 3D electron gas, the main deviation from this occurs close to the Fermi energy where a peak with an energy width of about 2 eV is present. This structure will play a relevant role in the superconducting properties.
Unlike the electronic structure, phonons are strongly pressure and material dependent. Clearly a key role is played by the occurrence of the II order R3m to Im3m phase transition. Far away from it (i.e. at very high pressure) we have three sets of well separated phonon modes: acoustic (below 60 meV), optical modes that are transverse with respect to the S/Se-H bond (between 100 to 200 meV) and, above 200 meV, stretching modes of the Se/S-H bond. These are clearly seen in Fig. 3b for H 3 Se at 200 GPa. As pressure reduces, the bond structure of the system tends to destabilize because, from a four-fold coordination in the Im3m structure, it goes to a threefold coordination in the R3m one. This means that one of the high-energy stretching mode slowly softens at Γ. This softening can be clearly seen in H 3 S at 200 GPa (Fig. 3a) where a H-S stretching mode went down to about 60 meV. Eventually, as pressure lowers this softens to zero energy, marking the occurrence of the phase transition, at about 150 GPa in H 3 S and slightly below 100 GPa in H 3 Se. In fact, at 100 GPa this mode has, in H 3 Se, almost zero energy (see Fig. 3c).
In spite of these important changes in the phononic energy dispersion, the overall coupling strength [58,59] λ does not show large variations over the pressure range, as we can see from Tab. I. Naturally the coupling increases near the phase transition due the optical mode softening, however, as this is restricted to a relatively small region near the Γ point, the effect is not dramatic. On the other hand there is definitively a difference in the coupling strength of the Se (λ ∼ 1.5) with respect to the S system (λ ∼ 2.5), indicating that selenium, due to its larger ionic size, provides a better electronic screening of the hydrogen vibrations.

SUPERCONDUCTING PROPERTIES
We have computed, by means of SCDFT, critical temperatures of H 3 S and H 3 Se in the pressure range of stability of the Im3m structure, these are collected in Tab. I. The predicted T C for the H 3 S system is 180 K at 200 GPa and 195 K at 180 GPa, in agreement with the measured value of 185 K at 177 GPa. On the other hand our predic-tion for the deuterium substituted system D 3 S is 141 K, at 200 GPa. That is much larger than the measured[1] T C of 90 K. This huge experimental isotope effect is therefore not consistent with our calculations. However the good agreement obtained with the T C of the H 3 S system seem to exclude an explanation in terms of anharmonic effects in the hydrogen vibrations, as suggested in Ref. 60. Nevertheless, the theoretical isotope coefficients α S = 0.05 and α H = 0.4 (defined as α A = − M A T c ∂T C ∂M A , and computed at 200 GPa with a three point numerical differentiation) clearly indicate and confirm [60] the dominant contribution of hydrogen phonon modes to the superconducting pairing.
Our prediction for H 3 Se at 200 GPa is of 131 K, this reduction of T C is clearly not an isotopic effect of the substitution S to Se. As mentioned in the previous section, it is caused, instead, by a different coupling strength of the hydrogen modes in the Se environment. In spite of the much lower coupling strength λ the reduction of T C is not very large with respect to the sulfur system, as expected from the fact that the critical temperature at high coupling increases with the square root of λ (while is exponential at low coupling) [58,59].
To compute the critical temperatures we have used SCDFT, this choice allow us to deal with the unusual superconducting properties of these systems from first principles without relying on conventional assumptions coming from low pressure experience. There are two aspects of these systems that are uncommon, that make the use of conventional [61] Eliashberg methods difficult to apply to these systems. First, the strong variation of the electronic density of states at the Fermi level, that is pinned to a rather sharp peak in the DOS, second the extremely large el-ph coupling and phonon frequencies that lead to a very broad region around E F where the interaction is dominated by phonons over Coulomb repulsion.
The effect of the energy dependence of the DOS can be appreciated by comparing Eliashberg results with SCDFT when neglecting the Coulomb interaction (see Tab. I). At 200 GPa the two theories [65] disagree by 54 K, SCDFT giving 284 K while Eliashberg gives 338 K. This difference comes from the energy dependence of the DOS, while Eliashberg assumes a flat DOS, in the SCDFT we can easily check this assumption by assuming a flat DOS, and for this case the SCDFT calculation would lead to 334 K, in agreement with the Eliashberg result.
Physically the reduction of T C , occurring when the real DOS is considered, arises from the fact that the phononic pairing extends in a rather large region around the Fermi level, over the DOS peak structure of these systems (see Fig. 2). Beyond the range of the phononic pairing the coupling is dominated by the Coulomb interaction. As, in the static limit, this is repulsive, a superconducting system compensates it by a phase shift in the gap (i.e. in the quasiparticle orbitals), therefore making this repul- sion contribute to the condensation (in unconventional superconductors exactly the same happens but directly at the Fermi level). This mechanism is called Coulomb renormalization [59] since it renormalizes the repulsive Coulomb scattering that occurs at low energy. The phase shift occurs at | k | ω log but the scattering processes become less and less important as | k | increases (going down as 1/ ). Therefore the most important energy region is where the DOS of the H 3 S and H 3 Se systems shows a dip, implying that the phase space available for this process is small and its effect weak. Note that in order to reproduce the T C coming from SCDFT within the Allen-Dynes (AD) formula one should assume a µ * of 0.16, that is actually much larger than the value of µ itself ( 0.1). Making clear that the the Morel-Anderson theory [66] can not even be applied.
The superconducting pairing is distributed over many phonon modes and over the Brillouin zone in q-space, despite the presence of several Fermi surface sheets and with different orbital character across the Fermi level, we obtained a isotropic (weakly k-dependent) gap at the Fermi level and the effect of anisotropy [67] on T C is negligible (< 1 K).

CONCLUSIONS
We have presented a theoretical investigation on the crystal structure and superconductivity of H 3 S. An extensive structural search confirms the H 3 S stoichiometry as the most stable configuration at high pressure. By   Table I. Calculated critical Temperatures and gaps. λ is the electron phonon coupling parameter [58,59]; ω log is the logarithmic average of the α 2 F function [58,59]; T SCDFT c is the critical temperature from SCDFT including RPA screened Coulomb repulsion; T SCDFT,ph c is the phonon only SCDFT critical temperature; T AD,µ * =0 c is the critical temperature from the Allen-Dynes modified McMillan formula [59,62,63], at µ * = 0 (i.e. with no Coulomb pairing); T AD,µ * =0.1 c is the same AD formula but with the conventional value of µ * = 0.1 ; ∆(T=0) is the superconducting gap computed [64] from the SCDFT calculations means of parameter-free SCDFT we have predicted a T C of 180 K at 200 GPa, in excellent agreement with experimental results. This confirms H 3 S to be the material with the highest known superconducting critical temperature. The mechanism of superconductivity is clearly the same that was predicted for metallic hydrogen [5,7,9]: the combinate effect of high characteristic frequency due to hydrogen light mass and strong coupling due to the lack of electronic core in hydrogen. Still the working pressures of this superconductor is too high for any technological application [68]. Nevertheless the discovery of metallic superconducting hydrogenic bands already at 150 GPa gives hope that further theoretical and experimental research in this direction may lead to even lower hydrogen metallization pressures and higher temperature superconductivity. Here we predict that H 3 S is stable in the cubic Im3m already at 100 GPa with a very high T C of 123 K, a value which is comparable to the cuprate superconductors. J.A.F.L. acknowledge fruitful discussion with Maximilian Amsler on crystal prediction and the financial support from EU's 7th framework Marie-Skłodowska-Curie scholarship program within the "ExMaMa" Project (329386)