Gapped Dirac cones and spin texture in thin film topological insulator

The protected surface states of topological insulators (TIs) form gapless Dirac cones corresponding non-degenerate eigenstates with helical spin polarization. The presence of a warping term deforms the isotropic cone of the most simple model into snowflake Fermi surfaces as in Bi2Se3 and Bi2Te3. Their features have been identified in STM quasiparticle interference (QPI) experiments on isolated surfaces. Here we investigate the QPI spectrum for the TI thin-film geometry with finite tunneling between the surface states. This leads to a dramatic change of spectrum due to gapping and a change in pseudo spin texture that should leave distinct signatures in the QPI pattern. We consider both normal and magnetic exchange scattering from the surface impurities and obtain the scattering t-matrix in Born approximation as well as the general closed solution. We show the expected systematic variation of QPI snowflake intensity features by varying film thickness and study, in particular, the influence on back scattering processes. We predict the variation of QPI spectrum for Bi2Se3 thin films using the observed gap dependence from ARPES results.


I. INTRODUCTION
The most dramatic manifestation in a topological insulator (TI) is the presence of protected surface states which are helical spin locked nondegenerate eigenstates with a gapless Dirac dispersion.This has been manifestly demonstrated e.g. by countless (spin-polarized) ARPES experiments [1][2][3] and STM-QPI [4][5][6][7][8][9][10] investigations on isolated surfaces, in particular in the canonical examples Bi 2 Se 3 and Bi 2 Te 3 .The theory of QPI spectra on single surfaces has been presented in numerous effective model investigations [10][11][12][13][14] and recently also with an ab-initio approach [15].Due to their helical nature these surface states have forbidden backsckattering by normal (scalar) impurities leading to the well-know weak anti-localization (WAL) effects in the surface magnetoconductance experiments [16,17].However the strict WAL due to the destructive interference caused by a Berry phase ∓π applies only to isolated surfaces.When we consider thin films there will be a tunneling between opposite surface states due to their finite decay lenghts into the bulk.The tunneling mixes the surfaces states leading to a finite gap at the Dirac point and it lifts their protection as witnessed by the suppression of helical polarization and modification of the Berry phase for states close to the gap energy [18].This has been directly observed in magnetotransport experiments in Bi 2 Se 3 where the weak anti-localization as indicator of the ungapped states rather suddenly breaks down below a thickness of five quintuple layers (QL) [16] (In Bi 2 Se 3 1QL =c/3 =9.55 Å; c= lattice constant).The gapping of the Dirac cones in thin films has also been observed directly with the ARPES technique [19,20].
Here we want to investigate theoretically what kind of influence the inter-surface tunneling in TI thin-films may have on the QPI pattern in scanning STM experiments.In the known QPI experiments on single surfaces the Fourier transformed scanning image is dominated by the scattering vectors that interconnect the characteristic points of the real 'snowflake' type Fermi surfaces, with those of the backstattering characteristically missing.Therefore our investigation demands that we take into account not only the inter-surface tunneling but also includes the higher order intra-surface warping terms that deform the isotropic Dirac cones with Fermi surface (FS) circles into the observed snowflake geometry.The QPI pattern of the isolated surfaces should be dramatically modified for small bias voltages such that ω = eV ≈ ±|t| where V is the STM bias voltage and t is the inter-surface tunneling leading to a gap energy 2|t|.Furthermore the QPI spectrum due to scattering from magnetic impurities should also be sensitive to the change of spin texture close to the gap threshold.It is induced by the tunneling consisting primarily in a suppression of helical polarization for small energies.In this work we intend to give a detailed theoretical analysis of the tunneling effects in TI thin films that one may expect to see as a guidance to QPI experiments with systematic variation of film thickness.We present the corresponding variation of the QPI spectrum using the thickness dependent surface state gap obtained from ARPES.As a theoretical tool we use t-matrix QPI theory [21,22], including both effect of warping and tunneling.We derive the QPI spectrum in simple Born approximation or closed fully analytical t-matrix result, using both normal and magnetic example of (momentum-isotropic) impurity scattering.
The paper is organized as follows: In Secs.II, III we define the model of warped and inter-surface coupled thin film TI states and discuss their spin texture in Sec.IV.In the main formal part of Sec.V we derive the full t-matrix for the assumed impurity scattering potential and discuss the simplified Born approximation.These results are then applied to calculate the QPI spectrum in Sec.VI in the normal and magnetic scattering cases with full t-matrix approach.An extended discussion of the numerical results as function of bias voltage and film thickness is presented in Sec.VII.Finally Sec.VIII gives the conclusions.

II. MODEL FOR WARPED TOPOLOGICAL INSULATOR DIRAC CONES
The protected surface states in topological insulators like e.g.Bi 2 Se 3 and Bi 2 Te 3 may be derived from a k • p-type Hamiltonian model for the 3D bulk states [23,24] using four strongly spin-orbit coupled basis states.The surface states are obtained by solving the 1D Dirac equation along the surface normal direction (ẑ) under appropriate boundary conditions [23,24].The effective low energy Hamiltonian for the 2D topological surface states can then be parametrized in the well-known form [11,[25][26][27] Here the first term describes the massless Dirac states forming an isotropic cone dispersion in the 2D surface Brillouin zone (SBZ) where the velocity v is the slope of the cone and its position is the SBZ projection of a bulk time reversal invariant (TRI) point of the bulk Brillouin zone, e.g., the Γ point for the above compounds.The second 'warping term' distorts the cone anisotropically in accordance with crystal symmetry in such a way that the 2D Fermi surfaces (cuts through the cone) are now six-pronged 'snowflakes'.The amount of distortion is determined by the strength λ of the warping term.The 'pseudo-spin' σ corresponds to the Kramers degeneracy of states at the projected TRI points which is lifted due to spin orbit coupling when moving away from them.Furthermore the surface wave vector k and the azimuthal angle θ k = tan −1 (k y /k x ).The above 2×2 Hamiltionian can be explicitly written as where k = E * k = vk is the isotropic (λ = 0) Dirac dispersion.An overall parabolic term 0k = k 2 /2m * that breaks particle hole-symmetry for larger energies [26] will be neglected.We also introduced a momentum scale k c = v/λ and energy scale E * = vk c to define dimensionless momentum variables k = k/k c and energies ˆ k = k /E * = k, see also Table I.The single surface Hamiltonian may be diagonalized by a unitary transformation according to with the warped Dirac cone energy now given by 1 2 , (4) or dimensionless Êk = E k /E * .The dispersion of Eq. ( 4) shown in Figs.1(i and j) results in a warped constant energy that surfaces have the shape of six-pronged 'snowflakes' (Fig. 1).The unitary transformation to the eigenstates is given (n = 0−5) and dents in between.Left column corresponds to isolated surfaces (t = 0) and center column to coupled surfaces (t = 0.8).In the latter the surfaces are absent for small ω and reappear for ω t, becoming similar to uncoupled case for ω |t|.(i,j) Gap opening in warped Dirac cone dispersion due to tunnelling t for momentum along dent (i) and tip (j) directions, where each curves corresponds to a different inter-surface tunneling values, as 0 ≤ t ≤ 1 in steps of 0.2, from blue to brown lines.Here, and in the rest of paper we set the energy unit as E * .by the matrix The columns of S k are the eigenvectors |ψ κ (in pseudospin ↑↓ basis) to the eigenvalues ±E k (κ = ±).These are called 'chiral' or 'helical' basis states because in the isotropic case λ = 0 they are also eigenstates to the chirality operator defined by κ z = (σ × k) • ẑ = ±1.The pseudo-spin mixing angle φ k is given by In the full circle 0 ≤ θ k < 2π, cos 3θ k changes sign six times at θ n = (2n + 1) π 6 ; (n = 0 − 5).To obtain a continuous variation of φ k (θ k ) at these boundaries and to have a well defined uniform convergence to the isotropic limit φ k → π 2 for vanishing warping (λ → 0) we define φ k = tan −1 (1/ k2 cos 3θ k ) for cos 3θ k > 0 and φ k = tan −1 (1/ k2 cos 3θ k ) + π for cos 3θ k < 0. This amounts to taking the second (upper) branch of tan −1 in the latter case.The k and θ k dependence of the pseudospin mixing angle φ k is shown in Fig. 2, and it indeed varies continuously with θ k around the (half-)circle centered at the isotropic limit φ k = π/2.
For calculation of the QPI spectrum one needs the Green's function in pseudo spin representation.For the single surface problem it is given by (7) Defining the warping energy by , the Green's function is obtained as In the isotropic Dirac cone case one simply has to set ∆ k = 0 and E k = k .

III. SPLITTING OF DIRAC CONES BY INTER-SURFACE TUNNELING IN TI THIN FILMS
The warped cones of single TI surfaces have been abundantly demonstrated by ARPES experiments [28,29].As outlined in the introduction in thin films the tunneling of helical Dirac states between the two surfaces introduces a gap in the spectrum signifyng the transition to topologically trivial surface states close to the gap threshold as function of thickness.This has been observed again directly with ARPES [19,20] in agreement with magnetotransport results [16].Theoretically it has been studied in detail [23,24,30] by starting from the bulk thin film states and introducing the approriate boundary conditions.A nonmonotonic, even oscillating behaviour of surface state tunneling and gap size as function of thickness d is possible although this has sofar not been observed in photoemission (Sec.VII).
For the purpose of investigating QPI signatures of the topological transition in thin films we employ a phenomenlogical model starting from the independent helical states on the two isolated surfaces.They are coupled by the thickness dependent tunneling matrix element t(d) of top (T) and bottom (B) surface states with equal helicities (Fig. 3).The 4×4 thin film Hamiltonian in the space with basis |κα (helicity κ = ±1 and surface index α = T,B) may be written as 2 × 2 block matrix in T,B space according to in the pseudo-spin space with basis |σα ( σ =↑, ↓) (left) and helicity space with basis |κα (κ = ±1) (right).In both representations α = T,B is the surface index.The minus sign in the last element appears because states of equal helicities on B,T k and k = k + q momenta leads to surface density oscillations δNT (q, ω) described by Eq. ( 32) and scanned by the tip at bias voltage eV = ω.
surfaces belong to opposite half-cones (see Fig. 3).We use the convention that σ, κ, α denote the vector of 2 × 2 Pauli matrices in pseudo-spin, helicity and surface layer (T,B) space, respectively.Furthermore the unit in each space is denoted by σ 0 , κ 0 , α 0 .The film Hamiltonian may be diagonalized by two transformations for the κ = ±1 helicities separately with the simple form where t = t/E * .The total transformation matrix from pseudo-spin single surface representation to helicity-film eigenstate representation that diagonalizes H k is then given by W k .The column vectors of this 4 × 4 matrix are the helicity eigenstates of the film (explicitly given in Appendix A).They have the energies resulting from the diagonalization Therefore at the Dirac point k = 0 a gap of size ∆ t = 2|t| opens between T,B hybridized film states (τ = 1, 2) each twofold degenerate due to time reversal symmetry with helicities κ = ±1.

IV. MIXED TOPOLOGICAL SURFACE STATES AND THEIR SPIN TEXTURE AND BERRY PHASE
The explicit form of the surface states is given by the columns of Eq. (A1).For the upper and lower split half cones with energies E (1,2) (k) = ± Ẽk the pairs of doubly degenerate (κ = ±) upper and lower (τ = 1, 2) band states | ψτκ can be written, respectively as  12) decouple to isolated surface states on T,B such that | ψk 1± become states with opposite chirality on opposite surfaces with positive energies and likewise | ψk 2± with negative energies.From these states we may now caluclate the change of the surface spin texture on T or B (they are time-reveresed copies) of the model as function of the tunneling matrix element t between the isolated helicity surface states.The spin textures of states 1,2 will be opposite on the same surface and identical on opposite surfaces.Therefore, restricting to the top surface we have for the spin expectation value of each pair corresponding to each half cone: ) Using the previous expressions for the mixing angles φ k , ψ k in terms of the polar angle θ k we get explicitly for in-and out-of-plane spin polarization: Here and therefore the total length of the pseudo-spin on T is When the tunneling between surfaces becomes negligible ψ k → 0 and σ tot 2 T → 1. Contour plots of σ T and σ z T are shown in Fig. 4 for isolated and coupled surfaces.In particular the in-plane component which is maximum (≤ 1) at the Γ point and along the tips (white region) in Fig. 4(c) is suppressed to zero by the tunneling at the Γ point in (d) (dark green area).We note that ab-initio calculations of the spin-textures for thick slabs have also been perfomed [31].The Dirac point of topological surface states may be viewed as a monopole in momentum space.It is connected with a nonvanishing Berry-phase when encircling the Dirac point on a closed path C containing the origin.First we consider the case of warped Dirac cone of a single (T or B) surface.Then the Berry phase is given by (κ = ±) [32,33] where the Berry connection described by depends on the warping through φ k .Because of the antiperiodic property cos[φ(θ k + π 3 )] = − cos[φ(θ k )], the contour integral over the cosine vanishes and we get γ ± = ∓π.This means that the topologically non-trivial Berry phase for the Dirac cone states is not influenced by the warping effect since the latter does not destroy the k = 0 singularity.Therefore when considering the effect of the inter-surface tunneling t on the Berry phase we may safely neglect the warping for simplicity.The tunneling leads to film states determined by the T,B inter-surface mixing angle ψ k (Eq.12).The Berry connection for the isotropic (φ k = π 2 ) film states (τ = 1, 2, κ = ±), | ψk τ κ with a given bare Dirac cone energy E = E k = vk may be calculated as (17) Therefore the Berry phase of states at the gap threshold which have parabolic dispersion vanishes while for energies much larger than the gap it approaches the values γ τ κ = (−1) τ κπ, (τ = 1, 2, κ = ±1) of the isolated surface Dirac cones.Here κ = ± correspond now to degenerate pairs (they become the T +, B− states for τ = 1 and B+, T − states for τ = 2 when t, ψ k → 0).The reduction of the Berry phase close to the mass gap of Dirac electrons leads to a violation of topological protection.Therefore WAL due to destructive interference caused by the Berry phase ∓π is suppressed and breaks down completely for large enough mass gap in ultrathin films with d ≤ 5QL [16].Furthermore the acquired degeneracy of film states opens the backscattering channel for the gapped states reducing the surface mobility [16] and also influencing the QPI signatures.

V. SELFCONSISTENT T-MATRIX THEORY FOR IMPURITY SCATTERING
The STM-QPI method measures the electronic density fluctuations on the surface caused the interference of scattered and ingoing waves at an impurity site.It must be stressed that this is a single impurity effect [21] although the amplitude of the density fluctuations will be proportional to the number of impurities.In the situation of a thin film another aspect is important.Due to the tunneling the film states are eigenstates composed of surface states on both top and bottom surface.Therefore, even if the density fluctuations are measured on the top surface they will also be influenced by the scattering on the bottom surface due to the tunnelling.This effect which is illustrated in Fig. 4 has to be included in the calculation.For the impurity scattering potential in pseudospin (σ) and FIG. 6. Momentum integrated Green's function g0(ω) (Eq.23) which determines the frequency dependence of the full t-matrix elements in Eq. ( 24).It becomes develops singular behaviour at the gap edge ω = t.
T,B surface space (α) we assume the generic form Here V c denotes charge and V m exchange scattering by normal and magnetic impurities, respectively, which does not depend on the momentum transfer q = k − k of the scattering process.In this case a closed solution for the scattering tmatrix necessary to compute the QPI spectrum is possible in pseudo-spin representation which then will have to be transformed to the helical eigenstate basis of the surfaces.The second exchange scattering term in Eq. ( 18) corresponds to a frozen impurity spin (along z-direction), neglecting any dynamical spin flip processes.This spin polarization can be achieved either by tiny magnetic field or, since the impurity concentration is finite, may be the result of long-range exchange interactions among the magnetic impurities that lead to quasi-static behaviour (long relaxation times) for each impurity spin.

A. Closed solution in pseudo-spin basis
For a momentum-independent scattering potential the tmatrix equation in pseudo-spin basis may be solved as We made the approximation that impurities on T,B do not scatter between the surfaces but only within each of them, i.e. we approximate R −1 = R −1 T ⊕ R −1 B and therefore t T B = t BT = 0. However the propagation, i.e. tunneling between surfaces is fully taken into account by the non-diagonal Green's function elements given below.Furthermore we assume that the type and scattering potential of impurities on T, B are the same and therefore set VT = VB = V .Then The spin flip amplitudes t ↑↓ , t ↑↓ that appear formally vanish identically (Appendix C) so that we only have to consider the diagonal elements t ↑↑ , t ↓↓ in the following.Ultimately this is due to the absence of spin-flip terms in the magnetic scattering term of Eq. ( 18).Furthermore the film Green's function Ĝk (iω n ) in pseudospin basis is given by With h k defined in Eqs.(1, and 2).The individual t-matrix elements t σσ may be evaluated (see Appendix C) from Eq. ( 19).
It is also useful to introduce (anti-)symmetric combinations by t s,a = 1 2 (t ↑↑ ± t ↓↓ ).Furthermore we need the determinants D = det[R T,B ] which are given by The momentum-integrated Green's function in the above expressions which depend only on frequency iω n are defined by (see Eq. ( 11)): the (anti-) symmetrized t-matrix elements are then obtained as Obviously the anti-symmetric amplitudes require a magnetic scattering and for normal scattering (V m = 0, t a = 0), we only obtain the symmetric scattering amplitude For this calculation we used the pseudo-spin basis because then the scattering potential is isotropic allowing for a closed solution.But for later calculation of the QPI spectrum we now have to transform back to the helical eigenbasis of the surfaces.The dynamics of the scattering matrix is determined by that of the momentum-integrated g0 (iω n ) which is shown in Fig. 6.

B. transformation to helical basis
The transformation to helical basis is accomplished by the unitary matrix S k in Eq. ( 5) which is identical for T,B.The is then explitcitly given in terms of the pseudo-spin based solution Eq. ( 24) by its elements for both T,B blocks as The (T,B-independent) form factors α ± ij (i, j = ±) that are functions of (φ k , θ k , φ k , θ k ) originating from the helical eigenstates are given explicitly in Appendix B. Before proceeding to the QPI spectrum it is useful to discuss the general solution for the t-matrix in the simple case of the Born approximation (BA) corresponding to single scattering from the impurity, i.e. to first order in V c and V m .

C. Born approximation and selection rules for scattering
In BA we have t s = V c and t a = V m .Treating normal (V c ) and magnetic (V m ) scattering separately, we get the simple complementary results, respectively [c.f.Eq. (B2)]: First we consider the limiting case of isotropic Dirac cones, i.e. vanishing warping term λ → 0. Then 2 and similar for k .Explicitly we have  40)) for normal scattering Vc ≡ 1 only (Vm = 0).The images are for for various bias voltages eV = ω (rows) and inter-surface tunneling strength t (columns).When the latter is increased the typical QPI structure appear only for ω > t due to the vanishing of constant energy surfaces for smaller frequencies (Fig. 1).The details of the QPI images are explained in Sec.VII.
while for purely magnetic scattering the complementary result is This demonstrates that backward scattering between same helicity states is forbidden for normal impurity scattering while it is allowed in the case of magnetic impurities.This will have direct influence on the QPI spectrum (Sec.VI).
It is also useful to derive the explicit expressions of form factors for forward (f ; θ k = θ k ) and backward (b; θ k = θ k +π) scattering including the effect of finite warping (λ = 0) when φ k = 0 in general.Using φ k = φ k in case (f ) and the identities sin φ k 2 = cos φ k 2 ; cos φ k 2 = sin φ k 2 in case (b) we nevertheless get for normal scattering V c (V m = 0) the identical result as in Eq. ( 29).Thus the pseudo-spin mixing angle φ k of the warped case does not enter in the f, b amplitudes and the result is identical to the isotropic Dirac case.In particular this means that the backscattering for equal helicity states remains forbidden even for the eigenstates in the warped Dirac cone.This is natural since the normal scattering does not react to the changed spin texture caused by the warping term.
On the other hand for exchange scattering V m (V c = 0) the change in spin texture will be important and therefore the result for the scattering matrix will depend on the pseudo spin mixing angle φ k caused by warping.In BA we obtain from the only nonvanishing second term in Eq. (B3) In the isotropic cone limit (λ → 0, φ k = π 2 ) we recover the result of Eq. (30).We note that in the warped case when φ k = π 2 in general the forward scattering for equal helicity states no longer vanishes as in Eq. ( 30) due to the presence of the perpendicular spin component ∼ cos φ k of the eigenstates [Eq.( 13)].

VI. QUASIPARTICLE INTERFERENCE SPECTRUM
The spectral Fourier component of the surface density modulation visible by STM-QPI is given by [21,34,35] Here ω = eV with V denoting the variable STM-tip bias voltage (Fig. 3).It is assumed that the tip is placed to the top (T) surface, therefore only the spin-trace over the T-block of t-matrix and Green's function product has to be performed.Explicitly we have two contributions given by As it stands everything is still written on pseudo-spin basis.Since the Green's functions are diagonal in the helical bases one should transform to the latter, using Eq. ( 25) and the similar transformation for the Green's function: where α, α = T, B. For the t-matrix we restricted to T,B diagonal elements only, neglecting the impurity scattering between the bottom and top surfaces.This cannot be done for the Green's functions because the inter-surface tunnelling or hybridization is essential for the low energy surface states of the TI film.Therefore the T,B nondiagonal Green's function elements have to be kept.This leads to two terms in the kernel Eq. ( 33) of the QPI spectrum involving scattering at the top (first term), and via hybridization of surface states also at the bottom (second term) surface.They are schematically shown in Fig. 3.
The T,B Green's function blocks are all diagonal in the helical eigenstate basis.Therefore it is economic to transform Eq. ( 32) to this basis using Eqs.(25, and 34): The latter leads to Where TT and BB correspond to upper and lower sign, respectively.Furthermore Using the diagonal Green's functions the kernel in Eq. ( 33) may be transformed to helical basis as T T .

A. Born approximation for QPI
To obtain a better insight in the general expression for the QPI spectrum we first simplify to the case of Born approximation with only single scattering events at the impurity included.In this case we have t T,B s = V c and t T,B a = V m .Again we treat normal (c) and magnetic (m) scattering separately.Inserting into Eq.( 38) and using the explicit form factors and expressions for Ã±k we obtain for normal scattering, after considerable algebra: , (39) with and likewise for magnetic scattering , (40) with Note that in these expression the gapped cone energies Ẽk of the film appear in the denominator while the ungapped energies E k of isolated surfaces remain in the numerator.It is instructive to consider the limit of the isotropic Dirac cone when λ → 0 and then φ k = π 2 .Inserting this into the form factors m ± c,m (k, k ) the above formulas then simplify to .
(41) For the decoupled case (t = 0, Ẽk = E k ) these expressions can also be directly obtained from Eq. ( 32) by using the single-surface Green's functions of Eq. ( 8) for the isotropic case.Due to the sign change of the numerator of Λ m under k ↔ k which results from the helical spin-locking we have Λ m ≡ 0 in Born approximation.For a finite result one has to use the full t-matrix theory.

B. Full t-matrix expressions QPI spectrum
For computational convenience we also give here the rather lengthy explicit expression for the QPI spectrum including the general t-matrix elements that will mostly be used in the numerical calculations.The total QPI spectrum, equal to the sum T (q, iω n ) from Eq. ( 38) is then obtained as (42) In each term the first part describes the QPI contribution on the top surface due to scattering on top surface while the second part proportional to t 2 represents the QPI contribution on the top surface due to scattering on the bottom surface.Obviously this term can only be present when there is tunneling (described by G T B , G BT in Fig. 3) between the T,B surface states, therefore it vanishes for decoupled surfaces.In these equations the expressions for the full t-matrix elements tκκ (iω n ) in helical presentation are obtained from Eq. (26).They are composed of the irreducible t s,a matrix elements in Eq. ( 24) and the form factors α ± κκ of Appendix B.

VII. NUMERICAL RESULTS AND DISCUSSION
We now discuss the numerical results for the QPI spectrum under systematic variation of bias voltage ω = eV , intersurface tunnelling t of the thin film and normal (V c ) and magnetic (V m ) impurity scattering potential (assumed identical on both surfaces).We mostly use the full t-matrix approximation as given in closed form by Eq. (42) in this section except when stated otherwise.In Fig. 7, we show an overview over the images for normal V c -type scattering where the intersurface tunneling t varies along the columns and the bias voltage eV = ω (measured from the Dirac point) along the rows.Here for once the Born approximation [Eq.( 40)] is employed for comparison.Generally in QPI images the Fermi surface for a given ω is reproduced , with a doubling of the FS radius.However, the intensities at wave vectors connecting special Fermi surface points like tips and dents may be strongly enhanced or depressed.The overall extension of the QPI images increases with ω, the distance from the Dirac point, according to the diameter of the cut ω = E k (or ω = Ẽk ) through the warped cone which gives the snowflake FS.The shape of the pattern and its increasing radius with ω is clearly seen in the first row representing the isolated surfaces (t = 0).When the tunneling is turned on (2 nd and 3 rd row) the low energy spectrum is gapped and therefore around ω t (i.f.we use ω, t ≥ 0) the QPI image will be strongly modified: The radius shrinks and the anisotropic 'snowflakes' character is reduced leading to a more isotropic image.This is completely in accordance with the evolution of the Fermi surfaces in Fig. 1.Therefore QPI investigation of thin films can give full information on the thickness (tunneling) dependence of low energy quasiparticles close to the gap threshold.For larger energies (voltages) ω t above the gap the QPI image approaches that of the isolated surfaces (c.f. the two top right figures).
A detailed comparison of Fermi surface shape (first row a,d) and according QPI image is presented in Fig. 8(a-f) for normal (second row b,e) as well as magnetic (third row c,f) scattering mechanism.Similar as in Refs.[5,11] we can identify characteristic wave vectors q i (i=1-7) in Fig. 8(a,d) connecting special points where the azimuthal group velocity vanishes (tips and dents of the snowflake) and should therefore figure prominently in the QPI image.Indeed most of the q i vectors (with the exception of q 3 ) can be identified in Figs.8(b,c) for the isolated surface corresponding either to large (dark red) intensity or low and vanishing (deep blue) intensity areas with arrows representing the characteristic wave vectors pointing to them.Some of the latter in Fig. 8 correspond to the forbidden backscattering as they are connected by q 1 , q 4 scattering vectors that pass through the origin such that k = −k.These regions are, however, quite narrow because of the vicinity of (equivalent) closeby allowed q 6 , q 7 vectors that are associated with large intensities.The q 2 , q 5 scattering vectors can also be identified though less prominently.The overall pattern for magnetic scattering in (c) looks quite similar although the pointwise intensities are largely different.In particular for magnetic impurities the backscattering vectors q 1 , q 4 are allowed (Secs.V C,VI A) and are also associated with finite QPI spectral intensity.There is another weak, though noteworthy peculiar feature of the QPI image: The single snowflake of the FS experiences a doubling in the QPI image.The outer one has the same orientation as the FS and is associated mostly with scattering vectors close to the group q 1 , q 4 , q 6 , q 7 and their equivalents while the inner smaller image is rotated by π/3 with respect to the FS and is mostly associated with scat- tering vectors close to q 2 , q 5 and equivalents.
When the inter-surface tunneling t is turned on and the gap is opened at the Dirac point, at first for ω |t| there is no drastic change in QPI image characteristics as seen from the upper right corner of Fig. 7, except that the characteristic scattering vector lenghts |q i | are shrinking.This becomes more dramatic when ω approaches the (half-) gap size |t| and the FS structure becomes more circular (i).Furthermore the tunneling opens backscattering channels close to the gap edge increasing the relative size of the high intensity areas (dark red).As a result the QPI snowflake image gradually 'melts' into a circular drop shape where only the most prominent scattering vectors q 6 , and q 7 can still be discerned.Their length is now close to 2k F of the nearly spherical Fermi surface image [Figs.7(d,i) and 8(e,f)] corresponding to the parabolic surface dispersion close to the gap edge (c.f.Fig. 1).This gradual transformation of the QPI image from the six-pronged (double) 'snowflakes' at zero or small tunneling (isolated surfaces) to the almost circular 'raindrop' structure at large tunneling (few quintuple layer films) should be worthwhile to investigate experimentally with STM-QPI technique.
Instead of looking at the QPI image of |Λ(q x , q y , ω)| in the 2D SBZ as before we may consider the QPI 'quasiparticle dispersion' defined by the complementary image |Λ(q = q, ω)| in the (q, ω)-plane for fixed momentum direction q Γ K or Γ M in the SBZ which should directly demonstrate the change of surface quasiparticle spectrum with tunneling strength.The dispersion results for normal scattering are shown in Fig. 9.For isolated surfaces (t = 0) the Dirac cone dispersion can clearly be identified by the envelope high intensity region.For T,B surfaces connected by the tunneling the gap |t| opens progressively and the parabolic dispersion for ω ≥ |t| is again seen in by a sharp prominent high intensity envelope.The destruction of the Dirac cone by the gap opening via inter-surface tunneling has sofar been observed directly only in ARPES experiments for Bi 2 Se 3 thin films between 1 − 6 QL thickness [19] and it would be highly desirable to investigate this also by the complementary QPI method as proposed here.The existence of the Dirac cone dispersion for the isolated surfaces (t = 0 in Fig. 9) has indeed been demonstrated before with the QPI method [6][7][8].
In an alternative constant-ω presentation of QPI images |Λ(q x , q y , ω = eV )|, we again keep the bias voltage fixed and change the tunneling t and scattering strengths V c , and V m .This is presented in Fig. 10, where we show the images for a constant ω = 1.25 and as function of t (columns) and relative normal (V c ) to magnetic (V m ) scattering strength (rows).While the overall image features are similar, the substructure of the large intensity regions (dents) depends in a subtle way on the relative size of V c , and V m .Finally in Fig. 11, we show an example of the evolution of QPI spectrum for normal scattering from very low (BA) to very strong scattering strength.The global structure of the image is preserved but details show again subtle changes: The high intensity spots for small scattering merge into a ring for large V c , but still keeping low intensity at the backscattering vectors.Furthermore the interior of the structure acquires more intensity presumably due to the increasing importance of multiple scattering processes with growing V c contained in the full t-matrix approach.
Sofar the tunneling matrix element t has been treated as an arbitrary but a fixed parameter.It is, however, an effective (film-) surface states parameter that derives from the solution of the true boundary value problem [23,24] starting from the bulk states of the bulk Hamiltoninan given in k • pparametrized form [30,36]. Using these parameters it was shown in Ref. 24 that the dependence of t(d) on film thickness d should be given by t(d) = t 0 exp(−d/d 0 ) sin(d/d 0 ) where the reference scales t 0 , d 0 , d 0 are given in terms of the k • p -parameters of the bulk [30,36].The t(d) dependence contains an overall exponential decay with increasing d but also an oscillatory term.The latter may lead to a closing of the gap for intermediate d-values [24].With respect to expected QPI pattern dependence on thickness this should have dramatic visible consequences.Let us assume for each thickness the bias voltage is kept at a constant value with respect to the Dirac point (which may itself depend on the thickness) in each film.Then, as the film thickness is varied the gap strongly varies and the according QPI pattern should change between relatively isotropic circular pattern for large |t| and strongly six-pronged snowflake patterns for small |t|.Observation of such sequence in STM-QPI experiments should be taken as direct proof for the gap oscillation.
However, sofar ARPES experiments in Bi 2 Se 3 have not shown any indications of a gap oscillation or closing [19] but rather show a monotonous decrease of |t(d)| in the range 2QL < d < 8QL, where it vanishes on the upper value, restoring the isolated Dirac surface states.This is shown by the red line in Fig. 12.On the top row we show the expected QPI image at ω = 0.5 for the five measured (half-) gap sizes on the red line in corresponding sequence.For increasing d one observes the monotonic appearance of the snow-flake features out of the isotropic QPI image for maximum gap at 2QL caused by the concommitant decrease in |t|.If this behaviour could be observed in STM-QPI it would support the suggestion of ARPES that the gap decays indeed monotonically with increasing d and the predicted oscillation, for some reason, is not present in Bi 2 Se 3 .

VIII. SUMMARY AND CONCLUSION
In this work we proposed a theory of quasiparticle interference for topological insulator thin films.We start from a model of single Dirac cone surface states with a warping term included to reproduce the realistic snowflake FS shape of the isolated TI surface in compounds like Bi 2 Se 3 and Bi 2 Te 3 .The main physical effect of thin film geometry is implemented within a model where equal helicity states can tunnel from one surface to the other, described by the tunneling energy t.This has two consequences: A gap opening is observed due to the (momentum-dependent) mixing of top/bottom states and the helical spin texture and Berry phase of Dirac states is strongly modified close to the gap edge at ω = |t|, leading to a vanishing of the in-plane spin component and Berry phase.
The warped surface cones have been previousy found for isolated surfaces both by ARPES and QPI methods.The gap opening of surface states in thin films as function of thickness has sofar only been directly seen with the former method.
Here we investigated what is to be expected for QPI experiments.To calculate the QPI signal we use the t-matrix approach with a full summation to infinite order in the scattering strength as well as the Born approximation.For impurity scattering we employ a generic potential consisting of isotropic scalar and exchange scattering whose strengths are model parameters.We find a closed solution of the t-matrix including the effect of inter-plane propagation.The impurity scattering and the tunneling to the STM tip is expressed in the pseudospin basis whereas the eigenstates of the warped cones have a complicated helical spin texture.Therefore the analytical espressions involve a summation over terms with different azimuthal form factors which also control the behaviour under backward scattering.
The QPI images can be interpreted in terms of characteristic scattering vectors connecting tips/dents in the snowflake FS for a given bias voltage.The allowed scattering and forbidden backscattering can be identified for normal impurities as long as the inter-plane tunneling is moderate.For constant ω the latter, when the gap size increases, shrinks the dimension of the QPI image and its anisotropic features.This turns the sixpronged snowflakes of QPI with strong intensity anisotropy into a more isotropic circular 2k F image of the Fermi surface of the gapped Dirac cones, also due to reopening of backscattering.Once ω < |t| falls below the (half-) gap size all specific QPI features are suppressed.It would be worthwhile to investigate this systematic variation of predicted QPI signal in experiment.In particular as previous surface state calculations in Bi 2 Se 3 and Bi 2 Te 3 have shown that t(d) behaves nonmonotonically or even oscillates with film thickness d.This would imply that as thickness is reduced continuously the QPI image of TI surface state will oscillate between anisotropic snowflake and nearly isotropic circular patterns.However, ARPES experiment sofar have not confirmed this prediction but rather find a monotonous reduction of the gap with increasing film thickness.This discrepancy would justify to perform complementary STM-QPI experiments to investigate this issue further.Finally the dispersion of QPI thin film images in (q, ω) plane for fixed direction of scattering vector q shows the breaking of the Dirac cone by gapping and appearance of quadratic low energy dispersion.Their observation would give a further instructive comparison to the dispersions observed by ARPES as function of film thickness.

FIG. 1 .
FIG. 1. (a-h) Constant-energy (ω) surfaces of warped Dirac cones with tips at θn = (2n+1) π 6 ;(n = 0−5) and dents in between.Left column corresponds to isolated surfaces (t = 0) and center column to coupled surfaces (t = 0.8).In the latter the surfaces are absent for small ω and reappear for ω t, becoming similar to uncoupled case for ω |t|.(i,j) Gap opening in warped Dirac cone dispersion due to tunnelling t for momentum along dent (i) and tip (j) directions, where each curves corresponds to a different inter-surface tunneling values, as 0 ≤ t ≤ 1 in steps of 0.2, from blue to brown lines.Here, and in the rest of paper we set the energy unit as E * .

FIG. 2 .
FIG. 2. (a) Dependence of pseudo spin mixing angle φ k due to warping term and inter-surface state mixing angle ψ k due to T,B tunneling as function of wave number.Here θ k =0, i.e. the wave vector points to the dents in Fig. 1.(b) Azimuthal θ k angle dependence of warped and T,B hybridized cone energies Ê = Ẽk /E * and mixing angles φ k , ψ k .Full lines for t = t/E * = 0.8, broken lines for isolated T,B surfaces ( t = 0).Here φ k does not depend on t because the latter is diagonal in the helicities.

FIG. 3 .
FIG. 3. Sketch of surface scattering processes in QPI for thin film geometry.Surface normals ẑB = −ẑT are opposite.Isolated T, B Dirac cones with helicities κ = ±1 are indicated.t(d) is thickness ddependent tunneling energy between equal helicity states on T,B.Impurity scattering of (top) cone states with momentum k to k = k + q is possible on both surfaces (VT , VB) due to effect of tunneling; simlar for the bottom states.Interference of waves with k and k = k + q momenta leads to surface density oscillations δNT (q, ω) described by Eq. (32) and scanned by the tip at bias voltage eV = ω.

2 t E k π 4
) in terms of the isolated surface chiral states |ψ ακ .These eigenstates are related by combined time reversal ± ↔ ∓ and interchange of surfaces T, B ↔ B, T .For states close to the gap with Ẽk ≥ t (k k c ) one has ψ k π 4 meaning these are true film states with equal amplitudes on B,T surfaces.On the other hand far from the gap with Ẽk t (k k c ), we have ψ k ≈ 1 and consequently the above states are mostly localized on either T or B surface with only small amplitude on the opposite B or T surfaces.Therefore for large thickness when t → 0 the state in Eq. (

FIG. 4 .
FIG. 4. Contour plot of spin texture.Left (a,c) for t = 0: In the cone center σz T vanishes whereas σ T 1 (white region).Right (b,d) for coupled surfaces t = 0.8 both σz T , σ T vanish around k ≈ 0. The sign of σz T alternates (green/red) when moving around the cone center.See also Fig. 5.
FIG.7.QPI spectrum in Born approximation (Eq.(40)) for normal scattering Vc ≡ 1 only (Vm = 0).The images are for for various bias voltages eV = ω (rows) and inter-surface tunneling strength t (columns).When the latter is increased the typical QPI structure appear only for ω > t due to the vanishing of constant energy surfaces for smaller frequencies (Fig.1).The details of the QPI images are explained in Sec.VII.

FIG. 8 .
FIG. 8. Comparison of characteristic wave vectors qi(i = 1, 7) conntecting Fermi surface tips and dents (a,d) with prominent features in the QPI spectrum for normal charge (Vc) (b,e) and magnetic (Vm) (c,f) impurity scattering.Most qi can be clearly associated in FS and QPI images, only q3 has too weak intensity.Note that intensity for back-scattering vectors q1, q4 passing through the origin is suppressed for normal impurity scattering (b,e) while it is finite for magnetic scattering (c,f).

FIG. 9 .
FIG.9.The dispersion characteristics of QPI spectrum along symmetry directions of the SBZ MΓ (−π, 0) and Γ K (0, π) as function of inter-surface tunneling strength.The gap opening of the Dirac cones can clearly be followed in the QPI dispersion.

FIG. 11 .
FIG. 11.The evolution of QPI pattern from small to large normal scattering potential Vc at the energy ω = 0.75 for the case of t = 0.4.

TABLE I .
Dirac cone parameters for Bi2Se3 (Ref.12 and 26).The intrinsic momentum and energy scales kc and E * are associated with the warping parameter λ.