Evolution of the topologically protected surface states in superconductor $\beta$-Bi$_{2}$Pd from the three-dimensional to the two-dimensional limit

TL;DR

This study uses first-principles calculations to explore how topologically protected surface states in $eta$-Bi$_{2}$Pd evolve from three-dimensional bulk to two-dimensional thin films, highlighting the role of spin-orbit coupling and film thickness.

Contribution

It provides a detailed analysis of the evolution of topological surface states in $eta$-Bi$_{2}$Pd as a function of film thickness using first-principles methods.

Findings

Topological surface states emerge at 9 triple-layers.

Rashba-type surface states appear at 11 triple-layers.

Van der Waals corrections are essential for accurate structural modeling.

Abstract

The recent discovery of topologically protected surface states in the noncentrosymmetric $\alpha$-BiPd and the centrosymmetric $\beta$-Bi$_{2}$Pd has renewed the interest in the Bi-Pd family of superconductors. Here, we employ first-principles calculations to investigate the structure, electronic, and topological features of $\beta$-Bi$_{2}$Pd, in bulk and in thin films of various thicknesses. We find that the Van der Waals dispersion corrections are important for reproducing the experimental structural parameters, while the spin-orbit interaction is critical for properly describing the appearance of topological electronic states. By increasing the thickness of the slab, the Dirac-cone surface states and the Rashba-type surface states gradually emerge at 9 and 11 triple-layers.

Tables

Table 1. Table 2: The band gap at DP and RP for 1-15 β 𝛽 \beta –Bi 2 Pd TLs calculated with inclusion of SOI.

Slabs	Gap at DP (eV)	Gap at RP (eV)
1 TL	0.703	1.273
3 TLs	0.084	0.320
5 TLs	0.012	0.000
7 TLs	0.002	0.000
9 TLs	0.000	0.000
11 TLs	0.000	0.000
13 TLs	0.000	0.000
15 TLs	0.000	0.000

Full text

I Introduction

The prediction and discovery of topological insulators (TIs) Fu07 ; Hsieh08 ; ZhangHJ ; XiaY ; Hsieh09 ; Chen09 ; Liu ; Hasan10 ; Sawai10 ; Eremeev10 ; Yang ; Kolmogorov16 ; Bansil16 have triggered the quest for other classes of materials that host exotic quantum states, such as topological superconductors XLQi2 ; Hasan15 ; Sato16 and Weyl Huang_NCOM15 ; Weng_PRX15 ; Lv_NPHYS15 ; Yang_NPHYS15 ; Xu_NPHYS15 , Dirac Wang_PRB12 ; Ali14 ; Borisenko_PRL14 ; Liu_SCI14 ; Neupane_NCOM14 , and nodal-line Bian_NCOM16 ; Wu_NPHYS16 semimetals. Similar to TIs, topological superconductors (TSCs) are characterized by a full paring gap in the bulk and topologically protected gapless states on the edge or surface that can support massless Majorana fermions XLQi2 ; Schnyder-PRB15 ; Ando15 ; Xu15 ; Hasan15 . These exotic states may find promising applications in spintronics and quantum computing Fu08 ; Nayak .

Although numerous materials have been identified as TIs, only a handful of systems have been found to exhibit the signatures of topological superconductivity. The proposed candidate TSCs have been realized by introducing carriers into topological insulators, such as Cu-intercalated Bi2Se3 Fu10 ; Hor ; Kriener ; Sasaki11 and In-doped SnTe Sasaki12 , or by applying pressure to topological parent compounds, such as Bi2Te3, Sb2Te3, and Bi2Se3 ZhangPNAS ; ZhuSR ; Kirshenbaum-PRL ; KongPP . More recently, two members of the Bi-Pd family of superconductors have been found to exhibit topologically protected surface states with Rashba-like spin splitting. While the superconducting state of noncentrosymmetric BiPd appears to be topologically trivial Mondal ; Sun ; Neupane , the pairing symmetry of centrosymmetric $\beta$-Bi2Pd remains to be elucidated Sakano ; Herrera ; Che ; Kacmarcik ; Biswas ; Margine ; LvYF . Despite the fact that spin- and angle-resolved photoemission spectroscopy (ARPES) has revealed the presence of notrivial surfaces states at the Fermi level Sakano , no Andreev bound states associated with Majorana fermions have been observed in point-contact spectroscopy Che . On the other hand, according to a recent scanning tunneling microscopy study, Majorana zero modes have been identified in $\beta$-Bi2Pd crystalline films LvYF . Along with the recently discovered PbTaSe2 superconductor Ali_PRB14 ; Bian_NCOM16 ; Chang_PRB16 , $\beta$-Bi2Pd constitutes a promising candidate for the long-sought stoichiometric TSC (i.e, a three-dimensional TSC realized in the absence of any external factors such as doping or pressure).

In the present work, first-principles calculations are carried out to systematically study the thickness-dependent evolution of the electronic structure and topological properties in $\beta$-Bi2Pd thin films. The bulk material has a layered structure consisting of stacked Bi–Pd–Bi triple layers (TLs). Owing to the weak coupling between two adjacent TLs, it has been possible to achieve TL-by-TL growth of $\beta$-Bi2Pd thin films using molecular beam epitaxy LvYF . This technique has also been used to grow high-quality thin films of Bi2Se3, Bi2Te3, and Sb2Te3 of various thicknesses by precisely controlling the growth conditions He ; Li_AdvMat10 ; WangG . In addition to the electronic properties of thin films, we also investigate the effect of van der Waals (vdW) and spin-orbit interaction (SOI) on the structure and the band dispersion of bulk $\beta$-Bi2Pd.

II Methods

The calculations are performed at the density functional theory level, employing the projector augmented wave method as implemented in the VASP package Kresse3 . The Perdew, Burke, and Ernzerhof (PBE) PBE form of the generalized gradient approximation is chosen to describe the exchange-correlation energy. Bi 5$d^{10}$6$s^{2}$6$p^{3}$ and Pd 4$p^{6}$4$d^{9}$5$s^{1}$ orbitals were included as valence electrons. To properly treat dispersion corrections, we use the non-local vdW density functional optB86b-vdW vdW1 ; vdW2 . Relativistic effects are included in the calculations in terms of the SOI SOI . The energy cutoff of the plane-wave basis-set is 500 eV. The reciprocal space integration is performed over 16$\times$16$\times$16 and 10$\times$10$\times$1 k-point grids in the Brillouin zone (BZ) for bulk and thin films, respectively. The structural parameters of the bulk structure are fully optimized until the Hellmann-Feynman force on each atom is less than 0.01 eV/Å. To model thin films, we use a slab supercell configuration with a 15 Å  thick vacuum layer in the (001) direction to eliminate the interslab interaction. The free-standing slab models are constructed with the bulk structural parameters optimized with optB86b-vdW.

III Results

III.1 Bulk properties

We begin with a review of the properties of bulk $\beta$–Bi2Pd which crystallizes in the centrosymmetric body-centered tetragonal (bct) crystal structure with space group I4/mmm (No.139) as shown in Figs. 1(a)-(b). The Wyckoff positions for the Bi and Pd atoms are 4e(0,0,$z_{Bi}$) and 2a(0,0,0), respectively. In the conventional unit cell, the atomic structure can be visualized as a superposition of triple layers (TLs) with a Bi–Pd–Bi sequence at the center of the unit cell along the $c$-axis. Within each triple layer every Pd atom occupies the center of the cube formed by 8 Bi atoms. The bulk BZ and the BZ projected onto the (001) surface are also shown in Fig. 1(d).

In Table I, we list our calculated equilibrium structural parameters found with PBE, PBE-vdW, and PBE+SOI, together with the corresponding values from previous theoretical Shein and experimental Zhuravlev ; Imai studies. PBE-vdW and PBE+SOI refer to calculations performed with vdW and SOI, respectively. Overall, our relaxed lattice parameters are in good agreement with previous calculations performed with the full-potential augmented plane wave method (FLAPW)-PBE/PBE+SOI Shein . When compared to the experimental values, the in-plane $a$ and out-of-plane $c$ lattice parameters are overestimated by approximately 1.4-1.6% (2.1-2.4%) and 0.9-1.1% (2.0-2.2%) in the case of PBE (PBE+SOI). The inclusion of the vdW dispersion corrections gives a much better agreement with the experimental values and, in this case, the two lattice parameters $a$ and $c$ are overestimated and underestimated by approximately 0.7%, respectively. The importance of the long-range dispersive interactions in predicting ground states and describing the interlayer distance in different classes of materials has long been recognized Kolmogorov ; Mishra ; SaPRL ; WangSb2Te3 ; SaNanoscale .

A simple picture of the chemical bonding can be obtained from the analysis of the charge density in the (110) plane of $\beta$–Bi2Pd. As shown in Fig. 1(c), the bond strength distribution predicts strong covalent bonding within the triple layer (0.035 $e$/au3 at the middle of Bi–Pd bonds) and weak van der Waals interaction between adjacent TLs (0.012 $e$/au3 at the middle of Bi–Bi bonds). These values are similar to the ones reported in the three-dimensional (3D) TIs Bi2Se3, Bi2Se3 and Sb2Te3 Mishra ; WangBi2Se3 ; WangSb2Te3 .

To examine the effect of the vdW and SOI on the electronic structure of bulk $\beta$–Bi2Pd, we compare in Figs. 2(a)–(d) the band structures obtained with PBE, PBE-vdW, PBE+SOI, and PBE-vdW+SOI PBE-vdW+SOI . As shown in Fig. 2(a), the vdW interaction has a limited effect, the bands at the Fermi level remain practically unchanged, while the deeper bands are only slightly pushed down. The picture is completely different when the SOI is included, the band structure undergoes significant changes over the whole BZ [see Fig. 2(b)]. First, one band at the $\Gamma$ point is downshifted by 1.46 eV and resides just above the Fermi level. Second, a 0.46 eV gap opens at the band crossings around the $Z$ point approximately 2 eV below the Fermi level [within the gray shaded area in Fig. 2(d)]. Third, two continuous bulk band gaps are formed across the whole BZ at the Fermi level. These gaps are highlighted in cyan and yellow shades in Fig. 2(d). Finally, similar to the results without SOI, we find that the reduction in the interlayer spacing due to vdW corrections has only a small effect on the band structure obtained in the presence of SOI [see Fig. 2(c)].

The inversion symmetry of $\beta$–Bi2Pd allows one to perform a parity analysis and estimate the Z2 invariant for the three SOI-induced gaps Fu08 . By checking the parity at the eight time-reversal invariant momenta (1$\Gamma$, 2X, 4N, and 1Z), we find that the middle gap is trivial, while the upper and lower gaps have a non-trivial topological character in agreement with previous calculations Sakano . The two topological gaps could harbor gapless surface states and ARPES measurements have indeed revealed that a Dirac-type band dispersion appears in the lower gap, 2.41 eV below the Fermi level Sakano . The absence of a gapless surface state in the upper gap in the ARPES spectra is due to its localization above the Fermi level. However, a spin polarization analysis has clearly showed that both topological and trivial surface states cross the Fermi level Sakano .

The above results demonstrate that while the vdW dispersion corrections are important for obtaining good agreement with the experimental reference values for structural parameters, the SOI is critical for properly describing the appearance of topological electronic states.

III.2 Thin films properties

We further study the thickness-dependent band structures of $\beta$–Bi2Pd thin films in order to establish the critical (minimum) thickness at which the topologically protected surface states develop. We consider several thin films made of 1, 3, $\cdots$, 15 TLs. Thin films with an even number of TL are excluded since they break the inversion-symmetry present in the bulk. All these slab models are constructed with the bulk structural parameters optimized with PBE-vdW. To check the effect of the lattice relaxation, we have fully optimized the 3 TLs slab with PBE-vdW. No visible changes were found between the band structures calculated with and without SOI in the case of 3 TLs slabs with optimized and unoptimized lattice parameters.

Figure 3 illustrates the evolution of the surface band dispersion as a function of slab thickness along the $\overline{\rm{M}}$-$\overline{\Gamma}$-$\overline{\rm{X}}$ direction. When the film thickness is small enough, the interactions between the states localized on the top and bottom surfaces lead to the opening of a gap at the Dirac point (DP). This area is indicated by the green rectangle in Fig. 3. Here, we label the surface states (SSs) that form the Dirac cone as SS1 and SS2. For 1 TL, a relatively large gap of about 0.70 eV is formed -2.45 eV below the Fermi level at the $\overline{\Gamma}$ point between a nearly parabolic band (SS2) and an almost flat band (SS1). For 3 TLs, the upper band moves down by 0.40 eV and the lower band moves up by 0.21 eV, such that the gap is reduced to 0.08 eV. For 5 TLs, the two bands are almost touching each other at an energy level of -2.19 eV and the gap is around 0.01 eV. Once the thickness of the slab reaches 9 TLs, the gap completely disappears and the two bands with nearly linear dispersion cross each other forming a surface Dirac cone (S-DC) Sakano . We find that the energy of the DP is ED=2.20 eV, very close to the ED=2.41 eV value extracted from the ARPES data Sakano . In Table 2, we list the gap values at the DP for all slabs considered. The thickness-dependent behavior of the gap at the DP found here is similar to the behavior observed in typical 3D topological insulators He ; Li_AdvMat10 ; WangG ; Jiang ; Bian12 ; Yan14 .

Besides the S-DC, an additional topological feature develops in thicker films at the $\Gamma$ point approximately 2.26 eV above the Fermi level. As shown in Fig. 3, with increasing thickness more bands emerge and two evolve into a Rashba-like crossing point above 3 TLs. These surface states are labeled as SS3 and SS4 in Fig. 3. The gap values at the Rashba point (RP), presented in Table 2, also support this feature. Our finding is consistent with a previous result reported for a slab of 11 TLs Sakano . We also show in Fig. 4 the orbital-resolved band dispersion without and with SOI for 11 TLs. At the DP, both SS1 and SS2 have Pd $d_{z^{2}}$ character, while at the RP, both SS3 and SS4 have Bi $p_{x+y}$ character. The SOI-induced crossing at the DP and at the RP can be seen more clearly in Figs. 4(c-f).

In order to characterize the spatial spread of the SSs (i.e., the penetration depth of the surface-state wave functions into the bulk), we calculate the real space charge distribution of the SSs at the $\overline{\Gamma}$ point. In Fig. 5, we show the charge density integrated over the $xy$–plane as a function of the $z$ coordinate. For the SS1 and SS2 states, starting from the 7 TL film, the charge density is concentrated mostly inside the top and bottom TLs and the wave functions penetrate only 15 Å  in depth from the surface. In contrast, the SS3 and SS4 states only become surface states starting from the 11 TLs film and have a much larger penetration depth of about 30 Å. Based on both the band structure and the charge density analysis, the critical thickness for emergence of the Dirac-cone SSs is 9 TLs, while that for the Rashba-type SSs is 11 TLs.

To further characterize the distribution of the surface states over atomic layers, we show in Fig. 6 the charge density integrated over the $xy$–plane for 11 TLs at several $k$-points in the vicinity of the $\overline{\Gamma}$ point [$k_{1}$–$k_{3}$ points indicated with blue vertical lines in Fig. 3(d)]. For the $k_{1}$ point closer to $\overline{\Gamma}$, the distribution of the wave functions corresponding to SS1, SS2, and SS3 remains localized to the surface region. On the contrary, the real-space distribution of SS4 displays some weight at the center of slab. Moving away from the $\overline{\Gamma}$ point, the surface states penetrate deeper into the slab and the wave functions acquire a bulk-like spatial distribution. Note that the SS3 still has the characteristic of surface state even at $k_{2}$. Using this charge distribution analysis, we determine the states away from the $\overline{\Gamma}$ point that retain a charge distribution localized in the surface region. These are shown as red circles in Fig. 3(d). The localization of the SSs found in our study is similar to the experimental results for $\alpha$–BiPd Benia and theoretical results for Bi2Se3 ZhangNJP .

IV Conclusion

In summary, we have investigated the structural, electronic, and topological properties of the centrosymmetric superconductor $\beta$–Bi2Pd in the bulk and in thin films using first-principles calculations. We have shown that the structural parameters obtained with the vdW interaction are in very good agreement with experiment, and that the inclusion of SOI is necessary for correctly describing the band structure. In addition, we have uncovered that thicknesses of 9 and 11 TLs are required for the appearance of the Dirac and Rashba surface states, respectively. The penetration depth and the charge distribution of these states have also been discussed. These findings may prove important for further exploration of the superconducting state in $\beta$–Bi2Pd thin films.

V Acknowledgments

B. Wang thanks B. Sa and H. Shi for helpful discussions.