Probing the wavefunction of the surface states in Bi$_2$Se$_3$ topological insulator: a realistic tight-binding approach

We report on microscopic tight-binding modeling of surface states in Bi$_2$Se$_3$ three-dimensional topological insulator, based on a sp$^3$ Slater-Koster Hamiltonian, with parameters calculated from density functional theory. The effect of spin-orbit interaction on the electronic structure of the bulk and of a slab with finite thickness is investigated. In particular, a phenomenological criterion of band inversion is formulated for both bulk and slab, based on the calculated atomic- and orbital-projections of the wavefunctions, associated with valence and conduction band extrema at the center of the Brillouin zone. We carry out a thorough analysis of the calculated bandstructures of slabs with varying thickness, where surface states are identified using a quantitative criterion according to their spatial distribution. The thickness-dependent energy gap, attributed to inter-surface interaction, and the emergence of gapless surface states for slabs above a critical thickness are investigated. We map out the transition to the infinite-thickness limit by calculating explicitly the modifications in the spatial distribution and spin-character of the surface states wavefunction with increasing the slab thickness. Our numerical analysis shows that the system must be approximately forty quintuple-layers thick to exhibit completely decoupled surface states, localized on the opposite surfaces. These results have implications on the effect of external perturbations on the surface states near the Dirac point.


I. INTRODUCTION
Topological insulator 1,2 (TI) materials host on their boundaries a novel type of topological states of quantum matter, which, unlike the quantum Hall state, exist without the breaking of time-reversal symmetry. 3,4 Theoretical prediction and subsequent experimental demonstration of these topological states in both two-5,6 (2D) and three-dimensional 7-17 (3D) systems have given rise to what is now one of the most rapidly developing fields in condensed matter physics. Apart from providing a test platform for fundamental concepts, the study of TIs holds promise for novel applications in materials science and chemistry, 18 spintronics 19 and quantum computation. 20,21 However, to be able to fully explore the potential of TIs, it is essential to have a detailed knowledge of the nature and properties of topological surface states in real TI materials, 22 as well as a quantitative understanding of how they respond to external perturbations. [23][24][25] Experimentally, these questions are being addressed with advanced surface-sensitive experimental probes, such as spin-and angle-resolved photoemission spectroscopy 16 [(SR)-ARPES] and scanning tunneling microscopy 26,27 (STM).
Along with experimental advances, there is a growing need for atomistic modeling of TIs that would enable quantitative predictions and direct comparison with experiment. Significant progress has been made in using ab initio methods to calculate electronic 13,28-30 and magnetic [31][32][33][34] properties of TIs. However, such methods suffer from severe computational limitations, particularly in the case of slab geometry as well as surface supercell calculations, which are employed in studies of impuritydoping effects. In addition, more accurate ab initio methods often lack the conceptual transparency and flexibility of the model Hamiltonian approaches, which have been of fundamental importance for driving progress in this research field. 3,4 Microscopic tight-binding (TB) models, which have already proved successful in quantitative description of electronic and magnetic properties of semiconductors, 35,36 may provide a convenient platform to address similar issues in TIs. Several studies have recently appeared in the literature, in which TB descriptions with different level of complexity have been introduced, ranging from models built on a simplified lattice structure 37 or a restricted orbital basis set inferred from symmetry arguments 38,39 to fully microscopic models, with parameters extracted from density functional theory (DFT). 28,[40][41][42] To date, the latter class of models is still the least represented among the model Hamiltonian approaches to TIs.
In this work we employ a microscopic TB model to study the properties of surface states in Bi 2 Se 3 , a prototypical 3D TI, which belongs, along with Bi 2 Te 3 and Sb 2 Te 3 , to the family of binary tetradymite semiconductors with layered structure. 43 Although these materials have been studied for decades due to their excellent thermoelectric properties, [44][45][46] they have recently attracted considerable attention as 3D TIs, e.g. materials that exhibit topologically protected conducting surface states with linear (Dirac) dispersion and helical spin-texture, traversing the bulk insulating gap. [14][15][16][17] Due to a relatively large band gap (0.3 eV for Bi 2 Se 3 ) and rather simple surface states, consisting of a single Dirac cone, 16 the Bi 2 Se 3 family of 3D TIs is the most studied both exper-imentally and theoretically.
Our treatment is based on the sp 3 Slater-Koster Hamiltonian. 47 We use the parametrization developed by Kobayashi,40 by fitting to DFT calculations. Throughout this work, our strategy has been to make use of the computational efficiency and simplifications, offered by the TB approach, in order to investigate key features of the surface states in Bi 2 Se 3 3D TI, which are inaccessible by ab initio methods. Importantly, we consider slabs with thicknesses ranging from 1 to 100 quintuple layers (QLs), which corresponds to length scales in the range of 1-100 nm. In contrast, thicknesses typically investigated in ab-initio-based studies do not exceed several quintuple layers. 30 In agreement with previous reports, 48-50 we find a gap due to interaction between opposite surfaces, which decreases with increasing the slab thickness. Starting from 5 QLs, the size of the gap becomes smaller than 10 −3 eV, and one can identify surface states with linear dispersion and helical spin-texture. For each slab thickness we determine the surface character of Bloch states using the procedure put forward in Ref. 29, i.e. based on the contribution of the real-space projected wavefunction onto the two surfaces of the slab.
Explicit calculations of the atomic-and orbitalprojections of the wavefunctions, associated with valence and conduction band extrema in both bulk and slab geometry, allowed us to construct a phenomenological picture of band inversion. The latter effect is induced by spin-orbit interaction and is responsible for the occurrence of topological surface states across the bulk insulating gap. 13 Furthermore, based on a similar analysis, we were able to track the changes in the spatial distribution and the spin character of the surface states wavefunctions at and in the vicinity of the Dirac point, for increasing slab thickness. Our calculations showed that the states corresponding to top and bottom surfaces become completely decoupled, i.e. spatially separated, only for very thick slabs containing 40 QLs. We also calculated the spin-orientation of the surface states in momentum space as a function of thickness. The disturbances in the helical spin-texture, expressed through in-plane and out-of-plane tilting angles of the spin, are shown to be significant for thin slabs up to 10 QLs and are manifestations of both inter-surface interaction and the proximity to bulk states.
The rest of the paper is organized as follows. In Section II we discuss the details of the TB model and some computational aspects. The results of the simulations are presented in Section III. We begin in Section III A by analyzing the effect of SO on the bulk bandstructures, where the bulk band inversion is interpreted as a characteristic change in the spatial distributions of the orbital-resolved projections of the wavefunctions at conductance and valence band extrema at theΓ point. The results of surface bandstructure calculations are presented in Section III B, in particular the spatial character of the states emerging within the bulk gap and the opening of the gap at finite thicknesses are described quantitatively in Section III B 1.
We also discuss the procedure to identify the inverted character on conduction and valence bands in the presence of SO for a slab geometry. In Section III B 2, based on the analysis of the wavefunctions of the states at the Dirac point, a quantitative criterion of an infinitely thick slab, with no interaction between the surfaces, is formulated. In addition, we comment on the orbital character of the wavefunctions at the Dirac point, especially on the non-negligible contribution of the in-plane p-orbitals. The spin-properties of the Dirac-cone states and deviations from the perfect spin-momentum locking due to inter-surface interaction and the presence of bulk states are calculated in Section III B 3. Finally, Section IV contains our conclusions.

II. METHODS
We begin with a brief description of the crystal structure and chemical bonding of materials belonging to the Bi 2 Se 3 family. Bi 2 Se 3 has a rhombohedral crystal structure with crystallographic R3m [D 5 3d ] space group, with five atoms in the unit cell [see Fig. 1(a)]. The crystal is formed by stacking of hexagonal monolayers of either Bi or Se in a close-packed fcc (ABC) fashion. The sequence of five atomic layers, i.e. Se1-Bi-Se2-Bi-Se1, forms a QL [see Fig. 1(b)], with Se1 and Se2 indicating two nonequivalent positions of Se atom. Hence it is convenient to describe the crystal structure of bulk Bi 2 Se 3 in terms of QLs stacked along the direction perpendicular to atomic planes (c-axis).
The chemical bonding between atoms in the atomic layers within a QL is of covalent-ionic type with dominant covalent character. Adjacent QLs are weakly bound through Se1-Se1 bonds by van der Waals forces, 44 allowing easy cleavage of Bi 2 Se 3 on the (111) Se surface plane. Based on the electronic configurations of Bi and Se, with the outermost orbitals being of p-character for both types of atoms, it is reasonable to expect that the bonding within QLs is mainly due to interactions between p-orbitals. In fact, a simple but useful physical picture of the chemical bonding present in QLs features ppσchains formed by strongly interacting p-orbitals of atoms in nearest-neighbors atomic layers. Evidence of this coupling can be indirectly inferred from STM topographies of Bi-antisites defects in Bi 2 Se 3 . 45 In Sections III A and III B, we will comment on the orbital character of the valence and conduction band extrema in bulk Bi 2 Se 3 and of surface states in Bi 2 Se 3 slab, respectively.
For modeling of the electronic structure of Bi 2 Se 3 we employ the sp 3 TB model with Slater-Koster parameters obtained by Kobayashi 40 by fitting to bulk bandstructures calculated with DFT. Interactions between atoms in the same atomic layer and between atoms in first and second nearest-neighbor layers are included. The spinorbit interaction (SO), which is the key element leading to a non-trivial bulk band structure and metallic surface states, is incorporated in the intra-atomic matrix elements. 51 Thus the Hamiltonian of the system readŝ where k is the reciprocal-lattice vector that spans the Brillouin zone, i(i ′ ) is the atomic index, α(α ′ ) labels atomic orbitals {s , p x , p y , p z }, and σ(σ ′ ) denotes the spin. Here i refers to the atomic positions in the unit cell, while i ′ =i runs over all neighbors of atom i, including atoms in the adject cells, with r ii ′ being the vector connecting the two atoms (r ii ′ =0 for i=i ′ ). The coefficients t αα ′ ii ′ are the Slater-Koster parameters (for i=i ′ t αα ′ ii ′ ≡ t αα ii give the on-site energies) andĉ σ † iα (ĉ σ iα ) is the creation(annihilation) operator for an electron with spin σ at the atomic orbital α of site i. The second term in Eq. (1) represents the on-site SO, where |i, α, σ are spinand orbital-resolved atomic orbitals,ˆ L is the orbital angular momentum operator andˆ S is the spin operator; λ i is the SO strength. We refer to Ref. 40 for the exact parametrization of t αα ′ ii ′ and λ i . In bulk-bandstructure calculations we use the rhombo-hedral unit cell with five nonequivalent atoms, with the cell repeated periodically in x-, y-and z-directions. In calculations involving the Bi 2 Se 3 (111) surface we consider a slab consisting of N quintuple layers, or, equivalently, of 5N atomic layers. Since each atomic layer is an equilateral triangular lattice, it is sufficient to assign one atom per each layer, which gives a total of 5N atoms in the slab unit cell. The slab is finite along the z-direction (QL-stacking axis), with the unit cell repeated periodically in the x-y plane. Using this TB model we were able to compute bandstructures of slabs with thicknesses up to 100 QLs with a reasonable computational cost.

A. Band inversion in bulk Bi2Se3
The presence of gapless edge or surface states, robust against time-reversal-invariant perturbations, distinguishes a topological insulating phase from a trivial one. However, in order to elucidate the origin of the topological order in existing TIs and to facilitate the search for new TI materials, it is necessary to have a set of criteria that allow us to differentiate between topologically trivial and non-trivial insulators, based on the information about their bulk properties. Within the framework of topological band theory, 4,7-9,11,52 time-reversal-invariant insulators are classified according to a Z 2 topological invariant, assigned to their bandstructure. In 2D there is a single Z 2 invariant, which distinguishes the quantum spin-Hall state (2D TI), as in HgTe/CdTe quantum wells, from an ordinary insulator. In 3D a set of four Z 2 invariants leads to a classification based of three classes, namely strong TIs, weak TIs and ordinary insulators. For systems with inversion symmetry, the Z 2 invariant is determined by the product of the parity eigenvalues of occupied bands at the time-reversal-invariant momenta in the Brillouin zone. 8 Using this scheme, a few strong TIs have been predicted, such as strained α-Sn and HgTe, Bi 1−x Sb x and Bi 2 Se 3 -like 3D TIs.
The emergence of topological order can be understood using the concept of band inversion. 8,13,53,54 Clearly, a necessary condition for a TI is the existence of a nontrivial bulk band gap. In fact, in the above mentioned materials, a topological phase transition, or alternatively a change in the value of the Z 2 invariant which in turn implies the existence of gapless states on the boundary, is accompanied by a visible and well-defined change in the bandstructure. This change occurs precisely in the insulating gap and can be observed as a function of an external or an intrinsic parameter. In HgTe, which is a zero-gap semiconductor but acquires a gap due to external potential in a quantum well structure, the bands with s-and p-character are inverted with respect to their usual sequence at the Γ point. 5 A similar mechanism is realized in strained α-Sn, which has recently been shown to exhibit a 3D TI phase, with Dirac-like surface states em-anating from the second lowest valence band across the strained-induced gap. 42 In Bi 1−x Sb x alloy the inverted bandstructure is characterized by a change in the order of the bands with even (L s ) and odd (L a ) parity at the L point of the Brillouin zone compared to pure Bi. This inversion is induced by Sb doping and is due to the strong TI character of the valence band of pure Sb. 8 In Bi 2 Se 3 3D TI a non-trivial bulk band gap owns its existence to SO, with the parity of the valence and conduction band inverted at the Γ point. 13 Without SO the material would be a trivial insulator. However, in the presence of SO the change in the parity of one of the occupied bands leads to a change in the value of the Z 2 invariant, signaling a topological phase transition.
An intuitive way of describing band inversion theoretically is to look at charge density distributions of the inverted bands at the point of the Brillouin zone, where the inversion is expected to occur, as functions of a characteristic parameter. 53 Here we apply this procedure to bulk Bi 2 Se 3 . The calculated bandstructures without and with SO are shown in Fig. 2(a) and (b), respectively. The effect of SO is clearly seen as an increase of the gap and a change in the curvature of the valence band at the Γ point. To further quantify the changes in the valence and conduction band character, we calculate the atomic-and orbital-projections of the wavefunctions associated with the eigenvalues at the valence-band minimum (VBM) [filled circles in Fig. 2], and the conduction-band maximum (CBM) [open circles in Fig. 2] at the Γ point. The orbital-resolved projection of the wavefunction |ϕ , corresponding to an eigenvalue ε n of the TB Hamiltonian in Eq. (1), onto atomic orbital |i, α, σ is calculated as |ϕ iα | 2 ≡ σ | i, α, σ| ϕ | 2 . The total weight of the wavefunction at the atomic site i is then given by The calculated spatial distribution and orbital character of the wavefunctions at VBM and CBM are plotted in Fig. 3 for all five atoms of the bulk unit cell.
The most noticeable feature of Fig. 3 is the inversion in the spatial distribution of the wavefunctions when SO is switched on. Without SO the CBM wavefunction has a node at Se2, which is in fact the inversion center of the bulk crystal, with the maximum weight on the two Bi atoms. In contrast, the maximum weight of the VBM wavefunction is distributed over Se2 and two Se1 atoms. With SO the situation is the opposite: now the VBM wavefunction has a node at Se2 while the CBM wavefunction is predominantly localized on Se atoms.
As expected, both VBM and CBM states are predominantly of p-character. The s-orbital contribution to the wavefunctions is approximately 10% for calculations with or without SO and the s-projections are not affected by the inversion. Note that in the presence of SO, the states at the Γ point are mostly originating from p z orbitals, however the p x(y) (in-plane) contribution is not negligible (40% for VBM and 20% for CBM).
In the next section, we will discuss the results of bandstructure calculations for Bi 2 Se 3 slabs, namely the emergence of conducting states across the inverted bulk gap in the presence of SO. In particular, we will show that the band inversion, characteristic of the bulk gap, can also be identified in the surface bandstructures.
B. Wavefunction-based analysis of surface states in Bi2Se3 slab

Surface bandstructures and thickness-dependent gap
The bandstructures of Bi 2 Se 3 slab of varying thickness, calculated using the TB model including SO, are shown in Fig. 4. A clear and sizable gap is found for 1 and 2 QLs (∆=0.84 eV and ∆=0.18 eV, respectively). Starting from 3 QLs, two surface bands, resembling the Dirac states, extend in the range of [-0.1;0.4] eV, with valence and conduction bands beginning to form below and above this range. Already for 3 and 4 QLs the two bands appear to be almost touching at theΓ point, however the calculated gap is not negligible and is found to be ∆=0.043 eV for 3QL and ∆=0.007 eV for 4 QLs. The presence of the gap due to finite thickness is a recently discovered feature of 3D TI thin films, which has been attributed to tunneling between surface states localized on the opposite surfaces of the film. 48 The size of the gap decreases with increasing the film thickness and is almost negligible for 6 QLs, which is quantitatively consistent with our numerical observations. Such gap-opening mechanism has been considered for possible applications in TI-based MOSFET devices. 55 The value of the gap for increasing slab thickness is plotted in Fig. 5. Interestingly, we find that the gap initially decreases and reaches the value of 1.7 · 10 −5 eV for 5 QLs but then increases up to 6.2 · 10 −4 eV for 6 QLs. After this sudden increase, the gap continues to decrease exponentially and we find ∆=2.6 · 10 −6 eV for 10 QLs. A similar non-monotonicity in the thickness-dependence of the gap was found in the range of 4-6 QLs in ab initio calculations for Bi 2 Te 3 thin films. 29 This result also resembles the oscillating behavior of the gap found in theoretical work based on effective models. 38,49,50 However, direct comparison with these effective-model results is not straightforward since, in contrast to our atomistic model with slab thickness measured in terms of elementary building blocks of the real material (QLs), there the gap is calculated as a continuous function of the size along the z-direction.
Since the gap decreases exponentially, we do not display the results for thicknesses greater that 10 QLs. Note, however, that for 40 QLs the size of the gap is smaller than the numerical accuracy of our calculations. For thicknesses larger than this critical value, we identify the crossing of the two branches of Dirac states with opposite group velocity as the Dirac point, which is found at 0.09 eV. In Section III B 2 we will investigate in more detail the effect of the coupling between the two surfaces of the slab on the states at the Dirac point. In particular, we will look at the spatial character of the wavefunctions associated with these states and their dependence on the slab thickness.
We will now comment on another feature of Fig. 4, namely on the procedure used to identify the surface states in the bandstructure calculations. The character of each Bloch state ε n (k), where n is the band index and k is momentum, is determined by the spatial distribution of the corresponding wavefunction: 29 if the relative weight of the atom-projected wavefunction i |ϕ i | 2 on the top and bottom QLs exceeds a critical value (critical percentage) γ, ε n (k) is identified as a surface state. Such criterion is reasonable for slabs with thickness greater than 3 QLs, since the typical penetration depth of the surface states in Bi 2 Se 3 is of the order of 1 QL, or 1 nm. 28 For ultra-thin slabs with thickness below 3QLs, the criterion has to be modified, namely we calculate the wavefunction projections onto top two and bottom two atomic layers. The value of the critical percentage is found empirically. Starting from a rough estimate of γ based on the ratio between the penetration depth and the slab thickness, the value should then be optimized in such a way that a small change around this value does not significantly modify the resulting distribution of the surface states. We find, however, that for 1-3 QLs the identification of such an optimized value of γ is problematic. To illustrate this point, we present the results of calculations for two different values of γ. In the case of 1 QL and γ=65%, for instance, all states in the range of energies considered can be identified as surface states. Hence, for such a thin slab it is reasonable to use a stronger criterion. As one can see from Fig. 4(a), for γ=90% only the upper band preserves the surface character while the lower band does not satisfy the criterion (we use terms "lower" and "upper" for the two bands defining the gap). Interestingly, this observation is consistent with scanning tunneling spectroscopy measurements on Sb 2 Te 3 ultra-thin films. 56 A similar situation is found for 2 QLs, i.e. for two different values of γ the lower band is clearly not a surface band. For a stronger criterion with γ=50% the upper band also partially looses its surface character. In the case of 3 QLs, varying γ in the range of 50-60% induces some visible changes in the distribution of the surface states, however, the surface character of the Dirac-like states is well captured for any value of γ in this range. For slabs with thicknesses greater than 3 QLs, the search for an optimized value of γ is significantly simplified. Already for 4 QLs, changing γ from 50% to 60% does not produce any significant difference. For 5 QLs-and 20 QLs-thick slabs we use γ=60%. For a very thick slab of 100 QLs, we focus particularly on the states in the vicinity of the Dirac point, e.g. in the range of [0.0;0.3] eV, and do not apply the procedure for the near-continuum of bulk states appearing above and below this range. In Fig. 4(g) we present the result for 100 QLs with γ=40%, where the surface character of the Dirac states is clearly confirmed.
Before concluding this section we repeat the procedure that was used in Section III A to illustrate the mechanism of band inversion in bulk Bi 2 Se 3 . For a thick slab [see Fig. 4(f) and (g)], one can identify a clear gap between valence and conduction bands, which is transversed by conducting surface states. As the thickness of the slab increases, the gap approaches its bulk value, which is found at the Γ point in the bulk bandstructure calculated with SO [ Fig. 2(b)]. VBM(CBM) at theΓ point can be defined as the first state below(above) the Dirac point. In the case without SO, no conducting surfaces states are present and a gap is found at theΓ point for all slab thicknesses (calculations are not shown here). Similarly to the SO case, as the thickness of the slab increases the value of the gap approaches that found in the bulk calculation without SO [ Fig. 2(a)]. The atomic-layer projected wavefunctions at VBM and CBM of a 20 QLs-thick slab without and with SO are shown in Fig. 6(a) and (b), respectively. One immediately notices that the spatial character of VBM and CBM is inverted when SO is switched on: similarly to the bulk case, with(without) SO the VBM(CBM) wavefunction has nodes on Se2 atoms. The effect becomes even more appreciable, if we sum the wavefunction projections on Se1, Se2 and Bi atoms over the entire slab [ Fig. 6(c)]. The result is clearly similar to what we found in the bulk case, indicating that the inverted character of the valence and conduction bands in the presence of SO is also an intrinsic property of the surface, with the difference that in the surface bandstructures one finds the topological surface states, dispersing linearly across the inverted band gap.

Probing the wavefunction at the Dirac point: effect of finite thickness
We will now consider the surface states and their corresponding wavefunctions, found at theΓ point in our surface bandstructure calculations. In a hypothetical situation, when the slab is thick enough so that the two surfaces do not interact with each other, one expects to find four degenerate states at theΓ point, namely one Kramers degenerate pair for each of the two surfaces. Away from theΓ point, the four-fold degeneracy is lifted and one expects a doubly degenerate state at each momentum k, with degeneracy guaranteed by inversion symmetry. In addition, for each state at k, there is an identical state with opposite spin at −k due to time-reversal symmetry. The analysis of the thickness-dependent gap, carried out in the previous section, shows that the limit of an infinitely thick slab with zero interaction between the two surfaces is most likely realized for 40 QLs in our system. Although the size of the gap already for 6 QLs is considerably smaller than the value found in ultra-thin slabs, for 40 QLs the gap becomes identically zero within the numerical precision of our calculations (note that we employ exact diagonalization to calculate the eigen-spectrum of the Hamiltonian at each k). By looking closely at theΓ point, we indeed find four degenerate states ε i (i=1,.., 4) in the case of 40 QLs, as shown in Fig. 7(c). Within each pair of degenerate states, ε 1,3 and ε 2,4 , the states have opposite spins and each pair is localized on either top or bottom surface. Note that the a-th Cartesian component of the spin of each Bloch state is calculated as s a εn (k)=Tr[ρ nn σ a ], where σ={σ a } is the set of Pauli matrices and ρ nn is the n-th diagonal element of the density matrix constructed from the eigenfunctions of the Hamiltonian [Eq. (1)] at momentum k. From Fig. 7(c) we can also estimate the decay length of the surface states to be approximately 10 atomic layers, or, equivalently, 2 QLs in agreement with previous reports. 28 For slabs with thicknesses below 40 QLs, we also find four states at theΓ point, however there is a finite gap between degenerate states ε 1,2 and ε 3,4 [see the inset in Fig. 7(b)]. Within each pair of states, the Bloch states with identical energy have opposite spins but, in contrast to the 40 QLs case, their wavefunctions are distributed over both top and bottom sides of the slab, signaling a non-negligible interaction between the the two surfaces. For 20 QLs, since the difference in energy between ε 1,2 and ε 3,4 is finite but small (less than 10 −6 eV), the spatial distributions of the corresponding wavefunctions are nearly identical and are therefore indistinguishable on the scale of the graph [ Fig. 7(b)]. In the case of 5 QLs, when the gap at theΓ point is more appreciable, there is a clear difference between the spatial profiles of the wavefunctions, associated with the two degenerate states [ Fig. 7(a)]. The difference is the largest in the middle of the slab, where the tails of the wavefunctions of the sur-face states residing on the opposite sides overlap, creating a mixed state.
Finally, we comment on the orbital character of surface states at the Dirac point. The orbital-resolved atomiclayer projections of the wavefunctions associated with one of the four quasi-degenerate states (ε 1 ) at theΓ point in a 20 QLs-thick slab are plotted in Fig. 8. Similarly to the VBM and CBM states at the Γ point of the bulk bandstructure with SO [ Fig. 2(b)], this state is predominantly of p-character with s-orbital contribution less than 10% percent. Although the contribution of p z -orbitals is the largest, the relative weight of in-plane (p x(y) ) orbitals to the orbital-resolved wavefunction is non-negligible and is of the order of 40%. The importance of the SO-induced in-plane orbital contribution to the wavefunction of the Dirac states in Bi 2 Se 3 has been recently demonstrated using orbital-selective SR-ARPES measurements. 22

Spin-resolved surface states on the Dirac cone
In this section we analyze the surface states away from theΓ point. We describe the spatial distribution and the spin-properties of the wavefunctions, associated with energy states found inside the inverted band gap at nonzero values of momentum k. In particular, we quantify the effect of the finite slab thickness on the helical spincharacter of the Dirac cone states.  (a) shows an equi-energy contour located slightly above the Dirac point, within the energy range containing the Dirac cone, for a 5QLs-thick slab. There is a doubly degenerate state at each k on the contour and a corresponding doubly degenerate state at −k with the same energy. We first consider a pair of doubly degenerate states ε 1,2 (k) and ε 1,2 (−k) with k = (0, k y ) as shown in Fig. 9(b). The atomic-layer projections of the wavefunction for all four states are presented is Fig. 9(c). Each of the degenerate states within a pair ε 1,2 (±k) is localized on either top or bottom surface. We will refer to the state whose wavefunction is localized predominantly on the top (bottom) surface as a top (bottom) state.
Furthermore, we calculate the spins of top and bottom states at several momenta on the equi-energy contour. Figures 10(a) and (b) show the projection of the spin on the x-y plane in momentum space for top and bottom states, respectively. The direction of the spin at each k appears to be tangential to the equi-energy contour, which is a manifestation of the spin-momentum locking intrinsic to 3D TIs. 16 The direction of the spin is exactly opposite for the states ε 1 (k) and ε 2 (k), hence the top and bottom surface states have opposite helicities. These observations indicate that despite the opening of the gap at theΓ point due to interaction between the two surfaces for a 5QLs-thick slab, the states appearing across the bulk band gap at this thickness are of topological character. In order to investigate the effect of the slab thickness on the helical spin-texture of the surface states, we determine the orientation of the spin s of the Bloch states ε 1,2 (k) for thicknesses in the range of 5 to 20 QLs. The orientation of vector s in momentum space at each k can be described by two angles: θ, or in-plane tilting angle, which is the angle between s and the normal to the contour at point k, and ψ, or out-of-plane tilting angle, which in the angle between s and the direction perpendicular to the plane [see Fig. 11(c)]. For a perfect spin-momentum locking θ=ψ=90 • . Note that the angles θ and ψ oscillate as functions of k along the equi-energy contour due to the three-fold rotational symmetry of the slab crystal structure. 30 In addition, at each value of momentum the in-plane and out-of-plane tilting angles of the top and bottom states are equal in magnitude and have opposite signs. Hence we compute the deviations of θ and ψ from 90 • only for the top state. The resulting absolute-value deviations, averaged over k on the equienergy contour, ∆θ and ∆ψ , are plotted as functions of the thickness in Figs. 11(a) and (b), respectively. Both in-plane and out-of-plane deviations decrease with increasing the slab thickness. The dependence of ∆θ shows a kink at 6 QLs, consistent with the increase of the gap at theΓ point (see Fig. 5), while ∆ψ decreases monotonically. Importantly, both ∆θ and ∆ψ do not decrease to zero with increasing the thickness but instead saturate to a constant value. The saturation value is small for ∆θ (2 · 10 −4 • ) and relatively large for ∆ψ (0.03 • ). The saturation starts at approximately 10 QLs for ∆θ and at 7 QLs for ∆ψ . The non-zero residual deviation from the perfect helical spin-texture at large thicknesses is a measure of the proximity of the equi-energy contour to bulk states. The further the contour is from the Dirac point, the stronger is the effect of the bulk continuum, with hexagonal wrapping in momentum space due to crystal symmetry, on the surface states located on the contour. For an equi-energy contour located at a higher energy (0.26 eV) we find considerably larger saturation values, namely 0.5 • for ∆ψ and 0.1 • for ∆θ . To summarize, both in-plane and out-of plane deviations decrease with increasing the slab thickness and with moving the equi-energy contour closer to the Dirac point. However in the vicinity of the Dirac point the inplane tilting angle appears to be more affected by the finite thickness while the out-of-plane tilting angle is more sensitive to the presence of bulk states.

IV. CONCLUSIONS
We have employed the sp 3 tight-binding model, with parameters extracted from ab initio calculations, to model the electronic structure of bulk and (111) surface of Bi 2 Se 3 3D TI. We presented a quantitative description of the band inversion mechanism for both bulk and slab geometry, which involves a detectable change in the spatial distribution of p-orbital projections of the wavefunctions of conduction and valence bands induced by spin-orbit interaction. The surface bandstructures, with the spatial character of Bloch states determined by quantitative criteria, were calculated for slabs with thicknesses up to 100 QLs. This thickness is well beyond what is accessible by ab initio approaches. Based of the calculated thickness-dependent gap due to inter-surface interaction and the atomic-projections of the wavefunctions for states at theΓ point, we found that the infinitethickness limit, characterized by zero gap and no interaction between the surfaces, is realized for 40 QLs. Furthermore, our calculations showed that the disturbances in the helical spin-texture of the Dirac-cone states, caused by the finite slab thickness and the proximity to bulk states, persist for thicknesses up to 10 QLs even in the vicinity of the Dirac point.
We would like to comment on our numerical observation that the top and bottom surface states become completely decoupled only for thick slabs containing 40 QLs. Strictly speaking only states with linear dispersion and identically zero gap at the Dirac point can be identified as Dirac states. However, in practice our observation does not contradict the commonly accepted results that 6 QLs is the critical thickness which determines the topological character of the material. 48 Indeed, we found that the gap for thicknesses greater that 5 QLs is very small (of the order of 10 −4 eV and smaller) and the surface states appearing inside the bulk insulating gap at these thicknesses have nearly linear dispersion and helical spintexture. Nevertheless, the presence of a weak but finite interaction between the opposite surfaces for slabs of less than 40 QLs thick might be crucial for the analysis of subtle effects in the vicinity of the Dirac point. Preliminary calculations have shown that even a small perturbation preserving the time-reversal symmetry can open a gap for thicknesses as large as 20 QLs. 57 This may further hinder the identification of the gap induced purely by the presence of time-reversal-breaking perturbations, such as magnetic impurities, 23 in calculations based on tight-binding models similar to ours and certainly in ab initio studies, where thicknesses are limited to only few QLs.
Effects due to finite thickness can be also expected in Landau levels spectroscopy of 3D TI thin films. 58 On one hand, it has been shown experimentally that in Bi 2 Se 3like materials the characteristic field-independent (zeroth) Landau level is absent for slabs with thicknesses smaller than 3 QLs but emerges already for 4 QLs. 56 On the other hand, a recent theoretical study suggests a splitting of the zero-th Landau level due to hybridization between top and bottom surface states. 59 Based on these considerations, one might expect that a non-negligible inter-surface interaction can lead to more subtle internal structure of the zero-th Landau level, persisting even for relatively large thicknesses.
We anticipate that microscopic tight-binding models, combined with input from ab initio calculations, will play an increasingly important role in practical calculations of various properties of TI materials. These include the detailed character of the surfaces states wavefunction, which to some extent can be already probed experimentally, 22 and the interplay between the surfaces states and external perturbations. 24 In fact, a finite-cluster tightbinding approach, based on the model used in the present work, has been already employed in the study of native defects in Bi 2 Se 3 , showing good agreement between the calculated local densities of states around the defects and experimental STM topographies. 60 In connection to this last point, we would like to mention that despite extensive studies of the effect of magnetic doping in 3D TIs, a consistent microscopic description of a single magnetic impurity in a TI environment appears to be incomplete, especially when compared to the progress that has been made in investigating similar questions in semiconductors, both theoretically 35,36 and experimentally. 61,62 A realistic tight-binding approach can be indispensable in providing such a microscopic description for single impurities in 3D TIs.