First Principles Studies on 3-Dimentional Strong Topological Insulators: Bi2Te3, Bi2Se3 and Sb2Te3

Bi2Se3, Bi2Te3 and Sb2Te3 compounds are recently predicted to be 3-dimentional (3D) strong topological insulators. In this paper, based on ab-initio calculations, we study in detail the topological nature and the surface states of this family compounds. The penetration depth and the spin-resolved Fermi surfaces of the surface states will be analyzed. We will also present an procedure, from which highly accurate effective Hamiltonian can be constructed, based on projected atomic Wannier functions (which keep the symmetries of the systems). Such Hamiltonian can be used to study the semi-infinite systems or slab type supercells efficiently. Finally, we discuss the 3D topological phase transition in Sb2(Te1-xSex)3 alloy system.


I. INTRODUCTION
The topological insulators (TIs) [1][2][3][4][5] , are interesting not only because of its fundamental importance but also because of its great potential for future applications [6][7][8] . TI is a new state of quantum matter, and is distinct from simple metal or insulator in the sense that its bulk is insulating (with a bulk band gap), while its surface (or edge) is metallic due to the presence of gapless surface (edge) states. Those surface states are spin splitted but with double degeneracy at the Dirac point, which is protected by the time-reversal symmetry. The surface states, consisting of odd number of massless Dirac cones, are robust against time-reversal-invariant perturbations and also are very different with Graphene where spin degeneracy is reserved (the concept of pseudo-spin describing the sub-lattice degeneracy is used for Graphene). At the early stages of studies for TIs, due to the absence of realistic materials, most of the discussions were based on the model Hamiltonians [1][2][3][4][5] . However, within last couple of years, this field is strongly motivated by the discoveries of several real compounds.
Two-dimension (2D) TI (or called quantum spin Hall insulator) has been proposed and realized in HgTe/CdTe quantum well structure through varying the well thickness [9][10][11] . Threedimension (3D) topological nature of Bi 1−x Sb x alloy has been demonstrated by fine tuning of the alloy concentration 12-14 . In both cases, the bulk band gap is of the order of tens meV, too small for the potential applications. In our recent studies, however, we predicted that Bi 2 Se 3 family compounds 15 (i.e., Bi 2 Te 3 , Bi 2 Se 3 and Sb 2 Te 3 ), which are well-known thermoelectric materials [16][17][18] , are strong topological insulators with surface states consisting of single Dirac cone at the Γ point. Particularly, for Bi 2 Se 3 , a bulk band gap of 0.3 eV is predicted, much larger than the energy scale of room temperature. The compounds are stoichiometric, chemically very stable, and easy to be synthesized, and yet its surface states are very simple. The existence of such surface states in this family compounds have been experimentally confirmed on samples ranging from bulk to thin film [19][20][21][22][23][24] .
More and more studies are now going on for Bi 2 Se 3 family compounds, quantitative studies and comparison with experiments become more and more important. On the other hand, fully self-consistent ab initio calculations are quite lim-ited due to the large system size. In order to study both the bulk and the surface behavior simultaneously, which is necessary to identify its topological nature, a large slab type supercell is needed, and heavy calculations are involved. To simplify future studies, it is therefore desirable to have a highly accurate effective Hamiltonian, which will be one of the main purposes of this paper. It's well-known that the maximally localized Wannier functions (MLWF) method introduced by Marzari and his coworkers 25,26 has played an important role in constructing the effective Hamiltonian. However, the algorithm involved does not hold the symmetry when the MLWF are constructed from Bloch functions. To solve this problem, we introduce a set of projected atomic Wannier functions (PAWF), whose symmetry can be easily controlled, yet keeping high accuracy. With this set of PAWF basis, highly accurate effective Hamiltonian can be constructed. Using this Hamiltonian, we will demonstrate that the topological nature as well as the details of surface states can be all well reproduced.
The organization of the present paper is as follows. In section II, we describe the crystal structure and symmetries of Bi 2 Se 3 family compounds. Section III is devoted to the construction of the effective Hamiltonian in the basis of PAWF. In section IV, this Hamiltonian is applied to Bi 2 Se 3 system, and the topological nature and surface properties are analyzed. In section V, we discuss the 3D topological phase transition by doping Se into Sb 2 Te 3 . Final remarks are summarized in section VI.

II. CRYSTAL STRUCTURE AND SYMMETRIES
The Bi 2 Se 3 family compounds have the rhombohedral crystal structure with space group D 5 3d (R3m), we take Bi 2 Se 3 as an example in the following. As shown in Fig. 1(a), the system has a layered structure with five atomic layers as a basic unit (cell), named a quintuple layer (QL). The inter-layer bonding within the QLs is strong because of the dominant covalent character, but the bonding between the QLs is much weaker due to the van der Walls type interaction. The binary (with twofold rotation symmetry), bisectrix (appearing in the reflection plane), and trigonal (with threefold rotation symmetry) (c) The first Brillouin zone (BZ). Four nonequivalent time-reversalinvariant-momentum (TRIM) points Γ (0, 0, 0), L(π, 0, 0), F(π, π, 0), and Z(π, π, π) are denoted in the 3D BZ. The corresponding surface 2D BZ is represented by the dashed blue hexagon, and Γ, M, K are the corresponding TRIM special k points in the surface BZ.

III. EFFECTIVE HAMILTONIAN FROM PAWF
The construction of the PAWF basis and the effective Hamiltonian is a post-process of ab initio calculations. The fully self-consistent ab initio calculations are performed for the bulk compounds in the framework of density functional theory (DFT) 28 using BSTATE (Beijing Simulation Tool for Atom Technology) package 29 with the plane-wave ultra-soft pseudo-potential method 30 . Generalized gradient approximation (GGA) of the Perdew-Burke-Ernzerhof (PBE) type is used for the exchange-correlation potential 31 . To guarantee the convergence, we use 340 eV for the cut-off energy of wave-function expansion, and 10×10×10 K-points mesh in the BZ. All structure parameters of Bi 2 Se 3 , Sb 2 Te 3 , Bi 2 Te 3 are chosen from experimental data listed in Table I.
After completing the self-consistent ab initio calculations, two steps are followed to construct the PAWF basis and the effective Hamiltonian. First we disentangle an isolated group of bands through minimizing the invariant part of the spread functional of wave functions (Ω I ) 26 . In the second step, the Hilbert space of this isolated group of bands is directly projected onto atomic p orbitals (keeping the angular symmetry). Finally a set of unitary rotation matrices U(k) can be obtained, which can be straightly used to construct the PAWF and the effective Hamiltonian. The details are elucidated in the following.

A. Extracting of isolated group of bands
For a simple band insulator, the disentangling procedure is typically not necessary, and it is straightforward to choose a set of Wannier basis expanding the identical space of occupied Bloch states. However, in our case, aiming to describe the electronic characters near the Fermi level, both occupied and unoccupied bands are necessary, which are not detached from higher bands. By partial density of states (PDOS) analysis, we know that the states near the Fermi level mainly come from the contribution of p orbitals of both Bi and Se atoms.
Since each atom has three p orbitals, the total number of orbitals of interest is N = 15, consisting of 9 valence bands and 6 unoccupied bands. Following the procedure described in MLWF 26 , we introduce a big enough outer energy window (-20 eV, 20 eV), and a small inner windows (-6 eV, 2 eV). A smooth subspace (of N bands) can then be constructed by optimizing the spread functional, namely Ω I , as suggested in Ref. [ 26 ].

B. PAWF basis and effective Hamiltonian
The Bloch functions ψ mk of the isolated group of N bands are directly projected onto the atomic p orbitals, and the unitary rotation matrices U(k) can be obtained by proper orthonormalization 25,32 . With these unitary rotation matrices U(k), we can obtain well-localized PAWF by rotating the original Bloch functions in the following way: where n and m=1,2 · · · N are band index, N k is the total number of k-points, W R n denotes the n-th PAWF orbital, which is centered at the lattice vector R.
Using the PAWF as basis, the effective Hamiltonian H W (k) can be obtained correspondingly by rotating the original Hamiltonian, Because our PAWF are constructed from atomic orbitals, which have atomic symmetry, we can further symmetrize the Hamiltonian using their atomic characters. Using Fourier transformation, the PAWF effective Hamiltonian in real space can be given by The Hamiltonians at any other k-point can be now obtained by transforming H W (R) back into k space, All the hopping matrix elements in the effective Hamiltonian are directly calculated via constructing the PAWF. This approach differs a lot from the conventional Slater-Koster TB method 33 where hopping integrals up to certain neighbors are obtained by a fitting method. The present method has several special advantages. First, it is not necessary to maximize the localization of Wannier functions (WF). Such step will usually break the local symmetry of WF. In our case, the WF could be constrained to satisfy local crystal symmetry, and the PAWF effective Hamiltonian can be symmetrized, while keeping the necessary accuracy. Second, the atomic spin-orbit coupling (SOC) can be easily implemented. Since the WF are constructed from atomic orbitals, they are already quite localized although not maximally localized by definition. The PAWF effective Hamiltonian is expected to be useful for large supercell calculations without extra errors from the symmetry problem.

C. The spin-orbit coupling
The SOC is important for the topological nature of Bi 2 Se 3 family compounds, and should be included in all analysis. Since SOC is mostly local atomic physics and has little kdependence, the simplest way to supplement SOC on top of our PAWF effective Hamiltonian (denoted as H 0 from now on) is to include an additional local term H soc in the total Hamiltonian 34 , where H soc comes from the SOC interaction, while V and σ represent total potential and Pauli spin matrices respectively. As discussed above, our PAWF are constructed from atomic p orbitals (having atomic angular symmetry), the SOC Hamiltonian can be straightforwardly written down, in terms of those PAWF basis, |p x , ↑>, |p y , ↑>, |p z , ↑>, |p x , ↓>, |p y , ↓>, |p z , ↓>, where ↑ (↓) indicates the spin. With this basis set for each atom, H soc part can be expressed as where λ denotes the SOC parameter. We take the SOC parameters of Bi, Se, Sb and Te atoms from Wittel's spectral data (λ Bi = 1.25 eV, λ S e = 0.22 eV, λ S b = 0.4 eV, λ T e = 0.49 eV) 35 .

D. Accuracy of effective Hamiltonian: results for bulk
To demonstrate the quality of the PAWF effective Hamiltonian, here we compare several important properties calculated from the PAWF Hamiltonian with that obtained from the ab initio calculations. The first property to be compared is the the calculated bulk band structure. As shown in Fig. 2(a)  The parity numbers of nine occupied bands and the first conduction band are shown (the corresponding band energy increases from left to right). The parity products of nine occupied bands are given in brackets.
Parity is another important quantity to distinguish the topological characters of Bi 2 Se 3 family compounds. Because the present systems posses the inversion symmetry, the method proposed by Fu and Kane 5 is used here. For all three systems, we calculate the parity eigenvalues of nine occupied bands and the first conduction band at the time-reversal invariant momentum (TRIM) points Γ, L, F, and Z. The results for the Γ point, obtained from the effective Hamiltonian, are listed in Table II. As can be clearly seen, the parity products of occupied bands are -1 for the Γ points. For the other TRIM points, L, F, and Z, the detailed numbers are not listed here, but our results give +1 for three of the compounds. From the parity numbers, we can therefore identify the corresponding Z 2 invariants 1,4,5 as 1 for all the three compounds. This means they have nontrivial topological nature, in good agreement with our previous studies 15 . From above comparisons, we conclude that the quality of the PAWF effective Hamiltonian is as good as that obtained from the ab initio calculations.

IV. THE PROPERTIES OF SURFACE STATES
Once the effective Hamiltonian has been constructed from above procedures, we will be able to calculate the surface states of semi-infinite systems from the iterative Green's function method 36 , as described in our previous studies 15 . The existence of gapless spin-filter surface states is the direct man- ifestation of the topological nature. The crystal structure of Bi 2 Se 3 family compounds can be understood as the stacking of QLs along the z direction. The inter-layer bonding between two QLs is much weaker than that inside the QL, it is natural to expect that the cleavage plane should be between two QLs. This fact has been well confirmed by the recent experiments on the layer-by-layer MBE growth of ultra-thin film [21][22][23][24] . We therefore focus our studies on this type of surface termination (with Se1 atomic layer as the top most layer).
Our previous studies have shown that gapless surface states with single Dirac point (located at Γ) exist for all three compounds, Bi 2 Se 3 , Bi 2 Te 3 and Sb 2 Te 3 , but not for Sb 2 Se 3 . The dispersion of surface states around Dirac point is highly linear, and these surface states can be described by a simple continuous model 15 . Here, we will not repeat those results, on the other hand however, will focus on some other details of surface states, including the penetration depth, the spin-resolved Fermi surface, and the chiral spin texture of the surface states.

A. The penetration depth of surface states
The spread of surface states, or in other words its spacial penetration depth into the bulk, is an essential quantity for the potential applications of these surface states if any. In order to clarify this point, a free standing slab model is constructed using our PAWF effective Hamiltonian. We use a slab consisting of 25 QLs, which is thick enough to avoid the direct coupling between the two surfaces, i.e, the top and the bottom surfaces of the slab. In order to make the calculations more realistic, we further take into account the surface correction. The bulk PAWF effective Hamiltonian is used to construct the Hamiltonian of the slab, however for those layers close to the surface, the on-site energies of PAWF should be modified from its bulk counterparts due to the presence of vacuum (i.e, the surface potential). By comparing with the fully self-consistent ab initio calculations of thin slab, all the correction of on-site energies for PAWF (∆E n = E sur f ace n − E bulk n ) can be obtained.
The calculated band structures of our 25-QLs Bi 2 Se 3 slab along Γ → K and Γ → M lines are shown in Fig. 4 and Fig. 5. Clearly, within the bulk energy gap, there exist two surface bands (marked by Σ 1 and Σ 2 ), which are degenerate at Γ point (where so called Dirac cone exists). The two surface bands (Σ 1 and Σ 2 ) are sampled by several k-points (with k i , i=1,2,3). The real space distribution of eigen wave functions for those sampling k-points are plotted in Fig. 4(b)-(c) and Fig. 5 the inset. Finally, we can conclude that surface states are very localized to the surface region, and the penetration depth is about 2 or 3 QLs.
To make further comparison, we have performed a fully self-consistent ab initio calculations for a 3-QLs slab. The real space charge distribution of surface states can thus be obtained, by integrating the states within an energy window of 10 meV around the Dirac point. As shown in Fig. 6, the charge distribution is really localized at the surface region. Further more the charges are not centered at the atom, but rather distributed mostly between the Bi and Se atoms, where strong covalent bonding is expected.

B. The spin-resolved Fermi surface
Probing the π Berry phase enclosed by the Fermi surface of surface states is one of the most direct methods to distinguish TI. Now we will analyze the spin-resolved Fermi surfaces around the Dirac point of the semi-infinite Bi 2 Se 3 system. With the Green's functions calculated from the effective Hamiltonian, the spin-filter surface states and the corresponding Fermi surfaces can be obtained directly.
As shown in Fig. 7, when the Fermi level is close to the Dirac point, the corresponding Fermi surface is nearly a perfect circle; while if the Fermi level is away from the Dirac point, the properties of the surface states are significantly affected by the bulk states and thus satisfy the crystal symmetry. For example, the Fermi surface for E f =0.1eV (shown in Fig. 7(a)) looks like a small circle, while the Fermi surface for E f =0.25eV (shown in Fig. 7(b)) looks like a hexagon satisfy- ing c 3z symmetry. The spin orientation for states around the Fermi surface is marked as well by arrows. The magnitude of the spin z component is very small, only about 1/40 of the in-plane component. Thus the spin almost completely lies in the plane. Moving around the Fermi surface, the spin orientation rotate simultaneously, forming a spin-orbit ring, which carries a π Berry phase. This signifies the topological nontrivial properties of Bi 2 Se 3 . In addition, as shown in Fig. 7, the spin orientation of the ring belongs to the left chirality (the normal direction of the semi-infinite system is defined as z direction) which is again one of the important manifestations of nontrivial topological characters.

C. The linear dispersion of surface states
The gapless surface states that connect the bulk valence and conduction bands have almost linear energy dependence near the Dirac point. Linear band structure in 2D should lead to linear density of states (DOS). To be specific, in the bulk energy gap, the energy bands are mostly from the surface states, so the DOS in the bulk energy gap are expected to be highly linear. Fig. 8(a) shows the band structures of a 3-QLs slab. Its corresponding DOS (presented in Fig. 8(b)) shows very nice linear energy dependence within the bulk energy gap as expected. This type of DOS can be easily measured by low temperature scanning tunneling spectroscopy (STS), which will provide an indirect method to probe the existence of linear surface states.

V. 3D TOPOLOGICAL PHASE TRANSITION
The above discussions presented in sections II and III concentrate on the effective Hamiltonian and the detailed understanding to the properties of surface states. By such a way, we can analyze the topological nature of specific compound. On the other hand, however, the topological nature can be also understood from simple bulk studies, where a gap closereopening transition (a topological phase transition) can be obtained by tuning some parameters 37 . We will show in this section that such a phase transition can be indeed obtained in Sb 2 (Te 1−x Se x ) 3 alloy system.
The virtual crystal approximation (VCA) method proposed by Bellaiche 38 is used here to simulate the doping process. The solid solution elements are treated as different atomic species with their own weights, rather than being treated as a whole virtual atom. With this method, we investigate the Se doping dependence of band structure of Sb 2 (Te 1−x Se x ) 3 . We assume that crystal parameters and positions of inner atoms change linearly with doping, and they are obtained by linear interpolation in calculation. The evolution of the gap energy of Sb 2 (Te 1−x Se x ) 3 as a function of Se concentration is calculated and illustrated in Fig. 9. With increasing Se doping, the valence band maximum and conduction band minimum get closer gradually (0 < x ≤ 0.94), attributed to the gradually decreased SOC strength. Consequently, the two bands cross at the critical point x = 0.94, resulting in a 3D topological quantum phase transition. An band-inverse appears when the doping concentration further increases. Therefore we can conclude that the occurrence of 3D topological phase transition is driven by the SOC interaction. In reality, the most stable crystal structure of Sb 2 Se 3 is slightly distorted from the rhombohedral structure of Sb 2 Te 3 , however, as long as the Sb 2 Se 3 is topologically trivial, such 3D topological phase transition will be expected. At the phase transition point (x=0.94 as predicted in our calculations), the 3D Dirac cone may be expected in the bulk band structure 37 .