Symmetry analysis of strain, electric and magnetic fields in the \text{Bi}_{2}\text{Se}_{3}-class of topological insulators

Symmetry analysis of strain, electric and magnetic fields in the -class of topological insulators

Mathias Rosdahl Jensen Technical University of Denmark, Department of Photonics Engineering, Kgs-Lyngby, 2800, Denmark    Jens Paaske Center for Quantum Devices, Niels Bohr Institute, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen, Denmark    Anders Mathias Lunde Center for Quantum Devices, Niels Bohr Institute, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen, Denmark    Morten Willatzen Beijing Institute of Nanoenergy and Nanosystems, Chinese Academy of Sciences, Beijing, China
and Technical University of Denmark, Department of Photonics Engineering, Kgs-Lyngby, 2800, Denmark

Based on group theoretical arguments we derive the most general Hamiltonian for the -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 . 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 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 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 BiSe class of materials as a function of strain and electric field.


I 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 -class of materials Zhang et al. (2009). This class, also known as the tetradymite group, contains compounds 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. For the simplest case of a 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 momentumFu (2009). For other surface terminations the Dirac cone of the surface states becomes anisotropic with elliptical curves of constant energy, as reported for the surface in Ref. Moon et al., 2011. Anisotropic Dirac fermions have interesting transport properties due to the different group velocity in different directions Trescher et al. (2015), and have attracted attention in other Dirac materials Feng et al. (2014); Lu et al. (2016).

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 Refs. Luo et al., 2012; Young et al., 2011 that strain strongly affects the band gap at the point of bulk , 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 where strain up to was induced by lattice mismatch Li et al. (2017).

In this paper, we use the method of invariants Lew Yan Voon and Willatzen (2009) 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 Liu et al. (2010a). Since this model is based solely on the symmetry of the crystal, it is valid for all materials in the class.

Figure 1: For unstrained 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 . 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.

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 Jensen et al. (). 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 surface states of a semi-infinite topological insulator. We analytically derive the spectrum of the surface states, 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 Fig. 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.

Ii Model Hamiltonian

In line with Ref. Liu et al., 2010a, we wish to derive a 4-band model Hamiltonian describing the class, which now includes all allowed terms to orders cubic in the wave vector , linear in the strain tensor , linear in magnetic field or electric field , and linear in both and .

In general any Hamiltonian can be written in terms of Dirac matrices defined in terms of the Pauli matrices by:


and their commutators . A general matrix can then be written:


where the coefficients must be real to ensure hermiticity. In the present case the coefficients are functions of , , and .

Here we use the basis , , , , derived from the orbitals of the Bi and Se atoms. The upper sign denotes the inversion eigenvalue and the the total angular momentum. For a more complete discussion see Ref. Liu et al., 2010a. In this basis the matrices do not represent spin, but are related to spin by:


as noted in Ref. Silvestrov et al., 2012.

In this basis, the transformation operators corresponding to the symmetries of the crystal are inversion , three-fold rotation around the -axis , two fold rotation around the -axis and time-reversal , where is the complex conjugation operator. Now the -matrices can be characterized by their irreducible representation under the group as well as their eigenvalues under and .

ii.1 Group theory

Next, we construct polynomials in wave vector, strain, electric and magnetic fields transforming according to the irreducible representations of the group . This group is a direct product group of with the group of the inversion operator. Here we will consider separately and simply add inversion eigenvalues as well as time-reversal eigenvalues. The group has three irreducible representations: , which is the trivial representation, , which is one-dimensional as well and the two-dimensional representation .

i 0
0 -i
0 0 0 1
1 0 0 0
Table 1: Multiplication table for the irreducible representations of the group , (a), and the coupling constants for the basis functions of the relevant products , (b) and , (c). Note that the multiplication table is completely general since it is in principle concerning equivalence classes of representations, while the coupling constants are only valid for a specific representation, here the representation with basis functions . Adopted from Ref. Koster et al., 1963.

The wave vector component transforms according to the one-dimensional representation , i.e. it only changes by multiplication of a constant under the transformations of . The in-plane components, and , transform according to the two-dimensional irreducible representation . Here we will use the basis , where . 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 LABEL:sub@table:coupling33, with , and we get:


where we have defined . Since the term is simply 0, and we have 3 independent second order terms in the in-plane momenta: , which is invariant, and the pair transforming according to the two-dimensional irreducible representation . The inversion and time-reversal eigenvalues follow the simple rule and . Proceeding like this we can go to any desired order, and include external fields as well. We note that the strain tensor transforms as Bir and Pikus (1974). This result follows from the definition of strain


where is the displacement vector. The resulting irreducible basis functions are summarized in Table 2.

dim=1 dim=1 dim=2
- -
+ +
+ -
- +
Table 2: Polynomials of , strain and magnetic fields and the matrices under the transformations of the group , inversion and time reversal. To simplify the notation we have introduced , , . In the 2-dimensional representations we have changed basis to . These have been constructed such that they are real, and each pair is hermitian conjugates of each other and the pairs transform like under the group .

ii.2 Model Hamiltonian

For the Hamiltonian to be invariant under the group it can only contain terms from the invariant representation, . As seen in the multiplication Table LABEL:sub@table:multiplicationtable, terms can only be formed by combining terms from the same irreducible representation. For invariance under inversion and time-reversal, we must combine terms and matrices with the same eigenvalues, i.e. from the same cell in Table 2. For the representation we can use Table LABEL:sub@table:coupling33 to construct an invariant term, and are 1 dimensional so we simply take the product.

For example, the terms belong to the irreducible representation and are odd under both TR and I. The pair of matrices , defined as , transform exactly the same way, and using Table LABEL:sub@table:coupling33 we can therefore combine them as follows:


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 . Proceeding this way we achieve the following Hamiltonian, for now excluding electric and magnetic fields:


where repeated indices are summed over and where: \cref@addtoresetequationparentequation


with , , . 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 , first order in strain and terms first order in both. Neglecting strain, this Hamiltonian reduces to the one derived in Ref. Liu et al., 2010a, except for three terms combining the in-plane and out-of-plane components of the wave vector, , and , which were not included in Ref. Liu et al., 2010a. The independent strain terms depend only on and , and are either proportional to the identity matrix or . Hence, the effects of these terms can be described within the unstrained model by making the parameters strain dependent and , where: \cref@addtoresetequationparentequation


Similarly the terms linear in and or , can be described by strain-dependent model parameters in the unstrained model. The shear strain terms and the terms with give new terms linear in , which cannot be described by making the parameters of the unstrained model strain-dependent. The effects of are similar to the effects of , since they contribute to the real parts of , and and the imaginary parts of and . For the same reason and have similar effects. As we demonstrate below, the new terms linear in the wave vector give rise to new physics at the 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 , and -independent strain terms may therefore be lumped into the strain dependent parameters and .

ii.3 Magnetic field

The Hamiltonian in Eq. (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 are real parameters and . This contribution was discussed in Ref. Liu et al., 2010a.

ii.4 Electric field

Inversion symmetry can be broken by an electric field, giving rise to a Hamiltonian of the form:


where are real parameters and .

Iii 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 of Eq. (10) and (14) can be diagonalized analytically giving:


In Fig. 2 an example of the band structure with and without strain is plotted. The terms in the Hamiltonian linear in both strain and wave vector split the bands in the and directions closer to the point compared to the case without strain, where the bands in the and directions are only split by the terms. The electric field breaks the inversion symmetry, and the degeneracies of the bulk conduction and valence bands are split. At the point the degeneracies are protected by TR symmetry which is not broken. Another important quantity of the electronic structure is the gap between the valence and conduction band at the point, given by:


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 Ref. Luo et al., 2012 both and were investigated. Here the lattice constant in the -plane or the -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 lattice constants was reported Luo et al. (2012). In Ref. Young et al., 2011 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 Ref. Luo et al., 2012. The results show good agreement and both predict a topological phase transition at approximately uni-axial strain. According to our model the topological phase transition will occur when:


It should be emphasized that strain can only change the gap by a term proportional to the matrix in the Hamiltonian, which is equivalent to simply changing the band gap parameter . In the following sections we will then simply use a strain-dependent band gap parameter , and use the parameters of Ref. Young et al., 2011 to get a quantitative strain dependence.

Figure 2: Band structure close to the gamma point, with (solid) and without strain (dashed) in the absence of electric and magnetic fields. 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 values. Here we have used the parameters of Ref. Liu et al., 2010a for the unstrained model, and . We note that gives the same dispersion for the parameters used here.

Band gap at point ()

Figure 3: The band gap at point as a function of strain according to the DFT calculations of Refs. Young et al., 2011; Luo et al., 2012. The curves from Ref. Young et al., 2011 are calculated using the approximate linear relationships between the in-plane and out-of-plane lattice constant reported in Ref. Luo et al., 2012. The fits to the data from Ref. Luo et al., 2012 show good agreement, especially for . Going to third order makes only a slight improvement compared to second order.

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 Refs. Zhang et al., 2009; Liu et al., 2010a; Nechaev and Krasovskii, 2016. In this paper,we use the values from Ref. Liu et al., 2010a for the parameters not related to strain.

iii.1 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. The inverse effective mass tensor is given by:


With we get the effective masses along the and axes:


for and along the axis


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 appendix A.

Figure 4: The bulk conduction and valence band effective masses of BiSe plotted as a function of the strain component and in the absence of electric and magnetic fields. In the calculations, we have used the DFT computed strain parameters from Ref. Luo et al., 2012 and Ref. Liu et al., 2010a for the parameters not related to strain. The solid (dashed) lines refer to the conduction (valence) band effective masses.

In Fig. 4 the bulk effective masses of BiSe are plotted as a function of the strain component and in the absence of electric and magnetic fields. Conduction (valence) band masses are shown as solid (dashed) lines and the effective mass superscript () refers to the conduction (valence) band. In this case where only and the electric field is zero the effective masses along the and 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.

Iv Surface state spectrum and spin structure

We shall now proceed to calculate the surface state (SS) spectrum and spin structure for the model Hamiltonian (10). In this and the following section we exclude terms, i.e. the last three terms in Eq. (10) and the dependent terms in and . First, we consider a semi-infinite topological insulator filling the half space, implemented via the boundary conditions and for . The translational invariance is broken in the direction and we make the substitution . Using the ansatz the stationary Schrödinger equation implies the following secular equation:


where . This equation has 4 complex solutions , each with two corresponding eigenspinors given by:


The general solution to the Schrödinger equation can then be written on the form:


where the coefficients are to be determined from the boundary conditions. For surface states we require that for , implying that . To satisfy the boundary condition at the surface we need two different solutions with positive real part. By complex conjugation of Eq. (21) we see that for a solution , will also be a solution. Hence, the solutions are either imaginary or occur in pairs related by . This allows for three distinct cases: (i) All solutions are imaginary. (ii) Two complex solutions related by and two imaginary solutions. (iii) Two pairs of complex solutions and . Existence of SSs is only possible in case (iii) and we therefore assume that and have positive real parts. Hence, we limit the sum in Eq. (24) to . Applying the boundary condition at we obtain the following secular equation for nontrivial solution to the coefficients :


and combining this with Eq. (21) we find the SS spectrum:


For small 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 , 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 and discarded terms of order . Both the directional dependence of the group velocity and the ellipticity of the contours of constant energy are shown in Fig. 5.

For , the linear term in Eq. (21) vanishes and and can be calculated analytically:


where: \cref@addtoresetequationparentequation


Imposing the boundary condition with Eq. (24) we find the exact wave functions of the two degenerate eigenstates at , which are related by TR: \cref@addtoresetequationparentequation


Note that the term couples the spin blocks and in contrast to the case without strain the eigenstates at the point are mixed up/down spinstates. The conditions for existence of surface states at are:


As expected, the surface states can only be destroyed by changing the sign of , i.e. closing and reopening the bulk band gap.




(a) Group velocity

(b) Unstrained