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

TL;DR

This paper uses group theory to derive a comprehensive Hamiltonian for Bi2Se3 topological insulators, analyzing how strain, electric, and magnetic fields influence their electronic and surface state properties, with implications for device applications.

Contribution

It provides the most general Hamiltonian including higher-order terms and analytically explores strain effects on electronic structure and surface states of Bi2Se3.

Findings

Strain causes anisotropic Dirac cones with elliptical constant energy contours.

Strain modifies spin-momentum locking and can induce a perpendicular spin component.

Surface state band gap varies oppositely to the bulk gap under strain.

Abstract

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.

Tables

Table 1. (a)

$\otimes$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(1\right)}\oplus {\Gamma }^{\left(2\right)}\oplus {\Gamma }^{\left(3\right)}$

Table 2. (a)

$\otimes$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(3\right)}$
${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(3\right)}$	${\Gamma }^{\left(1\right)}\oplus {\Gamma }^{\left(2\right)}\oplus {\Gamma }^{\left(3\right)}$

Table 3. (b)

	$u^{\left(2\right)}{v}_1^{\left(3\right)}$	$u^{\left(2\right)}{v}_2^{\left(3\right)}$
${\psi }_1^{\left(3\right)}$	i	0
${\psi }_2^{\left(3\right)}$	0	-i

Table 4. (c)

	$u_1^{\left(3\right)}{v}_1^{\left(3\right)}$	$u_1^{\left(3\right)}{v}_2^{\left(3\right)}$	$u_2^{\left(3\right)}{v}_1^{\left(3\right)}$	$u_2^{\left(3\right)}{v}_2^{\left(3\right)}$
${\psi }^{\left(1\right)}$	$0$	$\frac{1}{\sqrt{2}}$	$\frac{1}{\sqrt{2}}$	$0$
${\psi }^{\left(2\right)}$	$0$	$\frac{i}{\sqrt{2}}$	$\frac{-i}{\sqrt{2}}$	0
${\psi }_1^{\left(3\right)}$	0	0	0	1
${\psi }_2^{\left(3\right)}$	1	0	0	0

Table 5. Table 2: Polynomials of 𝐤 𝐤 \mathbf{k} , strain and magnetic fields and the Γ Γ \Gamma matrices under the transformations of the group D 3 subscript 𝐷 3 D_{3} , inversion and time reversal. To simplify the notation we have introduced ϵ | | = ϵ x ​ x + ϵ y ​ y \epsilon_{||}=\epsilon_{xx}+\epsilon_{yy} , ϵ ± = ϵ x ​ x − ϵ y ​ y ± 2 ​ i ​ ϵ x ​ y subscript italic-ϵ plus-or-minus plus-or-minus subscript italic-ϵ 𝑥 𝑥 subscript italic-ϵ 𝑦 𝑦 2 𝑖 subscript italic-ϵ 𝑥 𝑦 \epsilon_{\pm}=\epsilon_{xx}-\epsilon_{yy}\pm 2i\epsilon_{xy} , ϵ z ± = ϵ z ​ x ± i ​ ϵ z ​ y subscript italic-ϵ limit-from 𝑧 plus-or-minus plus-or-minus subscript italic-ϵ 𝑧 𝑥 𝑖 subscript italic-ϵ 𝑧 𝑦 \epsilon_{z\pm}=\epsilon_{zx}\pm i\epsilon_{zy} . In the 2-dimensional representations we have changed basis to Γ ± i , j = Γ i ± i ​ Γ j subscript superscript Γ 𝑖 𝑗 plus-or-minus plus-or-minus subscript Γ 𝑖 𝑖 subscript Γ 𝑗 \Gamma^{i,j}_{\pm}=\Gamma_{i}\pm i\Gamma_{j} . 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 − } subscript 𝑘 subscript 𝑘 \{k_{+},k_{-}\} under the group D 3 subscript 𝐷 3 D_{3} .

$T$	$I$	${\Gamma }^{\left(1\right)}$	${\Gamma }^{\left(2\right)}$	${\Gamma }^{\left(3\right)}$
		dim=1	dim=1	dim=2
-	-	$\begin{array}{c}\mathrm{Re}\left(k_{+}^3\right)\\ \mathrm{Im}\left(k_{+}{ϵ}_{z-}\right)\\ \mathrm{Re}\left(k_{+}{ϵ}_{+}\right)\\ {\Gamma }_4\end{array}$	$\begin{array}{c}{k}_z\\ k_z^3\\ k_{\left|\right|}^2{k}_z\\ \mathrm{Im}\left(k_{+}^3\right)\\ {ϵ}_{zz}{k}_z\\ {ϵ}_{\left|\right|}{k}_z\\ \mathrm{Re}\left(k_{+}{ϵ}_{z-}\right)\\ \mathrm{Im}\left(k_{+}{ϵ}_{+}\right)\\ {\Gamma }_3\end{array}$	$\begin{array}{c}\{k_{+},k_{-}\}\\ k_{\left|\right|}^2\{k_{+},k_{-}\}\\ k_z^2\{k_{+},k_{-}\}\\ \{i{k}_z{k}_{-}^2,-i{k}_z{k}_{+}^2\}\\ {ϵ}_{zz}\{k_{+},k_{-}\}\\ {ϵ}_{\left|\right|}\{k_{+},k_{-}\}\\ k_z\{{ϵ}_{z+},{ϵ}_{z-}\}\\ k_z\{i{ϵ}_{-},-i{ϵ}_{+}\}\\ \{i{k}_{-}{ϵ}_{z_{-}},-i{k}_{+}{ϵ}_{z+}\}\\ \{k_{-}{ϵ}_{+},k_{+}{ϵ}_{-}\}\\ \{{\Gamma }_{+}^{1,2},{\Gamma }_{-}^{1,2}\}\end{array}$
+	+	$\begin{array}{c}{k}_z^2\\ k_{\left|\right|}^2\\ {ϵ}_{zz}\\ {ϵ}_{\left|\right|}\\ {\Gamma }_5\end{array}$		$\begin{array}{c}\{i{k}_z{k}_{+},-i{k}_z{k}_{-}\}\\ \{k_{-}^2,k_{+}^2\}\\ \{i{ϵ}_{z+},-i{ϵ}_{z-}\}\\ \{{ϵ}_{-},{ϵ}_{+}\}\end{array}$
+	-	$\begin{array}{c}{\Gamma }_{45}\end{array}$	$\begin{array}{c}{E}_z\\ {\Gamma }_{35}\end{array}$	$\begin{array}{c}\{E_{+},E_{-}\}\\ \{i{\Gamma }_{+}^{25,51},-i{\Gamma }_{-}^{25,51}\}\end{array}$
-	+		$\begin{array}{c}{B}_z\\ {\Gamma }_{12}\\ {\Gamma }_{34}\end{array}$	$\begin{array}{c}\{B_{+},B_{-}\}\\ \{{\Gamma }_{+}^{13,23},{\Gamma }_{-}^{23,31}\}\\ \{{\Gamma }_{+}^{14,24},{\Gamma }_{-}^{14,24}\}\end{array}$

Table 6. Table 3: Experimental measurements of the surface state gap, induced by coupling between opposite surfaces. First row is from Ref.  Zhang  et al. , 2010 , second row is from Ref.  Sakamoto  et al. , 2010 . In the second row the gap increases from 2 to 3 QL, signifying an oscillatory dependence of the gap on the thickness.

$L$	2 QL	3 QL	4 QL	5 QL
$\Delta \left(\mathrm{eV}\right)$	$0.25$	$0.138$	$0.07$	$0.041$
$\Delta \left(\mathrm{eV}\right)$	$0.28$	$0.34$

Equations

$${\Gamma }_1={\sigma }_1\otimes {\tau }_1,{\Gamma }_2={\sigma }_2\otimes {\tau }_1,{\Gamma }_3={\sigma }_3\otimes {\tau }_1,$$

$${\Gamma }_4={\sigma }_0\otimes {\tau }_2,{\Gamma }_5={\sigma }_0\otimes {\tau }_3,$$

$$H_m=\epsilon \mathbf{I}+\sum _{i=1}^5{d}_i{\Gamma }_i+\sum _{i,j}{d}_{ij}{\Gamma }_{ij},$$

$$s_x=\frac{{\sigma }_1\otimes {\tau }_3}{2}{s}_y=\frac{{\sigma }_2\otimes {\tau }_3}{2}{s}_z=\frac{{\sigma }_3\otimes {\tau }_0}{2},$$

$${\Gamma }^{(3)}\otimes {\Gamma }^{(3)}={\Gamma }^{(1)}\oplus {\Gamma }^{(2)}\oplus {\Gamma }^{(3)},$$

$${\Gamma }^{(1)}$$

$${\Gamma }^{(2)}$$

$${\Gamma }^{(3)}$$

$${ϵ}_{ij}$$

$${\Gamma }^{(1)}\text{mathchar}58{\psi }^{(1)}$$

$$\propto k_z(i{k}_{-}^2{\Gamma }_{-}^{1,2}-i{k}_{+}^2{\Gamma }_{+}^{1,2}).$$

$$H_m={\mathcal{E}}_0$$

$$-Z_1\mathrm{Re}(k_{+}^3){\Gamma }_4+Z_2\mathrm{Im}(k_{+}^3){\Gamma }_3$$

$$+Z_3{k}_z(i{k}_{-}^2{\Gamma }_{-}^{1,2}-i{k}_{+}^2{\Gamma }_{+}^{1,2}),$$

$${\mathcal{E}}_0$$

$$\mathcal{M}$$

$${\alpha }_1$$

$$+i{Y}_3{ϵ}_{z-}+Y_4{ϵ}_{+},$$

$${\alpha }_2$$

$$+Y_3{ϵ}_{z-}-i{Y}_4{ϵ}_{+},$$

$${\alpha }_3$$

$${\beta }_1$$

$${\beta }_2$$

$${\beta }_3$$

$$\tilde{M}(ϵ)$$

$$\tilde{C}(ϵ)$$

$$H_B$$

$$H_E$$

$$=\left(\begin{array}{cccc}0& i{W}_z{E}_z& 0& i{W}_{||}{E}_{-}\\ -i{W}_z{E}_z& 0& -i{W}_{||}{E}_{-}& 0\\ 0& i{W}_{||}{E}_{+}& 0& -i{W}_z{E}_z\\ -i{W}_{||}{E}_{+}& 0& i{W}_z{E}_z& 0\end{array}\right),$$

$$E$$

$$+{\left|{\beta }_i{k}_i-i{Z}_1\mathrm{Re}(k_{+}^3)+Z_2\mathrm{Im}(k_{+}^3)\right|}^2+W_z^2{E}_z^2+W_{||}^2{E}_{||}^2$$

$$\displaystyle\pm\bigg{(}2\Big{|}\left(\alpha_{i}k_{i}+iZ_{3}k_{z}k_{-}^{2}\right)W_{||}E_{+}$$

$$\displaystyle+\left(\beta_{i}k_{i}-iZ_{1}\operatorname{Re}(k_{+}^{3})+Z_{2}\operatorname{Im}(k_{+}^{3})\right)W_{z}E_{z}\Big{|}^{2}$$

$$+4W_{||}^2{E}_{||}^2{\left|{\beta }_i{k}_i-i{Z}_1\mathrm{Re}(k_{+}^3)+Z_2\mathrm{Im}(k_{+}^3)\right|}^2$$

$$+4W_z^2{E}_z^2{\left|{\alpha }_i{k}_i+i{Z}_3{k}_z{k}_{-}^2\right|}^2$$

$$\displaystyle-2\operatorname{Re}\Big{(}\left(\alpha_{i}k_{i}+iZ_{3}k_{z}k_{-}^{2}\right)W_{||}E_{+}$$

$$\displaystyle-\left(\beta_{i}k_{i}-iZ_{1}\operatorname{Re}(k_{+}^{3})+Z_{2}\operatorname{Im}(k_{+}^{3})\right)W_{z}E_{z}\Big{)}\bigg{)}^{\frac{1}{2}}\Bigg{)}^{\frac{1}{2}}.$$

$${\Delta }_{\Gamma }=2\sqrt{(M+M_1{ϵ}_{zz}+M_2{ϵ}_{||}{)}^2+W_z^2{E}_z^2+W_{||}^2{E}_{||}^2}.$$

$${ϵ}_{zz}=-\frac{M}{M_1}-\frac{M_2}{M_1}{ϵ}_{||}.$$

$${\left[\frac{1}{m^{\ast }}\right]}_{ij}=\frac{1}{{\hslash }^2}\frac{{\partial }^{2}E}{\partial k_i\partial k_j}.$$

$${\left[\frac{{\hslash }^2}{m}\right]}_{ii}$$

$${\left[\frac{{\hslash }^2}{m}\right]}_{zz}$$

Full text

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 $\text{Bi}_{2}\text{Se}_{3}$-class of materials Zhang et al. (2009). This class, also known as the tetradymite group, contains compounds $\text{M}_{2}\text{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. 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 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 $(221)$ 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 $\Gamma$ point of bulk $\text{Bi}_{2}\text{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 $\text{Bi}_{2}\text{Se}_{3}$ where strain up to $3\%$ 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 $\Gamma$ 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 $\text{Bi}_{2}\text{Se}_{3}$ class.

Using this model, we determine analytical expressions for the modified bulk band structure and the effective mass tensor near the $\Gamma$-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 $\text{Bi}_{2}\text{Se}_{3}$ class, which now includes all allowed terms to orders cubic in the wave vector $\mathbf{k}$, linear in the strain tensor $\epsilon_{ij}$, linear in magnetic field $\mathbf{B}$ or electric field $\mathbf{E}$, and linear in both $\mathbf{k}$ and $\epsilon_{ij}$.

In general any $4\times 4$ Hamiltonian can be written in terms of Dirac $\Gamma$ matrices defined in terms of the Pauli matrices by:

[TABLE]

and their commutators $\Gamma_{ij}=\tfrac{1}{2i}[\Gamma_{i},\Gamma_{j}]$. A general $4\times 4$ matrix can then be written:

[TABLE]

where the coefficients $\varepsilon,d_{i},d_{ij}$ must be real to ensure hermiticity. In the present case the coefficients are functions of $\mathbf{k}$, $\epsilon_{ij}$, $\mathbf{B}$ and $\mathbf{E}$.

Here we use the basis $|P1_{-}^{+},\tfrac{1}{2}\rangle$, $-i|P2_{+}^{-},\tfrac{1}{2}\rangle$, $|P1_{-}^{+},-\tfrac{1}{2}\rangle$, $i|P2_{+}^{-},-\tfrac{1}{2}\rangle$, derived from the $p$ orbitals of the Bi and Se atoms. The upper sign denotes the inversion eigenvalue and the $\pm\tfrac{1}{2}$ the total angular momentum. For a more complete discussion see Ref. Liu et al., 2010a. In this basis the $\sigma$ matrices do not represent spin, but are related to spin by:

[TABLE]

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

In this basis, the transformation operators corresponding to the symmetries of the crystal are inversion $I=\sigma_{0}\otimes\tau_{3}$, three-fold rotation around the $z$-axis $R_{3}=e^{i\pi\sigma_{3}/3}\otimes\tau_{0}$, two fold rotation around the $x$-axis $R_{2}=i\sigma_{1}\otimes\tau_{3}$ and time-reversal $T=i(\sigma_{2}\otimes\tau_{0})K$, where $K$ is the complex conjugation operator. Now the $\Gamma$-matrices can be characterized by their irreducible representation under the group $D_{3}$ as well as their eigenvalues under $I$ and $T$.

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 $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 time-reversal eigenvalues. The group $D_{3}$ has three irreducible representations: $\Gamma^{(1)}$, which is the trivial representation, $\Gamma^{(2)}$, which is one-dimensional as well and the two-dimensional representation $\Gamma^{(3)}$.

The wave vector component $k_{z}$ transforms according to the one-dimensional representation $\Gamma^{(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 $\Gamma^{(3)}$. Here we will use the basis $\{k_{+},k_{-}\}$, where $k_{\pm}=k_{x}\pm 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

[TABLE]

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 $\{u_{1},u_{2}\}=\{v_{1},v_{2}\}=\{k_{+},k_{-}\}$, and we get:

[TABLE]

where we have defined $k_{||}^{2}=k_{x}^{2}+k_{y}^{2}$. Since $\{u_{1},u_{2}\}=\{v_{1},v_{2}\}$ the $\Gamma^{(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_{-}^{2},k_{+}^{2}\}$ transforming according to the two-dimensional irreducible representation $\Gamma^{(3)}$. The inversion and time-reversal eigenvalues follow the simple rule $\pm\otimes\pm=+$ and $\pm\otimes\mp=-$. Proceeding like this we can go to any desired order, and include external fields as well. We note that the strain tensor $\epsilon_{ij}$ transforms as $k_{i}k_{j}$Bir and Pikus (1974). This result follows from the definition of strain

[TABLE]

where ${\bf u}$ is the displacement vector. The resulting irreducible basis functions are summarized in Table 2.

II.2 Model Hamiltonian

For the Hamiltonian to be invariant under the group $D_{3}$ it can only contain terms from the invariant representation, $\Gamma^{(1)}$. As seen in the multiplication Table LABEL:sub@table:multiplicationtable, $\Gamma^{(1)}$ 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 $\Gamma^{(3)}$ representation we can use Table LABEL:sub@table:coupling33 to construct an invariant term, $\Gamma^{(1)}$ and $\Gamma^{(2)}$ are 1 dimensional so we simply take the product.

For example, the terms $\{ik_{z}k_{-}^{2},-ik_{z}k_{+}^{2}\}$ belong to the irreducible representation $\Gamma^{(3)}$ and are odd under both TR and I. The pair of matrices $\{\Gamma^{1,2}_{+},\Gamma^{1,2}_{-}\}$, defined as $\Gamma^{i,j}_{\pm}=\Gamma_{i}\pm i\Gamma_{j}$, transform exactly the same way, and using Table LABEL:sub@table:coupling33 we can therefore combine them as follows:

[TABLE]

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:

[TABLE]

where repeated indices are summed over and where:

[TABLE]

with $\epsilon_{||}=\epsilon_{xx}+\epsilon_{yy}$, $\epsilon_{\pm}=\epsilon_{xx}-\epsilon_{yy}\pm 2i\epsilon_{xy}$, $\epsilon_{z\pm}=\epsilon_{zx}\pm i\epsilon_{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 $\mathbf{k}$, 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, $k_{z}^{2}(k_{+}\Gamma^{1,2}_{-}+k_{-}\Gamma^{1,2}_{+})$, $k_{||}^{2}k_{z}\Gamma_{3}$ and $k_{z}(ik_{-}^{2}\Gamma^{1,2}_{-}-ik_{+}^{2}\Gamma^{1,2}_{+})$, which were not included in Ref. Liu et al., 2010a. The $\mathbf{k}$ independent strain terms depend only on $\epsilon_{||}$ and $\epsilon_{zz}$, and are either proportional to the identity matrix or $\Gamma_{5}$. Hence, the effects of these terms can be described within the unstrained model by making the parameters strain dependent $M\rightarrow\tilde{M}(\mathbf{\epsilon})$ and $C\rightarrow\tilde{C}(\mathbf{\epsilon})$, where:

[TABLE]

Similarly the terms linear in $\mathbf{k}$ and $\epsilon_{||}$ or $\epsilon_{zz}$, can be described by strain-dependent model parameters in the unstrained model. The shear strain terms and the terms with $\epsilon_{xx}-\epsilon_{yy}$ give new terms linear in $\mathbf{k}$, which cannot be described by making the parameters of the unstrained model strain-dependent. The effects of $\epsilon_{zx}$ are similar to the effects of $\epsilon_{xy}$, since they contribute to the real parts of $\beta_{1}$, $\alpha_{2}$ and $\alpha_{3}$ and the imaginary parts of $\alpha_{1}$ and $\beta_{2}$. For the same reason $\epsilon_{xx}-\epsilon_{yy}$ and $\epsilon_{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 $\Gamma_{5}$, and $k$-independent strain terms may therefore be lumped into the strain dependent parameters $\tilde{M}(\mathbf{\epsilon})$ and $\tilde{C}(\mathbf{\epsilon})$.

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:

[TABLE]

where $g_{1p,2p,1z,2z}$ are real parameters and $B_{\pm}=B_{x}\pm iB_{y}$. 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:

[TABLE]

where $W_{||,z}$ are real parameters and $E_{\pm}=E_{x}\pm iE_{y}$.

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 $\Gamma$ point. The Hamiltonian $H=H_{m}+H_{E}$ of Eq. (10) and (14) can be diagonalized analytically giving:

[TABLE]

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 $x$ and $y$ directions closer to the $\Gamma$ point compared to the case without strain, where the bands in the $x$ and $y$ directions are only split by the $\mathbf{k}^{3}$ terms. The electric field breaks the inversion symmetry, and the degeneracies of the bulk conduction and valence bands are split. At the $\Gamma$ 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 $\Gamma$ point, given by:

[TABLE]

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 $\text{Bi}_{2}\text{Se}_{3}$ and $\text{Bi}_{2}\text{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 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 $6\%$ uni-axial strain. According to our model the topological phase transition will occur when:

[TABLE]

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

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:

[TABLE]

With $\mathbf{E}=0$ we get the effective masses along the $x$ and $y$ axes:

[TABLE]

for $i=x,y$ and along the $z$ axis

[TABLE]

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.

In Fig. 4 the bulk effective masses of Bi2Se3 are plotted as a function of the strain component $\epsilon_{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 $\epsilon_{zz}\neq 0$ 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 $\Gamma$ 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 $\mathbf{k}^{3}$ terms, i.e. the last three terms in Eq. (10) and the $\mathbf{k}$ dependent terms in $\alpha_{i}$ and $\beta_{i}$. First, we consider a semi-infinite topological insulator filling the $z<0$ half space, implemented via the boundary conditions $\Psi(0)=0$ and $\Psi(z)\rightarrow 0$ for $z\rightarrow-\infty$. The translational invariance is broken in the $z$ direction and we make the substitution $k_{z}\rightarrow-i\partial_{z}$. Using the ansatz $\psi_{\lambda}e^{z\lambda}$ the stationary Schrödinger equation implies the following secular equation:

[TABLE]

where $D_{\pm}=D_{1}\pm B_{1}$. This equation has 4 complex solutions $\lambda_{i}$, each with two corresponding eigenspinors $\psi_{j}(\lambda_{i})$ given by:

[TABLE]

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

[TABLE]

where the coefficients $C_{ij}$ are to be determined from the boundary conditions. For surface states we require that $\Psi(z)\rightarrow 0$ for $z\rightarrow-\infty$, implying that $\operatorname{Re}(\lambda_{i})>0$. 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 $\lambda_{i}$, $-\lambda_{i}^{*}$ will also be a solution. Hence, the solutions are either imaginary or occur in pairs related by $\lambda_{i}=-\lambda_{j}^{*}$. This allows for three distinct cases: (i) All solutions are imaginary. (ii) Two complex solutions related by $\lambda_{1}=-\lambda_{2}^{*}$ and two imaginary solutions. (iii) Two pairs of complex solutions $\lambda_{1}=-\lambda_{3}^{*}$ and $\lambda_{2}=-\lambda_{4}^{*}$. Existence of SSs is only possible in case (iii) and we therefore assume that $\lambda_{1}$ and $\lambda_{2}$ have positive real parts. Hence, we limit the sum in Eq. (24) to $i\in\{1,2\}$. Applying the boundary condition at $z=0$ we obtain the following secular equation for nontrivial solution to the coefficients $C_{ij}$:

[TABLE]

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

[TABLE]

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 $\tilde{M}(\mathbf{\epsilon})$, 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 $\epsilon^{1}k^{1}$ and discarded terms of order $\epsilon^{1}k^{2}$. Both the directional dependence of the group velocity and the ellipticity of the contours of constant energy are shown in Fig. 5.

For $k_{||}=0$, the linear term in Eq. (21) vanishes and $\lambda_{1}$ and $\lambda_{2}$ can be calculated analytically:

[TABLE]

where:

[TABLE]

Imposing the boundary condition $\Psi(0)=0$ with Eq. (24) we find the exact wave functions of the two degenerate eigenstates at $k_{||}=0$, which are related by TR:

[TABLE]

Note that the $\alpha_{3}$ term couples the spin blocks and in contrast to the case without strain the eigenstates at the $k_{||}=0$ point are mixed up/down spinstates. The conditions for existence of surface states at $k_{||}=0$ are:

[TABLE]

As expected, the surface states can only be destroyed by changing the sign of $\tilde{M}(\mathbf{\epsilon})$, i.e. closing and reopening the bulk band gap.

Due to TR-symmetry the spin of the SS at $(k_{x},k_{y})$ must be opposite of the spin of the SS at $(-k_{x},-k_{y})$. For a $(111)$ surface without strain, the model Hamiltonian is invariant under rotations in the $xy$-plane of any angle. Hence $\langle S_{z}\rangle=0$ since $(k_{x},k_{y})$ and $(-k_{x},-k_{y})$ are related by a rotation around the $z$-axis which does not change $S_{z}$. The terms linear in both wave vector and strain break the full rotation symmetry, and non-zero $\langle S_{z}\rangle$ values are allowed. At finite $k_{||}$, Eq. (21) was solved numerically, and using $\Psi(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 $\langle S_{z}\rangle$ occur. Such out-of-plane spin components also show up to third order in $k$ in the case of hexagonal warpingFu (2009). The actual values of $S_{z}$ shown in Fig. LABEL:sub@fig:spin_strainxy 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}$.

IV.1 Effective 2D model

Using the surface states at $k_{||}=0$ as basis we can derive an effective 2D model for the surface states. We use the wave functions in Eq. (29), but with $\epsilon_{ij}=0$ to have a basis independent of strain. We note that in the absence of strain, $\psi_{1}$ and $\psi_{2}$ represent spin up/down states respectively. The matrix elements of the $\Gamma$-matrices in this basis read:

[TABLE]

The effective 2D model including $\mathbf{k}$ up to second order becomes:

[TABLE]

where $\mathbf{d}(\mathbf{k})=(d_{1}(\mathbf{k}),d_{2}(\mathbf{k}),d_{3}(\mathbf{k}))$ is given by:

[TABLE]

Diagonalizing the effective Hamiltonian gives the spectrum:

[TABLE]

This spectrum is different from Eq. (26), obtained from the 3D model. The two spectra become equal, however, if we take $\alpha_{3}=0$. We note that the Berry curvature of the SS band vanishes identically, i.e. $\Omega(\mathbf{k})=\hat{d}\cdot[(\partial_{k_{x}}\hat{d})\times(\partial_{k_{y}}\hat{d})]=0$.

V Localization of surface states and finite size effect

In the case of a finite slab, a gap will open in the SS spectrum at $k_{||}=0$. This is due to the overlap of SS on opposite surfaces. For a given thickness the overlap will be highly dependent on the penetration depth of the SS. From the wave functions the expectation value of the distance from the surface $d=-\langle z\rangle$ can be calculated to:

[TABLE]

and we see that $d$ depends on strain through the band gap parameter, and the quantity $\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}$ involving the coefficients of the terms first order in $k_{z}$.

The eigenenergies in the finite system can be calculated by imposing the boundary condition $\Psi(\pm L/2)=0$, where $L$ is the slab thickness. Now the wave function can be written as Eq. (24), but using all 4 solutions to Eq. (21). At the $\Gamma$ point, the solutions are $\pm\lambda_{1,2}$ with $\lambda_{1,2}$ given by Eq. (27). The secular equation of the nontrivial solution for the coefficients gives:

[TABLE]

for even/odd states respectively. This is the same form as found using the model without strain Shan et al. (2010) , but here the parameters $\tilde{M}(\mathbf{\epsilon})$ and $\tilde{C}(\mathbf{\epsilon})$ as well as $\lambda_{1,2}$ from Eq. (27) depend on strain. The parameter $\tilde{C}(\mathbf{\epsilon})$ is not important, since it simply shifts all energies. Again the important quantities are $\tilde{M}(\mathbf{\epsilon})$ and $\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}$.

To find the gap of the SSs, $\Delta$, Eq. (36) has been solved numerically for $L=n\cdot$9.547\text{,}\mathrm{\SIUnitSymbolAngstrom} for integers $n$ between 2 and 6 where $L=$9.547\text{\,}\mathrm{\SIUnitSymbolAngstrom} is the thickness of 1 QL. First we have analyzed the effect of the bulk band gap changed by strain, and taking $\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}=A_{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 III, the bulk gap is most sensitive to $\epsilon_{zz}$. Using the relation

[TABLE]

from Ref. Young et al., 2011 with $2|\tilde{M}(\mathbf{\epsilon})|=E_{g}$ we have plotted the SS gap as function of $\epsilon_{zz}$ in Fig. 6a, using the value $M=$0.28\text{,}\mathrm{eV}$$ from Ref. Liu et al., 2010a for the band gap in unstrained $\text{Bi}_{2}\text{Se}_{3}$ and setting $A_{11}=0$.

In Fig. 6b, we have plotted the dependence of the SS gap on $\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}$ assuming a constant bulk band gap parameter of $M=$0.28\text{,}\mathrm{eV} from Ref. Liu *et al.*, [2010a](#bib.bib11) and setting $M_{1}=M_{2}=0$. Without strain $\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}=A_{1}=$2.26\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom} according to Ref. Liu et al., 2010a, 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 $|\alpha_{3}|^{2}+|\beta_{3}|^{2}<\frac{4\tilde{M}|D_{+}D_{-}|}{B_{1}}$ the solutions to Eq. (21), $\lambda_{i}$, are complex giving an oscillatory spatial dependence of the SS. Oscillatory surface states give oscillations in the SS gap as a function of $L$, which has been theorized to signify topological phase transitions between a trivial 2D insulator and a quantum spin hall phaseLu et al. (2010); Liu et al. (2010b); Linder et al. (2009). We note here that these oscillations are present only using the parameters of Ref. Zhang et al., 2009, not the parameters of Ref. Liu et al., 2010a as shown in Ref. Jensen, 2014. Experimental evidence of oscillations in the gap remains inconclusive, since both a simple decay of the gap for 2-5 QL has been reportedZhang et al. (2010) as well as an increase in the gap from 2 to 3 QLSakamoto et al. (2010) signifying a topological phase transition between 2 and 3 QL. The observed gaps are summarized in Table 3.

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 $\sqrt{|\alpha_{3}|^{2}+|\beta_{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 $\text{Bi}_{2}\text{Se}_{3}$ makes it easy to fabricate integer numbers of quintuple layers. 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 QLZhang et al. (2010). 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.

VI Conclusion

We have derived 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 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 surface states 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 surface states. The possibility of 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.

Acknowledgements.

MRJ and MW gratefully acknowledge financial support from the Danish Council of Independent Research (Natural Sciences) grant no.: DFF-4181-00182. AML gratefully acknowledges financial support from the Carlsberg Foundation. The Center for Quantum Devices is funded by the Danish National Research Foundation.

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

[TABLE]

[TABLE]

[TABLE]