Superconductivity of Light-Elements Doped H$ {}_{3} $S

Pressurized hydrogen-rich compounds, which could be viewed as precompressed metallic hydrogen, exhibit high superconductivity, thereby providing a viable route toward the discovery of high-temperature superconductors. Of particular interest is to search for high-temperature superconductors with low stable pressure in terms of pressure-stabilized hydrides. In this work, with the aim of obtaining high-temperature superconducting compounds at low pressure, we attempt to study the doping effects for high-temperature superconductive $ \mathrm{H_3S} $ with supercells up to 64 atoms using first principle electronic structure simulations. As a result of various doping, we found that Na doping for $ \mathrm{H_3S} $ could lower the dynamically stable pressure by 40 GPa. The results also indicate P doping could enhance the superconductivity of $ \mathrm{H_3S} $ system, which is in agreement with previous calculations. Moreover, our work proposed an approach that could reasonably estimate the superconducting critical temperature ($ T_{c} $) of a compound containing a large number of atoms, saving the computational cost significantly for large-scale elements-doping superconductivity simulations.

Structure searches below 200 GPa were performed for the C-S-H system [19-21, 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 Im3m H 3 S [21,26], 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 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 * hanyuliu@jlu.edu.cn elements in Im3m H 3 S with supercell approach at the pressure range of 150-250 GPa by first-principle calculations. Our calculations mainly focus on H 24 S 7 X and H 48 S 15 X, where 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 H 3 S system. Furthermore, the P doping could enhance the superconductivity of the parent H 3 S system, which is ascribed to octahedra units [SH 6 ] and [PH 6 ]. 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.

II. COMPUTATIONAL DETAILS
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 projectoraugmented wave method [30]. Monkhorst k meshes [31] spacing 2π×0.1Å −1 was used to sample the first Brillouin Zones. The electronic density of states was also computed by VASP with 20×20×20 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 k mesh of 16×16×16 and q mesh of 4×4×4 and tested the convergence with k mesh of 20 × 20 × 20 and q mesh of 5 × 5 × 5 for H 24 S 7 X. For H 48 S 15 X, k mesh was used as 12 × 12 × 12 and q mesh was used as 3 × 3 × 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 where are the correction factors. µ * , λ and ω 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] to compare with that of MAD equation, which is realized by the Elk code [40]. T c could be obtained once the superconducting gap ∆(iω j ) becomes zero in numerically solving ME equation.

III. RESULTS AND DISCUSSIONS
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 (H 24 S 8 ) for 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 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 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 H 24 S 7 X and 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 H 24 S 7 X compounds, for example, imaginary frequency was found for Γ point with X=H, indicating dynamical instability. The dynamical stability of 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 H 3 S has a lower dynamical stable pressure (140 GPa) compared to 180 GPa of the parent H 3 S [2]. The absence of imaginary phonon frequency in the simulated phonon dispersion implies the dynamical stability of 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 λ of 2.31. As a result, T c of H 24 S 7 Na can reach 191 K at 150 GPa by using ME equation with µ * = 0.10.
As shown in previous studies [42,43], the P doped H 3 S has the high superconductivity due to the enhanced density of states at the Fermi level of parent H 3 S. We have also performed the simulations for elucidating the physical mechanism of this compound by using a supercell of 64 atoms (H 48 S 15 P). As is shown in Fig. 3(a) and Fig. S4, there are several flat bands along the k path F → Q → Z at the Fermi surface, with the derivative ∂E n /∂k almost zero. This indicates that the doping of P alters the two van Hove singularities [42,44] 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)   To explore the superconductivity of the C-doped H 3 S, we have computed the electron-phonon coupling strength of H 32 S 7 C, H 36 S 11 C and H 48 S 15 C. For larger supercells corresponding to lower-proportion C-doped 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 H 3 S without performing actual electron-phonon coupling sim- ulations for a large supercell. Given that similar structures of the low-proportion C-doped H 3 S shares similar magnitude of average mass and electron-phonon interaction at the same pressure, we could thus estimate λ of the MgB 2 type superconductors by the Hopfield expression [45,46] Then, the T c can be estimated by [46] T c = ω exp − 1 λ 1+λ − µ * 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 3 S up to 256 atoms as shown in Table S3 and Fig. 4. The T c of H 3 S 1−x C x with x=0.1250-0.0156 could reach 148-192 K with µ * =0.10 at 260 GPa. This suggests that the estimated superconductivity of C-doped H 3 S is much lower than room temperature as predicted in previous studies [22,23]. 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 I.

IV. CONCLUSION
In summary, we have computed the superconductivity of light-elements doped 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 H 3 S while the doping of P can increase the density of states at the Fermi level as well as the superconductivity of H 3 S. Remarkably, we found that the existence of octahedra [PH 6 ] and [SH 6 ] units with squeezed P-H and S-H bonds in 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 H 3 S. Furthermore, we have proposed an estimation approach to reasonably estimate the superconductivity of the low proportion C-doping 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.

A. Mcmillan-Allen-Dynes Equation
The Mcmillan-Allen-Dynes equation can be expressed as [1] T λ and ω log are isotropic electron-phonon coupling parameter and logarithmically averaged phonon frequency, respectively. λ depends on a critical quantity Where Ω is the volume of the first Brillouin zone, and g ijν (k, q) is the electron phonon matrix element [1,2]. To analyze the impact of each quantity on T c , we first rewrite Eq. S1 in a simplified form notice correction factor f = f 1 f 2 ∼ 1 Since ω log is nearly irrelevant to the magnitude of α 2 F , thus the change of magnitude of α 2 F only affect λ. Thus setting the screened Coulomb parameter µ * fixed, we have Furthermore, the larger λ enables f larger as well. Thus the increase of magnitude of α 2 F will enhance T c . Note that Thus regardless of anisotropy, λ ∼ |g| 2 N (0)N ph (ω). By further inspecting Eq. S4, we found that if at specific area (in the direct product of two reciprocal spaces V k ⊗ V q ) the electron density of states, phonon density of states and electron phonon coupling are simultaneously large, the λ can be greatly enhanced while keeping ω log at large values. If we only focus on q, we can integrate by k and define We can obtain below expression We found that λ ∝ N ph γ.

B. Estimation Approach of C-Doped H 3 S and Its Validity
For large supercells (H 3 S 1−x C x ), we computed the density of states at Fermi surface N (0) x , then given that T c have a roughly "ω 2 free" expression [1,3] T c = 0.18 √ λω 2 = 0.18 N (0)I 2 /M (S11) in Allen-Dynes limit and we assume similar structures share similar scattering matrix element I 2 and similar average inonic mass, we can define a pseudo "electron-phonon coupling constant" λ x by the λ from a known structure: λ x ≡ λN (0) x /N (0). Here N (0) is the density of states of known structure at the Fermi surface. For the reason of self-consistency, we use the density of states calculated without considering electron-phonon interpolation to estimate all the N (0) and N (0) x . Then with T c of the known structure, the transition temperature T cx is estimated to be We estimate the T cx of H 3 S 1−x C x by both the 64 atoms supercell H 3 S 0.9375 C 0.0625 and H 3 S, to obtain the values T cx,0.0625 and T cx,0 resprectively and estimate our final result as To illustrate the validity of this approach, we estimate the T c of H 3 S 0.9375 X 0.0625 by H 3 S 0.875 X 0.125 and H 3 S, in this scenario, T cx = x 0.125 T cx,0.125 + 1 − x 0.125 T cx,0 . All the results from the estimated approach and directly computed by Quantum Espresso are shown in Table S2, where their difference is only around 5% at maximum.