Superconductivity of Light-Elements Doped H₃S

The search for the high-temperature superconducting hydrides at high pressures has attracted attention in condensed matter physics field. In this regard, many hydrides with relatively high $T_{c}$ were identified under high pressure [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Among these high-temperature superconducting hydrides, the theoretically predicted $\mathrm{H_{3}S}$ with $T_{c}$ of 203 K [1, 2, 3] and $\mathrm{LaH_{10}}$ with $T_{c}$ of 250-260 K in [5, 6, 8, 7], as well as CaH₆ being the first example of clathrate hydrides ever predicted [16, 17, 18], were experimentally synthesized. Recently, the carbonaceous sulfur hydride was found to possess extremely high superconductivity with a $T_{c}$ as high as 288 K at 267 GPa [12], however, the actual crystal structure and the mechanisms of such extremely high superconductivity remains unclear and missing [19, 20, 21, 22, 23].

Structure searches below 200 GPa were performed for the C-S-H system [21, 19, 20, 24, 25], while no high superconductive structure is yet identified. Later, it was reported that the high superconductivity in the C-S-H compound could be explained by the doping of C into the $Im\bar{3}m$ $\mathrm{H_{3}S}$ [26, 21], whereas these simulations are based on the virtual crystal approximation (VCA) [27]. This approximation is employed with linearly mixed pseudopotentials and couldnot take the details of symmetry breaking and local distortions into account. These issues, in principle, could be addressed by performing the electronic simulations of doping carbon into a sufficiently large supercell of $\mathrm{H_{3}S}$, if the computational power is allowed.

In this paper, we investigated the structures and physical properties of partially substituting sulfur by the light elements in $Im\bar{3}m$ $\mathrm{H_{3}S}$ with supercell approach at the pressure range of 150-250 GPa by first-principle calculations. Our calculations mainly focus on $\mathrm{H_{24}S_{7}X}$ and $\mathrm{H_{48}S_{15}X}$, where $\mathrm{X}$ denotes the doping elements from H to Cl without He and Ne in the periodic table. As a result, we found that the Na doing could lower the dynamically stable pressure of 40 GPa compared to the parent $\mathrm{H_{3}S}$ system. Furthermore, the P doping could enhance the superconductivity of the parent $\mathrm{H_{3}S}$ system, which is ascribed to octahedra units [SH₆] and [PH₆]. In addition, we proposed an estimation approach to investigate the low proportion C-doping effects at 260 GPa. The results suggest the estimated $T_{c}$ is much lower than the room temperature.

The structural optimization was done by the Vienna Ab initio Simulation Package (VASP) [28], with pseudopotentials employing generalized gradient approximation (GGA) based Perdew-Burke-Ernzerhof (PBE) type exchange correlation functional [29] and projector-augmented wave method [30]. Monkhorst $\bm{k}$ meshes [31] spacing $2\pi\times 0.1$ Å^-1 was used to sample the first Brillouin Zones. The electronic density of states was also computed by VASP with $20\times 20\times 20$ $\bm{k}$ mesh and was analyzed by VASPKIT [32]. The phonon properties and superconducting properties were computed by the Quantum-Espresso (QE) package, with vanderbilt ultra-soft pseudopotentials [33]. We have adopted $\bm{k}$ mesh of $16\times 16\times 16$ and $\bm{q}$ mesh of $4\times 4\times 4$ and tested the convergence with $\bm{k}$ mesh of $20\times 20\times 20$ and $\bm{q}$ mesh of $5\times 5\times 5$ for $\mathrm{H_{24}S_{7}X}$. For $\mathrm{H_{48}S_{15}X}$, $\bm{k}$ mesh was used as $12\times 12\times 12$ and $\bm{q}$ mesh was used as $3\times 3\times 3$. The smearing method was Methfessel-Paxton first-order spreading [34] of 0.03 Ry. The cutoff energy for basis of plane waves was employed to be 100 Ry. Then, the transition temperatures were estimated by McMillan-Allen-Dynes (MAD) formula [35]

$$T_{c}=f_{1}f_{2}\dfrac{\omega_{\mathrm{log}}}{1.20}\mathrm{exp}\left[-\dfrac{1.04(1+\lambda)}{\lambda-\mu^{*}(1+0.62\lambda)}\right]$$

(1)

where

	$\displaystyle f_{1}$	$\displaystyle=\left\{1+\left[\dfrac{\lambda}{2.46(1+3.8\mu^{*})}\right]^{3/2}\right\}^{1/3}$		(2)
	$\displaystyle f_{2}$	$\displaystyle=1+\dfrac{\lambda^{2}(\omega_{2}/\omega_{\text{log}}-1)}{\lambda^{2}+[1.82(1+6.3\mu^{*})(\omega_{2}/\omega_{\text{log}})]^{2}}$

are the correction factors. $\mu^{*}$, $\lambda$ and $\omega_{\text{log}}$ indicate the screened Coulomb parameter, electron-phonon coupling constant and the logarithm average over phonon frequency, respectively.

We have also computed the results by Migdal-Eliashberg (ME) theory [36, 37, 38, 39]

	$$\displaystyle Z(i\omega_{j})=1+\dfrac{\pi T}{\omega_{j}}\sum\limits_{j^{\prime}}\dfrac{\omega_{j^{\prime}}}{\sqrt{\omega_{j^{\prime}}^{2}+\Delta^{2}(i\omega_{j^{\prime}})}}\lambda(i\omega_{j}-i\omega_{j^{\prime}})$$		(3)
	$$\displaystyle\begin{aligned} Z(i\omega_{j})\Delta(i\omega_{j})&=\pi T\sum\limits_{j^{\prime}}\dfrac{\Delta(i\omega_{j^{\prime}})}{\sqrt{\omega_{j^{\prime}}^{2}+\Delta^{2}(i\omega_{j^{\prime}})}}\\ &\times[\lambda(i\omega_{j}-i\omega_{j^{\prime}})-\mu^{*}]\end{aligned}$$		(4)
	$$\displaystyle\lambda(i\omega_{j}-i\omega_{j^{\prime}})=\int\mathrm{d}\omega\dfrac{2\omega\alpha^{2}F(\omega)}{\omega^{2}+(\omega_{j}-\omega_{j^{\prime}})^{2}}$$		(5)

to compare with that of MAD equation, which is realized by the Elk code [40]. $T_{c}$ could be obtained once the superconducting gap $\Delta(i\omega_{j})$ becomes zero in numerically solving ME equation.

We began our simulations on investigating the validity of the supercell approach by computing the electronic properties and phonon properties of primitive cell and a supercell of 32 atoms ($\mathrm{H_{24}S_{8}}$) for $\mathrm{H_{3}S}$ at 200 GPa, as shown in Fig. S1 [41]. As shown in Fig. S2 [41], the superconductivity using the supercell is well consistent with that simulated from primitive cell of $\mathrm{H_{3}S}$, which is also in agreement with previous results [2, 3]. The relevant information is listed in Table S1 [41].

Furthermore, we have systematically investigated the doping of $\mathrm{H_{3}S}$ by the elements from H to Cl without He and Ne in the periodic table using supercells of 32 and 64 atoms. The structures of $\mathrm{H_{24}S_{7}X}$ and $\mathrm{H_{48}S_{15}X}$ are provided in Fig. S3 [41]. The results indicate the doping could destabilize the structure for several compounds. For the $\mathrm{H_{24}S_{7}X}$ compounds, for example, imaginary frequency was found for $\Gamma$ point with X=H, indicating dynamical instability. The dynamical stability of $\mathrm{H_{24}S_{7}X}$ compounds within the range of 150-250 GPa is summarized in Fig. 1. Moreover, we found that 12.5% doping of Na into $\mathrm{H_{3}S}$ has a lower dynamical stable pressure (140 GPa) compared to 180 GPa of the parent $\mathrm{H_{3}S}$ [2]. The absence of imaginary phonon frequency in the simulated phonon dispersion implies the dynamical stability of $\mathrm{H_{24}S_{7}Na}$, as shown in Fig. 2. It is clearly seen that the strongest electron-phonon interaction mainly emerges around $P$ point, which may lead to a large $\lambda$ of 2.31. As a result, $T_{c}$ of $\mathrm{H_{24}S_{7}Na}$ can reach 191 K at 150 GPa by using ME equation with $\mu^{*}=0.10$.

As shown in previous studies [42, 43], the P doped $\mathrm{H_{3}S}$ has the high superconductivity due to the enhanced density of states at the Fermi level of parent $\mathrm{H_{3}S}$. We have also performed the simulations for elucidating the physical mechanism of this compound by using a supercell of 64 atoms ($\mathrm{H_{48}S_{15}P}$). As is shown in Fig. 3(a) and Fig. S4, there are several flat bands along the $\bm{k}$ path $F\to Q\to Z$ at the Fermi surface, with the derivative $\partial E_{n}/\partial k$ almost zero. This indicates that the doping of P alters the two van Hove singularities [44, 42] and can result in a peak of density of states right at the Fermi surface. Moreover, we found that the symmetry is preserved well after doping due to the existence of the near degeneracies in electronic band structures as shown in Fig. 3(a), which can also be a critical factor to induce Van Hove singularities [23]. The high electronic density of states at the Fermi level could contribute to the large magnitude of phonon linewidth, as shown in Fig. 3(b) and contrasted with Fig. S5 [41]. We have also found the hybridization of $s$ orbitals of H with $p$ orbitals of S and P near the Fermi surface (Fig. 3(c)). The negligible curvature of the flat bands indicates well localization of corresponding $s$ and $p$ electrons, revealing the strong P-H and S-H chemical bondings. Fig. 3(d) shows the existence of [$\mathrm{PH_{6}}$] and [$\mathrm{SH_{6}}$] units, where the distances between P-H and S-H are compressed to 1.470 Å and 1.474 Å, compared with 1.493 Å in $\mathrm{H_{3}S}$ and 1.494 Å-1.511 Å for other S-H bonds in $\mathrm{H_{48}S_{15}P}$. Finally, we found that the 6.25% P doping could enhance the $T_{c}$ of $\mathrm{H_{3}S}$ around 20-30 K at 200 GPa and 250 GPa.

To explore the superconductivity of the C-doped $\mathrm{H_{3}S}$, we have computed the electron-phonon coupling strength of $\mathrm{H_{32}S_{7}C}$, $\mathrm{H_{36}S_{11}C}$ and $\mathrm{H_{48}S_{15}C}$. For larger supercells corresponding to lower-proportion C-doped $\mathrm{H_{3}S}$, the simulations could not be afforded due to computational demanding of these simulations as well as the limition of our computational power. Therefore, we attempt to estimate the superconductivity of low-proportion C-doped $\mathrm{H_{3}S}$ without performing actual electron-phonon coupling simulations for a large supercell. Given that similar structures of the low-proportion C-doped $\mathrm{H_{3}S}$ shares similar magnitude of average mass and electron-phonon interaction at the same pressure, we could thus estimate $\lambda$ of the $\mathrm{MgB_{2}}$ type superconductors by the Hopfield expression [45, 46]

$$\lambda=\dfrac{N(0)I^{2}}{M\omega^{2}}$$

(6)

Then, the $T_{c}$ can be estimated by [46]

$$T_{c}=\omega\exp\left(-\dfrac{1}{\frac{\lambda}{1+\lambda}-\mu^{*}}\right)$$

(7)

Further details of our estimation approach is also provided in the supplemental material [41]. As is shown in Table S2, values of $T_{c}$ computed from our approach vary only about 5% at maximum compared with that directly computed by QE, indicating the validity of our computational scheme. We thus investigated the superconductivity of the doping system of H₃S up to 256 atoms as shown in Table S3 and Fig. 4. The $T_{c}$ of H₃S_1-xC_x with $x$=0.1250-0.0156 could reach 148-192 K with $\mu^{*}$=0.10 at 260 GPa. This suggests that the estimated superconductivity of C-doped $\mathrm{H_{3}S}$ is much lower than room temperature as predicted in previous studies [23, 22]. The reason for this discrepancy is possibly because different approches were employed. As for other predicted compounds, the superconductivity-related information is included in Table 1.

In summary, we have computed the superconductivity of light-elements doped $\mathrm{H_{3}S}$ using supercell approach within the framework of first-principle electronic structure. Our simulations indicate that the doping of Na can lower the dynamically stable pressure of $\mathrm{H_{3}S}$ while the doping of P can increase the density of states at the Fermi level as well as the superconductivity of $\mathrm{H_{3}S}$. Remarkably, we found that the existence of octahedra [$\mathrm{PH_{6}}$] and [$\mathrm{SH_{6}}$] units with squeezed P-H and S-H bonds in $\mathrm{H_{48}S_{15}P}$ which are likely to be related to the high density of states at the Fermi level, with the higher $T_{c}$ of 20-30 K compared with $\mathrm{H_{3}S}$. Furthermore, we have proposed an estimation approach to reasonably estimate the superconductivity of the low proportion C-doping $\mathrm{H_{3}S}$ without performing electron-phonon calculations on a quite large supercell. Our current work may inspire future work toward searching for high-temperature superconductivity in light-elements doping systems.

This work was supported by the Major Program of the National Natural Science Foundation of China (Grant No. 52090024), National Natural Science Foundation of China (Grant No. 12074138, 11874175, and 11874176), Fundamental Research Funds for the Central Universities (Jilin University, JLU), Program for JLU Science and Technology Innovative Research Team (JLUSTIRT), and The Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB33000000). This work used computing facilities at the High-Performance Computing Centre of Jilin University.