Hybridization of topological surface states with a flat band

TL;DR

This paper investigates how topological surface states interact with a flat bulk band, revealing transformations into different surface states and the emergence of Fermi arcs in Dirac semimetals through a combined theoretical model.

Contribution

It introduces a hybrid model combining 3D topological and 2D pseudospin Hamiltonians to explicitly analyze surface state evolution with flat band hybridization.

Findings

Hybridization converts Dirac-like states into branches below and above the flat band.

Upper surface state branch is identified as Dyakonov-Khaetskii states.

Fermi arcs emerge when the flat band crosses the conduction band.

Abstract

We address the problem of hybridization between topological surface states and a non-topological flat bulk band. Our model, being a mixture of three-dimensional Bernevig-Hughes-Zhang and two-dimensional pseudospin-1 Hamiltonian, allows explicit treatment of the topological surface state evolution by continuously changing the hybridization between the inverted bands and an additional "parasitic" flat band in the bulk. We show that the hybridization with a flat band lying below the edge of conduction band converts the initial Dirac-like surface states into a branch below and one above the flat band. Our results univocally demonstrate that the upper branch of the topological surface states is formed by Dyakonov-Khaetskii surface states known for HgTe since the 1980s. Additionally we explore an evolution of the surface states and the arising of Fermi arcs in Dirac semimetals when the flat band crosses the conduction band.

Equations

$$\hat{H}=\left(\begin{array}{cc}{\hat{H}}_0({\hat{k}}_x,{\hat{k}}_y,{\hat{k}}_z)& {\hat{H}}_z({\hat{k}}_z)\\ {\hat{H}}_z({\hat{k}}_z{)}^{†}& {\hat{H}}_0^{\ast }(-{\hat{k}}_x,-{\hat{k}}_y,-{\hat{k}}_z)\end{array}\right),$$

$${\hat{H}}_0=\left(\begin{array}{ccc}{C}_0+M_0& \hslash v_{\Vert }{\hat{k}}_{+}\sin \alpha & \hslash v_{\Vert }{\hat{k}}_{-}\cos \alpha \\ \hslash v_{\Vert }{\hat{k}}_{-}\sin \alpha & C_0+S_0& 0\\ \hslash v_{\Vert }{\hat{k}}_{+}\cos \alpha & 0& C_0-M_0\end{array}\right)+\mathcal{O}({\text{k}}^2)$$

$${\hat{H}}_z=\left(\begin{array}{ccc}0& 0& \hslash v_{\perp }{\hat{k}}_z\\ 0& 0& 0\\ \hslash v_{\perp }{\hat{k}}_z& 0& 0\end{array}\right)+\mathcal{O}({\text{k}}^2),$$

$$\begin{gathered} {\hslash }^2{v}^2(k_x^2+k_y^2)(E_{hh}{\cos }^2\alpha +E_{lh}{\sin }^2\alpha )+ \\ +{\hslash }^2{v}^2{k}_z^2{E}_{hh}=E_c{E}_{hh}{E}_{lh}, \end{gathered}$$

$$\left(\begin{array}{cc}\hat{A}+i\hat{B}-E& \hat{K}\\ {\hat{K}}^{†}& \hat{A}-i\hat{B}-E\end{array}\right)\left(\begin{array}{c}{\Phi }_1(x,y,z)\\ {\Phi }_2(x,y,z)\end{array}\right)=0,$$

$$\hat{A}=E_c-{\hslash }^2{v}^2{\hat{k}}_z\frac{1}{E_{lh}}{\hat{k}}_z-{\hslash }^2{v}^2{\hat{k}}_x\left(\frac{{\cos }^2\alpha }{E_{lh}}+\frac{{\sin }^2\alpha }{E_{hh}}\right){\hat{k}}_x-$$

$$-{\hslash }^2{v}^2{\hat{k}}_y\left(\frac{{\cos }^2\alpha }{E_{lh}}+\frac{{\sin }^2\alpha }{E_{hh}}\right){\hat{k}}_y,$$

$$\hat{B}={\hslash }^2{v}^2{\hat{k}}_y\left(\frac{{\cos }^2\alpha }{E_{lh}}-\frac{{\sin }^2\alpha }{E_{hh}}\right){\hat{k}}_x-$$

$$-{\hslash }^2{v}^2{\hat{k}}_x\left(\frac{{\cos }^2\alpha }{E_{lh}}-\frac{{\sin }^2\alpha }{E_{hh}}\right){\hat{k}}_y,$$

$$\hat{K}={\hslash }^2{v}^2{\hat{k}}_z\frac{\cos \alpha }{E_{lh}}{\hat{k}}_{-}-{\hslash }^2{v}^2{\hat{k}}_{-}\frac{\cos \alpha }{E_{lh}}{\hat{k}}_z.$$

$$E_c=C_0(x,y,z)+M_0(x,y,z)-E,$$

$$E_{hh}=C_0(x,y,z)+S_0(x,y,z)-E,$$

$$E_{lh}=C_0(x,y,z)-M_0(x,y,z)-E.$$

$${\hat{H}}_{2x2}(\mathbf{k})=d_0(\mathbf{k}){\sigma }_0+d_x(\mathbf{k}){\sigma }_x+d_y(\mathbf{k}){\sigma }_y+i{d}_z(\mathbf{k}){\sigma }_z,$$

$$M_0(z)=M_0^{(I)}+\left(M_0^{(II)}-M_0^{(I)}\right)\theta (z),$$

$$C_0(z)=C_0^{(I)}+\left(C_0^{(II)}-C_0^{(I)}\right)\theta (z),$$

$$S_0(z)=S_0^{(I)}+\left(S_0^{(II)}-S_0^{(I)}\right)\theta (z),$$

$${\Phi }_{1,2}^{(II)}\sim \exp \left(-{\lambda }_z^{(II)}z\right)\exp \left(i{k}_{x}x+i{k}_{y}y\right)\text{for}z>0,$$

$${\Phi }_{1,2}^{(I)}\sim \exp \left({\lambda }_z^{(I)}z\right)\exp \left(i{k}_{x}x+i{k}_{y}y\right),\text{for}z<0,$$

$${\lambda }_z^{(n)}=\sqrt{k^2\left(1+\frac{E_{lh}^{(n)}-E_{hh}^{(n)}}{E_{hh}^{(n)}}{\sin }^2\alpha \right)-\frac{E_c^{(n)}{E}_{lh}^{(n)}}{{\hslash }^2{v}^2}}.$$

$$\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\end{array}\right),\left(\begin{array}{cc}\frac{1}{E_{lh}(E)}\frac{\partial }{\partial z}& -i{k}_{-}\frac{\cos \alpha }{E_{lh}(E)}\\ i{k}_{+}\frac{\cos \alpha }{E_{lh}(E)}& \frac{1}{E_{lh}(E)}\frac{\partial }{\partial z}\end{array}\right)\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\end{array}\right).$$

$${\left({\lambda }_z^{(I)}{E}_{lh}^{(II)}+{\lambda }_z^{(II)}{E}_{lh}^{(I)}\right)}^2={\left(E_{lh}^{(II)}-E_{lh}^{(I)}\right)}^2{k}^2{\cos }^2\alpha .$$

$$\begin{array}{r}E=C_0\pm \sqrt{M_0^2+{\hslash }^2{v}^2(k_x^2+k_y^2)+{\hslash }^2{v}^2{k}_z^2},\\ E=C_0-M_0,\end{array}$$

$$A{\left({\hslash }^2{v}^2{k}^2\right)}^2-B{\hslash }^2{v}^2{k}^2+E^2(E+M_0^{(I)}{)}^2(E+M_0^{(II)}{)}^2=0,$$

$$\begin{gathered} A=\frac{(M_0^{(I)}-M_0^{(II)}{)}^2}{4}{\sin }^4\alpha + \\ +(E+M_0^{(I)})(E+M_0^{(II)}){\sin }^2\alpha , \\ B=(E^2-M_0^{(I)}{M}_0^{(II)})(E+M_0^{(I)})(E+M_0^{(II)}){\sin }^2\alpha + \\ +(E+M_0^{(I)}{)}^2(E+M_0^{(II)}{)}^2. \end{gathered}$$

$$E_{\mathrm{as}\pm }=-\frac{M_0^{(I)}+M_0^{(II)}}{2}\pm \frac{M_0^{(I)}-M_0^{(II)}}{2}\cos \alpha ,$$

$$M_0(y)=M_0^{(I)}+\left(M_0^{(II)}-M_0^{(I)}\right)\theta (y),$$

$$C_0(y)=C_0^{(I)}+\left(C_0^{(II)}-C_0^{(I)}\right)\theta (y),$$

$$S_0(y)=S_0^{(I)}+\left(S_0^{(II)}-S_0^{(I)}\right)\theta (y).$$

$${\Phi }_{1,2}^{(II)}\sim \exp \left(-{\lambda }_y^{(II)}y\right)\exp \left(i{k}_{x}x+i{k}_{z}z\right)\text{for}y>0,$$

$${\Phi }_{1,2}^{(I)}\sim \exp \left({\lambda }_y^{(I)}y\right)\exp \left(i{k}_{x}x+i{k}_{z}z\right),\text{for}y<0,$$

$${\lambda }_y^{(n)}=\sqrt{k_x^2+\frac{k_z^2-\frac{E_c^{(n)}{E}_{lh}^{(n)}}{{\hslash }^2{v}^2}}{1+\frac{E_{lh}^{(n)}-E_{hh}^{(n)}}{E_{hh}^{(n)}}{\sin }^2\alpha }},$$

$$\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\end{array}\right),\left(\begin{array}{cc}{R}_{+}\frac{\partial }{\partial y}+R_{-}{k}_x& k_z\frac{\cos \alpha }{E_{lh}}\\ k_z\frac{\cos \alpha }{E_{lh}}& R_{+}\frac{\partial }{\partial y}-R_{-}{k}_x\end{array}\right)\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\end{array}\right),$$

$$R_{\pm }(E)=\frac{{\cos }^2\alpha }{E_{lh}(E)}\pm \frac{{\sin }^2\alpha }{E_{hh}(E)}.$$

$$\begin{gathered} {\left({\lambda }_y^{(I)}{R}_{+}^{(I)}(E)+{\lambda }_y^{(II)}{R}_{+}^{(II)}(E)\right)}^2= \\ =k_z^2{\left(\frac{1}{E_{lh}^{(II)}}-\frac{1}{E_{lh}^{(I)}}\right)}^2{\cos }^2\alpha + \\ +k_x^2{\left(R_{-}^{(II)}(E)-R_{-}^{(I)}(E)\right)}^2. \end{gathered}$$

Full text

I Introduction

The research on topological materials constitutes one of the most active areas in modern condensed matter physics. Initiated by Kane and Mele, who introduced topology as a new property of two-dimensional (2D) insulators Kane and Mele (2005), their idea has subsequently been generalized to three-dimensional (3D) materials Fu et al. (2007), yielding a whole family of novel topological insulators (TIs). In general, the nontrivial topology of TIs arises from the inversion between two bands with opposite parity, resulting in the appearance of gapless states at the boundaries, which are insensitive to impurities and disorder Moore and Balents (2007); Roy (2009); Hasan and Kane (2010); Qi and Zhang (2011). In 2D TIs counter-propagating one-dimensional (1D) gapless states with opposite spin arise at the edges while 3D TIs feature helical surface states consisting of a single Dirac cone, where the spin points perpendicular to the momentum.

To date multiple 2D and 3D TIs have been experimentally verified (2D TIs Bernevig et al. (2006); König et al. (2007); Liu et al. (2008); Krishtopenko and Teppe (2018a); Krishtopenko et al. (2018); Wu et al. (2018); Krishtopenko and Teppe (2018b); 3D TIs Xia et al. (2009); Zhang et al. (2009); Chen et al. (2009); Hsieh et al. (2009); Fu (2011); Hsieh et al. (2012); Dziawa et al. (2012); Okada et al. (2013)), yet many materials with inverted band structure are bulk metals with additional helical surface states. The most prominent example for this case is the HgTe-class. Here the presence of the heavy-hole $|\Gamma_{8},\pm 3/2\rangle$ band in combination with an inverted pair of the electron $|\Gamma_{6},\pm 1/2\rangle$ and light-hole $|\Gamma_{8},\pm 1/2\rangle$ bands transform HgTe into a bulk semimetal. In order to observe the topologically nontrivial surface states, as predicted by the theoretical analysis based on the $Z_{2}$ topological invariant Fu and Kane (2007), it is necessary to open a gap between conduction $|\Gamma_{8},\pm 1/2\rangle$ and valence $|\Gamma_{8},\pm 3/2\rangle$ band. In the 3D case this can be achieved by applying tensile biaxial strain to the HgTe bulk film Fu and Kane (2007); Dai et al. (2008) which has been demonstrated experimentally in Ref. Brüne et al. (2011).

Theoretically the appearance of surface states in 2D and 3D TIs can be understood within the Bernevig-Hughes-Zhang (BHZ) model describing the band inversion in 2D Bernevig et al. (2006) and 3D systems Zhang et al. (2009). However, in addition to the Dirac-like surface states arising from the band inversion, gapless HgTe should also host parabolic Dyakonov-Khaetskii (DK) surface states, theoretically predicted in the 1980s Dyakonov and Khaetskii (1981). As the DK states are caused by the coupling between the light-hole $|\Gamma_{8},\pm 1/2\rangle$ and heavy-hole $|\Gamma_{8},\pm 3/2\rangle$ band, they should remain even in the presence of the strain-induced band gap Kibis et al. (2019). Thus, the complete picture of the surface states in strained HgTe should differ significantly from the predictions based on the 3D BHZ model Zhang et al. (2009). The picture of the surface states becomes even more complex in the case of compressive biaxial strain, when the $|\Gamma_{8},\pm 1/2\rangle$ and $|\Gamma_{8},\pm 3/2\rangle$ bands touch at certain points of the Brillouin zone Mahler et al. (2019), which is particularly similar to unstrained Cd3As2 crystals Bodnar ; Akrap et al. (2016); Desrat et al. (2018) known to be 3D Dirac semimetals.

In order to gain better insight into this particular query we investigate analytically the transformation of Dirac-like surface states induced by the hybridization between the inverted $|\Gamma_{8},\pm 1/2\rangle$ and $|\Gamma_{6},\pm 1/2\rangle$ bands with an additional heavy-hole $|\Gamma_{8},\pm 3/2\rangle$ band. To include an additional band, we have combined the 3D BHZ Hamiltonian Zhang et al. (2009) with the 2D pseudospin-1 Dirac-Weyl Hamiltonian (cf. Refs Malcolm and Nicol (2015); Raoux et al. (2014)) by introducing an effective hybridization strength. This allows us to explore the evolution of topological surface states at different position of the ”parasitic” band by varying the hybridization strength. At a specific value of hybridization strength, our linear model qualitatively represents the picture of the surface states in the vicinity of the $\Gamma$ point known from the tight-binding calculations for Cd3As2 Wang et al. (2013) and HgTe Chu et al. (2011). The results demonstrate that the parabolic Dyakonov-Khaetskii surface states known from HgTe Dyakonov and Khaetskii (1981) stem from the modification of the Dirac-like surface states by a hybridization with an additional $|\Gamma_{8},\pm 3/2\rangle$ band in the bulk.

The paper is structured as follows. We introduce a general analytical model based on a combination of the 3D BHZ and 2D pseudospin-1 Dirac-Weyl Hamiltonian in Section II. Subsequently we discuss the topological surface states at different strengths of the hybridization with the flat band as well as for different boundaries at different position of the flat band, including the 3D TI and Dirac semimetal cases. Finally, the main results are summarized in Section IV.

II Theoretical model

For the analytical model we first consider a ”modified” 6-band Hamiltonian including a variable hybridization with a ”parasitic” bulk band:

[TABLE]

Here the asterisk stands for complex conjugation and ”${\dagger}$” corresponds to Hermitian conjugation. The elements $\hat{H}_{0}(\hat{k}_{x},\hat{k}_{y},\hat{k}_{z})$ and $\hat{H}_{z}(\hat{k}_{z})$ in Eq. (1) are written as

[TABLE]

and

[TABLE]

where $\hat{k}_{\pm}=\hat{k}_{x}\pm i\hat{k}_{y}$ with $\hat{k}_{x}$, $\hat{k}_{y}$, $\hat{k}_{z}$ being momentum operators. We note that $C_{0}$ corresponds to a set of zero energies, $S_{0}$ describes the position of the ”parasitic” band and $v_{\|}$ as well as $v_{\perp}$ are the values of velocity for the massless particles. In our case, the $x$, $y$ and $z$ axes are oriented along the (100), (010) and (001) crystallographic directions. For simplicity, we further assume $v_{\|}=v_{\perp}=v$, which can be found in HgCdTe crystals Orlita et al. (2014); Teppe et al. (2016). The mass parameter $M$ describes the inversion of bands with opposite parities at which $M_{0}>0$ correspond to the normal band ordering and $M_{0}<0$ to an inverted one Bernevig et al. (2006); Zhang et al. (2009).

An important quantity of $\hat{H}_{0}$ is the parameter $\alpha$, which describes the hybridization of the topological surface states with the ”parasitic” bulk flat band. An exact $\alpha$ value for a given system can by obtained by k$\cdot$p perturbation theory up to linear-in-$k$ order developed in the vicinity of critical points of the Brillouin zone considering all point group symmetries of the bulk crystal. For instance, the Hamiltonian in Eq. (1) at $\alpha=\pi/3$ is essentially the 6-band Kane Hamiltonian regarding the $\Gamma_{6}$ and $\Gamma_{8}$ bands, which describes the band structure in the vicinity of the $\Gamma$ point of zinc-blende crystals (see Appendix). This means that, depending on $\alpha$, the Hamiltonian in Eq. (1) interpolates between 3D BHZ Hamiltonian with decoupled flat band at $\alpha=0$ for Bi2Se3-class materials Zhang et al. (2009) and the Kane Hamiltonian at $\alpha=\pi/3$.

Additionally, one can see that $\hat{H}$ in Eq. (1) with $M_{0}=S_{0}=0$, $\hat{k}_{z}=0$ and $\alpha=\pi/4$ corresponds, up to a simple unitary transformation, to the 2D pseudospin-1 Dirac-Weyl Hamiltonian for massless fermions (cf. Refs Malcolm and Nicol (2015); Raoux et al. (2014)). We note that $\hat{H}(\alpha)$, $\hat{H}(-\alpha)$ and $\hat{H}(\pi/2\pm\alpha)$ are all related by unitary transformation. Also there are no quadratic terms considered in $\hat{H}_{z}$ and $\hat{H}_{0}$ for Eqs. (2) and (3) as their form strongly depends on the crystalline symmetry and may differ for two crystals with different values of $\alpha$. Therefore, the universality of the model cannot be preserved beyond the linear approximation but it is in good agreement with magnetooptical experiments for real crystals with $\alpha=0$ Assaf et al. (2016); Krizman et al. (2018) and $\alpha=\pi/3$ Orlita et al. (2014); Teppe et al. (2016); Akrap et al. (2016); Desrat et al. (2018).

Supplementary one can see that $\hat{H}$ in Eq. (1) is invariant under inversion symmetry, but a real crystal may not necessarily feature an inversion center in the unit cell, which can results in additional terms in the Hamiltonian. An explicit form of these terms also depends on the crystalline symmetry and may differ for two crystals with the same value of $\alpha$. For instance, breaking the inversion symmetry in compressively strained HgTe and unstrained Cd3As2, both represented by $\alpha=\pi/3$, results in a transition from Dirac- into Weyl-semimetal in HgTe Ruan et al. (2016) while for Cd3As2 it retains a fourfold degenerate Dirac node Wang et al. (2013). In both crystals, the strength of these terms extracted from experimental data is small Orlita et al. (2014); Teppe et al. (2016); Akrap et al. (2016); Desrat et al. (2018), and can be neglected. Having said that, we will retain the Hamiltonian in a general form for the reasons of universality and explicitly consider the inversion symmetrical case.

Under these assumptions the Hamiltonian in Eq. (1) has three eigenvalues, each double degenerate due to the time-reversal symmetry. The eigenvalues $E$ follow the equation:

[TABLE]

where $k_{x}$, $k_{y}$, $k_{z}$ are the quantum numbers of the momentum operators, $E_{c}(E)=C_{0}+M_{0}-E$, $E_{hh}(E)=C_{0}+S_{0}-E$ and $E_{lh}(E)=C_{0}-M_{0}-E$.

III Topological surface states

For the dispersion of the topological surface states, we further consider an interface between two semiconductors with conventional ($M_{0}^{(I)}>0$, CdTe) and inverted ($M_{0}^{(II)}<0$, HgTe) band structure. Since the position of all the bands can be different on both sides of the interface it is justified to make the parameters $C_{0}$ and $S_{0}$, in addition to $M_{0}$, dependent on the coordinates. For coordinates far away from the interface the values of these parameters naturally tend to the values inherent to the bulk materials. However, the concrete form of $M_{0}(x,y,z)$, $C_{0}(x,y,z)$ and $S_{0}(x,y,z)$ depends on the smoothness of the junction and the crystallographic orientation.

We will consider two different cases, corresponding to the abrupt junction oriented along different crystallographic directions. As for the other parameters, we consider $v$ and $\alpha$ to be independent of coordinates with the same values from both side of the junction, like it is sufficient for the boundary between CdTe and HgTe Teppe et al. (2016). Under these assumptions it can be seen from Eq. (1), that the Hamiltonian remains Hermitian even in the presence of the junction, presuming the operators $\hat{k}_{x}$, $\hat{k}_{y}$ and $\hat{k}_{z}$ do not commute.

The initial Schrödinger equation with the $6\times 6$ Hamiltonian $\hat{H}$ can be considered as a set of differential equations for the $(\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4},\varphi_{5},\varphi_{6})^{T}$ envelope functions, resulting in an 6$\times$6 differential matrix. First, we note that each pair of $(\varphi_{1},\varphi_{4})$, $(\varphi_{2},\varphi_{5})$ and $(\varphi_{3},\varphi_{6})$ represents the electron states with opposite spin orientation in the given band. Namely, for $\alpha=\pi/3$ in HgTe-class materials, they correspond to the $|\Gamma_{6},\pm 1/2\rangle$, $|\Gamma_{8},\pm 3/2\rangle$ and $|\Gamma_{8},\pm 1/2\rangle$ band, respectively (see Appendix). Second, the absence of k-dependent terms in the diagonal elements allows one to express four of the six envelope functions by the two other functions for the same band. Such procedure, which is known as Gaussian elimination, is often used for multi-band Hamiltonians Pfeffer and Zawadzki (2003); Bernardes et al. (2007); Gavrilenko et al. (2011); Krishtopenko et al. .

Considering the pair of $(\varphi_{1},\varphi_{4})$ envelope functions we can express the four functions $\varphi_{2}$, $\varphi_{3}$, $\varphi_{5}$, $\varphi_{6}$ in terms of $\varphi_{1}$ and $\varphi_{4}$ by keeping the right order of non-commuting operators. Then, substituting the expressions for $\varphi_{2}$, $\varphi_{3}$, $\varphi_{5}$, $\varphi_{6}$ into the other two equations, we obtain the $2\times 2$ energy-dependent Hamiltonian $\hat{H}_{E}$, which describes the evolution of the vector $\mathbf{\Phi}=(\Phi_{1},\Phi_{2})^{T}=(\varphi_{1},\varphi_{4})^{T}$ for two spin states:

[TABLE]

where $E$ is the eigenvalue, and

[TABLE]

[TABLE]

[TABLE]

Here, we kept the previous notations, and $E_{c}$, $E_{hh}$ and $E_{lh}$ are the same as for Eq. (4). This changes in the following where they also become dependent on the coordinates:

[TABLE]

We note that the eigenvalue problem in Eq. (5) is very similar to conventional Schrödinger equation with the single-band Hamiltonian including non-diagonal spin-orbit interaction Pfeffer and Zawadzki (2003); Bernardes et al. (2007); Gavrilenko et al. (2011). Moreover, the procedure described above can be also performed for the other pairs $(\varphi_{2},\varphi_{5})$ and $(\varphi_{3},\varphi_{6})$ of the envelope functions.

Although the $2\times 2$ energy-dependent Hamiltonian $\hat{H}_{E}$ in Eq. (5) is non-Hermitian, its eigenvalues are real and are the same as those of $\hat{H}$ in Eq. (1). Recent progress achieved over the last two decades proves that a consistent quantum mechanics can be also built on non-Hermitian Hamiltonians with $\mathcal{PT}$ symmetry Bender and Boettcher (1998); Bender et al. (2002, 2007), where $\mathcal{P}$ and $\mathcal{T}$ are parity and time-reversal operators, respectively. One can see that $\hat{H}_{E}$ in Eq. (5) can be presented in the form:

[TABLE]

where $\sigma_{0}$ is a 2$\times$2 unity matrix, $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$ are Pauli matrices, $d_{0}(\mathbf{k})$, $d_{x}(\mathbf{k})$, $d_{y}(\mathbf{k})$, $d_{z}(\mathbf{k})$ are the real functions of k. The Hamiltonian $\hat{H}_{\mathrm{2x2}}(\mathbf{k})$ always has the real eigenvalues guaranteed by $\mathcal{PT}$ symmetry (where $\mathcal{PT}=\sigma_{x}\mathcal{K}$ and $\mathcal{K}$ is the complex conjugation operator) Bender and Boettcher (1998); Bender et al. (2002). In general, non-Hermiticity arises from the presence of energy or particle exchanges with its environment Konotop et al. (2016); El-Ganainy et al. (2018). The presence of heterojunction in the crystal induces additional interaction between all the bands from the opposite sides of the boundary, which results in the non-Hermiticity of $\hat{H}_{E}$ describing the pair $(\varphi_{1},\varphi_{4})$.

III.1 Surface states for the boundary parallel to (001) crystallographic plane

Let us now consider an abrupt semi-infinite boundary parallel to the (001) crystallographic plane, placed at $z=0$. In this case, $M_{0}(x,y,z)$, $C_{0}(x,y,z)$ and $S_{0}(x,y,z)$ only depend on $z$:

[TABLE]

where $\theta(z)$ is a step-like function defined as $\theta(z)=0$ at $z<0$ and $\theta(z)=1$ at $z\geq 0$. As mentioned above, the region I corresponds to CdTe with $M_{0}^{(I)}>0$, while the region II represents HgTe with $M_{0}^{(II)}<0$. We demonstrate in the following, that the step-like interface allows us to calculate the dispersion of surface states analytically at arbitrary values of $\alpha$. However, the analytical solution can be found at $\alpha=0$ for several smooth interfaces as well Tchoumakov et al. (2017).

The step-like form of $M_{0}(z)$, $C_{0}(z)$ and $S_{0}(z)$ in Eq. (III.1) does not only preserve translation symmetry along the $x$ and $y$ directions but also facilitates the reduction of the eigenvalue problem in Eq. (5) to a set of homogeneous differential equations for the regions $z<0$ and $z>0$. Then the dispersion of the surface states can be found by applying the boundary conditions at $z=0$, which were obtained after integrating Eq. (5) across the small region in the vicinity of $z=0$. Based on the above arguments, $k_{x}$ and $k_{y}$ are the good quantum numbers and the wave-function of the surface states localized in the vicinity of $z=0$ has the form:

[TABLE]

where $\lambda_{z}^{(I)}$ and $\lambda_{z}^{(II)}$ are written as:

[TABLE]

Here, $n=I,II$ and $k^{2}=k_{x}^{2}+k_{y}^{2}$. Note that Eq. (11) can be also derived from Eq. (4) by formally substituting $-k_{z}^{2}\rightarrow\lambda_{z}^{2}$. One can show that the following functions should be continuous across the junction:

[TABLE]

Applying the boundary conditions to $\Phi_{1,2}^{(I)}$ and $\Phi_{1,2}^{(II)}$, the secular equation for the non-trivial solution leads to

[TABLE]

Equations (11) and (13) give the energy dispersion relations of the surface states. It is seen that the surface states at $k=0$ exist if $E_{lh}^{(I)}$ and $E_{lh}^{(II)}$ are of different sign, which requires different signs of $M_{0}^{(I)}$ and $M_{0}^{(II)}$. Substitution of Eq. (11) into Eq. (13) results in a biquadratic equation, which can be solved analytically.

For simplicity, we set $C_{0}^{(n)}=0$ and analyze the case $S_{0}^{(n)}=-M_{0}^{(n)}$. The latter for instance corresponds to the Kane fermions in unstrained HgCdTe crystals Orlita et al. (2014); Teppe et al. (2016). In this case, Eq. (4) provides an energy dispersion for the bulk states, which is independent of $\alpha$:

[TABLE]

while Eq. (13) for the surface states is reduced to

[TABLE]

where

[TABLE]

One can see that for $\alpha=0$, i.e in the absence of hybridization with a flat band, in addition to $E=-M_{0}^{(I)}$ and $E=-M_{0}^{(II)}$, Eqs. (15) and (16) also give $E=\pm\hbar vk$. The latter coincides with the results obtained within 3D BHZ model with the open boundary conditions Shan et al. (2010).

The presence of hybridization splits initial Dirac-like surface states $E=\pm\hbar vk$ into several branches below and above the flat bands in the materials at both sides of the boundary. Particularly, the surface states lying between the flat bands $E=-M_{0}^{(I)}$ and $E=-M_{0}^{(II)}$ are pushed away from the edges of the flat bands at $\alpha\neq 0$. At large values of $k$, their asymptotic energies are written as

[TABLE]

which can be found from Eqs (15) and (16) with $A=0$. We note that $E_{\mathrm{as\pm}}$ is independent of $k$ only in the linear approximation of $\hat{H}$ in Eq. (1), used for the analytical investigation of the surface states at different values of $\alpha$. Including the quadratic terms, whose explicit form depend on $\alpha$, results in non-zero curvature for $E_{\mathrm{as\pm}}$. The latter has been shown by numerical tight-binding calculations on cubic lattices for the case of $\alpha=\pi/3$, corresponding to real HgTe crystals Chu et al. (2011).

Figure 1 shows the dispersion of the bulk and surface states for different values of $\alpha$ in the range $(0,\pi/2]$ for an unstrained film with parameters of HgTe ($\hbar v=850$ meV$\cdot$nm Teppe et al. (2016) and $M_{0}^{(II)}=-150$ meV) sandwiched between CdTe barriers ($M_{0}^{(I)}=450$ meV). We note that the energy range $-M_{0}^{(I)}<E<M_{0}^{(I)}$ corresponds to the band gap in the barriers. Although the bulk dispersion in the film remains the same for any values of $\alpha$ (due to $S_{0}^{(n)}=-M_{0}^{(n)}$, see Eq. (14)), the dispersion of the surface states strongly depends on hybridization with the flat bands in both materials. At small values of $\alpha$ (see Fig. 1(a)), it consists of four branches $E=E_{\mathrm{as+}}$, $E=\pm\hbar vk$ and $E=E_{\mathrm{as-}}$ anticrossing in the vicinity of the crossing points. In this case, the values of $E_{\mathrm{as+}}$ and $E_{\mathrm{as-}}$ are very close to $-M_{0}^{(II)}$ and $-M_{0}^{(I)}$, respectively, since $\cos\alpha\approx 1$ in Eq. (17). This picture can be also treated within the conventional degenerate perturbation approach.

With increasing of $\alpha$, the surface states for $-M_{0}^{(I)}<E<-M_{0}^{(II)}$ are pushed away from the energies of the flat bulk bands in the materials at both sides of the boundary toward the regions where these bulk states are absent. This is clearly represented by the evolution of the asymptotic energies $E_{\mathrm{as+}}$ and $E_{\mathrm{as-}}$, which are getting closer to each other when $\alpha$ increases as shown in Fig. 1(a-c). Note that dispersion of the surface states remains linear in the vicinity of the $\Gamma$ point of the Brillouin zone. At the specific value of $\alpha=\pi/2$, the asymptotic energies coincide both being equal to $-(M_{0}^{(I)}+M_{0}^{(II)})/2$, and the surface states become degenerate, see Fig. 1(d). The latter means the absence of odd-in-$k$ terms in their dispersion. We note that the surface states, similar to those provided in Fig. 1(b) for $\alpha=\pi/3$, were also obtained by more sophisticated numerical calculations based on tight-binding extension of the 6-band Kane Hamiltonian with square terms Chu et al. (2011). Although, the results of Ref. Chu et al. (2011) depend on the constant of artificial cubic lattice used in the calculations, they qualitatively reproduce the dispersion of the surface states at small quasimomentum obtained from our analytical model at $\alpha=\pi/3$.

In addition to the modification of the surface states in the range of $-M_{0}^{(I)}<E<-M_{0}^{(II)}$, the hybridization with the flat bulk bands also yields new ”massive” branches in the regions above and below the flat bands. We refer to the upper ”massive” surface states above the flat band as the Dyakonov-Khaetskii (DK) branch. Dyakonov and Khaetskii Dyakonov and Khaetskii (1981) were the first, who predicted the massive states at the surface of HgTe crystal. They derived analytically this branch in 1981 by using Luttinger Hamiltonian for the $\Gamma_{8}$ bands Luttinger (1956) with an open boundary conditions. In 1985, existence of the localized states at the HgTe/CdTe interface was also predicted for the quantum wells Lin-Liu and Sham (1985) and superlattices Chang et al. (1985).

Although the Luttinger Hamiltonian, used in Refs Dyakonov and Khaetskii (1981); Lin-Liu and Sham (1985), does not formally consider the inverted $|\Gamma_{6},\pm 1/2\rangle$ band, this Hamiltonian can be obtained from the 6-band Kane Hamiltonian with the HgTe/CdTe interface by assuming $M_{0}^{(II)}\rightarrow-\infty$ and $M_{0}^{(I)}\rightarrow\infty$. Therefore, the upper ”massive” surface states in Fig. 1 and solution of Dyakonov and Khaetskii Dyakonov and Khaetskii (1981), obtained for particular case of $\alpha=\pi/3$, have the same origin. As seen from Fig. 1, such DK branch is caused by the band inversion in the presence of hybridization with the flat bulk band.

Now we consider a bulk crystal, in which the flat band does not coincide with the bottom of the conduction band, i.e. $S_{0}^{(II)}\neq-M_{0}^{(II)}$. The case of $S_{0}^{(II)}<-M_{0}^{(II)}$ corresponds to an external tensile biaxial strain, which opens a band gap, yielding 3D TI state Fu and Kane (2007); Brüne et al. (2011). The opposite case of $S_{0}^{(II)}>-M_{0}^{(II)}$ is realized in the compressively strained HgTe films Mahler et al. (2019) or in unstrained Cd3As2 crystals Akrap et al. (2016); Desrat et al. (2018). As in previous case, we set $C_{0}^{(n)}=0$ and assume $S_{0}^{(I)}=-M_{0}^{(I)}$ in CdTe layer.

Fig. 2 shows a picture of the bulk and surface states at different values of $\alpha$ in the $(0,\pi/2]$ range for tensile strained film the with parameters of HgTe ($S_{0}^{(II)}=110$ meV) sandwiched between CdTe barriers. Note that now the energy dispersions are calculated numerically on the basis of Eq. (4) and Eqs. (11), (13) for the bulk and surface states, respectively. It is seen that energy of the bulk states and the value of a band-gap between the flat and conduction bands strongly depend on $\alpha$. The maximum gap is achieved in the absence of hybridization, while increasing of $\alpha$ leads to a band-gap vanishing. The value of $\alpha=\pi/2$ corresponds to a semimetal with circular nodal line at $k_{z}=0$ and $k=k_{N}$, where $k_{N}^{2}=2M_{0}^{(II)}(M_{0}^{(II)}+S_{0}^{(II)})/(\hbar^{2}v^{2})$.

As seen from Fig. 2, the surface states in a tensile strained film at different $\alpha$ values remain qualitatively the same, as in Fig. 1 for the unstrained film. The main difference is seen in the DK branch, which now exists in the band-gap for the bulk states for all $\alpha$. This is consistent with the general topological arguments claiming that tensile strained HgTe is a 3D TI with gapless surface states Fu and Kane (2007). However, these surface states can not be represented by massless Dirac fermions, as it is stated in some experimental works on HgTe strained by a CdTe substrate (see, for instance, Refs. Kozlov et al. (2016); Thomas et al. (2017); Noel et al. (2018)). Fig. 2 clearly shows that the surface states at the HgTe/CdTe boundary of strained HgTe-based 3D TI are ”massive” due to the hybridization with heavy-hole $|\Gamma_{8},\pm 3/2\rangle$ band and represented by DK branch Dyakonov and Khaetskii (1981).

In the opposite case of $S_{0}^{(II)}>-M_{0}^{(II)}$, the flat band crosses the conduction band at certain points of the Brillouin zone at $\alpha\neq 0$ yielding a 3D Dirac semimetal. At these points, the conduction and flat valence bands can be considered as two highly anisotropic and tilted cones Wang et al. (2013); Akrap et al. (2016); Desrat et al. (2018), whose nodes lie at $k_{x}=k_{y}=0$ and $k_{z}=\pm k_{D}$ with $k_{D}^{2}=(S_{0}^{2}-M_{0}^{2})/(\hbar^{2}v^{2})$, see Eq. (4). Note that at $\alpha=0$, the crossing points are located at the sphere defined by $\hbar^{2}v^{2}(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})=(S_{0}^{2}-M_{0}^{2})$, and the 3D Dirac semimetal state does not arise.

Fig. 3 presents the bulk and surface states in a film with $S_{0}^{(II)}>-M_{0}^{(II)}$ (where $S_{0}^{(II)}=190$ meV) at different strengths of hybridization with a flat band. As in the previous case of $S_{0}^{(II)}<-M_{0}^{(II)}$, the bulk and surface states both depend on the values of $\alpha$. Interestingly, the dispersion of the surface states in Fig. 3 for all $\alpha\neq 0$ starts from projection of the bulk Dirac nodes. The bulk band dispersion as a function of $k$ at $k_{z}=k_{D}$ is represented by the blue curves. We note that the particular case $\alpha=\pi/3$ is in a good qualitative agreement with the picture of the surface states obtained from the tight-binding calculations for Cd3As2 on a tetragonal lattice Wang et al. (2013) (see Fig. 3(a,b) therein).

III.2 Surface states for the boundary parallel to (010) plane

One of the inherent characteristics of the surface states in Dirac semimetals is the existence of a pair of surface Fermi arcs connecting bulk Dirac nodes projected on the surface boundary. The arcs meet at a sharp corner or ”kink” at the projected nodes. Such a kink is not allowed in a purely 2D metal, it is a special feature of the crystal symmetry-protected Weyl structure of the Dirac semimetals Potter et al. (2014). As the Dirac nodes are located along the (001) crystallographic direction, the surface boundary parallel to the (001) plane has no Fermi arcs.

Let us now briefly consider the surface states for the boundary containing two projections of the bulk Dirac nodes at $S_{0}^{(II)}>-M_{0}^{(II)}$. For the boundary plane parallel to the x-z plane and placed at $y=0$, $M_{0}$, $C_{0}$ and $S_{0}$ have a step-like dependence on $y$:

[TABLE]

Basing on the arguments similar to the case of (001) interface, the wave-function of the surface states for the (010) boundary has the form:

[TABLE]

where $\lambda_{y}^{(I)}$ and $\lambda_{y}^{(II)}$ are written as:

[TABLE]

where the index $n=I,II$ corresponds to $y<0$ and $y>0$, respectively.

Integration of Eq. (5) across the small region of $y=0$ gives the continuity function across the junction:

[TABLE]

where $R_{+}$ and $R_{-}$ are function of $E$ and $\alpha$:

[TABLE]

Applying these boundary conditions to $\Phi_{1,2}^{(I)}$ and $\Phi_{1,2}^{(II)}$, the secular equation for the non-trivial solution leads to

[TABLE]

Substitution of $\lambda_{y}^{(I)}$ and $\lambda_{y}^{(II)}$ into Eq. (13) also results in a quadratic equation for $k_{z}^{2}$ and $k_{x}^{2}$, which can be solved analytically. Such a quadratic equation gives the energy dispersion and energy contours of the surface states for the boundary parallel to the (010) crystallographic plane.

Fig. 4 shows the dispersions of the bulk and surface states as a function of $k_{z}$ for the bulk film with $S_{0}^{(II)}>-M_{0}^{(II)}$ and (010) surface boundary. Here, we set $C_{0}^{(n)}=0$, $S_{0}^{(I)}=-M_{0}^{(I)}$, $S_{0}^{(II)}=190$ meV and assume $k_{x}=0$. As for the (001) boundary, the picture of surface states for all $\alpha$ values consists of two branches above and below the bottom of the conduction band at $E=-M_{0}^{(II)}$. As it is seen, the upper DK branch for all $\alpha\neq 0$ crosses the bulk dispersion precisely at the Dirac nodes. This stems from the fact that two separated Dirac nodes are connected by the topological surface states Potter et al. (2014). A current picture of the surface states for $\alpha=\pi/3$ is also consistent with the tight-binding calculations for Cd3As2 Wang et al. (2013).

Fig. 5 provides energy contours for the surface states at $E=S_{0}^{(II)}$ at different strengths of hybridization with the flat band. In contrast to the (001) boundary, for which the energy contour of the surface states at the Dirac nodes reduces to a point, the nontrivial surface states at the (010) boundary are clearly visible. Its Fermi surface at $E=S_{0}^{(II)}$ is composed of two Fermi arcs with the kinks at the projected bulk Dirac nodes. As seen from Fig. 5, the length of Fermi arcs depends on the values of $\alpha$. This means that the period of quantum oscillations originated from cyclotron orbits weaving together Fermi arcs and chiral bulk states Potter et al. (2014), should also depend on the hybridization strength.

IV Summary

In conclusion, we have performed an analytical study of the hybridization between topological surface states and the non-topological flat band in the bulk. It was shown that the hybridization with the flat band divides the initially Dirac-like surface states, derived from the 3D BHZ model, into two branches, one below the flat band and another one above the edge of conduction band. The upper branch at $\alpha=\pi/3$ is formed by Dyakonov-Khaetskii surface states Dyakonov and Khaetskii (1981) known for HgTe since the 1980s. Adjusting the hybridization strength, we have explored the evolution of topological surface states in 3D TIs and 3D Dirac semimetals arising at different positions of the flat band. Our results show that the surface states lying inside the band gap of 3D TIs, as well as the Fermi arcs of 3D Dirac semimetals in HgTe and Cd3As2 are represented by the DK branch of the surface states. Although, we have applied the linear approximation for the bulk Hamiltonian, our model qualitatively represents the picture of surface states at small values of k, known for $\alpha=\pi/3$ from numerical tight-binding calculations on tetragonal Wang et al. (2013) and cubic lattices Chu et al. (2011). This work paves the way for further analytical investigations of different characteristics of the surface states hybridized with non-topological bands.

Acknowledgements.

The authors gratefully thank S. Gebert (Institut d’Electronique et des Systemes, Montpellier) for the helpful discussions and critical comments. This work was supported by MIPS department of Montpellier University through the ”Occitanie Terahertz Platform”, by the CNRS through LIA ”TeraMIR” and the French Agence Nationale pour la Recherche (Dirac3D project).

Appendix: 6-band Kane Hamiltonian

In order to demonstrate that the Hamiltonian in Eq. (1) for $\alpha=\pi/3$ corresponds to the 6-band Kane Hamiltonian, we consider the 8-band Kane Hamiltonian Bodnar , whose form is dependent on the choice of the basis set of the Bloch amplitudes for the $\Gamma_{6}$, $\Gamma_{8}$ and $\Gamma_{7}$ bands. In the given basis set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

the 8-band Kane Hamiltonian $\hat{H}_{\textbf{k}}^{8\times 8}$ in the presence of only the linear terms takes the form

[TABLE]

Here, $P$ is the Kane momentum matrix element, $\Delta$ is the spin orbit energy and $E_{c}$ as well as $E_{v}$ are the conduction and valence band edges, respectively. We note that the values of $E_{v}$, $E_{v}^{\prime}$ and $E_{v}^{\prime\prime}$ differ if there is biaxial strain in the (001) crystallographic plane Krishtopenko et al. (2016).

In the limit of large $\Delta$, the Hamiltonian $\hat{H}_{\textbf{k}}^{8\times 8}$ can be easily projected on the subspace, orthogonal to the the split-off $\Gamma_{7}$ band. As we are not interested in terms quadratic in $\mathbf{k}$,the projection is done by simply eliminating the fourth and the eight row and column of the matrix in Eq. (24):

[TABLE]

As it is seen, Eq. (25) coincides with Eq. (1) at $\alpha=\pi/3$ and $\hbar v_{\|}=\hbar v_{\perp}=\sqrt{6}P/3$.