Symmetry analysis of strain, electric and magnetic fields in the $\text{Bi}_2\text{Se}_3$-class of topological insulators

Based on group theoretical arguments we derive the most general Hamiltonian for the $\text{Bi}_2\text{Se}_3$-class of materials including terms to third order in the wave vector, first order in electric and magnetic fields, first order in strain and first order in both strain and wave vector. We determine analytically the effects of strain on the electronic structure of $\text{Bi}_2\text{Se}_3$. For the most experimentally relevant surface termination we analytically derive the surface state spectrum, revealing an anisotropic Dirac cone with elliptical constant energy counturs giving rise to different velocities in different in-plane directions. The spin-momentum locking of strained $\text{Bi}_2\text{Se}_3$ is shown to be modified and for some strain configurations we see a non-zero spin component perpendicular to the surface. Hence, strain control can be used to manipulate the spin degree of freedom via the spin-orbit coupling. We show that for a thin film of $\text{Bi}_2\text{Se}_3$ the surface state band gap induced by coupling between the opposite surfaces changes opposite to the bulk band gap under strain. Tuning the surface state band gap by strain, gives new possibilities for the experimental investigation of the thickness dependent gap and optimization of optical properties relevant for, e.g., photodetector and energy harvesting applications. We finally derive analytical expressions for the effective mass tensor of the Bi$_2$Se$_3$ class of materials as a function of strain and electric field.


Introduction
Topological insulators have an inverted band gap which engenders topologically protected surface states (SS). Exhibiting linear dispersion, the electrons at the surface resemble massless helical Dirac fermions, with spin locked to the momentum. The prime examples of three-dimensional topological insulators are among Bi 2 Se 3 -class of materials [1]. This class, also known as the tetradymite group, contains compounds M 2 X 3 where M is either Bi or Sb and X is a combination of Se, S and Te. The crystal structure consists of unit layers of five atomic layers, so-called quintuple layers (QL). For the simplest case of a (111) surface termination, where the surface is parallel to the quintuple layer, the Dirac cone is fully isotropic for small in-plane momentum, perturbed only by a hexagonal warping to third order in momentum [2]. For other surface terminations the Dirac cone of the SS becomes anisotropic with elliptical curves of constant energy, as reported for the (221) surface in [3]. Anisotropic Dirac fermions have interesting transport properties due to the different group velocity in different directions [4], and have attracted attention in other Dirac materials [5,6].
For both fundamental research and applications it is interesting to be able to tune the physical properties of these materials. One way to do this is by the application of strain. It was reported in [7][8][9] that strain strongly affects the band gap at the Γ point of bulk Bi 2 Se 3 , and could even close the gap at large strain values, i.e. induce a topological phase transition. Whether this is possible to achieve in real materials is not known, but a recent study points to the feasibility of strain induced effects in Bi 2 Se 3 where strain up to 3% was induced by lattice mismatch [10].
In this paper, we use the method of invariants [11] to derive the most general Hamiltonian at the Γ point to third order in wave vector, first order in strain, first order in electric and magnetic fields as well as terms to first order in both wave vector and strain. The allowed third order terms in the wave vector include three terms that were neglected in a previous analysis [12]. Since this model is based solely on the symmetry of the crystal, it is valid for all materials in the Bi 2 Se 3 class.
Using this model, we determine analytical expressions for the modified bulk band structure and the effective mass tensor near the Γ-point. Changes in the effective mass tensor with strain and electric field are important for transport and optical properties [13]. Comparing the strain dependence of the band gap to recent density functional theory (DFT) calculations, allows us to determine some of the strain related parameters in the band structure. From this bulk spectrum, we go on to investigate the effects of strain on the SS of a semi-infinite topological insulator. We analytically derive the spectrum of the SS, applying hard-wall boundary conditions at the surface for a semi-infinite topological insulator with a (111) surface termination. In contrast to the unstrained case, the full in-plane rotational symmetry to linear order in wave vector is broken by terms linear in both strain and wave vector, leading to an anisotropic Dirac spectrum. The spin expectation values of the SS are calculated, and it is shown that the spin-momentum locking is affected by strain, revealing a non-zero spin component perpendicular to the surface. The modified Dirac cone and spin structure of the SS is shown in figure 1. Finally, we show that strain affects the localization of the SS wave functions and thereby the band gap arising in thin films due to the coupling between SS on opposite surfaces. We show that the SS band gap changes oppositely to the bulk band gap under strain: increasing the bulk band gap, localizes the SS and decreases the SS band gap, and vice versa. This shows that not only the bulk, but also the SS band gap can be tuned via strain.

Model Hamiltonian
In line with [12], we wish to derive a 4-band model Hamiltonian describing the Bi 2 Se 3 class, which now includes all allowed terms to orders cubic in the wave vector k, linear in the strain tensor ò ij , linear in magnetic field B or electric field E, and linear in both k and ò ij .
In general any 4×4 Hamiltonian can be written in terms of Dirac Γ matrices defined in terms of the Pauli matrices by: where the coefficients ε, d i , d ij must be real to ensure hermiticity. In the present case the coefficients are functions of k, ò ij , B and E.
Here we use the basis ñ -+ |P1 ,  For unstrained Bi 2 Se 3 the linear dispersion of the helical surface states is isotropic in the plane to linear order in the wave vector. In the case of strain, the in-plane rotational symmetry is broken by terms to first order in both strain and wave vector leading to a Dirac cone with elliptical contours of constant energy, here shown for ò xy =5%. Hence, the group velocity of the surface electrons becomes dependent on the direction. The broken rotational symmetry also allows for a spin component perpendicular to the surface, as shown by the red arrows. complete discussion see [12]. In this basis the σ matrices do not represent spin, but are related to spin by: as noted in [14].
In this basis, the transformation operators corresponding to the symmetries of the crystal are inversion I=σ 0 ⊗τ 3 , three-fold rotation around the z-axis , two fold rotation around the x-axis R 2 =i σ 1 ⊗τ 3 and time-reversal T=i(σ 2 ⊗τ 0 )K, where K is the complex conjugation operator. Now the Γ-matrices can be characterized by their irreducible representation under the group D 3 as well as their eigenvalues under I and T.

Group theory
Next, we construct polynomials in wave vector, strain, electric and magnetic fields transforming according to the irreducible representations of the group D 3d . This group is a direct product group of D 3 with the group of the inversion operator. Here we will consider D 3 separately and simply add inversion eigenvalues as well as timereversal eigenvalues. The group D 3 has three irreducible representations: Γ (1) , which is the trivial representation, Γ (2) , which is one-dimensional as well and the two-dimensional representation Γ (3) .
The wave vector component k z transforms according to the one-dimensional representation Γ (2) , i.e. it only changes by multiplication of a constant under the transformations of D 3 . The in-plane components, k x and k y , transform according to the two-dimensional irreducible representation Γ (3) . Here we will use the basis {k + , k − }, where k ± =k x ±ik y . As an example we will show how to construct the second order terms transforming according to irreducible representations including only in-plane momentum. First we note that i.e. the second order terms of in-plane momenta transform according to a representation equivalent to the sum of these three irreducible representations. The basis functions for the irreducible representations can be formed using table 1(c), with {u 1 , u 2 }={v 1 , v 2 }={k + , k − }, and we get:   (a) 2 . Since {u 1 , u 2 }={v 1 , v 2 } the Γ (2) term is simply 0, and we have 3 independent second order terms in the in-plane momenta: || k 2 , which is invariant, and the pair -+ { } k k , 2 2 transforming according to the two-dimensional irreducible representation Γ (3) . The inversion and time-reversal eigenvalues follow the simple rule ±⊗±=+ and ±⊗m=−. Proceeding like this we can go to any desired order, and include external fields as well. We note that the strain tensor ò ij transforms as k i k j [15]. This result follows from the definition of strain where u is the displacement vector. The resulting irreducible basis functions are summarized in table 2.

Model Hamiltonian
For the Hamiltonian to be invariant under the group D 3 it can only contain terms from the invariant representation, Γ (1) . As seen in the multiplication table 1(a), Γ (1) terms can only be formed by combining terms from the same irreducible representation. For invariance under inversion (I) and time-reversal (TR), we must combine terms and matrices with the same eigenvalues, i.e. from the same cell in table 2. For the Γ (3) representation we can use table 1(c) to construct an invariant term, Γ (1) and Γ (2) are 1 dimensional so we simply take the product. For example, the terms - , transform exactly the same way, and using table 1(c) we can therefore combine them as follows: z z z 1 1 2 1,2 2 1,2 2 1,2 2 1,2 Table 2. Polynomials of k, strain and magnetic fields and the Γ matrices under the transformations of the group D 3 , inversion and time reversal.
To simplify the notation we have introduced    = + || xx yy , ò ± =ò xx −ò yy ±2iò xy , ò z± =ò zx ±iò zy . In the 2-dimensional representations we have changed basis to G = G  G . These have been constructed such that they are real, and each pair is hermitian conjugates of each other and the pairs transform like {k + , k − } under the group D 3 .
This forms an invariant, and thus it is an allowed term in the Hamiltonian. Any multiple of this term is of course also an invariant, and we get a free parameter in front, which we denote Z 3 . Proceeding this way we achieve the following Hamiltonian, for now excluding electric and magnetic fields: where repeated indices are summed over and where: yy , ò ± =ò xx −ò yy ±2iò xy , ò z± =ò zx ±iò zy . From the method of invariants it follows that all model parameters designated by capital roman letters must be real parameters, which are not determined by the symmetries of the crystal. The Hamiltonian above contains all allowed terms to third order in k, first order in strain and terms first order in both. Neglecting strain, this Hamiltonian reduces to the one derived in [12], except for three terms combining the in-plane and out-of-plane components of the wave vector, 2 1,2 , which were not included in [12]. The k independent strain terms depend only on  || and ò zz , and are either proportional to the identity matrix or Γ 5 . Hence, the effects of these terms can be described within the unstrained model by making the parameters strain-dependent Similarly the terms linear in k and  || or ò zz , can be described by strain-dependent model parameters in the unstrained model. The shear strain terms and the terms with ò xx −ò yy give new terms linear in k, which cannot be described by making the parameters of the unstrained model strain-dependent. The effects of ò zx are similar to the effects of ò xy , since they contribute to the real parts of β 1 , α 2 and α 3 and the imaginary parts of α 1 and β 2 .
For the same reason ò xx −ò yy and ò zy have similar effects. As we demonstrate below, the new terms linear in the wave vector give rise to new physics at the (111) surface of a semi-infinite topological insulator. Higher order terms in strain will always be both inversion and time-reversal symmetric, and in table 2 we see that the only matrices symmetric under both I and TR are the identity matrix and Γ 5 , and k-independent strain terms may therefore be lumped into the strain-dependent parameters  ( ) M and ( ) C .

Magnetic field
The Hamiltonian in equation (10) is invariant under both TR and inversion, thus all bands are doubly degenerate. The TR symmetry can be broken by a magnetic field, which must be included in the following form: where g 1z,2z,1z,2z are real parameters and B ± =B x ±iB y . This contribution was discussed in [12]. The Hamiltonian in equation (13) is responsible for the Zeeman effect. Orbital effects can be included by invoking minimal coupling,   + k k A e , with vector potential A, in the strained Hamiltonian equations (10)-(a), and the effects of strain on the Landau levels can be calculated along the lines of Liu et al [12]. We shall not pursue these effects further in this work.

Electric field
Inversion symmetry can be broken by an electric field, giving rise to a Hamiltonian of the form: where || W z , are real parameters and E ± =E x ±iE y .

Bulk spectrum and band-gap
In this section we will analyze the effects of strain and electric field on the bulk band structure close to the Γ point. The Hamiltonian H=H m +H E of equation (10) and (14) can be diagonalized analytically giving:   Note that the gap can only be increased by applying an electric field, whereas strain can increase or decrease the gap depending on the sign of the strain. The effects of strain have been investigated earlier by DFT calculations. In [7] both Bi 2 Se 3 and Bi 2 Te 3 were investigated. Here the lattice constant in the xy-plane or the z-direction was fixed relative to the known unstrained lattice constant, and the lattice constant in the other direction was relaxed before the calculations were performed. An approximate linear relationship between the in-plane, and out-of-plane Without strain the in-plane rotational symmetry is broken only by terms to third order in the wave vector, and only a small splitting between different in-plane direction is seen. Including strain allows terms to first order in the wave vector that breaks the in-plane rotational symmetry, and we see that a splitting of the in-plane directions occurs at lower k values. Here we have used the parameters of [12] for the unstrained model, and X i =Y i =10 eV Å. We note that ò xx −ò yy =10% gives the same dispersion for the parameters used here.
lattice constants was reported [7]. In [8] the band gap was calculated for many different strain configurations and then fitted to a second order polynomial of the strain components. In figure 3 the results of these calculations are shown together with second and third order fits to the results of [7]. The results show good agreement and both predict a topological phase transition at approximately 6% uniaxial strain. Similar results were reported in [9]. According to our model the topological phase transition will occur when: In general, the parameters of this symmetry-adapted model can be determined by using DFT calculations. The values of the parameters are highly dependent on the computational details as can be seen by comparing [1,12,17]. In this paper, we use the values from [12] for the parameters not related to strain. It should be emphasized that strain can only change the gap by a term proportional to the matrix Γ 5 in the Hamiltonian, which is equivalent to simply changing the band gap parameter M. In the following sections we focus on the on the response of the band gap under uniaxial strain, and include a second order contribution as well since this can easily be done using the result of [8]: with the parameters given in table 3.

Effective masses
An important quantity for transport and optical properties of a solid is the effective mass tensor. Here we calculate the diagonal components of the inverse effective mass tensor including both strain and electric field.  [7,8]. The curves from [8] are calculated using the approximate linear relationships between the in-plane and out-of-plane lattice constant reported in [7]. The fits to the data from [7] show good agreement, especially for Bi 2 Te 3 . Going to third order makes only a slight improvement compared to second order. The inverse effective mass tensor is given by: With E=0 we get the effective masses along the x and y axes: The plus (minus) sign corresponds to the conduction (valence) band. For a non-zero electric field each of the two bands is split into two subbands having the same effective mass, but differing by linear and cubic terms in the wave vector. The effective masses including both strain and electric field are given in the appendix.
In figure 4 the bulk effective masses of Bi 2 Se 3 are plotted as a function of the strain component ò zz and in the absence of electric and magnetic fields. Conduction (valence) band masses are shown as solid (dashed) lines and the effective mass superscript c (v) refers to the conduction (valence) band. In this case where only  ¹ 0 zz and the electric field is zero the effective masses along the x and y directions are the same. Note that since the maximum of the valence band shifts away from the Γ point at a sufficiently large strain value, the valence band effective masses change as a function of strain from positive to negative and reach infinite values in-between.

SS spectrum and spin structure
We shall now proceed to calculate the SS spectrum and spin structure for the model Hamiltonian (10). In this and the following section we exclude k 3 terms, i.e. the last three terms in equation (10) and the k dependent terms in α i and β i . First, we consider a semi-infinite topological insulator filling the z<0 half space, implemented via the boundary conditions Ψ(0)=0 and Y  ( ) z 0 for  -¥ z . The translational invariance is broken in the z-direction and we make the substitution  - ¶ k i  where D ± =D 1 ±B 1 . This equation has 4 complex solutions λ i , each with two corresponding eigenspinors ψ j (λ i ) given by: The general solution to the Schrödinger equation can then be written on the form: where the coefficients C ij are to be determined from the boundary conditions. For SS we require that Y  To satisfy the boundary condition at the surface we need two different solutions with positive real part. By complex conjugation of equation (22) we see that for a solution λ i , −λ i * will also be a solution. Hence, the solutions are either imaginary or occur in pairs related by * l l = -  . Existence of SSs is only possible in case (iii) and we therefore assume that λ 1 and λ 2 have positive real parts. Hence, we limit the sum in equation (25) to iä{1, 2}. Applying the boundary condition at z=0 we obtain the following secular equation for nontrivial solution to the coefficients C ij : 2 6   1  2  2  3  2  3  2 and combining this with equation (22) we find the SS spectrum: Re .

7
x y x y x y x y For small || k we see that the spectrum is linear but due to the strain, the group velocity has now become anisotropic and the contours of constant energy have become elliptical. Note that the position of the Dirac-point can be shifted by changing  ( ) M , however the relative position within the bulk band gap is not changed. Notice also that the second order term is not changed by strain, since we have systematically retained terms of order ò 1 k 1 and discarded terms of order ò 1 k 2 . Both the directional dependence of the group velocity and the ellipticity of the contours of constant energy are shown in figure 5. For = || k 0, the linear term in equation (22) vanishes and λ 1 and λ 2 can be calculated analytically: where: Imposing the boundary condition Ψ(0)=0 with equation (  Spin-momentum locking is still present, but spin and momentum are no longer perpendicular. Here we have used the parameters of [12] for the parameters of the unstrained model, and S z values are allowed. At finite || k , equation (22) was solved numerically, and using Ψ(0)=0 a numerical expression for the SS spinor was determined. Finally, the expectation values of the spin operators were calculated and are shown in figure 5. As expected, non-zero values of á ñ S z occur. Such out-of-plane spin components also show up in the absence of strain to third order in k in the case of hexagonal warping [2] and has been experimentally verified in [18]. The actual values of S z shown in figure 5(a) are of course strongly dependent on the chosen parameters, but the non-zero S z -component is generally present for any values of X i and Y i .

Effective 2D model
Using the SS at = || k 0 as basis we can derive an effective 2D model for the SS. We use the wave functions in equation (30), but with ò ij =0 to have a basis independent of strain. We note that in the absence of strain, ψ 1 and ψ 2 represent spin up/down states respectively. The matrix elements of the Γ-matrices in this basis read: The effective 2D model including k up to second order becomes: 3 4 x y zy xx yy x zx xy y Diagonalizing the effective Hamiltonian gives the spectrum: 3 5 x y x y This spectrum is different from equation (27), obtained from the 3D model. The two spectra become equal, however, if we take α 3 =0. We note that the Berry curvature of the SS band vanishes identically, i.e.

Localization of SS and finite size effect
In the case of a finite slab, a gap will open in the SS spectrum at =     l l l l l l for even/odd states respectively. This is the same form as found using the model without strain [19] , but here the parameters  ( ) M and ( ) C as well as λ 1,2 from equation (28) depend on strain. The parameter ( ) C is not important, since it simply shifts all energies. Again the important quantities are 2 . To find the gap of the SSs, Δ, equation (37) has been solved numerically for L=n·9.547 Å for integers n between 2 and 6 where L=9.547 Å is the thickness of 1 QL. First we have analyzed the effect of the bulk band gap changed by strain, and taking a b 1 the value without strain. We see that increasing the bulk band gap, decreases the SS band gap. Close to the phase transition this behavior is evident, since the SS must approach bulk states extending further into the material as the bulk gap approaches zero, giving a large coupling between opposite surfaces. As we saw in section 3, the bulk gap is most sensitive to ò zz . Using equation (18) we have plotted the SS gap as function of ò zz in figure 6(a). A similar result was reported in [9] for a slab of 8 QL using DFT, where a gap opens when tensile strain is applied. Our results points to the physical origin of this effect being the modification of the spatial distribution of the SS wave functions.
In figure 6( according to [12], but the quantitative dependence on strain has not been determined. We see that the SS gap closes at low values, which is related to the fact that for Table 4. Experimental measurements of the surface state gap, induced by coupling between opposite surfaces. First row is from [24], second row is from [25]. In the second row the gap increases from 2 to 3 QL, signifying an oscillatory dependence of the gap on the thickness. We have demonstrated that the gap of the SS can be tuned by strain. However, quantitative predictions require access to the parameters relating the strain tensor to a b + | | | | 3 2 3 2 . Tuning the band gap of the SS is important for both applications and fundamental research. One particular experimental challenge is that the layered structure of Bi 2 Se 3 makes it easy to fabricate integer numbers of QL. Thus the gap dependence on the layer thickness has only been experimentally studied with few data points, since the gap becomes too small to be measured already at 6 QL [24]. With strain, it is be possible to increase the SS gap and investigate the thickness dependence of the gap further. Another intriguing prospect would be tuning the gap via strain, so as to pass through the topological transition in a controlled manner.

Conclusion
We have derived the most general Hamiltonian for the Bi 2 Se 3 -class of materials including terms to third order in the wave vector, first order in electric and magnetic fields, first order in strain and first order in both strain and wave vector. We show that this model provides a description of a range of different effects of strain on the electronic structure of these materials. Specifically we have analytically derived the spectrum of the SS for a semi-infinite topological insulator, showing qualitatively the effects of strain on both the spectrum and the spin structure. The terms first order in both wave vector and strain break the full rotational symmetry close to = || k 0, leading to an anisotropy in the Dirac cone otherwise absent for the simple (111) surface termination. The spin structure is altered as well, and for some strain configurations a spin component perpendicular to the surface arises. In an analysis of the finite size effect we show that increasing the bulk band gap by virtue of strain, decreases the induced gap in the SS. Strain-tuning the SS gap gives new possibilities for experimental investigation of the thickness dependence of the SS gap and control of transport and optical properties.

Appendix. Effective masses
Here we give the expressions for the bulk effective masses including both contributions from strain and electric fields. Upper (lower) signs refer to the conduction (valence) band.