Full-bandwidth anisotropic Migdal-Eliashberg theory and its application to superhydrides

Migdal-Eliashberg theory is one of the state-of-the-art methods for describing conventional superconductors from first principles. However, widely used implementations assume a constant density of states around the Fermi level, which hinders a proper description of materials with distinct features in its vicinity. Here, we present an implementation of the Migdal-Eliashberg theory within the EPW code that considers the full electronic structure and accommodates scattering processes beyond the Fermi surface. To significantly reduce computational costs, we introduce a non-uniform sampling scheme along the imaginary axis. We demonstrate the power of our implementation by applying it to the sodalite-like clathrates YH$_6$ and CaH$_6$, and to the covalently-bonded H$_3$S and D$_3$S. Furthermore, we investigate the effect of maximizing the density of states at the Fermi level in doped H$_3$S and BaSiH$_8$ within the full-bandwidth treatment compared to the constant-density-of-states approximation. Our findings highlight the importance of this advanced treatment in such complex materials.


INTRODUCTION
Discovering and designing new and technologically relevant superconductors is one of the grand challenges of modern science [1].Conventional superconductivity arises from an intricate interplay between the electrons and the vibrational modes of the lattice, which can be condensed into a single parameter known as the electron-phonon (el-ph) coupling strength λ.This interaction leads to pairing electrons with opposing spins below the critical temperature T c , creating an energy gap at the Fermi surface and resulting in a zeroresistance superconducting condensate.Since the pioneering work of Bardeen, Cooper, and Schrieffer (BCS) [2], advancements in computational and theoretical techniques have allowed accurate calculations of λ and fully ab-initio predictions of T c [3].The density-functional theory for superconductors [4][5][6] and the anisotropic Migdal-Eliashberg theory (AME) [7,8] are state-of-the-art examples of such techniques that have contributed significantly in unraveling the properties of the superconducting states of seminal materials like MgB 2 [9][10][11] and NbS 2 [12] in unprecedented detail, and in predicting entirely novel classes of superconductors from first principles [13].One of the most topical examples is the class of the high-pressure superhydrides [14][15][16], which have revolutionized the search for high-T c superconductivity by demonstrating that detailed calculations of the electronic structure, phonon dispersion, and el-ph coupling can guide experiments in the search for new superconducting materials.Prominent examples would be LaH 10 , theoretically predicted in 2017 [17,18] and experimentally confirmed two years later [19], YH 6 [20][21][22][23] and CaH 6 [24][25][26], or most recently, LaBeH 8 , the first successfully synthesized ternary superhydride [27,28].
The AME formalism is particularly useful in describing the order parameter of weak and strong coupling su-perconductors.However, computing the el-ph matrix elements and numerically solving the Eliashberg equations requires extremely dense electron and phonon meshes in the Brillouin zone (BZ) to overcome the strong sensitivity to the sampling of the el-ph scattering processes involving states around the Fermi level [29].The AME implementation of the EPW code [30,31], which was developed by some of the present authors [32], enables the interpolation of a small number of el-ph matrix elements to arbitrary electron and phonon wave vectors in the Bloch representation using maximally localized Wannier functions [33].This has helped to bridge the gap between experiments and theory and has been widely used in the last few years to determine, among other superconducting properties, the momentum-and band-resolved superconducting order parameter of various anisotropic bulk materials [34][35][36][37][38], layered compounds [39][40][41], and two-dimensional systems [42][43][44][45][46].
Despite the extraordinary success of the AME implementation in EPW, its main shortage comes from the assumption that the density of states (DOS) is constant for a finite energy window around the Fermi level ε F (of the order of the Debye energies) where the superconducting coupling occurs [32].This approximation, widely employed in literature, is valid for a broad range of compounds but will break down for materials with narrow bands or critical points in the vicinity of ε F [47,48], such as van Hove singularities (VHSs) and Lifshitz transitions.
With the present work, we remedy this shortcoming.Our implementation goes beyond the limitations of the previous approach, by explicitly incorporating scattering processes of electrons with energies and momenta beyond the confines of the Fermi surface.This is made possible through the selfconsistent determination of the mass renormalization function, energy shift, and order parameter at every temperature while ensuring the system's charge neutrality (see Supplemen-tary Method 1).The corresponding theoretical considerations and equations are detailed in the methods section and in Supplementary Method 1.As this leads to an increased computational workload, we have also implemented a sparse, nonuniform sampling scheme over the imaginary axis, considerably lowering the number of Matsubara frequencies needed compared to the uniform sampling scheme, which, in practice, highly decreases the computational costs (see methods section).
We apply this implementation to two different classes of topical superhydrides, the sodalite-like clathrates YH 6 and CaH 6 , and the covalently-bonded H 3 S and D 3 S. Results and discussion are provided under General applications and benchmarking and material-specific computational details can be found in the methods section.Moreover, we present compelling evidence that the commonly used approach of computationally optimizing T c by maximizing the DOS at ε F (N(ε F )) is often ineffective, as the reported enhancements in T c are, in fact, artifacts resulting from the constant-DOS approximation.To support this claim, we conduct a detailed study of electron-and hole-doping effects in H 3 S and BaSiH 8 using our full-bandwidth implementation (see Application to doped hydrides).The results shed light on the limitations of the maximizing-N(ε F ) strategy, emphasizing the need for a more comprehensive and accurate approach in predicting superconducting properties, as provided with our implementation.

General applications and benchmarking
The emergence of superconductivity at record-breaking temperatures has reignited the hope of achieving superconductivity at ambient conditions [1].Indeed, the discoveries of near-room temperature superconductivity in H 3 S [49,50] and LaH 10 [18,19,51] at megabar pressures constitute a new landmark for superconductivity and for the prediction of entirely new materials with advanced functionalities fully ab-initio.
Among the numerous already predicted and experimentally confirmed superhydrides, H 3 S and D 3 S have received particular attention from the scientific community.Multiple independent experimental groups employing different characterization techniques have confirmed the existence of the cubic Im3m-H 3 S structure at pressures around 100-200 GPa and its superconducting state near room temperature [52][53][54][55][56][57][58][59].Furthermore, after the prediction of stable body-centered cubic structures of hydrogen that form sodalite-like cages containing Ca [20] and Y [24] atoms above 150 GPa, YH 6 and CaH 6 have also been comprehensively studied by many independent theoretical [17,18,60,61] and experimental groups [21-23, 25, 26], paving the way for the search for the holy grail of superconductivity [27,37,[62][63][64].Many hydrides exhibit interesting features, such as van Hove singularities near ε F and metallic hydrogen states with strong el-ph coupling [1], which make them a unique condensed matter platform to study superconductivity.
Due to their topical relevance and the amount of experimental data available, we have employed the full-bandwidth method to the sodalite-like clathrates YH 6 and CaH 6 , and to the covalently-bonded H 3 S and D 3 S.In the following, we will consider two levels of approximation when solving the AME.The first is the Fermi-surface-restricted approximation (FSR) [30], which, as discussed in detail in the methods section, assumes that the DOS around ε F is constant.The second, and main object of interest in this work, is the full-bandwidth method (FBW), which takes into account the full energy dependence of the DOS and thus allows for the inclusion of el-ph scattering processes away from ε F [31].
Furthermore, our implementation of FBW comes in two different flavors: (i) updating the chemical potential, µ, while solving self-consistently the AME equations to maintain the charge neutrality of the system (referred to as FBW+µ henceforth); and (ii) keeping the chemical potential fixed (referred to simply as FBW) to lighten the computational load.We will also compare our results for T c to the values given by the commonly employed semi-empirical modified McMillan equation [65,66] and the recently proposed machine-learned SISSO model [67].To benchmark our implementation, we first take a look at YH 6 and CaH 6 , two hydrides whose variation in the DOS around ε F is rather small, i.e., they are materials for which the FSR approach should be reasonably accurate.
In the following, we summarize the important physical properties to understand the emergence of a high-T c in these materials: The electronic band structures and DOS for YH 6 and CaH 6 at a pressure of 200 GPa are presented in Figs.1(a) and (d).The corresponding phonon dispersions and phonon DOS along with the isotropic Eliashberg spectral function α 2 F(ω) and the cumulative el-ph coupling strength λ(ω) are reported in Figs.1(b) and (e).The relatively high el-ph coupling is associated primarily with the Kohn anomalies observed in the phonon dispersion along the Γ-H direction for YH 6 and the H-N direction for CaH 6 [24].For YH 6 , numerous modes between 90 meV and 220 meV significantly contribute to the total el-ph coupling strength.Conversely, the primary contribution for CaH 6 is localized in the energy range of 100-160 meV, derived from the T 2g and the E g modes at the Γ-zone center belonging to the vibrations of the H 4 units [20].By integrating α 2 F(ω), the total el-ph coupling strengths for both YH 6 and CaH 6 are close to 2.0.This value and the α 2 F(ω) functions are in excellent agreement with the calculations presented in Refs.[20,60].
Figures 1(c) and (f) depict the anisotropic superconducting gap ∆ nk as a function of temperature for YH 6 and CaH 6 at 200 GPa using the FSR, FBW, and FBW+µ implementations, with µ * = 0.16.YH 6 exhibits two well-defined super- conducting gaps on the Fermi surface; a larger, broad energy gap ranging from 35 meV to 56 meV at low temperatures, and a smaller gap at approximately 28 meV, the latter originated from the small zone-centered Fermi surface pockets [12].The gaps of YH 6 close at T c ≈ 250 K within the FSR approximation and at T c ≈ 238 K within the FBW approach, independent of whether the chemical potential is updated selfconsistently or not.CaH 6 is a single-gap superconductor featuring a well-defined gap energy with a maximum of approximately 45 meV and broadness of about 5 meV; the T c value is around 200 K in all three implementations, FSR, FBW, and FBW+µ.
Experimentally, the critical temperature of YH 6 varies depending on the synthesis route.Maximum values for T c ranging from 220 K at 183 GPa [21] to 224 K at 166 GPa [23] have been reported.In contrast, CaH 6 exhibits a wider transition width of approximately 25 K, with an onset T c of 195 K at 185 GPa [25].A maximum T c of 215 K for CaH 6 has been observed at 172 GPa [26].At 200 GPa experiments report a T c of about 211 K for YH 6 and about 204 K for CaH 6 .Considering the considerable variation of T c with respect to different samples and the uncertainties for a particular measurement, the agreement between experimental values and our numerically determined ones is quite satisfactory.We expect that by in-cluding corrections from quantum anharmonic effects, which can be sizeable in hydrides [68][69][70], the difference between experiment and theory could be even further reduced, which will be a topic of future investigation.
As mentioned before, both YH 6 and CaH 6 exhibit a slowly varying DOS around ε F and in such cases, the FSR approximation is reasonable for describing the el-ph scattering process around the Fermi level.This is also evident from the fact that the chemical potential remains almost constant throughout the self-consistent solution within FBW+µ.Nevertheless, it is important to note that the FBW implementation offers a better agreement with experiments even for these simple cases.In contrast to the XH 6 materials described above, the behavior of the DOS around ε F is quite different in the covalently bonded H 3 S and D 3 S hydrides.As indicated in Fig. 2(a), the Fermi level is located right at the shoulder of a marked peak in the DOS, which will give rise to considerable differences between the FSR and FBW approaches.and (d) report the phonon dispersion, the phonon DOS, α 2 F(ω), and the cumulative λ(ω) for H 3 S and D 3 S at a pressure of 200 GPa.The corresponding α 2 F(ω) functions possess two main peaks in both compounds.For H 3 S, the dominant one is centered around 120 meV, and the second one, less intense, around 190 meV.For D 3 S, the peaks are shifted to lower frequencies, as expected due to the greater mass of deuterium atoms, with maxima centered around 90 meV and 130 meV, respectively.In both hydrides, the whole optical spectra from 30-40 meV to the Debye frequency contribute to the total electron-phonon coupling, which is found to be λ = 2.3 for H 3 S and λ = 2.2 for D 3 S, respectively.These values and the corresponding α 2 F(ω) functions are in excellent agreement with those reported in Refs.[49,71,72].
The solution of the AME equations reveals that H 3 S and D 3 S are single-gap superconductors with a broad energy gap distribution.Compared to the FSR treatment, the FBW calculation with fixed µ lowers the gap energy [see Figs.2(c) and (e)].This effect is even more pronounced when updating µ self-consistently (FBW+µ), as the chemical potential is shifted to higher energies, moving the Fermi level away from the peak of the van Hove singularity.These results emphasize the critical role of the VHS on the superconducting properties of H 3 S, as has been pointed out in other works as well [48,71,[73][74][75].
The highest measured T c in the study of Drozdov et al. [50] is 203 K at 155 GPa for H 3 S and 152 K at 173 GPa for D 3 S, with a variation of T c for different samples of up to 15 K.Samples with better crystallinity for H 3 S were later obtained by Mozaffari et al. [57] with a T c of 201 K at 155 GPa and a small transition width of 5.5 K; and also by Nakao et al. [56], where a sharp drop of the resistance was measured at T c = 200 K at 150 GPa.Minkov et al. [59] have used the same direct in-situ synthesis from elemental S and excess H 2 as in Refs.[56,57] to obtain better homogeneous samples for D 3 S, revealing that D 3 S reached a maximum T c of 166 K at 157 GPa, significantly higher (≈ 10 % difference) than previously reported values [50].
As can be appreciated in Tab.I, FBW+µ performs best in approaching the experimental critical temperatures T c exp among the different implementations for solving the AME equations.For H 3 S, our calculations provide T c = 232 K at 200 GPa, a percentage difference of 23-30 % compared to the experimental values of 172-184 K at 200 GPa from Refs.[50,53,59].The rather large differences originate from anharmonicity and the quantum motion of the nuclei, which are known to play a crucial role in H 3 S [68].Incorporating these effects is, in principle, possible within EPW, as demonstrated in Refs.[70,76] for example, but beyond the scope of the current work.Here, it is important to point out that FBW+µ provides a much better estimate for T c than FSR.The performance of FBW+µ in reproducing the experimental values is notably better for D 3 S, where anharmonic and quantum ionic effects are smaller.The full-bandwidth treatment only Table I.Summary of the obtained superconducting properties for YH 6 , CaH 6 , H 3 S, and D 3 S: The table lists the DOS N(ε F ) and the H partial DOS N H (ε F ) at the Fermi level, the electron-phonon coupling parameter λ, the logarithmic average of the phonon frequencies ω log , and the superconducting critical temperature T c according to the modified McMillan formula (mMc), the University of Florida machine-learning model (UF), the Fermi-surface-restricted approximation (FSR), and the full-bandwidth implementation without (FBW) and with (FBW+µ) the self-consistent chemical potential update scheme, and all experimental critical temperatures T c exp at 200 GPa available in the literature.The T c exp of Ref. [21] and Ref. [26] are used as reference for YH 6 and CaH 6 , respectively; and the T c exp of Refs.[50,53,59] are used as reference for H 3 S and D 3 S. slightly modifies the structure of the superconducting gap ∆ nk , i.e., the gap distributions are shifted to lower energies while retaining the overall shape at each temperature.The isotope effect coefficient, according to the BCS theory, is given by where M H and M D are the atomic mass of hydrogen and deuterium.The experimental values obtained for α are around 0.47 at 150 GPa [52].Within FSR, we obtain α = 0.54 at 200 GPa; for FBW we have α = 0.48; and for FBW+µ we have α = 0.45.These results demonstrate that the full-bandwidth method is imperative for accurately describing these systems within the Eliashberg formalism.As already pointed out by other authors [48,71,73], the conventional Eliashberg formalism within the FSR framework partially fails to accurately describe the T c behavior in H 3 S due to a substantial variation of the DOS near ε F originating from the van Hove-type singularity present there.The strong el-ph coupling of λ = 2.3 for H 3 S and λ = 2.2 for D 3 S, and the broad distribution of α 2 F over the vibrational spectra makes this scenario even more dramatic since the region around ε F , where the phonon scattering dominates over the Coulomb repulsion, is strongly enhanced by the el-ph interactions.Lastly, these results highlight the importance of updating the chemical potential while solving the AME equations within the FBW treatment.
At the end of this section, we also want to shortly discuss the modified McMillan [65,66] formula and a recent machine-learning approach to improve upon it [67], both of which can be used as an almost instant way to predict T c , but do not offer insight into other properties of the superconducting state.The modified McMillan formula (mMc, see Supplementary Note 1) is obtained from the Eliashberg theory by defining moments of the α 2 F spectral function and fitting an equation to match the experimental T c , taking into account a bit more than 200 data points.As the data set consisted of the (low-T c ) superconductors known at the time and a few pure model calculations, it does not reproduce the experimental critical temperatures of the extreme cases of contemporary highly-compressed high-T c materials, as is also evident for the materials chosen in this study [66].
The machine-learned equation proposed by a group of the University of Florida (UF), on the other hand, has been trained to match the solutions of the Migdal-Eliashberg equations with a dataset of thousands of real and artificially generated α 2 F functions, including high-pressure, high-T c hydrides as well.It thus performs much better than the mMc formula and matches the results obtained with AME fairly well.It is therefore a viable tool to determine an accurate value for T c quickly but does not offer insight into the superconducting gap function and its energy distribution, the superconducting DOS, and so on.We also observe shortcomings when trying to simulate the effects of doping, as detailed in the next section.

Application to doped hydrides
Doping can be a powerful method to tailor and fine-tune specific material properties, especially in electronic applications.It has been used to metallize semiconducting phases, induce superconductivity, and optimize specific properties of the superconducting phase [41,42,[77][78][79][80].In particular, this approach has been employed to increase the T c in known superconducting systems and is claimed to be a route to obtain (or to be the source of) room-temperature superconductivity in recently reported hydrides [64,[81][82][83][84][85][86][87].
The doping of hydrogen-rich superconductors and its tentatively beneficial effect on T c has been extensively studied in various works [83,84,[88][89][90][91][92][93][94][95][96][97][98], however, most of the theoretical predictions for an enhancement of T c rely on calculations based on the McMillan/Allen-Dynes formulas, or on the ME equations within the constant-DOS approximation, and hence focus mainly on maximizing the value of the DOS at the Fermi level and thereby that of λ.In particular, in systems with VHS-like peaks close to the Fermi level, this approach might severely overestimate the contribution of the electronic states available for the superconducting pairing and thus also the critical temperature.Furthermore, the exact value of the DOS close to a VHS-like shape is subject to large variations, which makes predictions even more error-prone [73].With the FBW implementation, we overcome these problems and limitations as demonstrated for H 3 S and BaSiH 8 , two materials exhibiting considerable variation in the DOS around the Fermi level.Applying doping via a rigid band shift, we explore the different regions of increased, maximal, and lowered DOS, and discuss the effects on T c within the mMc and MD formulas, and the FSR and FBW methods.We want to note here that for these calculations we do not update the chemical potential, as an explicit shift in the Fermi level can also be interpreted as a shift in the chemical potential, and hence such a calculation can be reproduced by one with a (slightly) different effective shift.
The well-studied hydride superconductor H 3 S is a perfect candidate material where doping to optimize T c appears very tempting due to the VHS-like peak in the DOS in close proximity to the Fermi level, as shown in Fig. 3(a)-(b).Many theoretical works have considered doping of H 3 S to increase its T c estimates, choosing dopants to bring the system's Fermi level closer to the maximum of the VHS-like peak.The stability of the Im3m H 3 S structure with incorporation of various elements has been systematically investigated using either the direct supercell approach with substitutional dop-ing or the virtual crystal approximation, followed by T c estimates based on the mMc formula or the isotropic FSR approach [83, 88-93, 95, 98-100].
While adding dopants may increase the T c by also enhancing the el-ph coupling, we want to specifically address the effect of a change in the number of available electronic states.To this end, we solved the AME equations in the FSR and FBW frameworks for shifts of the Fermi level ∆ε F between −0.5 eV and 0.1 eV in steps of 0.1 eV, corresponding to changes in the electron number ∆ elec of -0.30, -0.26, -0.20, -0.13, -0.06, and +0.05.
In Fig. 3(a), we present the T c values obtained within the FSR and FBW treatments of the AME equations and the results obtained using the mMc and UF semi-empirical formulas mentioned earlier.The T c values are plotted as a function of ∆ε F demonstrating that all approaches show a clear correlation between T c and DOS with the maximum T c occurring at around ∆ε F ∼ −0.1 eV to −0.2 eV, as highlighted in Fig. 3(b).
However, within the FSR approximation, there is a notable and unphysical increase in T c values around the maximum of the DOS.This behavior can be attributed to the limitation of the constant-DOS assumption in the FSR approach.Similarly, the UF model exhibits a strong dependence on doping.In contrast, the AME solutions in the FBW approach, which consider the full energy dependence of the DOS, exhibit a much less pronounced effect of doping on T c .Finally, the semi-empirical mMc formula consistently underestimates T c for strongly coupled systems when compared to results obtained from the AME equations.
In the context of our analysis of a doped H 3 S system, we would like to address the topic of speculated room temperature superconductivity within this system [101].Considering that the experimental T c values for pure H 3 S are around or below 200 K, and that the FBW calculations only show a maximum increase in T c of about 5-10% upon doping (see Fig. 4), achieving a conventional superconducting state at room temperature in the H 3 S parent phase remains elusive; at least within the assumption that the slight doping leaves the electronic structure unaltered.
On a more technical note, we want to add that VHS-like features also pose other problems when trying to arrive at robust numerical results.For example, the exact value of the DOS at the VHS-like peak is difficult to converge and requires extremely dense k-grids and small smearing values.To demonstrate this, we examined the VHS region by performing additional calculations with two different smearing values σ for the energy-conserving δ-functions.The results are shown in Figure 4.The critical temperature obtained with the FSR approximation is very sensitive to the smearing value chosen, the largest deviation in T c reaching up to 70 K, whereas the FBW implementation is considerably more robust, with at most a 10 K difference between the different smearing values tested.The reason for that is that the only quantity impacted by smearing in the FBW implementation is the screened Coulomb interaction term in Eq. ( 15), whereas smearing impacts all FSR equations via N(ε F ) and δ(ε mk+q − ε F ) (see Eqs. 18 and 19).As a result, besides being more rigorous, the FBW implementation also considerably improves convergence behavior for materials exhibiting a strongly varying DOS.

BaSiH 8
The Fm3m phase of BaSiH 8 was recently predicted and described in detail by some of the current authors [37,70].This ternary hydride has the crystal structure of the XYH 8 template first introduced for LaBH 8 [63], also assumed by other high-T c superhydrides that are stable at moderate pressures [27,102], like the recently synthesized LaBecompound [28].Up to now, BaSiH 8 is the XYH 8 compound with the lowest predicted pressure of dynamical stability of 5 GPa within the harmonic phonon theory.A more realistic estimate for the critical pressure of stability and synthesizability is provided by considering the kinetic stability of the compound, which places stability at a pressure above 30 GPa [37,70].
Independent of the actual critical pressure, this material exhibits a step-like feature around the Fermi level, with an almost constant region of high DOS below the Fermi level and a sharp drop to a region of very low DOS above ε F , as can be appreciated in Fig. 3(c), making it another perfect test bed to compare the FSR and FBW approaches.
We solved the AME equations within FSR and FBW approaches for shifts of the Fermi level ∆ε F between −0.5 eV and 0.5 eV in steps of 0.1 eV, corresponding to changes in the electron number ∆ elec of -0.52, -0.42, -0.31, -0.21, -0.10, +0.09, +0.15, +0.18, +0.20, and +0.22.In Figs.3(c) and (d), we show the T c values obtained for the two levels of approximations as a function of ∆ε F .As before, T c roughly follows the shape of the DOS, but, in contrast to H 3 S, we observe a considerable increase in T c for BaSiH 8 when employing the more rigorous FBW approach.Already in the undoped case the T c is raised to 87 K, which can be further increased by doping to about 92 K when shifting the Fermi level by -0.1 eV.This would place the critical superconducting temperature of BaSiH 8 above the technologically extremely important threshold set by the boiling temperature of nitrogen.
The behaviour of T c with respect to Fermi level shifts for the different methods is quite complex: For a shift of -0.1 eV or below the AME results are similar and also agree with the UF model, while mMc gives considerably lower values for T c .In a region around 0.1 eV shift, the FBW implementation predicts a larger value for T c than FSR, and the simple mMc and UF formulas provide an even smaller estimate.For ε Fshifts larger than 0.2 eV, mMc and UF actually give the largest values for T c and FBW the smallest.In other words, while the dependence of T c with respect to Fermi level shifts was smallest in FBW for H 3 S, it is actually varying the strongest within FBW for BaSiH 8 .
These intricate results underscore the importance of taking into account the energy dependence of the electronic DOS around the Fermi level, in particular for systems where the DOS is either strongly peaked close to ε F [84,86,[103][104][105], as for H 3 S, or highly asymmetric [76,[106][107][108][109], as for BaSiH 8 .In that light, we believe that the agreement between experimental measurements and theoretical predictions can be considerably improved by employing the FBW method not only for the class of superhydrides but also for other material systems showing similar features in the DOS [42,45,[110][111][112][113][114][115][116][117][118].

DISCUSSION
In summary, we have employed the anisotropic Migdal-Eliashberg formalism within the full-bandwidth formulation utilizing maximally-localized Wannier functions as implemented in the EPW suite.This approach enables us to calculate the momentum-and band-resolved superconducting gap more accurately, taking into account the electron-phonon scattering processes around the Fermi level, and not only restricted to the Fermi surface.In addition, we introduced a sparse, non-uniform sampling scheme over the imaginary Matsubara frequencies, which shows similar accuracy and much-improved efficiency compared to the uniform sampling scheme.
To validate the robustness of our methodology, we conducted comprehensive tests on two representative classes of superhydrides: the sodalite-like clathrates YH 6 and CaH 6 , as well as the covalent hydrides H 3 S and D 3 S.To assess the accuracy of our approach, we compared our results with previous ab-initio calculations and experimental data.Our results unequivocally demonstrate the indispensable role of the full-bandwidth formulation, particularly for compounds characterized by narrow bands or critical points in proximity to the Fermi level.A noteworthy illustration of the importance of employing the FBW formulation is evident in the case of H 3 S, which possesses a van Hove singularity.Our methodology effectively captures the intricate behavior of such systems, highlighting the superiority of the FBW approach in these critical scenarios.Furthermore, we emphasize the crucial impact of the chemical potential updating scheme within the FBW formulation, which substantially contributes to accurately describing the superconducting phase in these challenging cases.This aspect proves to be a vital component in achieving a comprehensive understanding of the superconducting properties of these complex materials.
In addition, we applied the FBW approach to investigate electron-and hole-doped hydride superconductors, namely H 3 S and the recently predicted low-pressure BaSiH 8 .These materials serve as prime examples of systems with distinct DOS features that deviate significantly from the constant-DOS assumption made in the FSR approximation.Previous studies have often focused on maximizing the DOS specifically at the Fermi level, aiming to design high(er)-T c superconductors by doping the system to shift ε F to the maximum of a VHS-like structure.Our calculations reveal significant pitfalls associated with such a simplistic FSR-based approach.We find instances of pronounced under-or overestimation of T c , highlighting the critical importance of adopting the FBW method, particularly in scenarios with strongly peaked DOS or closely adjacent DOS regions exhibiting extremely high and low values.By employing the FBW approach, we can accurately capture the intricate interplay of DOS features and better predict the behavior of superconductors in these complex cases.This sheds light on the limitations of the FSR approach and underscores the significance of our advanced methodology in studying and engineering novel superconducting materials with tailored properties for real-world applications.

Anisotropic Migdal-Eliashberg theory
The Eliashberg theory is a powerful many-body perturbation approach for describing conventional superconductors, where the Cooper pairing between two electrons stems from the interplay between the attractive el-ph coupling and the repulsive screened Coulomb interaction [8,119].The Nambu-Gor'kov's formalism [120,121] with a generalized matrix Green's function can be used to formulate the Eliashberg theory.The on-and off-diagonal elements of the 2 × 2 Green's function matrix describe the single-particle excitations in the normal state and Cooper-pair amplitudes in the superconducting state, corresponding to the standard and anomalous Green's functions, respectively [11,32,[120][121][122][123][124].The transition from normal to superconducting state manifests in anomalous Green's functions becoming nonzero below a material specific critical temperature.
The matrix Green's function is obtained from the Dyson equation where Ĝ0 nk (iω j ) is the non-interacting Green's function in the normal state with band index n and wavevector k, Σpa nk (iω j ) is the pairing self-energy, and iω j = i(2 j + 1)πT is the fermionic Matsubara frequency with T being the absolute temperature and j an integer.An expression for the selfenergy in terms of the electron Green's function can be obtained with the el-ph and electron-electron contributions given by the Migdal [7,123] and GW [125,126] approximations, respectively: Using the Pauli matrices, τ0 the el-ph and Coulomb contributions to the self-energy can be expressed as and where qν ] is the dressed phonon propagator for phonons with wavevector q and branch index ν, iω l = i2lπT is the bosonic Matsubara frequency with l an integer, g mnν (k, q) is the screened el-ph matrix element for the scattering between the electronic states nk and mk + q through a phonon of frequency ω qν , and V nk,mk+q is the static screened Coulomb interaction between electrons [123,127].In Eq. ( 5), only the off-diagonal components of the Green's function Ĝod nk (iω j ) are retained in order to avoid double counting the Coulomb effects that are already included in Ĝ0 nk (iω j ) [123].The anisotropic el-ph coupling strength is described as with N(ε F ) the DOS per spin at the Fermi level.Eq. ( 6) can be used to rewrite the el-ph self-energy in Eq. ( 4), which can then be taken together with Eq. ( 5) and inserted into Eq.( 3).The pairing self-energy then becomes: Ĝod mk+q (iω j ′ )τ 3 (7) To replace Ĝnk (iω j ) in Eq. ( 7), we expand the two components of the Dyson equation ( 2) in terms of the Pauli matrices as Ĝ0 nk (iω j ) where ε nk are the Kohn-Sham eigenenergies, and In Eq. ( 9), we introduced the mass renormalization function Z nk (iω j ), the energy shift χ nk (iω j ), and the order parameter ϕ nk (iω j ).Inserting Eqs. ( 8) and ( 9) into Eq.( 2), and inverting the resulting matrix leads to the following expression for Ĝnk (iω j ): Ĝnk iω j Z nk (iω j )τ 0 + ε nk − µ + χ nk (iω j ) τ3 + ϕ nk (iω j )τ 1 + φnk (iω j )τ 2 (10) It can be easily verified that ϕ nk (iω j ) and φnk (iω j ) are proportional within an arbitrary phase, and without loss of generality, one can choose the relative phase such that φnk (iω j ) = 0 [30,123].Eq. (10) with φnk (iω j ) = 0 can be used to rewrite Eq. ( 7) for the pairing self-energy as: where Equating the different components of the Pauli matrix elements in Eqs. ( 9) and ( 11) leads to a system of three coupled non-linear equations: This set of equations is supplemented with an equation for the electron number, which determines the chemical potential µ [127,128]: where N e is the number of electrons per unit cell (see Supplementary Method 1 for a discussion about the electron number equation).Equations ( 13)-( 16) involve electronic states that are not restricted to the Fermi surface or its immediate vicinity; hence, labeled as anisotropic full-bandwidth (FBW) Migdal-Eliashberg equations [129].We have recently implemented the above-described anisotropic FBW approach in the EPW code [31].
The set of coupled equations can be solved self-consistently at various temperatures for the temperature-dependent superconducting gap ∆ nk , given by: The Padé approximation [130,131] can then be used to obtain the continuation of ∆ nk (iω j ) from the imaginary to the real axis.The superconducting temperature T c is the highest temperature at which ϕ nk (iω j ) 0 has a nontrivial solution.
The contribution of the Coulomb interaction to the Eliashberg equation, through matrix elements V nk,mk+q , can be evaluated at the same level as the el-ph interaction for the simpler versions of the Eliashberg formalism [48,69,72,132,133].To reduce the computational cost, however, the common approach is to replace the N(ε F )V nk,mk+q terms with the semiempirical Morel-Anderson pseudopotential µ * [134].In practice, the numerical value of this parameter is connected with the cutoff frequency ω max of the Matsubara frequencies.With a typical choice of ω max being ten times the maximum phonon frequency ω ph , a value of µ * = 0.1-0.2results in a satisfactory agreement with experiment for many applications.The µ * can also be calculated from first-principles, using the double Fermi surface average of V nk,mk+q [12,43,60,63].
The numerical solution of the anisotropic FBW Migdal-Eliashberg equations is computationally very demanding.A common simplification of these equations consists in restricting the energy range close to the Fermi level by introducing the unity factor +∞ −∞ dε δ(ε nk − ε), and to assume that the DOS within this energy window is constant [11,32,123,127,[135][136][137][138].It can be shown that, within these approximations, the energy shift χ nk vanishes and the requirement in Eq. (( 16)) is automatically satisfied.As a result, only two equations for Z nk (iω j ) and ϕ nk (iω j ) need to be solved self-consistently: Equations ( 18) and ( 19) are the anisotropic FSR Migdal-Eliashberg equations.They have been the basis for the superconductivity calculations in the EPW code prior to the recent developments [30,32].

Sparse sampling of Matsubara frequencies
All calculated quantities in the Migdal-Eliashberg equations depend on Matsubara frequencies, which are proportional to the absolute temperature.This results in a computational challenge since solving these equations at low temperatures (e.g., necessary for low-T c superconductors) requires a larger number of frequencies within the same energy range.As discussed in the previous section, for the FBW+µ method, one needs a Matsubara frequency cutoff 6-12 times the Fermi energy window in order to converge the chemical potential and, consequently, the superconducting gap energy, leading to a considerable increase in computational cost.While a frequency cutoff of 4 eV has been found to be sufficient for calculations with FSR and FBW with fixed chemical potential at the Fermi level in the case of H 3 S, adopting FBW+µ required a Matsubara frequency 2-3 times larger as shown in Supplementary Figure 1.
The computational cost can be reduced by pruning the Matsubara frequencies [132].We implemented a sparse sampling scheme, described in the following, in the EPW code as an alternative to uniform sampling over the Matsubara frequencies.Denoting with integer N j the numerical index for the j th Matsubara frequency, iω j = i(2N j + 1)πT , the uniform and sparse grids can be obtained as: and where INT[ ] is the rounding to the closest integer, W is an adjustable weight factor, and N max is the maximum Matsubara index for a given energy cutoff ω max and temperature T : Eq. (( 21)) can be used to generate a Matsubara frequency grid with indices N j ≤ N max .For negative indices, the corresponding frequency can be easily obtained as ω −( j+1) = −ω j (with j > 0).The resulting mesh is uniform at lower frequencies, contributing the most to the summation in the Migdal-Eliashberg equations, and becomes logarithmically sparser with increasing the Matsubara frequency index.Increasing (decreasing) W results in a denser (sparser) grid sampling.With the default setting of W = 1.0, the sparse sampling scheme produces approximately 30% fewer Matsubara frequencies than the uniform one, while the first ∼40% points of the grid are still uniformly distributed.Numerical tests show that this approach maintains the accuracy of more expensive calculations that use the full uniform grid.
We systematically evaluated this sampling approach by computing the superconducting gap function within the FBW+µ method using both the uniform and sparse sampling scheme over Matsubara frequencies up to 6 eV for H 3 S, D 3 S, YH 6 , and CaH 6 .The primary outcomes are presented in Supplementary Figure 2, where the temperaturedependent behavior of ∆ nk is plotted for FBW+µ (blue lines) and FBW+µ+sparse (red lines).The results indicate that the sparse sampling method can accurately reproduce the superconducting gap structure obtained using the uniform sampling scheme.In summary, no significant difference was observed across all the compounds analyzed.Among all the hydrides, the largest deviation was observed for D 3 S, with a difference of only 3 K compared to the uniform sampling, corresponding to a percentage difference of 1.8 %.For all other hydrides, the difference between the sparse and uniform schemes was only 1 K. Thus, the sparse sampling can quantitatively reproduce the results obtained with the uniform sampling, but with a computational cost of approximately 40 % lower.

Computational details
The electronic structure calculations are performed within the Kohn-Sham scheme [139] of the density functional theory [140] as implemented in the Quantum ESPRESSO suite [141][142][143].The exchange and correlation effects are treated within the Perdew-Burke-Ernzerhof parametrization [144] using scalar-relativistic optimized norm-conserving Vanderbilt pseudopotentials [145,146].The Kohn-Sham orbitals are expanded in a plane-wave basis set with a kineticenergy cutoff of 100 Ry for H 3 S, D 3 S, and CaH 6 , and 80 Ry for YH 6 and BaSiH 8 .The charge density is computed using Γ-centered Monkhorst-Pack k-meshes [147] of 24 3 k-points for H 3 S, D 3 S, and CaH 6 , 16 3 k-points for YH 6 , and 12 3 kpoints for BaSiH 8 .The Brillouin-zone integration employs a Methfessel-Paxton smearing [148] of 0.01 Ry for H 3 S, D 3 S, and BaSiH 8 , and 0.04 Ry for YH 6 and CaH 6 .All lattice parameters and internal degrees of freedom were fully relaxed to reach a ground-state convergence of 10 −7 Ry in the total energy and 10 −6 Ry/a 0 for forces acting on the nuclei.BaSiH 8 is relaxed to a pressure of 30 GPa, all other compounds to 200 GPa.
The dynamical matrices and the linear variation of the self-consistent potential are calculated within the densityfunctional perturbation theory [149] on a regular phonon grid of 4 3 q-points for H 3 S, D 3 S, and 6 3 q-points for YH 6 , CaH 6 , and BaSiH 8 .The threshold for self-consistency is set to 10 −14 or lower.
The maximally localized Wannier functions (MLWFs) are constructed using the Wannier90 code [33,150].In the case of H 3 S and D 3 S, 10 Wannier functions are used to describe the electronic states near the Fermi level.The Wannier orbitals are three H-s-like functions and seven functions with s, p, d xy , d xz , and d yz angular momentum states associated with the S site, with a spatial spread ranging from 0.77 Å 2 to 1.53 Å 2 .For YH 6 , six H-s-like projections and five Y-d-like functions are used to construct the initial guess, resulting in a spatial spread between 1.74 Å 2 and 1.86 Å 2 .For CaH 6 , besides the six H-s-like projections, we also use sp and sp 3 hybrid orbital functions associated with the Ca site to construct the initial guess, yielding a spatial spread between 0.80 Å 2 and 1.20 Å 2 .For BaSiH 8 , we use H s, Si s, p, d z 2 , d x 2 −y 2 , and Ba p, d orbitals, giving a total of 22 Wannier functions, resulting in spatial spreads between 1.04 Å 2 and 1.79 Å 2 .
The fully anisotropic Migdal-Eliashberg equations [32] are solved using the EPW code [30,31].Electron energies, phonon frequencies, and electron-phonon matrix elements are computed on fine grids containing 48 3 kand q-points for H 3 S, D 3 S, YH 6 , and CaH 6 , and 30 3 kand q-points for BaSiH 8 .The lower boundary for the phonon frequency is set to 5 cm −1 for H 3 S and D 3 S, and 15 cm −1 for YH 6 , CaH 6 , and BaSiH 8 .The width of the Fermi window is set to 2 eV for H 3 S, D 3 S, YH 6 , CaH 6 , and BaSiH 8 .We set the Matsubara frequency cutoff to ω max = 6 eV.The smearing values for the energy-conserving δ-function and for the sum over q-space in the el-ph coupling are set to 25 meV and 0.05 meV, respectively, for H 3 S and D 3 S, 150 meV and 0.15 meV for YH 6 and CaH 6 , and 100 meV and 0.1 meV for BaSiH 8 .We solved the equations adopting a Coulomb pseudopotential of µ * = 0.16 for all materials except BaSiH 8 , where a value of 0.1 was chosen to be consistent with our previous works [37,70].The continuation of the superconducting gap along the imaginary axis to the real energy axis is determined by applying the approximate analytic continuation using Padé functions [130,131].
The doping calculations within rigid band model are performed by shifting the Fermi level in EPW before interpolating the el-ph matrix elements.All anisotropic Migdal-Eliashberg calculations are performed for electronic energies of ±1 eV around the (shifted) Fermi level and a Matsubara cutoff of ω max = 4 eV, if not stated otherwise.

Figure 1 .
Figure 1.Electron, phonon, and superconducting properties for sodalite-like clathrates: Panel (a) shows the calculated electronic band structure and DOS with respect to the Fermi energy ε F for YH 6 at 200 GPa.The solid blue lines represent the DFT bands, the dashed red lines the Wannier bands, the solid black line the total DOS, the shaded green area the projected DOS for hydrogen s states (H-s), and the dashed black line indicates ε F .Panel (b) shows the phonon dispersion (solid blue), the phonon density of states (PDOS, solid black) and its elemental contributions (shaded green and purple), the isotropic Eliashberg spectral function α 2 F(ω) (shaded ochre), and the cumulative electron-phonon coupling parameter λ(ω) (solid black).Panel (c) displays the distribution of the values of the anisotropic superconducting gap ∆ nk on the Fermi surface according to the FSR (blue), FBW (red), and FBW+µ (green) implementations for the Migdal-Eliashberg equations.Panels (d)-(f) show the corresponding results for CaH 6 at 200 GPa.

H 3 S
and D 3 S

Figure 2 .
Figure 2. Electron, phonon, and superconducting properties for the covalently-bonded materials: Panel (a) shows the calculated electronic band structure and DOS with respect to the Fermi energy ε F for H 3 S and D 3 S at 200 GPa.The solid blue lines represent the DFT bands, the dashed red lines the Wannier bands, the solid black line the total DOS, the shaded coloured areas the projected DOS for hydrogen s (H-s) and sulfur s and p states (S-s, S-p), and the dashed black line indicates ε F .Panel (b) shows the phonon dispersion (solid blue), the phonon density of states (PDOS, solid black) and its elemental contributions (shaded green and purple), the isotropic Eliashberg spectral function α 2 F(ω) (shaded ochre), and the cumulative electron-phonon coupling parameter λ(ω) (solid black).Panel (c) displays the distribution of the values of the anisotropic superconducting gap ∆ nk on the Fermi surface according to the FSR (blue), FBW (red), and FBW+µ (green) implementations for the Migdal-Eliashberg equations.Panels (d) and (e) show the corresponding results for D 3 S at 200 GPa.

Figures
Figures 2(b) and (d) report the phonon dispersion, the phonon DOS, α 2 F(ω), and the cumulative λ(ω) for H 3 S and D 3 S at a pressure of 200 GPa.The corresponding α 2 F(ω) functions possess two main peaks in both compounds.For H 3 S, the dominant one is centered around 120 meV, and the second one, less intense, around 190 meV.For D 3 S, the peaks are shifted to lower frequencies, as expected due to the greater mass of deuterium atoms, with maxima centered around 90 meV and 130 meV, respectively.In both hydrides, the whole optical spectra from 30-40 meV to the Debye frequency contribute to the total electron-phonon coupling, which is found to be λ = 2.3 for H 3 S and λ = 2.2 for D 3 S, respectively.These values and the corresponding α 2 F(ω) functions are in excellent agreement with those reported in Refs.[49,71,72].
Figures 2(b) and (d) report the phonon dispersion, the phonon DOS, α 2 F(ω), and the cumulative λ(ω) for H 3 S and D 3 S at a pressure of 200 GPa.The corresponding α 2 F(ω) functions possess two main peaks in both compounds.For H 3 S, the dominant one is centered around 120 meV, and the second one, less intense, around 190 meV.For D 3 S, the peaks are shifted to lower frequencies, as expected due to the greater mass of deuterium atoms, with maxima centered around 90 meV and 130 meV, respectively.In both hydrides, the whole optical spectra from 30-40 meV to the Debye frequency contribute to the total electron-phonon coupling, which is found to be λ = 2.3 for H 3 S and λ = 2.2 for D 3 S, respectively.These values and the corresponding α 2 F(ω) functions are in excellent agreement with those reported in Refs.[49,71,72].

Figure 3 .
Figure 3. Doping effects on the critical temperature in different approaches: Effects of doping (shifting the Fermi level) in H 3 S at 200 GPa (panels (a) and (b)) and BaSiH 8 at 30 GPa (panels (c) and (d)).Panel (a) shows the superconducting critical temperature T c as a function of the Fermi level shift ∆ε F , obtained within the mMc formula (green), the UF equation (orange), the FSR approximation (blue), and the FBW approach (red).The smaller subpanels on top show the corresponding DOS in a range of ±2 eV around the unshifted Fermi energy ε F , where the dashed lines mark the position of ε F + ∆ε F and the shaded red areas highlight the included electronic energy range of ε F + ∆ε F ± 1 eV.Panel (b) displays the distribution of the values of the anisotropic superconducting gap ∆ nk on the Fermi surface within the FSR (blue) and FBW (red) approach for H 3 S with ∆ε F = -0.1 eV, where we find the maximum absolute difference between the FSR and FBW T c (see blue and red lines in (a)).The inset shows the corresponding DOS subpanel from (a).Panels (c) and (d) show the corresponding results for BaSiH 8 at 30 GPa, where we find the maximum T c difference for ∆ε F = +0.1 eV.

Figure 4 .
Figure 4. Dependence of the doped T c results on the electronic smearing: Influence of the smearing value σ for the energyconserving δ-functions on the T c of doped H 3 S within the FSR (blue) and FBW (red) approach as a function of the Fermi energy shift ∆ε F .The solid lines represent the corresponding results shown in Figure 3(a) for a smearing value of σ = 25 meV, the dashed (dotted) lines represent the results for σ = 50 (100) meV.