Topological superconductor candidates PdBi₂Te₄ and PdBi₂Te₅ from a generic ab initio strategy

As one of the most important systems in both fundamental physics and topological quantum computation, topological superconductors (TSCs) have attracted increasing interest for their ability to support Majorana fermions and anyons with non-Abelian statistics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Currently, the search for TSCs candidates has been focused on two experimental schemes. One is the architectures by the combination of conventional superconductors with topological insulators (TIs) [18, 19, 20] or 1D nanowires [21, 22], but this approach brings high requirements for sample fabrication and interface engineering. The other route is to achieve TSCs in superconducting topological metals (SCTMs) that host both topological electronic structures at the Fermi level and superconductivity in one compound [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36], in which the topological surface states are gapped by the “self-proximity effect” of bulk superconductivity, thus avoiding the complications of interface engineering. This approach has successfully predicted the SCTM FeTe_0.55Se_0.45 [24, 25, 26] and similar compounds of iron-based superconductors [27, 28, 29, 30], owing to the favorable SC gap and non-trivial band topology. Beside the MZMs, 1D helical/chiral Majorana states have also been reported in domain walls of FeTe_0.55Se_0.45 [37] and the magnetism-superconductor heterostructures [38, 39, 40, 41, 42, 43, 44, 45]. It is also proposed that the propagating chiral Majorana states can be applied to realize non-Abelian quantum gate operations, which could be $10^{3}$ faster than the currently existing quantum computation schemes [46].

Encouraged by the success of Fe(Se,Te) [24, 25, 26], many topological materials that host both superconductivity and topological electronic structures are proposed  [47, 48, 49, 50, 51, 52]. However, very rare experimental progress of TSC has been made in such SCTMs. This is because, on the one hand, all of them are not the ideal SCTMs, whose band structures are too complicated, the topological surface states are usually buried in the bulk states and difficult to form the pairing required by TSC. On the other hand, lacking a direct characterization of the TSC properties from ab initio calculations also hinders the effective experimental search in such materials. Therefore, a general program that could calculate the TSC invariant from first-principles calculations is highly desirable.

In this work, we develop a program to characterize the superconducting topological invariant of 2D system from ab initio calculations. Besides, we also propose a new strategy to design ideal SCTMs by intercalating superconducting units into topological insulators. Following this strategy, PdBi₂Te₅ and PdBi₂Te₄ are found to be ideal SCTMs that host topological surface states at the Fermi level and superconductivity at $0.57$ K and $3.11$ K respectively. By performing the superconducting energy spectrum and topological invariant calculations, we identify that chiral TSC could be realized in the slab of such SCTMs by introducing considerable Zeeman splitting on the topological surface states, which can be realized by an external magnetic field or in proximity to FM insulators. Our strategy provides a new framework to enrich SCTMs and TSC candidates, and the program makes it possible to design and modulate the TSC system from ab initio calculations, which can also be extended to study the TSC properties in other system, such as magnetic TI/SC heterostructure, SC/FM heterostructure and SC/TI/SC heterostructure.

Inspired by the construction of magnetic TI MnBi₂Te₄ [53, 54], we propose that the SCTMs can be designed by intercalating the SC units into the TIs, as illustrated by the schematic of Fig. 1(a). As an ideal SCTM, the target crystal should be relatively stable in both energy and structure. More importantly, it must inherit the topological electronic structures of the parent TI near the Fermi level, and also the superconductivity of the parent SC as shown in Fig. 1(b). However, the combination of topological electronic structures and SC does not result in TSC eventually. The realization of TSC generally requires a delicate modulation of many parameters, such as SC pairing, Zeeman splitting and chemical potential, et al  [18, 19, 20, 21, 22, 23, 24, 25, 26, 43, 44, 45, 38, 39, 40, 41, 42]. Thus, the ability to characterize the TSC invariant and determine the required parameters in real materials from the ab initio calculations is not only of theoretical significance, but also highly desirable in experiment.

Here we develop a program to simulate the superconducting properties and characterize its topological invariant in 2D slab system from ab initio calculations, in which the necessary ingredients to realize chiral TSC based on SCTM are included, such as bulk band structures, SC pairing, Zeeman splitting, Rashba spin-orbit coupling and chemical potential. The workflow of this program is shown in Fig. 2. First, one should calculate the electronic structures of SCTM materials, and construct the localized Wannier functions that capture all electronic features from the first-principles calculations, referred as $\hat{H}_{\mathrm{bulk}}$. The next step is to construct the slab Hamiltonian $\hat{H}_{\mathrm{slab}}$ with open boundary condition along a certain direction [55]. In general, the spin-orbit coupling (SOC) and surface effect can be included automatically in $\hat{H}_{\mathrm{slab}}$ through the first-principles calculations with SOC. So that the topological properties, such as the surface states and spin-texture, can be directly studied by using $\hat{H}_{\mathrm{slab}}$. On the other hand, one can also construct a slab Hamiltonian $\hat{H}^{\mathrm{nsoc}}_{\mathrm{slab}}$ that excluded SOC from the non-SOC first-principles calculations, and add $\hat{H}_{\mathrm{SOC}}$ and $\hat{H}_{\mathrm{surf}}$ manually to simulate the variable SOC and surface effect in the topological electronic states and TSC. In this work, we will adopt the former type of $\hat{H}_{\mathrm{slab}}$, in which only the intrinsic SOC of the real material is included. With adopting particle-hole transformation, the $\hat{H}_{\mathrm{slab}}$ can be extended to BdG Hamiltonian $\hat{H}^{\mathrm{BdG}}_{\mathrm{slab}}$ by adding SC pairing $\hat{H}_{sc}$ and Zeeman splitting $\hat{H}_{\mathrm{z}}$. In the Nambu basis $\Phi_{\mathbf{k}}=(c_{\mathbf{k},j,\alpha,\uparrow},c_{\mathbf{k},j,\alpha,\downarrow},c^{{\dagger}}_{-\mathbf{k},j,\alpha,\uparrow},c^{{\dagger}}_{-\mathbf{k},j,\alpha,\downarrow})$, where the $c_{j,\alpha,\sigma}$ is the fermion operator denotes an electron at $j$ layer with orbital $\alpha$ and spin $\sigma(\uparrow,\downarrow)$, the BdG Hamiltonian is formulated as:

In Eq. 1, $\mu$ is the chemical potential, which can be used to simulate the carriers doping. $\Delta(\mathbf{k})$ denotes the SC pairing matrices, which could be both singlet and triplet pairing form. For the conventional $s$-wave SC, $\Delta(\mathbf{k})$ is expressed as:

where $\Delta_{s}$ is the magnitude of intrinsic bulk $s$-wave pairing, $\sigma_{y}$ is the Pauli matrix in spin space, $I_{slab}$ ($I_{orb}$) is an $\mathcal{N}_{slab}\times\mathcal{N}_{slab}$ ($\mathcal{N}_{orb}\times\mathcal{N}_{orb}$) identity matrix that represents the number of slab layers (Wannier orbitals). $\tau_{z}$ is the Pauli matrix in particle-hole space, $M_{z}$ is the Zeeman splitting energy, and $H_{z}=M_{z}\tau_{z}$ is used to simulate the influence of the external magnetic field or the proximity effect of the FM insulator. Thus, $H_{z}$ can be chosen to be applied for the whole slab or just few surface layers, depending on the slab thickness, strength of magnetic field, the type of the SC et al. In principle, chiral TSC can be achieved by modulating the SC pairing, Zeeman splitting and chemical potential [43, 44, 45, 38, 39, 40, 41, 42], which can be further revealed by calculating the superconducting energy spectrum and the superconducting topological invariant.

In the gaped 2D superconducting system, the topological superconductors are classified by BdG Chern number in the absence of time-reversal symmetry [3]. Such superconducting topological invariants can be characterized by the evolution of Wilson loop [56, 57, 58]. For the occupied quasiparticle states $|u^{\textrm{BdG}}_{n,k_{1},k_{2}}\rangle$, where $k_{1}$ and $k_{2}$ are momenta along two primitive vectors of the Brillouin zone (BZ), the Berry phase of the Wilson loop along $k_{2}$ at a fixed $k_{1}$ can be expressed as:

with the overlap matrix $M^{(i)}_{k_{1},mn}=\langle u^{\textrm{BdG}}_{m,k_{1},k^{(i)}_{2}}|u^{\textrm{BdG}}_{n,k_{1},k^{(i+1)}_{2}}\rangle$, where $k^{(i)}_{2}$ is the $i$-th discretized momenta along $k_{2}$ direction. The winding number of $\mathcal{W}(k_{1})$ with respect to $k_{1}$ is equal to the superconducting BdG Chern number $C_{\textrm{BdG}}$.

Next, we take TI Bi₂Te₃ [59, 60], SC PdTe [61, 62] and SC PdTe₂ [50, 51, 52] as parent compounds to demonstrate that our SCTMs strategy is feasible. Experimentally, Bi₂Te₃ (space group $R\bar{\mathrm{3}}m$, $a=4.35$ Å, $c=30.36$ Å), PdTe (space group space group $P\rm{6_{3}}/mmc$, a = 4.152 Å, c = 5.671 Å, T_c = 2.3 K) and PdTe₂ (space group $P\bar{\mathrm{3}}m$1, $a=4.03$ Å, $c=5.12$ Å, T${}_{c}=1.64$ K) all adopt the triangle lattice and have very similar in-plane lattice constants, which makes it much easier to integrate them together to form a new compound. According to our calculations, the stable unit of PdBi₂Te₅ and PdBi₂Te₄ adopt octuple-layer (OL) structure and septuple-layer (SL) structure respectively, as shown in Fig. 3(a)(also Fig. S1) and Fig. S2 of Supplementary Material(SM) [63]. They both favor the ABC stacking along c-direction, and form the rhombohedral unit cell as shown in Fig. 3(a), which is 73 meV/f.u. (73 meV/f.u. for PdBi₂Te₄) and 46 meV/f.u. (12 meV/f.u. PdBi₂Te₄) lower than the AA and AB stacking structures. The detailed crystal parameters and total energy of different stacking PbBi₂Te₅ and PbBi₂Te₄ are tabulated in the Table. S1 and Table. S2, respectively [63].

The formation energy of PdBi₂Te₅ and PdBi₂Te₄ are calculated to study their thermodynamic stability by using $E_{f}^{\mathrm{Pd_{m}Bi_{n}Te_{l}}}=E^{\mathrm{Pd_{m}Bi_{n}Te_{l}}}-mE^{\mathrm{Pd}}-nE^{\mathrm{Bi}}-lE^{\mathrm{Te}}$, with $E^{i}$($i$=$\mathrm{Pd_{m}Bi_{n}Te_{l}},\mathrm{Pd},\mathrm{Bi}~{}\mathrm{and}~{}\mathrm{Te}$) means the calculated total energy per formula in the ground state. The calculated $E^{PdBi_{2}Te_{5}}_{f}$ and $E^{PdBi_{2}Te_{4}}_{f}$ are $-3.184$ eV/f.u. and $-2.476$ eV/f.u., which means that 3.184 eV and 2.476 eV can be released during their synthesis processes from the constituent elements. To further manifest their thermodynamic stability, we construct the convex hull diagram in Fig. 3(b) with all of the synthesized Pd-Bi-Te compounds, whose crystal parameters and the calculated formation energy have been tabulated in Table. S3 and Table. S4, respectively [63]. Fig. 3(b) shows that PdBi₂Te₅ and PdBi₂Te₄ are 13 meV/atom and 61 meV/atom above the convex hull respectively. Moreover, considering that metastable PdBi₂Te₃, 52 meV and 3 meV higher than PdBi₂Te₅ and PdBi₂Te₄ as shown in Fig. 3(b), has been synthesized in experiments [64, 65], we thus conclude that PdBi₂Te₅ and PdBi₂Te₄ could be synthesized in experiments. For PdBi₂Te₅, we propose a synthetic route through the growth of Bi₂Te₃ and PdTe₂ layer by layer. Our calculated results reveal that bulk PdBi₂Te₅ is 59 meV/f.u. lower than the total energy of free standing Bi₂Te₃ and PdTe₂ layers, which strongly suggest that PdTe₂ layer tends to deposit on Bi₂Te₃ to form new PdBi₂Te₅ crystal. To investigate their dynamical stability, we calculate the phonon dispersion of PdBi₂Te₅ and PdBi₂Te₄, and plot them in Fig. 3(c) and Fig. S3(a) [63]. There are 24 (21) phonon modes with fully real positive frequencies for PdBi₂Te₅ (PdBi₂Te₄), which indicates that the rhombohedral unit cells are dynamically stable. Based on these results, we conclude that PdBi₂Te₅ and PdBi₂Te₄ are relatively thermodynamically and dynamically stability in the rhombohedral structure, and further experimental investigation is called for.

Then we study the electronic structures and topological properties of PdBi₂Te₅ and PdBi₂Te₄. Since PdBi₂Te₅ and PbBi₂Te₄ exhibit similar electronic structures and non-trivial band topology, we only show the detailed density of states (DOS), band structures, and topological surface states of PdBi₂Te₅ as an example in the main text, one can check the results of PdBi₂Te₄ in Section III and Figs. S3 of the SM [63]. In Fig. 3(d), we plot the total and projected DOS of PdBi₂Te₅, which gives rise to DOS(0 eV) = 1.91 states/eV at Fermi level, indicating its metallic nature and the possibility of superconductivity. The projected DOS demonstrates that the states between $-1$ eV and $1$ eV are dominated by the $p$-orbitals of Te hybridized with $d$-orbitals from Pd and $p$-orbitals from Bi. The hybridization is also manifested by the projected band structures shown in Fig. 3(e), which shows that two bands with $p$-orbital components from Te or Bi cross the Fermi level and form several Fermi surfaces. Further detailed orbital components analysis demonstrates that a continuous band gap (yellow region in Fig. 3(e)) and band inversion exists between the nominal valence band and conduction band around the Fermi level, which implies that PdBi₂Te₅ inherits the topological electronic nature of Bi₂Te₃ successfully. The nontrivial band topology can be confirmed by calculating the $Z_{2}$ topological invariant of time-reversal invariant insulators [66]. Given that rhombohedral PdBi₂Te₅ possesses inversion symmetry and a continuous band gap, the $Z_{2}$ topological invariant $\nu_{\text{TI}}=(1-P)/2$ is determined by the product $P$ of the parity of the wave function at the TRIM points [66]. Our calculated results give $Z_{2}$ index $\nu_{\text{TI}}=1$, confirming PdBi₂Te₅ is a $Z_{2}$ topological metal. To visualize the bulk–boundary correspondence, we calculate and plot the topological surface states on the (001) surface in Fig. 3(f). The surface states are similar to that of Bi₂Te₃ [59, 60], the Dirac cone at the $\Gamma$ point manifest approximately $-6.3$ meV below the Fermi level (the dashed line in Fig. 3(f)).

To investigate the superconducting property of PdBi₂Te₅, we perform the electron-phonon calculations based on density functional perturbation theory [67]. The calculated electron-phonon coupling constant $\lambda=0.43$ and logarithmic average phonon frequency $\omega_{\log}=97~{}cm^{-1}$, as tabulated in Table. S5 [63]. Furthermore, the superconducting transition temperature ($T_{c}$) is estimated by using the reduced Allen-Dynes formula [68, 69]:

where $\mu^{*}$ is the effective Coulomb potential. By adopting a typical $\mu^{*}$ = 0.1, the $T_{c}$ of PdBi₂Te₅ is estimated as 0.57 K. As comparison, the calculated $\lambda$ and $\omega_{\log}$ in PdTe₂ is 0.52 and 112 $cm^{-1}$, respectively. Accordingly, the estimated $T_{c}$ in PdTe₂ is 1.59 K, which agrees well with the experimental $T_{c}$ of 1.64 K [50, 51, 52]. These results clearly demonstrate that the SC in PdTe₂ is well inherited into the PdBi₂Te₅.

We now study the TSC property of the PdBi₂Te₅ slab by introducing the SC pairing and Zeeman splitting into the topological surface states. Usually, the Zeeman splitting is applied by external magnetic field or in proximity to a FM insulator, as illustrated in Fig. 4(a). As a concrete example, we use a 2D slab consisting of 10-OL PdBi₂Te₅, which is thick enough to avoid the hybridization between top layer and bottom layer (Fig. 3(f)). Since PdBi₂Te₅ is an intrinsic SC, the estimated $s$-wave superconducting gap $\Delta_{s}=1.0$ meV is introduced globally for all 10-OLs. The out-of-plane Zeeman splitting is applied only in the bottom layer consisting of one Bi₂Te₃ and one PdTe₂, by assuming PdBi₂Te₅ is the conventional SC from the parent type-I SC PdTe₂ [50]. The chemical potential $\mu$ is set at the energy of surface Dirac cone at the $\Gamma$ point (about $-6.3$ meV below the Fermi level). In Fig. 4(b), we show the low energy spectrum of $H_{\mathrm{BdG}}$ at $\Gamma$ point as a function of Zeeman splitting energy $M_{z}$, which manifest that the superconducting spectrum is fully gaped with an energy gap of $\Delta$ at $M_{z}=0$. As $M_{z}$ increases, the superconducting gap at the $\Gamma$ point closes and reopens. This behavior indicates that a topological phase transition happens at critical point $M_{z}/\Delta=3.1$, and this 2D slab enters chiral TSC phase characterized by a nonzero BdG Chern number and chiral Majorana edge states according to previous model simulations [38, 39, 40].

To firmly verify its topological property and visualize the low energy physics in the TSC phase, we calculate the superconducting energy spectrum at $M_{z}$=5 meV and $\Delta$=1 meV in Fig. 4(c). The corresponding Wilson loop evolutions for the occupied states are ploted in Fig. 4(d). The zoom-in image of Fig. 4(c) reveals that a full superconducting gap is opened in the whole BZ, indicating that the system is a well defined chiral TSC. The Wilson loop evolution exhibits a nontrivial chiral winding number 1, which directly confirms the superconducting BdG Chern number $C_{\textrm{BdG}}=1$. Given that the experimental accessible magnetization energy usually reaches a few tens of meV, our results provide a feasible guideline for discovery the chiral TSC phase in PdBi₂Te₅.

Finally, we would like to point out that the chiral TSC phase could also be realized in PdBi₂Te₄ as shown in Fig. S4 [63], which exhibits a similar superconducting spectrum gap closing behavior with respect to $M_{z}/\Delta$ as in PdBi₂Te₅. In addition, we emphasize that our material design strategy can also be applied to search for other SCTM candidates. For example, our calculated results demonstrate that AuBi₂Te₅ formed by SC AuTe₂ interacting into Bi₂Te₃ is also an ideal SCTM, whose detailed crystal structures, dynamic stability, electronic structures, and topological surface states are discussed in Section V and Fig. S5 of SM [63]. Therefore, we expect that SCTM AuBi₂Te₅ could also be a TSC candidate. Last, we would like to point out that the program can be extended to study many 2D topological superconducting heterostructure systems, such as magnetic TI/SC heterostructure, SC/FM heterostructure and SC/TI/SC heterostructure. This will make it possible to determine the accurate parameters of the TSC phase and simulate their TSC property in such systems from first-principles calculations. We expect our program to be also useful for optimizing the experimental setup, stimulating the field of TSC study.

This work is supported by the National Key Research and Development Program of China (2018YFA0307000), and the National Natural Science Foundation of China (12274154, 11874022). B.L. is supported by the Alfred P. Sloan Foundation, the National Science Foundation through Princeton University’s Materials Research Science and Engineering Center DMR-2011750, and the National Science Foundation under award DMR-2141966.

The first-principles calculations based on density functional theory are performed by the Vienna ab initio simulation package [70, 71] with treating Perdew–Burke–Ernzerhof type of generalized gradient approximation as the exchange-correlation potential [72]. The cutoff energy for wave function expansion is set as 450 eV, $k$-points grid 13$\times$13$\times$13 is used for sampling the first BZ. All crystal structures are fully optimized until the force on each atom is less than 0.01 eV/Å, and the SOC is included self-consistently. The electron-phonon coupling calculations with van der Waals correction [73] are carried out in Quantum Espresso [74] based on the perturbation theory. A Hermite-Gaussian smearing of 0.0025 Ryd is used for the electronic integration. The 8$\times$8$\times$8 $k$-mesh is used for the electron-phonon coupling strength $\lambda$ calculations, and the dynamical matrices are calculated on a 4$\times$4$\times$4 phonon-momentum grid. Besides, a 2$\times$2$\times$2 supercell is built to calculate the phonon dispersion by using PHONOPY [75]. For the surface calculation, the Wannier functions of Pd-$d$, Bi-$p$ and Te-$p$ orbitals are constructed by using WANNIER90 [76]. A slab consisting of 10-OL PdBi₂Te₅ layers with a bottom surface terminated as the Bi₂Te₃ layer is implemented in WannierTools [55], which is further used to calculate the electronic surface states, the superconducting spectrum, and the superconducting BdG Chern number.