---
abstract: 'Using a density functional based electronic structure method, we study the effect of perturbations on the surface state Dirac cone of a strong topological insulator Bi$_2$Se$_3$ from both the intrinsic and extrinsic sources. We consider atomic relaxations, and film thickness as intrinsic and interfacial thin dielectric films as an extrinsic source of perturbation to the surface states. We find that atomic relaxations has no effect on the degeneracy of the Dirac cone whereas film thickness has considerable effect on the surface states inducing a gap which increases monotonically with decrease in film thickness. We consider two insulating substrates BN and quartz as dielectric films and show that surface terminations of quartz with or without passivation plays critical role in preserving Dirac cone degeneracy whereas BN is more inert to the TI surface states. The relative orbital contribution with respect to bulk is mapped out using a simple algorithm, and with the help of it we demonstrate the bulk band inversion when spin-orbit coupling is switched on. The layer projected charge density distributions of the surface states shows that these states are not strictly confined to the surface. The spatial confinement of these states extends upto two to three quintuple layers, a quintuple layer consists of five atomic layers of Bi and Se.'
author:
- Jiwon Chang
- Priyamvada Jadaun
- 'Leonard F. Register'
- 'Sanjay K. Banerjee'
- Bhagawan Sahu
title: 'Intrinsic and extrinsic perturbations on the topological insulator Bi$_2$Se$_3$ surface states'
---

Introduction
============

Three dimensional (3D) topological band insulators (TI), Bi$_2$X$_3$, X=Se, Te and Sb$_2$Te$_3$ have attracted considerable attention from the condensed matter physics community because of the relatively simple crystal structure that hosts its novel surface states[@zhang1]. Many more 3D TI materials have been predicted[@new3d] by now and quest for studying their novel surface state properties has increased dramatically. The surface state of 3D TI are time-reversal symmetric (TRS) in some regions of momentum space and are therefore protected against perturbations which cannot break TRS. These considerations also suggest the possibility of dissipationless transport within these surface states where it would be relatively insensitive to non-magnetic disorder or any perturbations that protect the TRS properties of these states[@zahid]. Such properties have caused excitement in the electron device community. The advances in undertanding of structural, electronic, magnetic and transport properties of 3D TI, made possible by both experimental and theoretical studies[@zhang2], have allowed it to gain attention in the scientific community[@star1]. It is, however, not clear whether perturbations, both intrinsic and extrinsic, can affect the surface state Dirac cone, by breaking the TRS. Therefore, it is necessary to study 3D TI and the novel surfaces states, using [*ab-initio*]{} methods. To understand the effects of these variables on the surface states, we report on the density functional based electronic structure studies of one such strong 3D TI Bi$_2$Se$_3$. Our study suggests a critical thickness of six quintuple layers (QLs) needed to maintain the degeneracy of the electron-hole bands at the Dirac point and the surface states have spatial extension to within 2$\sim$3QLs. Atomic relaxations have no effect on the linear spectrum but thin dielectric films (crystalline BN and SiO$_2$) are found to play a critical role in breaking the surface state degeneracy, depending upon the nature of surface terminations and passivations. Our studies will be useful for designing electronic devices using a 3D TI and will have implications in interpreting experiments.

We begin by describing the bulk crystal structure of Bi$_2$Se$_3$ and the computational method used for this study in section II. In section III, we present the electronic structure of bulk TI and a simple algorithm to map out the surface states from the bulk band structure, and we use this method to explain band inversion when spin-orbit coupling is switched on. We then discuss the thin-film structure and the thickness dependent surface band structure of Bi$_2$Se$_3$ in section IV. The spatial extension of these surface states are then determined using charge densities projected onto each atomic layer and each atom. The effect of atomic relaxations and thin dielectric films on the surface electronic structure will be discussed in section V. Finally we present our summary and conclusions.

Bulk Structure and Computational Method
=======================================

The bulk crystal structure of Bi$_2$Se$_3$ is trigonal (or rhombohedral) with lattice constant of [*a*]{}= 0.984 nm and $\alpha$$\sim$25$^o$. The primitive cell consists of five atoms (two Bi and three Se atoms) arranged in a order Se2-Se1-Bi-Bi-Se1 where Se1 and Se2 are two inequivalent sites of Se atoms in the cell. Figure 1(a) shows a cartoon of bulk rhombohedral structure along with three primitive lattice vectors that span it (the tracing inside the large cell). The three-fold rotation axis is taken along the [*z*]{}-axis. Alternatively, one can construct a hexagonal cell from the trigonal cell so that the number of atoms is tripled with the lattice parameters [*a*]{}=0.4138 nm and [*c*]{}=2.8633 nm (Fig. 1(a), the larger cell). The building block of the hexagonal bulk Bi$_2$Se$_3$ crystal consists of five atomic layers, each layer containing only one atom. One such building block is referred to as quintuple layer (QL) (square shaded region in Fig. 1(a)). The hexagonal unit cell contains three such QLs, i.e., 15 atomic layers stacked along the [*z*]{}-direction. The atomic planes are arranged in a sequence Se1-Bi-Se2-Bi-Se1, which is different from the trigonal case, with the following interlayer distances: [*d*]{}(Se1-Bi)= 0.156344 nm, [*d*]{}(Bi-Se2)=0.190665 nm, [*d*]{}(Se2-Bi)=0.190665 nm, [*d*]{}(Bi-Se1)=0.156353 nm and [*d*]{}(Se1-Se1’)=0.259304 nm where Se1’ is the Se1 site in the next unit cell. Because of the relatively larger distance between two QLs, their anchoring along the [*z*]{}-direction is considered to be rather weak. Indeed, mechanical exfoliation of Bi$_2$Se$_3$ demonstrated recently[@balan1] hints at van der Waal’s type of bonding between two neighboring unit cells, much like a graphene exfoliation.

We use a DFT-based electronic structure method implemented with projector-augmented wave basis and pseudopotentials[@vasp1]. Since the topological features of the surface states arise from spin-orbit coupled states of the bulk, we invoke its spin-orbit coupling (SoC) feature[@vasp2] and we choose exchange-correlation as Perdew-Burke-Ernzerof (PBE) type[@perdew]. The choice of the PBE over local-density approximation is guided by its use in previous structural studies on bulk Bi$_2$Se$_3$ which reproduced the photoemission spectra very well[@zhang2]. These studies find the optimized lattice constant of Bi$_2$Se$_3$ close to the experimental value, therefore we did not optimize lattice constants of the bulk structure in our calculations. A kinetic energy cut-off of 25 Ry and a [**k**]{}-point mesh of 10 $\times$ 10 $\times$ 10 in the first Brillouin Zone (BZ) (Fig. 1(b)) are used for self-consistent steps as well as generating the bulk band structures. Larger energy cut-offs and [**k**]{}-mesh size were tested for total energy convergence. The first BZ of the bulk structure contains four time-reversal invariant (TRI) points, namely, $\bf \Gamma$, [**L**]{}, [**Z**]{} and [**F**]{} (Fig.1(b)).

Bulk band structure
===================

In this section, we address the band structure of bulk Bi$_2$Se$_3$, with and without SoC. We demonstrate the band inversion with switching on of SoC with the help of an algorithm which maps out orbital contribution to the electronic states from bulk bands. Figure 2(a) and (b), respectively, shows energy dispersion curves for the bulk Bi$_2$Se$_3$ without and with SoC, along high-symmetry directions connecting the four TRI points (Fig. 1(b). The direct band gap, in both cases, forms at the $\Gamma$-point and is 200 meV(285 meV) without(with) SoC. These gaps compare well with those obtained in other DFT simulations[@zhang2]. It was suggested theoretically, using a model Hamiltonian based on atomic orbitals and without SoC, that the valence band maximum (VBM) and conduction band minimum (CBM) mainly consists of Se and Bi $p_z$ orbitals[@zhang3], respectively (Fig. 2(a)). Switching on the SoC inverts the bands at the $\bf \Gamma$-point (Fig. 2(b)) and consequently the contributions from Bi and Se ions associated with these bands.

\

To understand the band inversion at the $\Gamma$-point visually, we adopt a simple algorithm to map out the most relevant orbital contribution to VBM and CBM[@park]. Since $p_z$ orbitals mainly contribute to these bands near the Fermi level, the algorithm projects the crystal wave-functions on to the angular momentum dependent spherical harmonics at a particular [**k**]{}-point with the associated band energies for each atom. The primitive unit cell of bulk Bi$_2$Se$_3$ contains three Se and two Bi atoms and using the above algorithm, we obtain [*s*]{} and three [*p*]{} orbitals contribution for each of the atom type at each [*k-point*]{} and band energy. We take the ratio of Se and Bi $p_z$ orbitals contributions compared to the total orbital contributions resulting from all the orbitals of five atoms at a particular [**k**]{}-point and the band energy. If this ratio is 40$\%$ (60$\%$) or greater for Bi (Se), we take the $p_z$ orbital contribution to that band from Bi or Se at that [**k**]{}-point with the associated band energy; otherwise, we neglect it. The choice of the percentage cut-off is guided by the number of each type of atoms in the unit cell: two for Bi versus three for Se. We find that changing the cut-off or critical percentage ratio makes the contributions to each band at each [**k**]{}-point different but the overall contribution to the VBM and CBM at the $\Gamma$ point remains the same. In this sense, the choice of critical percentage is arbitrary. We use the above algorithm both without and with SoC calculations, and as a result we could see the band inversion as seen in Fig. 2(b). In the literature, this inversion process is termed as phase transition from a trivial insulator (Bi$_2$Se$_3$ without SoC) to a topological insulator.

\

Thin film band structures and surface state localization
========================================================

The most interesting feature emerges when thin-film structures are considered. A gapless linear spectrum with continuum in the bulk is evident at the TRI point $\Gamma$ (for example Figure 4(e) for a thin film consisting of 6QLs). These states are contributed by both the surfaces of the thin film and are spin-polarized. The thin film structure is constructed by stacking up the quintuple layers along the [*z*]{}-direction with a vacuum region that forms the supercell in the DFT calculation (Fig. 3(a)). We consider several thicknesses of the thin film structure: 2QLs, 3QLs, 4QLs, 5QLs and 6QLs to study the robustness of the linear energy dispersion and for each thickness we used a vacuum size of about 3 nm. We did not relax the atomic positions. The relaxations and their effect will be considered later. Kinetic energy cut-off of 25 Ry and 10 $\times$ 10 $\times$ 1 [**k**]{}-point mesh in the surface BZ of the hexagonal Bi$_2$Se$_3$ (Fig. 3(b)) is used. On decreasing the thickness from 6QLs to 5QLs, a gap seems to be induced at the degenerate point of the linear spectrum at $\Gamma$ (Fig. 4(d)). The magnitude of the induced gap increases monotonically going from 5QLs to 4QLs to 3QLs and finally to 2QLs (Figs. 4(a-c)). These gaps are tabulated in Table I. This suggests a critical thickness of at least 6 QLs in experiments in order to access the novel surface states of Bi$_2$Se$_3$. We note that thickness dependent gaps in Bi$_2$Te$_3$ were realized experimentally[@Li] and predicted theoretically[@park]. In Bi$_2$Se$_3$, we are aware of only one theoretical work which suggests thickness dependent electronic structure[@liu].

To understand the origin of the induced gap in the linear spectrum, we first map out all the orbital contributions to the bulk band structure from all the atoms in first few layers of each thickness, from the top and bottom region of the thin film (Figs. 4(a)-(e)). For thickness of 3QLs-5QLs, the contributions from the top or bottom one quintuple layer is chosen whereas for 2QLs system, we consider the top or bottom two atomic layers for estimating the contributions. The contributions are calculated using the same procedure as explained in Section III with two exceptions: 1) The sum of contributions from all the orbitals and atom types are considered (not just from the $p_z$ orbitals) and 2) the critical percentages, for accepting it as a contribution, are different for different thicknesses. The critical percentage for 2QLs is fixed at 30 $\%$ whereas for 3QLs-6QLs, percentages are respectivley, 60$\%$, 55$\%$, 50$\%$ and 45$\%$. The choice of these critical percentages are arbitrary. We find that choosing different percentages do not change the overall features in the band structures. We note that thinner structures have larger critical percentage for accepting the contributions with reduction in percentages for thicker structures. This is motivated by the fact that thicker structures have more QLs than the thinner ones and we fixed the contributions to only 1 QL for all thicknesses, except for 2 QL structures where two atomic layers are considered to be a cut-off limit. So keeping the same critical percentage for all thicknesses will offset the number of quintuple layers to be considered for acceptance as a contribution and different thickness will need different number of quintuple layers. Alternatively, one can fix the critical percentage for all the thicknessess to the same value and let the number of quintuple layers that contribute to the bulk band structure differ from one thickness to another. Either consideration leads to the same conclusion, namely the contribution to the Dirac cone near $\Gamma$ originates from the first few layers of the thin film.

To locate the critical number of quintuple layers upto which surface states contributing to the Dirac cone is spread in real space, we use a method in which atomic layer dependent charge density from the surface wave functions is calculated. We choose an energy window of 50 meV around the Fermi level in 4QLs-6QLs bandstructures (Figs 4(c)-(e)) with the states in the neighborhood of the $\Gamma$-point and plot the layer projected relative charge density (Figs. 5(a)-(c)). The relative charge density is defined to be the ratio of maximum charge density in a given layer to the maximum among the maximal charge densities of all the layers in that structure. Our calculations hint at 2-3 QLs spatial extension of the surface states in 4QLs-6QLs. The charge density distribution for 2QLs and 3QLs is not shown because within the energy window of 50 meV considered for this calculations, bulk as well surface states both contribute and it is not a pure surface state contribution (Figs. 4(a) and (b)). The surface state localization length to within 2QLs-3QLs in real space hint at possible interactions due to their overlap in thinner structures. This surface state interactions result in opening of the band gap. However, the time-reversal symmetry is not broken at the $\Gamma$-point; the four-fold degenerate band splits into two-fold degenerate bands. The overlap in thinner structures results in finite charge density in a region where surface states for thicker structures had no densities (Figs. 5(a) and (b)). As the degree of overlap increases with decreasing thickness, the induced gap increases monotonically, as seen from Table I. The surface state localization deep inside the bulk region has implications for experiments where measured charge mobility can only be assigned to both suface and bulk states and it will be a challenge in experiments to separate these two contributions.

   Thickness (QL$\sim$1nm)     2        3        4        5        6
  ------------------------- -------- -------- -------- -------- --------
          Band gap           0.1281   0.0430   0.0155   0.0051   0.0000
      at $\bar{\Gamma}$                                         

  :  Induced band gaps (in eV) for thin-films with various thicknesses at the time-reversal invariant point $\bar{\Gamma}$.

\
\

Effects of Atomic Relaxations and the Dielectric films
======================================================

As the Bi$_2$Se$_3$ thin film is contructed from the bulk crystal, atomic relaxations are inevitable. As a result, both Dirac cone degeneracy and surface state localization length may be affected. Since the 6QLs electronic spectrum has no gap, it is natural to consider this structure for the relaxation study. We considered the same computational parameter as the unrelaxed case but allowed only the [*z*]{}-component of the atom to relax. Since a vacuum along [*z*]{}-direction is present, atomic motions are more likely near the vacuum region than within the rigid planar structure. The total energy is assumed to have converged when the [*z*]{}-component of the Hellman-Feynman force is smaller than 0.025 eV/Å. We find that the Dirac cone degneracy is not perturbed with atomic relaxations (Fig. 6(a)) and overall band structure features are retained except for small rearrangements of bands throughout the spectrum. This result suggests the robustness of the topologically protected surface states with respect to atomic relaxations. The extraction of surface states for each band and each [**k**]{}-point is performed as explained in Section III, and the results are denoted by circles (Fig. 6(b)). The charge density distribution for each layer is determined in a similar way to that in section IV. Figure 6(c) shows the relative charge density of the relaxed structure compared to the non-relaxed case. The overall trends exhibited for the non-relaxed case is maintained, but the charge density seems to be enhanced in the second QL from each side of the surface. This hints at strong spread of the surface states in the second QL as compared to the non-relaxed case.

We considered thin dielectric films of crystalline boron nitride (BN) and quartz (or SiO$_2$) to study their effect on surface states. We note here that amorphous SiO$_2$ is used in practical devices so our results for quartz are at best qualitative. We first discuss the interface structure of BN and Bi$_2$Se$_3$ and its effect on surface band structure of Bi$_2$Se$_3$. BN crystallizes in hexagonal structure with [*a*]{} = 0.2494 nm and [*c*]{}= 0.333 nm with d(B-N) =0.144 nm[@CrystalBN]. The bulk direct band gap is 5.97 eV[@kanda]. We use GGA for Bi$_2$Se$_3$ studies. However, use of GGA to optimize the bulk positions and lattice parameters for BN results in overestimated value of lattice constants and it underestimates the bandgap by 33$\%$ [@brink]. Therefore, for the interface studies, we used experimental values of bulk lattice parameters and positions. With these experimental values, our calculated bulk band gap value is same as fully optimized DFT calculations suggest[@openmx].

Our analysis suggests that to keep our computational burden in DFT-based calculations minimal, 3 $\times$ d(B-N) lattice structure can be matched in-plane with 1 $\times$ 1 Bi$_2$Se$_3$ cell or any multiple of this combination that can maintains 3:1 lattice size ratio resulting always in 4$\%$ natural strain. Any other combinations result in either difficulty in forming the periodic structure or large supercell-size with increase of number of atoms by orders of magnitude to allow strain less than 4$\%$. Along the stacking direction, 6 QLs of Bi$_2$Se$_3$ ($\sim$ 6 nm), which has no band gap for surface states, is put on six layers of BN ($\sim$ 1.7 nm). The choice of 6 BN layer is somewhat arbitrary and guided by the fact the size of vacuum and BN layers should be thick enough to avoid interactions of periodically repeated Bi$_2$Se$_3$ surface layers. Since the atomic displacements in Bi$_2$Se$_3$ play no role in breaking the surface state degeneracy, we strained the BN lattice, instead of the TI lattice, to fit with Bi$_2$Se$_3$ surface resulting in compressive strain. We did not relax the interfacial atomic positions to check whether the strained BN lattice can affect the Dirac cone at the $\Gamma$ point. We considered four configurations of Se positions with respect to boron and nitrogen positions in BN: Se on the top of B, on the top of N, on the hexagonal hole and in the interstitial region between B and N. We note here that BN layers cannot be put on both sides of Bi$_2$Se$_3$ film for Se on the top of B configuration because the same orientations are not maintained between the positions of Se and B (Fig. 7). This occurs because of need to satisfy simultaneously the symmetry of Bi$_2$Se$_3$ layers and BN layers in forming the supercell structure. The same holds true for Se on the top of N atom. However, for Se on the top of the hexagonal hole or at the interstitial region between B and N, the orientations are same on both sides (Fig. 8). For the sake of consistency, we calculated the band structures of the interface structures with BN only on one side of Bi$_2$Se$_3$ for all configurations. For the interface simulations, we used 7 $\times$ 7 $\times$ 1 [**k**]{}-point mesh and 2 nm thick vacuum layer. We adjusted the separation between the Se layer and the BN layer. We found the optimal separation of 0.3 nm, and all the studied configurations of Se with respect to BN layer are energetically close and give quite similar band structures. We chose one of these configurations namely Se on B for further discussions. Figure 9 shows the band structure of Bi$_2$Se$_3$ on BN at the optimal interfacial separation of 0.3 nm. It is evident that BN has negligible effect on the surface state dispersion of Bi$_2$Se$_3$ and the Dirac cone degeneracy at the $\Gamma$-point is preserved with no change in the band gap of BN layers. However, the bulk band gap value of Bi$_2$Se$_3$ in the slab structure changes slightly.

We now discuss our results for surface state dispersion in presence of crystalline quartz (or $\alpha$-SiO$_2$). We put quartz on both sides of the TI film. The hexagonal crystal structure of quartz contains four-fold coordinated oxygens, forming a layered structure with Si. The lattice parameters are [*a*]{}=0.4914 nm and [*c*]{}=0.5408 nm[@quartz]. Our DFT calculations of optimized lattice paramters with PBE potentials are found to be close to the experimental values, accurate to within 0.1$\%$. Therefore, in this section we take experimental lattice parameters for building our interface structure of TI and quartz. The bulk amorphous quartz has direct band gap of 8.9 eV[@sio2gap] and our DFT calculation on crystalline SiO$_2$ results in a bandgap value of 6 eV. The interface supercell structure consists of six quintuple layers of Bi$_2$Se$_3$ and two unit cells of SiO$_2$ stacked along the [*z*]{}-axis. The size of the SiO$_2$ and the Bi$_2$Se$_3$ cell along [*x*]{}-[*y*]{} direction chosen is, respectively, 1 $\times$ 1 and 2 $\times$ 2 to keep the computational burden minimal. This produces about 2.75$\%$ compressive strain on both sides of the TI surface. Consideration of larger sizes results in lower strain but the total number of atoms in the cell increases atleast an order of magnitude (250 versus 2500). We did not relax the interfacial atomic positions to check whether the resulting strain can affect the TI Dirac cone. We use 7 $\times$ 7 $\times$ 1 [**k**]{}-point mesh for the BZ integration. Periodicity is assumed along [*x*]{}- [*y*]{} direction. We performed electronic structure calculations at various interplaner distances between the quartz and Bi$_2$Se$_3$ to fix the optimal distance at which the dielectric oxide effects can be assessed. We consider two surface terminations for quartz: Si- and O-termination each with and without hydrogen passivation of dangling states (Fig. 10). Without hydrogen saturation in O-terminated quartz, the Dirac cone of TI is destroyed whereas in Si-terminated quartz, it is preserved. The interplaner separation at the interface region is found to be 0.3 nm for Si-terimated quartz without passivation. Hydrogen saturation reduces the optimal separation to 0.25 nm in O-termination but the optimal distance is insenstitive to the hydrogen saturation in Si-termination. Figures 11 show band structures of the TI in presence of oxygen terminated quartz with hydrogen passivation. Similar band structure is obtained with Si-termination (Figure not shown). The Dirac cone is seen to be insensitive to both the terminations with hydrogen saturation. The band gaps of thin films of quartz with hydrogen saturated oxygen and with Si terminations are calculated to be 9.23 and 8 eV, respectively. In presence of TI, gaps of each of these terminations reduce but the bulk gaps of TI remain unchanged.

Summary and Conclusions
=======================

We use a density functional based electronic structure method to study the effect of intrinsic and extrinsic perturbations on the linear spectrum of a strong topological insulator Bi$_2$Se$_3$. Narrow Bi$_2$Se$_3$ film thickness, atomic relaxation, and applied dielectric films of BN and quartz are considered as perturbations. The thickness of the Bi$_2$Se$_3$ film has a considerable effect on the surface state Dirac cone, inducing a gap due to overlap of surface states orginating from two sides of the film. We estimated surface state localization length, by determining the layer projected valence charge density from the states with energies spanning few tens of meV around the Fermi level. The localization length is found to be within 2-3QLs from each side of the film. With increasing overlap between the surface states in thinner structures, the induced gap increases monotonically. We map out the atom and orbital projected bands from the crystal wave-functions. We then extract the surface state contributions to the thin-film band structure and explain the band inversion in bulk TI when spin-orbit interaction is included. Our study hints at insensitivity of the Dirac cone to both type of perturbations. We hope that our studies will promote studies of interplay between surface states and other perturbations.

The authors acknowledge financial support from SWAN-NRI. We thank Texas advanced computing center (TACC) for computational support.

[999]{} For a popular review see X.-L. Qi and S.-C. Zhang, Physics Today 63, 33 (2010); X.-L. Qi and S. -C. Zhang, arXiv:1008.2026; M. Z. Hasan and C. L. Kane, arXiv:1002.3895. S. -Y. Xu, L.A. Wray, Y. Xia, R. Shankar, S. Jia, A. Fedorov, J.H. Dil, F. Meier, B. Slomski, J. Osterwalder, R.J. Cava, M.Z. Hasan, arXiv:1008.3557; J. Wang, R. Li, S.-C. Zhang, X.-L. Qi, arXiv:1008.2666; B. Yan, H. -J. Zhang, C. -X. Liu, X. -L. Qi, T. Frauenheim, S. -C. Zhang, Phys. Rev. B [**82**]{}, 161108 (2010); D. Xiao, Y. Yao, W. Feng, J. Wen, W. Zhu, X. Q. Chen, G. M. Stocks, Z. Zhang, Phys. Rev. Lett. [**105**]{}, 096404 (2010); W. Feng, J. Ding, D. Xiao, Y. Yao, arXiv:1008.0056; H. Jin, J. -H. Song, A. J. Freeman, M. G. Kanatzidis, arXiv:1007.5480; S. -Y. Xu, L. A. Wray, Y. Xia, R. Shankar, A. Petersen, A. Fedorov, H. Lin, A. Bansil, Y. S. Hor, D. Grauer, R. J. Cava, M. Z. Hasan; arXiv:1007.5111; H. Lin, R.S. Markiewicz, L.A. Wray, L. Fu, M.Z. Hasan, A. Bansil, Phys. Rev. Lett. 105, 036404 (2010); Y. Chen, Z. Liu, J. G. Analytis, J. -H. Chu, H. Zhang, S. -K. Mo, R. G. Moore, D. Lu, I. Fisher, S. -C. Zhang, Z. Hussain, Z. -X. Shen, arXiv:1006.3843; H. Lin, L.A. Wray, Y. Xia, S.-Y. Xu, S. Jia, R.J. Cava, A. Bansil, M.Z. Hasan, arXiv:1004.0999; S. Chadov, X.-L. Qi, J Kübler, G. H. Fecher, C. Felser, S.-C. Zhang, Nature Mater. [**9**]{}, 541 (2010) ; H. Lin, L.A. Wray, Y. Xia, S. Jia, R.J. Cava, A. Bansil, M.Z. Hasan, Nature Materials (2010); B. Yan, C. -X. Liu, H.-J. Zhang, C.-Y. Yam, X. -L. Qi, T. Frauenheim, S. -C. Zhang, Euro. Phys. Lett. [**90**]{}, 37002 (2010). O. E. Ayala-Valenzuela, J .G. Analytis, J. -H. Chu, M. M. Altarawneh, I. R. Fisher, R. D. McDonald, arXiv:1004.2311, Z. Alpichshev, J. G. Analytis, J.-H. Chu, I.R. Fisher, A. Kapitulnik, arXiv:1003.2233; D. Hsieh, Y. Xia, D. Qian, L. Wray, F. Meier, J. H. Dil, J. Osterwalder, L. Patthey, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, M. Z. Hasan, Phys. Rev. Lett. 103, 146401 (2009), D. Hsieh, D. Qian, L. Wray, Y. Xia, Y.S. Hor, R.J. Cava, M.Z. Hasan, Nature 452, 970-974 (2008). H. Zhang, C. -X. Liu, X. -L. Qi, X. Dai, Z. Fang and S.-C. Zhang, Nature Physics [**5**]{}, 438 (2009).A.A. Burkov, D.G. Hawthorn,PRL 105, 066802 (2010); Y. S. Hor, J. G. Checkelsky, D. Qu, N. P. Ong, R. J. Cava, arXiv:1006.0317; A. S. Núñez, J. Fernández-Rossier, rXiv:1003.5931; R. Zitko, Phys. Rev. B 81, 241414(R) (2010); J. G. Checkelsky, Y. S. Hor, R. J. Cava, N. P. Ong, arXiv:1003.3883; N. P. Butch, K. Kirshenbaum, P. Syers, A. B. Sushkov, G. S. Jenkins, H. D. Drew, J. Paglione, Phys. Rev. B 81, 241301(R)(2010); Y. S. Hor, P. Roushan, H. Beidenkopf, J. Seo, D. Qu, J. G. Checkelsky, L. A. Wray, D. Hsieh, Y. Xia, S.-Y. Xu, D. Qian, M. Z. Hasan, N. P. Ong, A. Yazdani, R. J. Cava, Phys. Rev. B 81, 195203 (2010); S. R. Park, W. S. Jung, Chul Kim, D. J. Song, C. Kim, S. Kimura, K. D. Lee, N. Hur, Phys. Rev. B [**81**]{}, 041405 (2010). L. Wray, S. Xu, J. Xiong, Y. Xia, D. Qian, H. Lin, A. Bansil, Y. Hor, R.J. Cava, M.Z. Hasan, accepted for publication in Phys. Rev. B(R). G. Brumfiel, Nature 466, 310-311 (2010). D. Teweldebrhan, V. Goyal, M. Rahman and A. A. Balandin, Appl. Phys. Lett. [**96**]{}, 053107 (2010); D. Teweldebrhan, V. Goyal and A. A. Balandin, Nano Lett. [**10**]{}, 1209 (2010); K. M. F. Shahil, M. Z. Hossain, D. Teweldebrhan and A. A. Balandin, Appl. Phys. Lett. [**96**]{}, 153103 (2010). G. Kresse and J. Furthmuller, Comp. Mater. Sci. [**6**]{}, 15 (1996); G. Kresse and J. Furthmuller Phys. Rev. B [**54**]{}, 11169 (1996); G. Kresse and D. Joubert Phys. Rev. B [**59**]{}, 1758 (1999). D. Hobbs, G. Kresse, and J. Hafner, Phys. Rev. B [**62**]{}, 11556 (2000); J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. [**77**]{}, 3865 (1996). C. -X. Liu, X.-L. Qi, H. J. Zhang, X. Dai, Z. Fang, S. -C. Zhang, Phys. Rev. B [**82**]{}, 045122 (2010). Y. -Y. Li, G. Wang, X. -G. Zhu, M. -H. Liu, C. Ye, X. Chen, Y. -Y. Wang, K. He, L. -L. Wang, X. -C. Ma, H. -J. Zhang, X. Dai, Z. Fang, X. -C. Xie, Y. Liu, X. -L. Qi, J. -F. Jia, S. -C. Zhang and Q. -K. Xue, Adv. Mater. [**22**]{}, 4002 (2010). A similar algorithm is adopted to project the surface bands in Bi$_2$Te$_3$ in K. Park, J. J. Heremans, V. W. Scarola, and D. Minic, Phys. Rev. Lett. (Accepted Tuesday Sep 21, 2010). C. -X. Liu, X.-L. Qi, H. J. Zhang, X. Dai, Z. Fang, and S. -C. Zhang, Phys. Rev. B [**82**]{}, 045122 (2010). http://cst-www.nrl.navy.mil/lattice/struk/b$_k$.html. K. Watanabe, T. Taniguchi, and H. Kanda, Nat. Mater. [**3**]{}, 404 (2004) G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J van den Brink, Phys. Rev. B [**76**]{}, 073103 (2007). For studying dielectric effects on TI surface states, we used numerical atomic orbital based DFT method (http://www.openmx-square.org). This is partly to accelerate the pace of interface calculations. We first benchmarked our numerical atomic orbital results for bulk SiO$_2$, bulk BN and bulk TI against the plane-wave based DFT calculations by choosing a suitable set of pseudo-atomic orbitals or basis sets, a good pseudopotential and the energy cut-off. These computational parameters were then used for the interface calculations. http://cst-www.nrl.navy.mil/lattice/struk/sio2a.html. R. B. Laughlin, Phys. Rev. B [**22**]{}, 3021 (1980)
