Proximity effect in a superconductor-topological insulator heterostructure based on first principles

Superconductor-topological insulator (SC-TI) heterostructures were proposed to be a possible platform to realize and control Majorana zero-modes. Despite experimental signatures indicating their existence, univocal interpretation of the observed features demands theories including realistic electronic structures. To achieve this, we solve the Kohn-Sham-Dirac-Bogoliubov-de Gennes equations for ultrathin Bi$_2$Se$_3$ films on superconductor PdTe, within the fully relativistic Korringa-Kohn-Rostoker method, and investigate quasiparticle spectra as a function of chemical potential and film thickness. We find a strongly momentum-dependent proximity-induced gap feature where the gap sizes highly depend on characteristics of the TI states. The interface TI Dirac state is relevant to the induced gap only when the chemical potential is close to the Dirac-point energy. Otherwise, at a given chemical potential, the largest induced gap arises from the highest-energy quantum-well states, whereas the smallest gap arises from the TI topological surface state with its gap size depending on the TI pairing potential.

In SC-TI heterostructures, it was predicted that the topological TI interface state creates proximityinduced topological superconductivity where the edge state or the vortex lattice can host Majorana zeromodes. Keeping this in mind, experiments were performed on Bi 2 Se 3 (111) or Bi 2 Te 3 (111) films (< 10 nm) grown on NbSe 2 or Nb substrates, using scanning tunneling miscroscopy/spectroscopy (STM/S) [5,7,9] and angle-resolved photoemission spectra (ARPES) [6,10]. A proximity-induced gap was measured when the Fermi level lies in the TI conduction band. The proximityinduced gap was shown to decrease with increasing TI film thickness. However, both characteristics of the gap and its dependence on film thickness remain elusive and inconsistent among different groups. The induced surface gap was reported to reach 15-40% of the bulk SC gap [5,7,10] even for thick TI films where the interface and top-surface Dirac states do not hybridize. Intriguingly, transport experiments on Bi 2 Se 3 with a Pb overlayer showed zero resistivity and a proximity-induced gap about 1 µm away from the interface [4].
Theoretical efforts have been made to understand the proximity effect in SC-TI heterostructures, using the Bogoliubov-de Gennes Hamiltonian for SC and effective models for TI [2,3,[20][21][22][23]. The TI model Hamiltonian ranges from the Fu-Kane-Mele model [24] to the k · p model [25] with parameter values from calculated band structures using density-functional theory (DFT). For chemical potential away from the Dirac point, triplet pairing was predicted [2,20], and proximity-induced pairing types were classified based on symmetries [22]. Despite this success, the current theoretical approaches have many limitations due to the lack of information on realistic band structures and interface effects as well as insufficient treatment of quasi-two-dimensional experimental systems. The reported theoretical studies either neglected quantum-well states (QWS) or included generic QWS, although the experimental Fermi level typically lies in the TI conduction band, and they oversimplified the Fermi surfaces despite their importance in the induced gap features.
We investigate quasiparticle spectra of heterostructures comprised of thin Bi 2 Se 3 overlayers on a s-wave SC PdTe substrate, by solving the Dirac-Bogoliubov-de Gennes (DBdG) equations [26] using the fully relativistic screened Korringa-Kohn-Rostoker (SKKR) Green's function method [27,28] within DFT. We calculate the band structure of Andreev bound states and determine proximity-induced gap features in the TI layers, consid-ering two different TI pairing potentials, as the overlayer thickness and chemical potential are varied. Several different types of induced gaps appear depending on the characteristics of the TI states for a given overlayer thickness and chemical potential. The induced gap arising from the interface TI Dirac state is the largest, while the gap from the top-surface TI state is smallest, independently of the thickness, chemical potential, and TI pairing potential. The induced gap associated with a given QWS type increases as the chemical potential decreases. The induced gap sizes do not have a spatial dependence except for the top-surface TI Dirac state for thick TI films. The induced gap size from the top-surface TI Dirac state is highly susceptible to the TI pairing potential, although that is not the case for the induced gap from the interface TI Dirac state.
We first discuss our systems of interest in Sec. II and present our first-principles calculations of the electronic structure of the heterostructures in normal state as a function of overlayer thickness in Sec.III.A. We then show our calculated spectral functions of the heterostructures in SC state with two different TI pairing potentials as a function of overlayer thickness and chemical potential in Sec.III.B. We make brief comparison of our results with relevant experimental data in Sec.III.C and make a conclusion in Sec.IV.

II. SYSTEMS OF INTEREST
To incorporate realistic electronic structures in the superconducting proximity effect, we simulate Bi 2 Se 3 films of 1-6 quintuple layers (QLs) overlaid on SC PdTe [ Fig. 1(a)] within the fully relativistic SKKR method [27,28]. We consider such heterostructures due to the following advantages: (i) good lattice match at the interface; (ii) reasonable SC transition temperature of PdTe; (iii) a single Dirac cone of Bi 2 Se 3 at a given surface [25]; (iv) experimental data on SC-Bi 2 Se 3 heterostructures (despite different SC substrates) [4][5][6][7]10]. Bulk Bi 2 Se 3 has a rhombohedral structure (space group 166, R3m) with experimental lattice constants a = 4.143 and c = 28.636Å [29]. As shown in Fig. 1, Se-Bi-Se-Bi-Se atomic layers form one QL (about 1 nm) along the [111] direction, and individual QLs are bonded via weak van der Waals interaction. The (111) surface has a hexagonal in-plane lattice with a lattice constant a = 4.143Å. Bulk PdTe has a NiAs-type hexagonal structure (space group 194, P 6 3 /mmc) with experimental lattice constants a = 4.152 and c = 5.671Å [30]. When a Se-terminated Bi 2 Se 3 (111) film is interfaced with a Te-terminated PdTe(001) substrate, the in-plane lattice mismatch between Bi 2 Se 3 and PdTe is about 0.2%. PdTe is a type-II SC with critical temperature of 4.5 K [30] and its SC gap ∆ PdTe is about 0.71 meV at zero temperature [30]. The electron-phonon coupling of PdTe is 1.4 [30], whereas the reported electron-phonon coupling of Bi 2 Se 3 has an upper bound of 0.43 [31].
As required in the SKKR formalism, the SC-TI heterostructures are divided into three regions I-III ( Fig. 1(a)). Region I and III are semi-infinite PdTe and vacuum layers, respectively. Region II consists of Pd-Te-Pd-Te atomic layers overlaid with Bi 2 Se 3 layers (1-6 QLs) and 4-6 vacuum layers on top. The heterostructures have two-dimensional translational symmetry with the lattice constant of Bi 2 Se 3 , 4.143Å. Our choice of the in-plane lattice constant is made to avoid an effect of strain on the topological surface states [32,33]. The outof-plane lattice constant for PdTe is slightly expanded to conserve the volume of PdTe. Otherwise, experimental lattice constants are used. We consider the band structure only along the Γ − K direction, k x axis, in this work.

III. RESULTS AND DISCUSSION
In multiple scattering theory the band structure is obtained from site-dependent Bloch spectral functions where i denotes the i-th atomic site and E n (k ) is the n-th band energy at in-plane momentum k . We calculate the BSF for the i-th site (located at r i ) from its retarded Green's function G + i (E, r i , k ): Density of states (DOS) within the SKKR method is obtained from an integral of the BSF over k , i.e., . A detailed description of solving the DBdG equations within the fully relativistic SKKR method can be found in Refs.27, 28. A brief method description with parameter values for our simulations is shown in the Supplemental Material (SM).

A. Electronic structure of normal state
First of all, we present calculated electronic structures of PdTe and Bi 2 Se 3 in the normal state, separately.  Figure 1(d) shows BSF contours of a semi-infinite Bi 2 Se 3 system. We find Dirac surface states within a bulk band gap of about 0.03 Ry as well as continuous conduction and valence band regions. For Bi 2 Se 3 slabs, all bands are doubly degenerate due to time-reversal and inversion symmetries, and several QWS appear in the conduction and valence band regions. For a N -QL slab, N − 1 (N − 2) QWS appear in the conduction (valence) band region [34], as shown in Fig. S1(b),(d) in the SM. For Bi 2 Se 3 slabs thinner than 5-6 QLs, the top and bottom surface states hybridize, opening an energy gap in the Dirac surface states [34,35]. This gap is referred to as a surface-hybridization gap. The SKKR-calculated band structures of PdTe and Bi 2 Se 3 agree with those using VASP code [36,37] (Fig. S1 in the SM).
Using the above converged potentials of PdTe and Bi 2 Se 3 , we perform fully relativistic SKKR calculations on the Bi 2 Se 3 -PdTe heterostructure in the normal state. Figure 2(a)-(f) show calculated normal-state BSF of the heterostructure with the TI film thickness varying from 1 to 6 QLs, where all layers are summed. The gray continuous spectrum in the BSF is similar to that of PdTe. Compare Fig. 2(a)-(f) with Fig. 1(b) or 2(i). We find that a shift of the Madelung potential lowers the TI Dirac point around −0.29 Ry and that the top-surface and interface Dirac states are shifted from each other with strong modification of the dispersion of the interface Dirac states. The slope of the dispersion near the Dirac point is substantially reduced and the states lose the interface-state character somewhat away from Γ. See Fig. 2(g) and (h). The top-surface (interface) Dirac states are identified as states with a large BSF weight onto the topmost (interface) QL. Strong hybridization of the interface states with the substrate causes the strong modification of their dispersion. Similar effects have been reported in various heterostructures involving Bi 2 Se 3 [38,39]. For thin films (< 5 QLs) we also observe an energy gap in the vicinity of the two Dirac points which decreases with increasing the TI overlayer thickness. This gap is induced by the hybridization between the interface and top-surface Dirac states. The shape and number of TI QWS, however, remain unchanged with the substrate, but they are quite broadened compared to the Dirac states, as shown in Fig. 2(b)-(f). This broadening may be caused by scattering of electrons from the substrate. Calculated BSF with finer resolution suggests that each broad QWS peak in Fig. 2  Let us now consider that PdTe is in the SC state with a bulk SC gap ∆ PdTe of 0.001 Ry and that the TI pairing potential ∆ TI eff is zero. We then calculate BSF of the heterostructures at a fixed chemical potential. See the SM for the detailed procedure. Figure 3 shows BSF contours of the 5-QL TI overlayer on SC PdTe as a function of energy and k x at chemical potential µ 4 = −0.254 Ry. Within the bulk SC gap, we find a highly momentumdependent proximity-induced gap feature, ∆ k . In partic-ular, we focus on seven induced gap sizes listed in Table I. Two large gaps ∆ k1 and ∆ k2 are clearly seen in Fig. 3(a)-(c). Gap sizes ∆ k3 , ∆ k4 , ∆ k5 , ∆ k6 , and ∆ k7 are zoomed in in Fig. 3(d)-(f). Except for ∆ k7 , all six gap sizes are observed in the BSF of all TI layers. For example, see Fig. 4(c) in the case of ∆ k1 . Thus, these gap sizes do not depend on z for a fixed TI-overlayer thickness. The induced gap size is larger with smaller k x . The dispersion near the six gap sizes arises from both electron and hole contributions. See Figs. 3(g),(h),(j),(k), S3, and S4. In the case of ∆ k7 , however, the dispersion appears only for the topmost QL and it exhibits a positive (negative) slope only from the electron (hole) contribution. See Fig. 3(i) and (l).
The dispersion near the first six induced gaps suggests that Cooper pairs of PdTe tunnel into the TI region, giving rise to the proximity-induced gaps and Andreev bound states. The momentum-dependence and gap sizes can be strongly modified with the interface type and the band structure of non-SC. By comparing the SC-state BSF to the normal-state BSF, we identify that ∆ k1 and ∆ k2 originate from Cooper-pair tunneling into QWS3 and that ∆ k3 and ∆ k4 from QWS2. Similarly, ∆ k5 and ∆ k6 arise from Cooper-pair tunneling into QWS1, and ∆ k7 from the top-surface Dirac state. See also Fig. S5. To elucidate the origin of the two different gap sizes from each QWS, we plot normalized SC-state and normal-state BSF at fixed energies and k x points near the gap size ∆ k1 . As shown in Fig. 4(a), the two SCstate BSF qualitatively differ from each other, while the two normal-state BSF are indistinguishable. The state with a larger spectral weight near the interface gives rise to a larger induced gap. Figure 4(c) also clearly shows the two induced gaps at k 1 . This two-gap feature is consistent with that of the normal-state BSF (Fig. S6). The induced gap size overall increases as Cooper pairs tunnel into higher-energy QWS because of stronger coupling with the substrate. See Fig. 4(b) and Table 1. The dispersion near ∆ k7 (Fig. 3(i) and (l)) is quite distinct from that near the other gaps. Cooper pairs do not seem to efficiently tunnel into the top-surface Dirac state which is strongly localized at the top surface, giving rise to zero induced gap.
We now investigate the effect of chemical potential on the proximity-induced gap at the 5-QL overlayer by considering µ 3 , µ 2 , and µ 1 , as indicated in Fig. 2(g). For µ 3 five distinct gap sizes are found. Similarly to the case of µ 4 , the smallest gap is zero from the top-surface Dirac state, while the rest four gaps are as significant as those for µ 4 . The dispersion near the four gaps has characteristics of Andreev states. As chemical potential decreases, the number of gaps decreases, whereas the gap size associated with a given QWS type (such as QWS1, QWS2) increases. See Fig. S7. For µ 2 and µ 1 , only two gap sizes are observed, and Andreev states exist only near the larger gap. Figure 4(d) and (e) show layer-resolved BSF with µ 2 , and Fig. 4(f) and (g) show BSF of the interface and topmost layers with µ 1 . The two gaps appear only closer to the interface QL or the topmost QL. The highly suppressed slope of the dispersion does not allow Cooper pairs to tunnel into the interface TI state until chemical potential reaches near the Dirac points. For µ 2 and µ 1 , only top-surface and interface TI states are involved. The induced gap from the interface TI state is largest among the gaps found.
To examine the effect of TI film thickness, we plot the smaller induced gap from each TI-state type (i.e., QWS1, QWS2, QWS3, interface or top-surface) vs thickness at µ 4 [ Fig. 4(j)]. We find that the induced gap from higherenergy QWS decays much slowly than that from lowerenergy QWS, and that the gap from the top-surface Dirac state becomes zero for thickness greater than 3 QLs. Figure 4(k) shows the "maximum" gap from the gap sizes shown in Fig. 4(j) as a function of overlayer thickness, at µ 3 and µ 4 . Interestingly, the "maximum" gap oscillates with thickness. For overlayers thinner than 5 QLs, the surface-hybridization effect is also seen in the induced gap. The gap size from the top-surface Dirac state becomes noticeable, and the induced gap sizes do not depend on z. Compare Fig. 4(h) and (i) with Fig. 4(f) and (g) for µ 1 . Tables in the SM list all induced gap sizes.

Induced spectral gap with ∆ TI eff = 0
We study an effect of the TI pairing potential ∆ TI eff by keeping ∆ PdTe same as before and considering that ∆ TI eff = 0. Although electron-phonon coupling of the TI cannot induce superconductivity itself, it can create a finite pairing potential induced by the superconductor [2]. Figure 5(a) shows effective pairing potential ∆ eff calculated self-consistently for the 5-QL overlayer on PdTe using semi-phenomenological parameters within the SKKR method. See the SM for the detailed procedure. As shown in the inset, the pairing potential decays very slowly in the TI region. With this calculated pairing potential, we obtain BSF for the 5-QL overlayer on SC PdTe at µ 4 . As listed in Table II, the induced gap from the top-surface Dirac state becomes noticeable with Andreev-state characteristics (Fig. 5(b)-(d)), whereas the gap from high-energy QWS does not change much compared to the case of zero TI pairing potential. We expect that topological edge states within the smallest induced gap can trap Majorana zero-modes with broken time-reversal symmetry [2,3] and that pairing type of the proximity-induced superconductivity may be identified from gap anisotropy [3] and/or unique spin-orbital textures of the Andreev states [40].

C. Comparison with experiment
Let us make brief qualitative comparison of our results with experimental and model-Hamiltonian studies, considering that our SC substrate differs from those used in experiments. Refs. 6, 10 showed momentum- is set to zero or not. The discrepancy may originate from the experimental difficulty in identifying the peak center of the spectral function corresponding to the bulk and surface states. We find that an accurate estimate of the induced gap from the TI surface state requires a high precision in momentum such as 10 −5 π/a. With a precision of 0.001 π/a, the induced gap can be overestimated by two orders of magnitude. However, the induced gap from bulk states is insensitive to the precision of momentum. Compared to STM experimental data [5,7], we obtain underestimated induced gaps even for ∆ TI eff = 0. Model-Hamiltonian studies predicted the zero gap from the top-surface TI state for thick TI overlayers [20,23], which is consistent with our result. However, the rich features of the induced gaps that we find were not obtained from the model-Hamiltonian approaches.

IV. CONCLUSION
In summary, we simulated Bi 2 Se 3 films of 1-6 nm overlaid on SC PdTe by solving the DBdG equations with two pairing potential profiles within the SKKR method, finding that the strong k -dependence of the proximityinduced gap arises from the unique TI band structure and its modifications under the SC substrate. The size of the induced gap strongly varies with characteristics of TI states which are partially occupied at the Fermi level. Cooper pairs tunnel into higher-energy QWS more efficiently and the induced gap from higher-energy QWS decreases more slowly with increasing TI film thickness. For thick TI films, the induced gap from the top-surface Dirac state becomes zero with zero TI pairing potential, whereas the gap can be substantial for finite TI pairing potential. For a given thick TI film, the induced gap from the interface Dirac state appears only near the of (c) near the induced gaps ∆ k3 , ∆ k4 , ∆ k5 , ∆ k6 , and ∆ k7 . In (a)-(f) both electron and hole contributions are included. (g),(j) Electron and hole BSF contours of (a). (h),(k) Electron and hole BSF contours of (c). (i),(l) Electron and hole BSF contours of (f). The smallest induced gap ∆ k7 is less than 0.1% of the bulk SC gap, as shown in (f).
interface QL, although the induced gap size from each QWS does not depend on z. Our findings demonstrate an importance of consideration of the realistic TI band structure for studies of the superconducting proximity effect, and they can be used for future larger-scale simulations and experimental studies of topological edge states or nanowires on SC substrates in pursuit of Majorana zero-modes.

Acknowledgments
We are grateful to Peter Rakyta for setting up initial Bi 2 Se 3 structures. The computational support was provided by San Diego Supercomputer Center under DMR060009N and Virginia Tech Advanced Research Computing. BU was supported by the Hungarian National Research, Development and Innovation Office under Contract No. K131938 and BME Nanotechnology FIKP grants. GC acknowledges support from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 754510.  4: (a) Normalized SC-state and normal-state BSF vs z at fixed energies and kx (near k1). The black (red) curve here corresponds to the black (red) BSF peak at −0.00032 (−0.00052) Ry in Fig. 3(a). (b) Normalized SC-state BSF vs z at fixed kx and energies with µ4 near ∆ k3 , ∆ k4 , ∆ k5 , ∆ k6 , and ∆ k7 . (c)-(e) Layer-resolved SC-state BSF vs energy at k1 with µ4 (Table 1), and at kx = 0.01369 and 0.03700 π/a with µ2 for 5-QL/PdTe, respectively. The BSF of each layer is shifted and scaled.