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

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.

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:

	$$\displaystyle\Gamma_{1}=\sigma_{1}\otimes\tau_{1},\quad\Gamma_{2}=\sigma_{2}\otimes\tau_{1},\quad\Gamma_{3}=\sigma_{3}\otimes\tau_{1},$$
	$$\displaystyle\Gamma_{4}=\sigma_{0}\otimes\tau_{2},\quad\Gamma_{5}=\sigma_{0}\otimes\tau_{3},$$		(1)

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

$\displaystyle H_{m}=\varepsilon\mathbf{I}+\sum_{i=1}^{5}d_{i}\Gamma_{i}+\sum_{i,j}d_{ij}\Gamma_{ij},$

(2)

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:

$\displaystyle s_{x}=\frac{\sigma_{1}\otimes\tau_{3}}{2}\quad s_{y}=\frac{\sigma_{2}\otimes\tau_{3}}{2}\quad s_{z}=\frac{\sigma_{3}\otimes\tau_{0}}{2},$

(3)

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$.

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:

	$\displaystyle E$	$\displaystyle=\mathcal{E}_{0}\pm\Bigg{(}\mathcal{M}^{2}+\mathinner{\!\left\lvert\alpha_{i}k_{i}+iZ_{3}k_{z}k_{-}^{2}\right\rvert}^{2}$
		$\displaystyle+\left|\beta_{i}k_{i}-iZ_{1}\operatorname{Re}(k_{+}^{3})+Z_{2}\operatorname{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}$
		$\displaystyle+4W_{||}^{2}E_{||}^{2}\mathinner{\!\left\lvert\beta_{i}k_{i}-iZ_{1}\operatorname{Re}(k_{+}^{3})+Z_{2}\operatorname{Im}(k_{+}^{3})\right\rvert}^{2}$
		$\displaystyle+4W_{z}^{2}E_{z}^{2}\mathinner{\!\left\lvert\alpha_{i}k_{i}+iZ_{3}k_{z}k_{-}^{2}\right\rvert}^{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}}.$		(15)

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:

$\displaystyle\Delta_{\Gamma}=2\sqrt{(M+M_{1}\epsilon_{zz}+M_{2}\epsilon_{||})^{2}+W_{z}^{2}E_{z}^{2}+W_{||}^{2}E_{||}^{2}}.$

(16)

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:

$\displaystyle\epsilon_{zz}=-\frac{M}{M_{1}}-\frac{M_{2}}{M_{1}}\epsilon_{||}.$

(17)

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.

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:

	$\displaystyle\Leftrightarrow 0$	$\displaystyle=D_{+}D_{-}\lambda^{4}+\bigg{(}2(E-\tilde{C}(\mathbf{\epsilon})-D_{2}k_{||}^{2})D_{1}$
		$\displaystyle-2(\tilde{M}(\mathbf{\epsilon})-B_{2}k_{||}^{2})B_{1}+|\alpha_{3}|^{2}+|\beta_{3}|^{2}\bigg{)}\lambda^{2}$
		$\displaystyle+i(2\operatorname{Re}(\alpha_{1}\alpha_{3}^{*}+\beta_{1}\beta_{3}^{*})k_{x}+2\operatorname{Re}(\alpha_{2}\alpha_{3}^{*}+\beta_{2}\beta_{3}^{*})k_{y})\lambda$
		$\displaystyle+(E-\tilde{C}(\mathbf{\epsilon})-D_{2}k_{||}^{2})^{2}-(\tilde{M}(\mathbf{\epsilon})-B_{2}k_{||}^{2})^{2}$
		$\displaystyle-|\alpha_{1}k_{x}+\alpha_{2}k_{y}|^{2}-|\beta_{1}k_{x}+\beta_{2}k_{y}|^{2},$		(21)

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:

	$\displaystyle\psi_{1}(\lambda_{i})$	$\displaystyle=\begin{pmatrix}E-\mathcal{E}_{0}+\mathcal{M}\\ \beta_{1}k_{x}+\beta_{2}k_{y}-i\beta_{3}\lambda_{i}\\ 0\\ \alpha_{1}k_{x}+\alpha_{2}k_{y}-i\alpha_{3}\lambda_{i}\end{pmatrix},$		(22)
	$\displaystyle\psi_{2}(\lambda_{i})$	$\displaystyle=\begin{pmatrix}\beta_{1}^{*}k_{x}+\beta_{2}^{*}k_{y}-i\beta_{3}^{*}\lambda_{i}\\ E-\mathcal{E}_{0}-\mathcal{M}\\ \alpha_{1}k_{x}+\alpha_{2}k_{y}-i\alpha_{3}\lambda_{i}\\ 0\end{pmatrix}.$		(23)

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

$\displaystyle\Psi(z)=\sum_{j=1}^{2}\sum_{i=1}^{4}C_{ij}\psi_{j}(\lambda_{i})e^{\lambda_{i}z},$

(24)

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}$:

$\displaystyle(\lambda_{1}+\lambda_{2})^{2}=\frac{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}{-D_{+}D_{-}},$

(25)

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

	$\displaystyle E$	$\displaystyle=\tilde{C}(\mathbf{\epsilon})+\frac{D_{1}}{B_{1}}\tilde{M}(\mathbf{\epsilon})$
		$\displaystyle\pm\sqrt{1-\frac{D_{1}^{2}}{B_{1}^{2}}}\Big{(}|\alpha_{1}k_{x}+\alpha_{2}k_{y}|^{2}+|\beta_{1}k_{x}+\beta_{2}k_{y}|^{2}$
		$\displaystyle-\frac{\operatorname{Re}((\alpha_{1}k_{x}+\alpha_{2}k_{y})\alpha_{3}^{*}+(\beta_{1}k_{x}+\beta_{2}k_{y})\beta_{3}^{*})^{2}}{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}\Big{)}^{\frac{1}{2}}$
		$\displaystyle+\left(D_{2}-\frac{D_{1}}{B_{1}}B_{2}\right)k_{||}^{2}.$		(26)

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:

$\displaystyle\lambda_{\alpha}=\sqrt{\frac{-F+(-1)^{\alpha}\sqrt{R}}{2D_{+}D_{-}}},$

(27)

where:

	$\displaystyle F$	$\displaystyle=|\beta_{3}|^{2}+|\alpha_{3}|^{2}-2\tilde{M}(\mathbf{\epsilon})B_{1}+2ED_{1}-2\tilde{C}(\mathbf{\epsilon})D_{1},$		(28a)
	$\displaystyle R$	$\displaystyle=F^{2}-4D_{+}D_{-}(E-\tilde{C}(\mathbf{\epsilon})+\tilde{M}(\mathbf{\epsilon}))(E-\tilde{C}(\mathbf{\epsilon})-\tilde{M}(\mathbf{\epsilon})).$		(28b)

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:

	$\displaystyle\Psi_{1}$	$\displaystyle=N\begin{pmatrix}i\operatorname{sgn}(D_{+})\sqrt{\frac{|D_{+}|}{|2B_{1}|}}\\ \frac{\beta_{3}}{\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}}\sqrt{\frac{|D_{-}|}{|2B_{1}|}}\\ 0\\ \frac{\alpha_{3}}{\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}}\sqrt{\frac{|D_{-}|}{|2B_{1}|}}\end{pmatrix}(e^{\lambda_{1}z}-e^{\lambda_{2}z}),$		(29a)
	$\displaystyle\Psi_{2}$	$\displaystyle=N\begin{pmatrix}0\\ \frac{\alpha_{3}^{*}}{\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}}\sqrt{\frac{|D_{-}|}{|2B_{1}|}}\\ i\operatorname{sgn}(D_{+})\sqrt{\frac{|D_{+}|}{|2B_{1}|}}\\ \frac{-\beta_{3}^{*}}{\sqrt{|\alpha_{3}|^{2}+|\beta_{3}|^{2}}}\sqrt{\frac{|D_{-}|}{|2B_{1}|}}\end{pmatrix}(e^{\lambda_{1}z}-e^{\lambda_{2}z}).$		(29b)

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:

$\displaystyle D_{+}D_{-}<0\quad\text{and}\quad\tilde{M}(\mathbf{\epsilon})B_{1}>0.$

(30)

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.