Prediction of intrinsic topological superconductivity in Mn-doped GeTe monolayer from first-principles

The recent discovery of topological superconductors (TSCs) has sparked enormous interest. The realization of TSC requires a delicate tuning of multiple microscopic parameters, which remains a great challenge. Here, we develop a first-principles approach to quantify realistic conditions of TSC by solving self-consistently Bogoliubov-de Gennes equation based on a Wannier function construction of band structure, in presence of Rashba spin-orbit coupling, Zeeman splitting and electron-phonon coupling. We further demonstrate the power of this method by predicting the Mn-doped GeTe (Ge1-xMnxTe) monolayer—a well-known dilute magnetic semiconductor showing superconductivity under hole doping—to be a Class D TSC with Chern number of −1 and chiral Majorana edge modes. By constructing a first-principles phase diagram in the parameter space of temperature and Mn concentration, we propose the TSC phase can be induced at a lower-limit transition temperature of ~40 mK and the Mn concentration of x~0.015%. Our approach can be generally applied to TSCs with a phonon-mediated pairing, providing useful guidance for future experiments.


INTRODUCTION
The topological phase of superconductors (SC) has recently received intense research interest as the superconducting quasiparticles residing in the non-trivial gapless/zero-energy boundary states are considered a form of Majorana fermions. Majorana fermions are their own anti-particles 1 and obey the non-Abelian exchange statistics 2 , which can be utilized for topological quantum computation 3 . Topological superconductors (TSC) exhibit various exotic phenomena, including zero modes on the magnetic vortex 4 , "fractional" Josephson effect 5 , non-local correlation 6 , and thermal responses 7 . By now, the theoretical aspects of TSCs are reasonably well understood, but the experimental confirmation remains a great challenge due to the requirement of tuning multiple microscopic parameters like the Fermi level, magnetic field, temperature, etc. Hence, it is highly desirable to predict more TSCs and quantify experimental conditions to advance the field.
Unlike the successful first-principles prediction of electronic and topological materials, theoretical predictions of TSCs are challenging because of the uncertainty in the parameters used to construct Bogoliubov-de Gennes (BdG) Hamiltonian. Usually, only the pre-conditions of TSC, e.g., Rashba splitting 8 or topological properties [9][10][11][12] in the normal state of known SC, were analyzed using first-principles method, but not the topology of superconducting quasi-particles. Instead, effective models of TSC states are constructed with empirical parameters, at the best partially fit to the first-principles results 13 . Meanwhile, there is a parallel development beyond the mean-field approximation employing more realistic number-conserving approach 14,15 , which is yet to be made material specific. Moreover, conventional first-principles approaches that estimate the superconducting transition temperature (T c ) by employing the empirically McMillan's formula 16 or solving the Migdal-Eliashberg formula 17 cannot be applied to the cases involving spin-orbit coupling (SOC) and magnetism (internal or external). Therefore, more versatile and accurate methods to predict T c for SC as well as TSC are highly desirable.
In this article, we attempt to further extend first-principles calculations to the field of TSCs, by developing a versatile approach to quantify realistic conditions of TSC. We construct and solve self-consistently a material-specific first-principles BdG Hamiltonian, based on Wannier function (WF) construction of band structure, in presence of Rashba SOC, Zeeman splitting and electron-phonon coupling (EPC). Furthermore, we demonstrate the usefulness of this method by predicting the Mn-doped GeTe (Ge 1-x Mn x Te) monolayer to be a TSC by constructing a firstprinciples phase diagram in the parameter space of temperature and Mn concentration.
Generally, TSC materials can be classified as intrinsic or extrinsic, depending on the experimental conditions of realizing the non-trivial phase. Intrinsic TSCs exhibit inherently a nontrivial superconducting gap without the need of applying an external field or constructing a heterostructure. They may be pwave SCs with natural spin-triplet pairing 18,19 , such as Sr 2 RuO 4 20 , Cu/Sr/Nb-doped Bi 2 Se 3 21 and non-centrosymmetric SCs 22 , or swave SCs with an effective spin-triplet pairing resulting from helical spin-polarized states, such as the two-dimensional (2D) topological electronic states 23,24 , and 1D 25,26 and 2D Rashba electronic states [27][28][29] which belong to the so-called Class D TSC without time-reversal symmetry (TRS). Extrinsic TSCs employ the same physical mechanisms, but realization of their non-trivial properties requires applying external fields or constructing heterojunctions. To the best of our knowledge, all the known Class D TSCs formed by the s-wave superconductivity are 1 extrinsic, such as the semiconductor nanowire with strong SOC 30 , the ferromagnetic atomic chains 31 , the nanoscale magnetic islands 32 , the ferromagnet 33 , and the topological surface 34 and edge states 35 proximitized with conventional SCs with/without applying external magnetic field. Notably, the signature of TSCs observed by applying an external magnetic field in a superconducting material, e.g., FeTe 0.55 Se 0. 45 36 , epitaxial GeTe 37 and β-Bi 2 Pd thin film 38 indicates the possible existence of intrinsic Class D TSC without needing the external magnetic field, which will further enrich the physics of TSC, in the same perspective as from quantum Hall effect (with magnetic field) to anomalous quantum Hall effect (without).
Given the necessary conditions for realizing Class D TSCs with 2D Rashba electrons [27][28][29] , i.e., inversion symmetry breaking, Zeeman gap opening and superconductivity, the IV-VI compound GeTe with Mn doping, a dilute magnetic semiconductors with a ferromagnetic Curie temperatures T FM c up to~200 K for epitaxial layers on BaF 2 (111) substrate [39][40][41][42][43][44] , caught our attention. The superconductivity of GeTe with p-type doping due to Ge vacancy was confirmed as early as the 1960s 45,46 . It is also known as a ferroelectric material with rhombohedral layered, non-centrosymmetric structure below the ferroelectric Curie temperature of~700 K 47 . Recently, a gradual opening of Zeeman gap in the Rashba bands of GeTe with Mn doping was observed, attributed to the entanglement of ferromagnetic and ferroelectric order 48 . Also, a recent experiment has reported possible signatures of extrinsic TSC in GeTe film under external magnetic field 37 .
Specifically, we focus on the recent experimentally exfoliated GeTe monolayer 49 , which was predicted to be useful in optoelectronic devices and may be a type-II Ising superconductor upon slight hole doping 50,51 . We first show that GeTe monolayer inherits all the key characteristics of its bulk phase by using conventional first-principles calculation. Then, the firstprinciples BdG Hamiltonian was constructed via a WF scheme, through which we found that the GeTe monolayer with the hole concentration of~7.4 × 10 13 cm −2 becomes superconducting below~120 mK and the Ge 1-x Mn x Te monolayer is a Class D TSC with T c~4 0 mK characterized by a non-zero Chern number and chiral Majorana edge modes. A phase diagram of Ge 1-x Mn x Te is constructed by employing the developed first-principles approach to guide experimental detection of the predicted SC and TSC phase. Since both the exfoliated GeTe monolayer 49 and epitaxial Ge 1-x Mn x Te thin film already exist [39][40][41][42][43][44] , our prediction should be readily testable experimentally. Our approach provides a benchmark to make material-specific predictions of TSCs by using first-principles calculations.

RESULTS AND DISCUSSION
Crystal and electronic band structure The crystal structure of GeTe monolayer is shown in Fig. 1a, which is a (111) layer fragment of its bulk phase. Each Ge(Te) atom is bonded with three Te(Ge) atoms, forming a buckled honeycomb lattice. The in-plane lattice constant a and buckling height h was optimized to be~3.955 Å and~1.565 Å, respectively, in good agreement with previous report 50 . Due to the absence of inversion symmetry, a large Rashba splitting arises in the electronic band structure (Fig. 1b). The electronic states are doubly degenerate at the Γ and M points, forming the so-called Kramers pairs, while the degeneracy was lifted away from these time-reversal invariant points. For the four valence bands near the Fermi level of our interest, hereafter we name the lower (upper) two bands as the Rashba (Ising) bands for clarity, referring to their respective electronic spin-texture near the Γ point ( Supplementary Fig. 1).
To predict the TSC formed by 2D Rashba electrons [27][28][29] , we focus on the Rashba bands with a significant Rashba splitting coefficient α R = 0.66-0.76 eV Å. It is comparable with that of heavy metals Au(111) and Bi(111) surface 52,53 , but slightly smaller than that of bulk GeTe 54 . A strong Rashba effect is desirable for the electrons to overcome the suppressing effect of Zeeman field on superconductivity. Doping 0.1 holes per primitive cell, corresponding to a hole concentration of~7.4 × 10 13 cm −2 , will move the Fermi level (E F ) to the Dirac point formed by the Rashba splitting (Fig. 1b). The electronic density of states (DOS) at Fermi level, i.e., N F , is thus increased from 0 to~1.4 states/eV/primitive-cell, which stems mainly from the p-orbitals of Te and Ge atoms (Supplementary Fig. 2). Figure 1c shows the spin-texture on the Fermi Having demonstrated the Rashba spin splitting in the GeTe monolayer, we now discuss the second ingredient, the Zeeman gap. It has been reported that a Zeeman gap can be opened in the bulk Ge 1-x Mn x Te with a ferromagnetic order parallel to the (111) direction 48 , which is the easy magnetization direction for small x 55 . By reproducing the experimental results of bulk Ge 1-x Mn x Te based on the virtual crystal approximation (VCA) 56 , the spin state of Mn dopants was determined to be S = 5/2 (Supplementary Note 1). Consequently, the out-of-plane high-spin state (S = 5/2) of Mn dopants is adopted in Ge 1-x Mn x Te monolayer under VCA. As expected, the Zeeman gaps δ z of Rashba and Ising bands opened at Γ increase monotonically with the increasing Mn concentration (Fig. 1d), and can be fit by the equation of δ z = 250 × x and δ z = 1550 × x meV, respectively. The different slopes are the result of different out-of-plane spin magnitude of the electronic states near the Dirac point versus near the valance band maximum ( Supplementary Fig. 1c).

Superconductivity with and without TRS
We next discuss the phonon-mediated superconductivity of the 0.1-hole-doped GeTe monolayer. From the calculated phonon spectra (Fig. 2a), we first confirm its dynamical stability by the absence of imaginary frequency. For the acoustic branch with the lowest vibration frequency, Kohn anomalies can be seen at certain q-points around Γ, which is favorable for enhancing EPC. Then the EPC strength is evaluated based on the conventional firstprinciples approach (see Methods sections). The calculated EPC strength λ qv of a specific phonon mode v at the wavevector q with the frequency ω qv shows two significant features (Fig. 2a). On one hand, all phonon modes can couple with electrons. This is further confirmed by the comparison between the frequency ω dependent phonon DOS F(ω) and isotropic Eliashberg spectral function where α 2 is the average of electron-phonon interaction (Fig. 2b). Meanwhile, the cumulative EPC strength λ(ω) increases quickly to 1.13 at the frequency of ω~10 meV, which is about 81% of the total EPC constant λ = 1.39. This indicates that the EPC stems mainly from the acoustic modes. The convergence of EPC calculation has been carefully checked (Supplementary Note 2). On the other hand, only the vibration modes with a finite wave vector can couple with electrons. This is because for all the FS contours surrounding Γ (Fig. 1c), only a finite length of phonon wave vectors can connect the initial and final scattering states. In addition, both α 2 F(ω) and λ qv illustrate that the soft modes associated with the Kohn anomalies help to enhance the EPC strength 57 .
To estimate the superconducting transition temperature T c , we construct a material-specific BdG Hamiltonian H BdG ðkÞ in the momentum space by employing the electronic Hamiltonian H WFs k ð Þ: Here the H WFs k ð Þ is obtained by the Fourier transform of the realspace Hamiltonian H ðRÞ WFs , and the latter can be constructed from fitting the first-principles band structure of specific materials using WANNIER90 code under the basis of WFs 58 . Each WF with the orbital index i contains two spin-components, leading to 2@ total WFs. The chemical potential E F in H BdG ðkÞ is the Fermi level where the superconducting gap Δ condenses under the basis vector of φ BdG ¼ ðφ WFs ; φ y WFs Þ T . Only the intra-orbital spin-singlet pairing is considered in the H BdG ðkÞ following the previous theoretical proposals [27][28][29] . Then we formulate the superconducting gap equation into the following form: Here E l;k are eigenvalues of the so constructed H BdG ðkÞ. V, l, k B , and T represent material volume, band index of quasi-particles, Boltzmann constant and temperature, respectively. The intraorbital spin-singlet pairing in the form of Eq. 2 ensures i = j, i.e., Δ ii Δ. The absolute pairing strength g ii is usually identical for the bands with similar orbital character in one specific material 57 , which is calculated as g ii ¼ λ À μ Ã ð Þ=N F g with μ * representing effective Coulomb repulsion 59 . This gap equation enables us to solve the superconducting gap self-consistently at different temperatures. Only the quasi-particle states within one Debye temperature θ D around zero energy, i.e., E l;k θ D , are summed over in the k > 0 half of Brillouin zone (BZ) considering the particle-hole symmetry of H BdG ðkÞ. Details of constructing the BdG Hamiltonian H BdG k ð Þ and formulating the gap equation can be found in the Methods sections.
We emphasize that this method is not only different from the conventional method employing the McMillan's formula 16 or solving the anisotropic Migdal-Eliashberg formula 17 in estimating T c , but also extend the first-principles approach to calculate the topological invariant of superconducting gap and the critical magnetic field/doping-concentration of superconductivity (see below). We check the correctness of Eq. 3 by reducing it to the well-known gap equation for a single-band s-wave SC 60 . Here ε k are eigenvalues of normal electronic state; details are given in the Methods sections. Its reliability is further confirmed by reproducing superconductivity of three representative known SCs, i.e., bulk lead ( Supplementary Fig.  4) 61 , bulk GeTe ( Supplementary Fig. 3d) 46 , and MoS 2 monolayer ( Supplementary Fig. 5)  For the 0.1-hole-doped GeTe monolayer, we assume the Debye temperature (~200 K) and Coulomb repulsion μ * to be same as that of bulk GeTe and extract the WFs using the p orbitals of Ge and Te. Also, we heuristically reduce the calculated total EPC constant λ from 1.39 to~0.76 by~45.5%, based on the benchmark of correlation effect in MoS 2 monolayer 66 . This should set a lower limit on EPC constant since the correlation effect of p-orbitals is usually weaker than that of d-orbitals. The resulting absolute pairing strength g (~0.4) is comparable to that of bulk GeTe (~0.49) 67 , which enables us to predict the superconducting gap Δ of the 0.1-hole-doped GeTe monolayer at different temperatures. From Fig. 2c, one can see the calculated Δ~18.6 μeV for both Rashba and Ising bands, which is gradually suppressed with the increasing temperature. The T c is around~120 mK, lower than that of GeTe film 46 . We anticipate that the predicted 2D superconductivity may be confirmed by growing GeTe monolayer on Si (111) wafers, as the epitaxial GeTe thin film was observed to be superconducting on this substrate 37 .
Next we simulate the superconductivity of Ge 1-x Mn x Te monolayer by adding an out-of-plane Zeeman energy B z in H WFs k ð Þ first: Here σ z is the Pauli matrix in spin space and the I @ @ is a @ @ identity matrix. Then the BdG Hamiltonian H z BdG k ð Þ can be reconstructed through the Eq. 1 and Eq. 2. The reliability of such treatment in simulating the SC without TRS is confirmed by reproducing the in-plane critical magnetic field of MoS 2 monolayer (Supplementary Note 3.2). 68 By diagonalizing the H z WFs k ð Þ with different B z in momentumspace, we obtain the Zeeman gap δ 0 z of Rashba and Ising bands opened at Γ point ( Supplementary Fig. 6a), which can be fit as δ 0 z = 0.122 × B z and δ 0 z = 2.0 × B z meV, respectively. Combining with the δ z fit to the first-principles results in Fig. 1d, one obtains the relationship between B z and Mn concentration, as B z = 2049 × x and B z = 775 × x meV for the Rashba and Ising bands, respectively. The self-consistently calculated T c ( Supplementary  Fig. 6b) and Δ (Supplementary Fig. 6c) domenstrate that they both decrease gradually with the increasing B z due to the pairing breaking effect of magnetism. The superconductivity of Rashba (Ising) bands is fully superssed when B z > 0.35 (0.23) meV, indicating a critical Mn doping concentration of x c~0 .017% (0.03%) (Fig. 2d). This value of x c is two orders of magnitude smaller than that (2%) of Mn doped MgB 2 69 , which is reasonable since the T c of GeTe monolayer is lower than MgB 2 by similar magnitude.
Topological superconductivity and phase diagram To realize TSC formed by 2D Rashba electrons, model analysis proposes that the half of the Zeeman gap opened at the Dirac point of Rashba bands, i.e., δ z /2, should be larger than the superconducting gap Δ 27-29 . In the following, the first-principles approach has been extended to characterize the TSC phase based on a material-specific BdG Hamiltonian H z BdG k ð Þ. Specifically, we take Δ = 0.2 meV and B z = 7.5 meV with δ z~0 .9 meV to construct the H z BdG k ð Þ of Ge 1-x Mn x Te monolayer via Eq. 1, Eq. 2 and Eq. 4. The relatively large B z and Δ are used to show the topological non-triviality more clearly. The H z BdG k ð Þ is analogous to the singleparticle Hamiltonian of electrons with an energy gap mathematically. By diagonalizing H z BdG k ð Þ in the momentum space, we obtain the dispersion relation of superconducting quasi-particles (Fig. 3a). One can clearly see that the superconducting gap is indeed opened, where the topological invariant, i.e., first Chern number (N c ), is well-defined.
For 2D systems, the Chern number of l-th band is calculated by integrating the Berry curvature Ω l ðkÞ ¼ ∇ A l ðkÞ over the first BZ: where A l ðkÞ is Berry connection. The total Chern number N c can be obtained by summing up the Chern numbers of all the states below the superconducting gap, which is quantized to −1. The Berry curvature resides mainly at the Γ point associated with the Zeeman gap opening (Fig. 3b), similar to the band inversion in the quantum anomalous Hall systems. Here we should emphasize that N c does not physically correspond to a quantized Hall conductance because charge is not conserved in the BdG Hamiltonian 24 . Two chiral Majorana edge modes localized at two different edges clearly exist in the continuous superconducting gap due to the bulk-boundary correspondence (Fig. 3c and 3d). The propagation of chiral Majorana fermions could lead to same unitary transformation as that in braiding Majorana zero modes 70 , and the deterministic creation and braiding of chiral edge vortices in hybrid structures were elaborated 71 .
We finally construct a phase diagram of the 0.1-hole-doped Ge 1-x Mn x Te monolayer in Fig. 4, to help guide future experimental detection of the predicted TSC phase formed by the superconducting Rashba bands. At the zero-temperature limit, the SC phase of the Rashba bands will be preserved for x < x c = 0.017% and the TSC phase will arise when x > x min = 0.014%, where the pre-condition of δ z /2 > Δ can be met ( Supplementary  Fig. 6d). At finite temperature, both the ferromagnetic and SC order should exist simultaneously for the formation of TSC phase. Referring to the ferromagnetic Curie temperatures T FM c of Ge 1-x Mn x Te that increases linearly with increasing Mn concentration up to x = 0.2 and can be fit by T FM c x ð Þ ¼ 333 x K ( Supplementary Fig. 7a) [39][40][41][42][43][44] , we estimate T FM c x ð Þ and xdependent T c will cross over at x 0 min ¼ 0:014% ¼ x min (Supplementary Fig. 7b), too. Consequently, the TSC phase can be X. Zhang et al. formed for x min < x < x c at the temperature where the SC order occurs. Given that the non-trivial phase can be readily realized in non-centrosymmetric SCs when p-wave component is stronger than s-wave 22 , our results indicate the TSC phase of Ge 1-x Mn x Te monolayer could be robust against parity mixing of Cooper pairs. We suggest preparing the desired Ge 1-x Mn x Te monolayer on BaF 2 (111) substrate 39-44 by molecular beam epitaxy since the growth is known to start in a 2D manner 41 . We anticipate that the chiral Majorana edge modes of Ge 1-x Mn x Te monolayer can be detected using Josephson effect 5 or charge transport 72 , and controlled by magnetic flux 73 . The effects of magnetic anisotropy and GeTe film thickness on the TSC phase are discussed in Supplementary Note 4 and Note 5.
Lastly, in addition to the monolayer Ge 1-x Mn x Te we demonstrated here, we suggest two more candidate materials for Class D TSC. Firstly, since the heterostructures of MnBi 2 Te 4 /Bi 2 Te 3 74,75 and Bi 2 Te 3 /NbSe 2 34,76 have already been fabricated, the MnBi 2 Te 4 /Bi 2 Te 3 /NbSe 2 hold high possibility to be synthesized. We demonstrate that this type of heterostructures with magnetized topological surface states are also a Class D TSC characterized with non-zero Chern number (Supplementary Note 6) 77 . Secondly, it was experimentally reported the desired Rashba-Zeeman splitting can be alternatively achieved by the magnetic order in the Si-terminated surface of HoRh 2 Si 2 78 . With a spin-singlet Cooper pairing being tunneled into this surface state by superconducting proximity effect, the Class D TSC will be readily emerging. By applying our developed first-principles BdG Hamiltonian approach, a complete phase diagram of these systems can be constructed in the near future.

Details of the first-principles calculations
The Vienna ab initio simulation pack 79,80 was utilized to calculate the electronic property of normal states based on the density-functional theory. The exchange-correlation of electrons was treated within the generalized gradient approximation in the form of Perdew-Burke-Ernzerhof 81 . The atomic structures of GeTe monolayer and thin film was set up by introducing a vacuum region of more than 15 Å to avoid the interactions between neighboring images. Structural relaxations and self-consistent calculations as well as the Zeeman gap calculations were performed on a uniform 30 × 30 × 1 (18 × 18 × 18) k-point sampling of the first BZ for monolayer (bulk) GeTe. The energy cutoff was set to 400 eV for plane-wave basis. The dipole correction was used to cancel the artificial electric field imposed by the periodic boundary condition of GeTe thin film.
The QUANTUM ESPRESSO package 82 was used to calculate the phonon spectra and EPC strength based on the density-functional perturbation theory 83 as well as fit the first-principles band structure by interfacing with the WANNIER90-2.1 code 58 . The Optimized Norm-Conserving Vanderbilt Pseudopotential 84 was employed and the kinetic energy cutoff was set to 100 Ry for wave functions. The hole doping was simulated by removing electrons from intrinsic GeTe monolayer and introducing the compensating jellium background to avoid divergence. The dynamic matrix and phonon frequency are computed on a 18 × 18 × 1 q-point mesh with a 18 × 18 × 1 k-point sampling, and a finer 36 × 36 × 1 k-point grid is used in the EPC calculations, where the DOS is converged (Supplementary Fig. 2b). Other q/k-point samplings (Supplementary Table 1) are also employed to check the convergence of EPC calculations. The phonon DOS F(ω) and the isotropic Eliashberg spectral function α 2 F(ω) as well as the cumulative frequency-dependent EPC strength λ(ω) are calculated using a 60 × 60 × 1 q-point sampling by means of the Fourier interpolation.  Specifically, the q-and v-resolved EPC strength λ qv is given by: where N F is the electronic DOS at the Fermi level, W k is the weight of wavevector k, ε nk is the eigenvalue for electronic wavefunction ψ nk with band index n and wavevector k, ω qv is the frequency of a phonon branch v at wavevector q, h is the reduced Planck constant, and M 0 is the ionic mass. g mn;v k; q ð Þ represents the scattering amplitude between the electronic states with wavefunction ψ nk and ψ mkþq , induced by derivative ∂ qv Ξ À Á of the self-consistent potential associated with phonon ω qv . δ is the Dirac delta function. The frequency ω dependent isotropic Eliashberg spectral function α 2 F(ω) and the cumulative EPC strength λ(ω) are then calculated from: λðωÞ ¼ 2 Here W q is the weight of wavevector q. The total EPC constant λ equals to λ(ω max ) with ω max being the maximum of phonon frequency.

Constructing the BdG Hamiltonian
To perform first-principles prediction of TSC, the main challenge is to construct a BdG Hamiltonian of superconducting quasi-particles from the electronic Hamiltonian of normal state. Here we propose a strategy to overcome this obstacle by employing the WFs φ WFs ¼ Here i is the orbital index, and the total number of WFs is 2@. and τ x is the Pauli matrix in particle-hole space.
For the materials with external/internal magnetism, we first add a Zeeman term B in H WFs ðkÞ using the vector of Pauli matrix σ in spin space: Then the first-principles BdG Hamiltonian H B BdG k ð Þ without TRS can be constructed for a specific material by using the above procedure. We should emphasize that the above construction procedure is practically also applicable to materials with revelatively strong SOC. It is noted that the WFs, rather than the maximally localized WFs, are obtained by the WANNIER90 code without minimization procedure, so that the resulting WFs can be approximately seperated into up (majority) and down (minority) pseudospin orbitals in the presence of SOC. This treatment has been widely adopted in the WF-based methods for investigation of topological materials, such as WannierTools 85 . Here we develop another route to quantifying realistic conditions of TSC by solving self-consistently BdG equation based on WFs.

Formulating the gap equation
Under the basis of Ψ k ¼ c 1;k ; c 2;k ; Á Á Á ; c i;k ; Á Á Á ; c y 1;Àk ; c y 2;Àk ; Á Á Á ; c y i;Àk ; Á Á Á T , the multi-band Hamiltonian with s-wave pairing can be written as: H BdG ðkÞ ¼ hðkÞ ÀΔ Δ Àh Ã ðÀkÞ With the relation of ∂H ∂Δij ¼ À P k c y i;k c y j;Àk þ c j;Àk c i;k , we can derive the gap equation of Δ ij as:  (12) Solving this gap equation self-consistently enables us to estimates the superconducting transition temperature T c and the critical magnetic field/ doping-concentration of specific materials based on the material-specific BdG Hamiltonian H BdG ðkÞ with/without TRS constructed from the state-ofart first-principles approach.
To confirm the above derivation, we apply the derived gap equation to the single-band Hamiltonian with s-wave pairing: Its eigenvalues are when the eigenvalues of normal electronic state satisfy ε k ¼ ε Àk . Substituting the four eigenvalues into the derived gap equation, we obtain the well-known gap equation of single-band s-wave superconductor: 60 1

DATA AVAILABILITY
All data needed to evaluate the conclusions of this paper are available within the paper and Supplementary Information.