Effective Hamiltonian for surface states of topological insulator thin films with hexagonal warping

The effective Hamiltonian of the surface states on semi-infinite slabs of the topological insulators (TI) Bi2Te3 and Bi2Se3 require the addition of a cubic momentum hexagonal warping term on top of the usual Dirac fermion Hamiltonian in order to reproduce the experimentally measured constant energy contours at intermediate values of Fermi energy. In this work, we derive the effective Hamiltonian for the surface states of a Bi2Se3 thin film incorporating the corresponding hexagonal warping terms. We then calculate the dispersion relation of the effective Hamiltonian and show that the hexagonal warping leads distorts the equal energy contours from the circular cross sections of the Dirac cones.


I. INTRODUCTION
Experimental measurements 1,2 of the Fermi surfaces of the surface states on semi-infinite slabs of the three-dimensional topological insulator (TI) Bi 2 Te 3 reveal deviations at moderate Fermi energies from the circular cross sections of the Dirac cones predicted by the Dirac fermion Hamiltonian H = v(⃗ p ×ẑ) · ⃗ σ. This discrepancy led Fu to introduce a hexagonal warping term to the Dirac Hamiltonian based on the symmetries of the underlying crystal structure in order to reproduce the experimental results. 3 The effective Hamiltonian now reads where the quintuple layers are on the x y plane. A hexagonal warping in the surface states of Bi 2 Se 3 which can also be described by Eq. (1) was experimentally discovered a year later. 4 We showed in an earlier work 5 that the additional hexagonal warping term in a semi-infinite Bi 2 Te 3 slab leads to the distortion of the Fermi surface and the opening of a bandgap by an in-plane magnetization, both of which would not have occurred in the absence of the hexagonal warping term.
In this work, we derive the effective Hamiltonian for the surface states of a TI thin film 6-8 of small, finite thickness, and infinite width and length which incorporates the hexagonal warping. A key feature that distinguishes a TI thin film from a TI slab of semi-infinite thickness is that there are now two surfaces, which we label as the top and bottom surfaces. Each of these surfaces admits states localized at the respective surfaces. The states localized on the top and bottom surfaces are not independent of each other but are correlated, for example by the boundary condition that their wavefunctions have to simultaneously vanish at both surfaces. The coupling between the localized states on the two surfaces leads to additional effects absent in semi-infinite slabs where only one surface is present.

II. FOUR-BAND HAMILTONIAN
We start from the effective four-band model Hamiltonian of Liu et al. 9 which describes both the bulk and surface states of the BiSe family of topological insulators where indicates that this is a four band Hamiltonian.) The t's can be interpreted as describing the orbital degree of freedom and the σ's the real spins. Fu's hexagonal warping term in Eq. (1) can be derived from Eq. (2) and originate from the R 1 and R 2 terms.
Following Lu et al., 8 we divide Eq. (2) into a 'perpendicular Hamiltonian' H (4B),⊥ containing constant terms and derivatives perpendicular to the two surfaces but excluding the R 1 and R 2 terms, and a 'parallel Hamiltonian' H (4B), ∥ containing the remaining terms. We solve for the eigenstates of the perpendicular Hamiltonian localized at the two film surfaces, and treat the H (4B), ∥ as a perturbation to H (4B),⊥ . The localized eigenstates of H (4B),⊥ are at least two-fold degenerate due to the spin degeneracy. Consistent with standard degenerate perturbation theory, the effective Hamiltonian for the surface states is then obtained obtained by projecting H (4B), ∥ onto the basis of the localized eigenstates of H (4B),⊥ .

III. EFFECTIVE HAMILTONIAN
We consider a thin film with normal parallel to the z crystal growth direction. Accordingly, the corresponding H (4B),⊥ containing constant terms and the z derivatives reads As H (4B),⊥ is diagonal in spin, we can consider the spin up and spin down states separately. We look for states localized at the surface boundaries with the form of exp(λz), so that k z → −iλ. For a given energy E f , diagonalizing Eq. (3) and equating the eigenenergies with E f gives This gives two values of λ 2 , whose square roots we denote as λ ± . We seek linear combinations of these exponentials which disappear simultaneously at the two surfaces at z = ±W/2 where W is the thickness of the film. Two such combinations are Using the ansatz that the spin up eigenstates are in the form of and putting, for example, χ ↑ into (H − E f ) χ ↑ = 0 gives a set of 2 equations which contain hyperbolic trigonometric functions of z, but which should nonetheless give 0 everywhere independent of the value of z. This indicates that the coefficients in front of the various hyperbolic trigonometric functions should vanish. The coefficients of cosh(λ + z) and cosh(λ − z) in the upper and lower components of (H − E f ) χ ↑ = 0 respectively yield two different equations for c − . Imposing the condition that these two equations for c − agree with each other gives This is a transcendental equation in E f due to the λ ± dependence on E f via Eq. (4), which can be solved numerically. Once we find an energy where the values of c − calculated from the 055706-3 Siu, Tan, and Jalil AIP Advances 6, 055706 (2016) equations resulting coefficients of cosh(λ ± λ z) agree, we can use either expression to obtain the value of c − . Applying a similar procedure on (H − E f )ϕ ↑ = 0 with ϕ ↑ as defined in Eq. (5) gives a corresponding equation to Eq. (6) for the ϕ states which can be solved numerically to give E ϕ and its corresponding eigenspinor. The Hamiltonian for the spin down states is identical to the one for spin up states. The spin down eigenstates and eigenenergies are thus given by the same equations as the spin up states. Now projecting H (4B), ∥ into the basis of ( χ ↑/↓ , ϕ ↑/↓ ) and introducing the shorthand In Eq. (7), theτ's can be interpreted as a degree of freedom indicative of the relative parities of the real and imaginary parts of the wavefunction. In particular, an examination of the eigenstates τ polarized along the ±τ i , (i = x, y, z), directions indicates that theτ x polarization is indicative of whether the states are localized near the top or bottom surface of the film.
In what follows, we use the material parameters for Bi 2 Se 3 from Ref. 9 for all the parameters in Eq. (2) except for R 1 and R 2 in our numerical calculations. Although the hexagonal warping in Bi 2 Te 3 is stronger than in Bi 2 Se 3 , the only set of numerical values for the material parameters in Eq. (2) available in published literature which we are aware of for Bi 2 Te 3 is in Ref. 9. This set of parameters does not admit a solution for real E f in Eq. (6), and hence does not support the existence of localized surface states. This contradicts experimental ARPES observations of localized surface states 10,11 in Bi 2 Te 3 thin films, and casts doubts on the validity of modeling Bi 2 Te 3 thin films with Eq. (2) and the set of parameters in Ref. 9. We therefore eschew Bi 2 Te 3 in favor of Bi 2 Se 3 where the available numerical parameters allow us to make predictions closer to actual experimental results. Deriving the effective Hamiltonian for a semi-infinite slab, we find that A 0 and R 1 in Eq. (2) correspond to v and λ respectively in Eq. (1). The value of A 0 obtained from first-principles calculations in Ref. 9 for Bi 2 Se 3 matches the experimentally determined value in Ref. 4 within 10%, but the value of R 1 is smaller than λ by a factor of 2.5. Accordingly we scale the values of R 1 and R 2 from Ref. 9 up by a factor of 2.5 in our calculations.
Using this set of parameters, we find numerically that M x is real for all values of thicknesses at which surface states exist. Dropping 0 terms and performing a rotation of theτ axis,τ → τ, so that the τ z polarization now indicates which surface the states are localized at, the effective Hamiltonian takes the form of where the explicit form of the parameters that appear above can be read off in comparison with Eq. (7). We have introduced -signs in front of some of the terms so that the parameters assume positive values within the thickness ranges considered.
Recalling that the τ z polarization is associated with whether the state is localized nearer the top or bottom surface of the film, the second line of Eq. (8) is the familiar result that the surface states on each of the two surfaces of a TI film are described by the usual Dirac fermion Hamiltonian v( ⃗ k × ⃗ σ) ·n wheren has opposite signs on the two surfaces. The terms proportional to τ x in the first line denote the coupling between the two surfaces of the film. The last line in Eq. (8) represents the hexagonal warping. The first term h 1 τ z (k 3 x − 3k x k 2 y )σ z on the last line has the familiar form of the hexagonal warping on a semi-infinite slab (the term proportional to λ in Eq. (1)) with an additional factor of τ z denoting that this term has opposite signs on the two surfaces. The second term h 2 (3k 2 x k y − k 3 y )τ z I σ is 0 on the semi-infinite slab as the corresponding expression for h 2 on a semi-infinite geometry is proportional to expectation values of terms which evaluate to 0. h 2 is however finite on the thin film geometry, and leads to an additional shift in energy cubic in momentum.

IV. EIGENSTATES OF EFFECTIVE HAMILTONIAN
We now exhibit some properties of the eigenstates of Eq. (8) by computing its dispersion relations numerically. Panels (a) and (b) of Fig. 2 show the dispersion relations for a 5 nm thick film and a 15 nm thick film respectively. These two thicknesses represent the two extremes of a thin film not much thicker than the minimum thickness needed to support the formation of localized surface states, and a thicker film at which the effective Hamiltonian parameters are almost at their infinite thickness asymptotic values. The biggest difference between the two thicknesses is that a bandgap at ⃗ k = 0 is present at the smaller thickness and absent at the larger one. This difference is due to the decrease of E δ from a finite value to 0 with increasing thickness. The increase of E 0 with thickness gives a slight upwards shift in the dispersion relations along the energy axis at different thicknesses. Since the dispersion relations at different thicknesses are otherwise qualitatively similar, we shall concentrate on the 15 nm thick film.
Compare now the cases where the hexagonal warping is present in the solid curves of panel (b) and switched off in the dotted curves. The hexagonal warping leads to an energy splitting of the otherwise two-fold degenerate bands. (The splitting is larger in thinner films.) Comparing between panels (c) and (d) of the figure, it is evident that the splitting occurs with different magnitudes along different directions in k-space. Panel (d) shows that this anisotropic splitting leads to a deviation of the equal energy contours from the circular cross section of the two overlapping Dirac cones to a profile which exhibits a six-fold rotational symmetry. This six-fold symmetry is a hallmark of the hexagonal warping and is reminiscent of that experimentally observed in semi-infinite bulk slabs of Bi 2 Te 3 and Bi 2 Se 3 . We finally consider the effects of applying an in-plane magnetization in the y direction modeled by the addition of a m y I τ σ y term to the Hamiltonian Eq. (8). Fig. 3 shows the dispersion relations and constant energy contours for a 15 nm thick film for m y = 0.005 eV. There is a significant difference between the effects of the magnetization in the thin film, and in a semi-infinite slab. In the latter, the hexagonal warping allows an in-plane magnetization opens up a bandgap. 5 In contrast, panel (a) shows that for a thin film there is no bandgap opening. Instead the coupling between the spin σ and surface τ degrees of freedom here leads causes the in-plane magnetization to generate a k space separation between the two initially overlapping Dirac cones centered at ⃗ k = 0. Comparing between the dispersion relations of panel (a) where the hexagonal warping is present (solid lines) and switched off (dashed lines), we see that the hexagonal warping shifts the band bottom of the E − k x curves at a given finite value of k y away from k x = 0. This shift occurs in opposite directions for the hole and particle bands, and distorts the constant energy contours away from the perfect circles of the Dirac cones. The distortion in the presence of an in-plane magnetization due to hexagonal warping has also been reported in semi-infinite Bi 2 Te 3 slabs. 5

V. CONCLUSION
In this work, we derived an effective Hamiltonian for the surface states of a topological insulator thin film incorporating the hexagonal warping terms and computed the numerical values of parameters within. We then calculated the eigenstates of the effective Hamiltonian at two representative thicknesses, and showed that the hexagonal warping leads to a lifting of the energy degeneracy of the energy bands as well as a six-fold symmetric distortion of the constant energy contours away from the usual circular cross sections of Dirac cones. In particular, the thin film with hexagonal warping differs from a semi-infinite thick TI slab in that whereas an in-plane magnetization opens up a bandgap in the latter, it does not in the former.