Possible Superconductivity Approaching Ice Point

Recently BCS superconductivity at 203 K has been discovery in a highly compressed hydrogen sulfide. We use first-principles calculations to systematically examine the effects of partially substituting the chalcogen atoms on the superconductivity of hydrogen chalcogenides under high pressures. We find detailed trends of how the critical temperature changes with increasing the V-, VI- or VII-substitution rate, which highlight the key roles played by low atomic mass and by strong covalent metallicity. In particular, a possible record high critical temperature of 280 K is predicted in a stable H3S0.925P0.075 with the Im-3m structure under 250 GPa.


I. INTRODUCTION
Since the discovery of superconductivity in mercury at 4 K in 1911 1 , scientists have been ardently pursuing new and better superconductors at higher temperatures. In a conventional superconductor, the vibrations of crystal lattice provide an attractive force that binds an electron with its time-reversal partner into a Cooper pair 2 . The Cooper pairs can Bose condense below a critical temperature (T c ), which had been believed to be no larger than 25 K prior to the discovery of MgB 2 3,4 . In this case, the Debye temperature could be as large as ∼ 10 3 K, and low T c is primarily related to small electronphonon couplings 5 . Since 1986 and 2008, respectively, the discoveries of copper-and iron-based superconductors have provided two new avenues for making high-T c superconductors [6][7][8][9][10][11][12] , while generating new excitements in fundamental physics. Although their unconventional mechanisms are still under hot debate, to date a record T c of 164 K has been experimentally realized in the cuprate family 9 and 56 K among the iron-pnictide compounds 12 .
One may thus wonder whether a T c higher than 164 K could be achieved and whether the conventional phononmediated mechanism could play a significant role in such a race. Both answers have become extremely positive and fascinating, given the recent theoretical prediction 13,14 and particularly the most recent experimental observation 15,16 of superconductivity at 203 K in a highly compressed hydrogen sulfide [13][14][15][16][17][18][19][20][21][22][23] . As pointed out by Bernstein et al. in a microscopic theory of H 3 S 17 , the substantial covalent metallicity leads to a large electron-phonon coupling, and the low atomic mass leads to high-frequency phonon modes. Both features play essential roles in increasing the T c of H 3 S. Notably, the former feature is analogous to that of MgB 2 25,26 whereas the latter effect is similar to those in the H-rich materials 27-33,35? -38 .
As a powerful method for the optimization of T c , the atomic substitution has been routine in superconductivity experiments 39,40,42? . However, there is yet any study on how the atomic substitution influences the T c of H 3 S and of the more general hydrogen chalcogenides. Here, we systematically examine the effects of (partially) substituting the chalcogenide atoms on the superconductivity and particularly the T c 's of hydrogen chalcogenides in the Im3m phase under high pressures, based on the first-principles calculations with virtual crystal approximation (VCA) [43][44][45][46][47] . In the V-and VII-substitution cases, we find that the significant changes of the DOS at Fermi surface and of the phonon linewidths, coming from the different number of valance electrons, are the principal factors affecting the electron-phonon coupling and T c . In the particular case of H 3 S 1−x P x at 200 GPa, we find that the DOS, the electron-phonon coupling constant, and the T c all first increase and then decrease as the P-substitution rate increases from zero to 0.2. In the optimized condition in which x = 0.075 and the pressure increases to 250 GPa, a possible record high T c of 280 K is predicted. In contrast, the T c decreases monotonously with increasing the VII-substitution degree. In the VI-substitution cases, the T c does not appear to increase substantially, because of the reduction of DOS with O-or Se-substitution, or due to the softening of logarithmically averaged phonon frequency with Te-substitution. These findings emphasize the importance of low atomic mass and strong covalent metallicity in conventional high-T c superconductors, pave the way for substantially enhancing T c by combining application of a high pressure and properly designed chemical substitution, and suggest that in principle it is not impossible to boost T c to ice point in optimized conditions.

II. RESULTS
Soon after the experimental report 15,16 on the record T c of a highly-compressed hydrogen sulfide 13,14 , further theoretical studies 17-23 have reached two consensuses: (i) the superconducting matter is most likely H 3 S with a strong covalent character in the Im3m phase, as sketched in Fig. 1a, and (ii) the superconductivity of H 3 S is phonon-mediated. Thus, the superconductivity of highly-compressed hydrogen sulfide can be accurately described by the celebrated Eliashberg theory 48 , which takes into account the renormalization of electronelectron repulsion by electron-phonon interactions. In this theory, the Allen-Dynes-modified McMillan formula 49,50 relates T c to the logarithmically averaged phonon frequency ω ln , the effective Coulomb repulsion µ * , and the electron-phonon coupling constant λ: where f 1 and f 2 are the strong coupling and the shape correction factors 49  enhance T c , i.e., to search for materials with high-frequency phonon modes and to increase electron-phonon couplings. As mentioned in the introduction, the H 3 S happens to constitute both advantages and hence has a record high T c . Indeed, within this framework recent theoretical calculations have well explain the high T c in the highly-compressed hydrogen sulfide [13][14][15][16][17][18][19][20][21][22][23][24] . Therefore, we will first reproduce the T c of H 3 S in the Im3m phase at 200 GPa and study the substitution of S atoms by other VI atoms. Comparing our results with previous ones will show the validity of our first-principles calculations with VCA.
In the VI-substitution systems at 200 GPa, we focus on studying Because when x takes a greater value in the case of O-or Te-substitution, we can find imaginary phonon modes at Γ point, indicating a lattice structure instability. The electronic band structures of these systems are shown in the supplementary materials, and our results for H 3 S and H 3 Se are in great agreement with previous reports based on similar first-principles calculations 14,51,52 . We find that the VIsubstitutions have little influence on the electronic band structure around the Fermi energy, except unimportant changes near Γ point (see supplementary materials). We further find that including spin-orbit couplings for the Te-substitution case does not introduce any significant correction, either. Notably, the main electronic effect is the reduction of DOS at Fermi surface in all of these VI-substitution cases, as summarized in Table I. The DOS decreases monotonically with increasing the VI-substitution rate. Moreover, the lighter the substitution element, the stronger the reduction of DOS.
Figures 1b-1e display the phonon spectra of four representative VI-substitution systems, and the magnitude of the phonon linewidth is indicated by the size of the red error bar. Our results for H 3 S and H 3 Se are again consistent with the earlier reports 19,52 . Overall there is a clear separation between H modes at high energy and S modes at low energy 19 . Evidently, with increasing the degree of heavy-element (Se or Te) substitution, all the three acoustic phonon modes decrease in frequency. In addition, the H-VI bond-bending modes (displacement of an H atom perpendicular to a H-VI bond) in the intermediate frequency region are softened, whereas the H-VI bond-stretching modes (displacement of an H atom parallel to a H-VI bond) in the high frequency region are hardened. Not surprisingly, the opposite trends occur for the case of Osubstitution. The above changes at low and intermediate frequency are expected by considering the relative atomic mass of the substitution elements, whereas the changes at high frequency might be explained by the dependence of the chem-    ical precompression effect 28,52,53 on the atomic radius of VI elements.
Further analyzing phonon linewidths reveals that they are larger for the H vibrational modes, particularly for the bondstretching modes at high frequency, similar to the results by Errea et al. 19 We note that the small magnitude of phonon linewidth in H 3 S 0.4 O 0.6 (Fig. 1d) is ascribed to the reduction of DOS. In order to expose more clearly the relative contributions of different phonon modes to the electron-phonon coupling constant λ, we define an integral function λ (ω) = 2 ω 0 ω −1 α 2 F (ω )dω , where λ (∞) is λ and α 2 F (ω) is the Eliashberg function [54][55][56] . As shown in Fig. 2a, the contribution from the intermediate frequency region (300 ∼ 1700 cm −1 ) has an increase in H 3 S 0.7 Te 0.3 , which results in a corresponding increase in the electron-phonon coupling constant. Such a variation is a consequence of the increase of phonon linewidths in the intermediate frequency region, which is present in Fig. 1e but absent in the O-and Sesubstitution systems (Figs. 1c-1d).
With the above results, we use equation (1) to estimate the T c for the investigated VI-substitution systems and the results are plotted in Fig. 2b. The effective Coulomb repulsion has been chosen to be µ * = 0.12, in the reasonable range of 0.1 ∼ 0.15 57 , such that the estimated T c for H 3 S and H 3 Se, 194 K and 160 K, are most close to the reported values 19,51 . In all of the VI-substitution cases, there is no substantial enhancement of T c compared with the 194 K in H 3 S, except a slightly higher T c of 198 K in H 3 S 0.9 O 0.1 . Notably, despite a higher λ in H 3 S 0.7 Te 0.3 , the T c is in fact reduced to 161 K. Because of the heavy mass of Te, the logarithmically averaged phonon frequency ω ln is only 90.8 meV in H 3 S 0.7 Te 0.3 , much lower than the 125.2 meV in H 3 S. In fact, the effect of VI-substitution at a fixed physical pressure can be viewed as a chemical pressure effect. Thus, our results that the T c would decrease in the various VI-substitution cases at 200 GPa are in agreement with the experimental observation 16 and an earlier theoretical study 21 , which have confirmed that the T c of H 3 S is peaked near 200 GPa.
Besides the substitution of S by elements in the same main group, we also further consider substitutions by the adjacent group elements, i.e., phosphorus and halogen groups. We mainly focus on the cases in which a small percentage (x = 0.0 ∼ 0.2) of S is replaced by P (Cl), and likewise Se replaced by As (Br). In sharp contrast to the VI-substitution cases, the averaged phonon frequency is not expected to have a significant drop after substitutions, because the atomic mass of after chemical substitutions are close to that of the prototypes. Thus, as the calculations will show below, the changes of T c follow closely the changes of λ in the V-and VIIsubstitution cases. Note that the absence of imaginary frequency in the phonon spectra, as exhibited in Fig. 3 (also see supplementary materials), ensures the lattice dynamic stabilities for all the cases of interest.
Due to the decrease of valence electrons after Vsubstitutions, the electronic energy bands shift up a little bit with respect to the Fermi energy (see supplementary materials). Such changes can also be reflected effectively by the changes in DOS at the Fermi surface, as summarized in Table I. It turns out that the DOS increase first and then decrease with increasing the V-substitution rates. Notably, the DOS can reach very large values, e.g., 10.55 in H 3 S 0.9 P 0.1 and 8.38 in H 3 Se 0.85 As 0.15 , in units of Hartree −1 per spin. Although the phonon spectra hardly change after V-substitutions, the phonon linewidths follow a close trend of the changes in DOS.  Fig. 2d. Analysis of λ (ω) plotted in Fig. 2c further confirms the enhancement of electron-phonon couplings in the intermediate and high frequency regions in H 3 S 0.9 P 0.1 . Very similar behaviors can be found as the As-substitution rate increases from zero to x = 0.2 in H 3 Se 1−x As x , as seen in Figs. 2c-2d (also see supplementary materials). We find that the maximal coupling constant values are 2.12 in H 3 S 0.9 P 0.1 and 1.97 in H 3 Se 0.85 As 0.15 . Based on the above results, we use equation (1) with µ * = 0.12 to estimate T c for the two V-substitution systems and plot them in Fig. 2d. The T c are found to be as high as 250 K in H 3 S 0.9 P 0.1 and 185 K in H 3 Se 0.85 As 0.15 , which are greatly enhanced from the 194 K (160 K) in the prototypical H 3 S (H 3 Se). Given the possible variation of the effective Coulomb repulsion µ * in the range of 0.1 ∼ 0.15, the T c in H 3 S 0.9 P 0.1 may also vary in the range of 227 ∼ 265 K.
We now take into account the influence of varying the high pressure on the superconductivity 15,16 in the case of H 3 S 1−x P x . Since we are particularly interested in optimizing the T c , we focus on two cases suggested by Fig. 2d: H 3 S 0.9 P 0.1 and H 3 S 0.925 P 0.075 . The main results are summarized in Table II and Fig. 4. In the case of H 3 S 0.9 P 0.1 , the T c reaches the highest value at 200 GPa and remains the same value up to 225 GPa, whereas in the case of H 3 S 0.925 P 0.075 , the T c increases monotonously as the pressure rises from 150 GPa to 250 GPa, beyond which a structure instability is found. Results in the latter case follow the fact that the DOS at the Fermi energy, the phonon linewidths, and the electronphonon coupling constants all increase gradually with increasing the pressure. For µ * = 0.1, the T c of H 3 S 0.925 P 0.075 at 250 GPa can reach 280 K, which is higher than the ice point. We note that the pressure of 250 GPa has already been feasible in the H 3 S experiment by Drozdov et al 16 . Given the range of µ * = 0.1 ∼ 0.15, the T c in such an optimized case is at least as high as 241 K, as indicated by the error bar in Fig. 4c. (For µ * = 0.12 that yields the T c of 194 K in H 3 S, we find T c = 264 K.) Finally, in the VII-substitution cases we have studied two cases: H 3 S 1−x Cl x and H 3 Se 1−x Br x . In general, the electronic energy bands shift down with respect to the Fermi energy, because of the increase of valence electrons; the DOS at the Fermi energy, the phonon linewidths, and hence the electron-phonon coupling constants all decrease monotonously as the degree of VII-substitutions increases. As a consequence, the T c drops monotonously with increasing the VII-substitution degree. These results are summarized in Table I and Figs. 2e-2f (also see supplementary materials).

III. DISCUSSION
We have investigated the influence of partial atomic substitution on the superconductivity of hydrogen chalcogenides  using first-principles calculations with VCA. Our study has highlighted the key roles of strong covalent metallicity and low atomic mass in boosting the T c of BCS superconductivity. The former can produce large electron-phonon couplings, whereas the latter can yield high-frequency phonon modes, and in fact the highly-compressed H 3 S constitutes both advantages. We now take H 3 S as the example to summarize our results. In the VI-substitution cases, the reduction of DOS at Fermi energy dilutes the covalent metallicity even though the coupling constant can be enhanced in the Te-substitution case, the T c is lowered because of the stronger atomic mass. This leads us to further study the cases of Cl-and P-substitutions, in which the atomic mass remains hardly changed. In the Psubstitution case, the DOS, the phonon linewidths, the coupling constant, and hence the T c all increase as the substitution rate increases from zero up to x = 0.1. In sharp contrast, the oppose trend occurs in the Cl-substitution case.
In particular, we have shown that in the optimized case of H 3 S 0.925 P 0.075 the T c may reach a record high value of 280 K at 250 GPa, which is a feasible pressure in current experiments and would not induce any structure instability. In light of our results, it might be possible that the silicon-substitution would also enhance the T c of high-pressure H 3 S superconductor. Because of the even less valance electrons, the optimum T c of H 3 S 1−x Si x might occur at even lower substitution concentration, which might be easier to realize in experiment. For ex-ample, we find T c = 274 K in the case of H 3 S 0.96 Si 0.04 (see supplementary materials). In the future, inclusion of anharmonic effects 19 may improve our predictions and thus better guide the experiments. Nevertheless, our finding is exciting. It not only suggests that partial atomic substitution may lead to possible superconductivity above the ice point in the highly compressed H 3 S, but also gives a hope that in principle low atomic mass and strong covalent metallicity may be designed in novel materials to realize high-T c BCS superconductivity under ambient pressures.

IV. METHODS
There are mainly two different procedures for disordered systems and partial atomic substitutions in first-principles calculations, i.e., an ordered supercell and the virtual crystal approximation (VCA) 43 . The former method is timeconsuming and technically difficult to deal with the case for small concentrations. Thus, we chose to use VCA in this work. In calculations based on VCA, the primitive periodicity is retained but composed of virtual atomic potentials interpolating between the behaviors of actual components. Even though this approach neglects the local deformations around atoms and cannot explore the disordered structures very accurately, it often produces acceptable and useful re-sults that have been verified in many research fields of condensed matter physics 42,44-47,59? -63 . The atomic substitution in the present work is simulated by the self-consistent VCA. For example, in H 3 S 1−x Se x the virtual pseudopotentials of the S 1−x Se x is represented by a pseudopotential operator V VCA = xV Se + (1 − x)V S , where V Se (V S ) is the pseudopotential of Se (S) atom.
The present studies, including the electronic structures, the phonon spectra, and the electron-phonon couplings, were carried out using the ABINIT package [64][65][66][67] with the local-density approximation (LDA). Hartwigsen-Goedecker-Hutter (HGH) pseudopotentials 68 were used in order to include spin-orbit couplings (SOC) for the heavy element, tellurium. The SOC for other elements were neglected since they are sufficiently light. By requiring convergence of results, the kinetic energy cutoff of 30 Hartree and the Monkhorst-Pack k-mesh of 40×40×40 were used in all calculations about the electronic ground-state properties. The phonon spectra and the electron-phonon couplings were calculated on a 8×8×8 qgrid. Since calculating electron-phonon couplings referred to integrals over the Brillouin zone, we also carefully checked convergence of the results on the aforementioned k-mesh and q-grid, by comparing them with results in denser samples (a 40×40×40 k-mesh and a 10×10×10 q-grid).