Surface States of a System of Dirac Fermions: A Minimal Model

TL;DR

This paper presents a minimal analytical model for surface and edge states in Dirac materials, applicable to both topological and non-topological systems, based on envelope wave functions and boundary conditions.

Contribution

It introduces a minimal formalism that describes surface and edge states of Dirac fermions across different systems, expanding theoretical understanding.

Findings

Analytical description of surface and edge states in Dirac systems.

Applicability conditions for the minimal model are discussed.

The model covers both topological and non-topological Dirac materials.

Abstract

A brief survey is given of theoretical works on surface states (SSs) in Dirac materials. Within the formalism of envelope wave functions and boundary conditions for these functions, a minimal model is formulated that analytically describes surface and edge states of various (topological and non-topological) types in several systems with Dirac fermions (DFs). The applicability conditions of this model are discussed.

Equations

$$E(p)=\pm \sqrt{m^2{c}^4+p^2{c}^2}$$

$$\{\left[E_n(0)+\frac{{\hat{\mathbf{p}}}^2}{2m_0}+V(z)\right]{\delta }_{n{n}^{\prime }}+\frac{\hat{\mathbf{p}}\cdot {\mathbf{p}}_{n{n}^{\prime }}}{m_0}$$

$$+\frac{\hslash }{4m_0^2{c}^2}{\left(\mathbf{p}\cdot \left[\mathbf{\sigma }\times \nabla V_0\right]\right)}_{n{n}^{\prime }}\}{\Phi }_{n^{\prime }}=E{\Phi }_n$$

$${\left({\Phi }_{\lambda }^{+}{\hat{v}}_z{\Phi }_{\nu }\right)|}_{z=0}=0$$

$$\hat{H}=\frac{{\hat{\mathbf{p}}}^2}{2m^{\ast }}+V(z)+{\hat{H}}_{BIA}+{\hat{H}}_{SIA},$$

$${\hat{H}}_{BIA}=\frac{{\gamma }_c}{{\hslash }^3}[{\sigma }_x{p}_x\left({\hat{p}}_y^2-{\hat{p}}_z^2\right)$$

$$+{\sigma }_y{p}_y\left({\hat{p}}_z^2-{\hat{p}}_x^2\right)+{\sigma }_z{p}_z\left({\hat{p}}_x^2-{\hat{p}}_y^2\right),$$

$${\hat{H}}_{SIA}=a_{SO}\left({\sigma }_x{p}_y-{\sigma }_y{p}_x\right){\partial }_{z}V(z).$$

$$\hat{T}=i{\sigma }_y\hat{K},$$

$$[{\sigma }_0-i\frac{R{\hat{p}}_z}{\hslash }-i\frac{2m^{\ast }{\gamma }_{c}R}{{\hslash }^4}\left({\sigma }_y{p}_y-{\sigma }_x{p}_x\right){\hat{p}}_z$$

$$-i\frac{m^{\ast }{\gamma }_{c}R}{{\hslash }^4}{\sigma }_z\left(p_x^2-p_y^2\right)+\frac{\left(\chi +{\chi }^{int}\right)R}{\hslash }\left({\sigma }_x{p}_y-{\sigma }_y{p}_x\right)$$

$${-\frac{2m^{\ast }{\gamma }_c^{int}}{{\hslash }^3}\left({\sigma }_y{p}_y-{\sigma }_x{p}_x\right)]\varphi |}_{z=0}=0.$$

$$\left(\begin{array}[]{cc}mc^{2}-E&c\bm{\sigma}\cdot\bm{k}\\ c\bm{\sigma}\cdot\bm{k}&-mc^{2}-E\end{array}\right)\left(\begin{array}[]{cc}\Psi_{c}\\ \Psi_{v}\end{array}\right)=0$$

$${\left[{\Psi }_v-i{a}_0\mathbf{\sigma }\cdot \mathbf{n}{\Psi }_c\right]}_S=0.$$

$$E_s({\mathbf{k}}_{||})=s\frac{2a_0}{1+a_0^2}|{\mathbf{k}}_{||}|+E_0.$$

$$E_0=m{c}^2\frac{1-a_0^2}{1+a_0^2}.$$

$$\frac{2a_0}{1+a_0^2}\frac{mc}{\hslash }-s|{\mathbf{k}}_{||}|\ge 0.$$

$$\left(\begin{array}[]{cc}-E&i\dfrac{E_{g}(z)}{2}+c\bm{\sigma}\cdot\bm{k}\\ -i\dfrac{E_{g}(z)}{2}+c\bm{\sigma}\cdot\bm{k}&-E\end{array}\right)\left(\begin{array}[]{cc}\Psi_{c}\\ \Psi_{v}\end{array}\right)=0$$

$$\Psi_{\pm}=\left(\begin{array}[]{c}0\\ e^{i\theta/2}\\ \pm e^{i\theta/2}\\ 0\end{array}\right)\exp\left\{\int_{0}^{z}\frac{E_{g}(z)}{2\hbar c}dz+i\bm{k}_{||}\bm{r}_{||}\right\},$$

$$E=-\frac{E_{g1}{\phi }_0}{E_{g0}}\pm |{\mathbf{k}}_{||}|\sqrt{1-{\left(\frac{2{\phi }_0}{E_{g0}}\right)}^2}$$

$$a_0=-\mathrm{sign}\left(E_{g0}\right)\sqrt{\frac{E_{g0}-{\phi }_0}{E_{g0}+{\phi }_0}}$$

$$E=\pm \frac{2a_{0}c\hslash }{1+a_0^2}\sqrt{k_z^2+\frac{(j+\Phi /{\Phi }_0{)}^2}{R^2}}+E_0,$$

$$\Psi \propto \exp \left(\frac{1}{\hslash v}{\int }_{z_0}^z\Delta (z)dz\right),$$

$$\left(\begin{array}[]{cccc}\varepsilon(k)+m(k)&A(k_{x}-ik_{y})&0&0\\ A(k_{x}+ik_{y})&\varepsilon(k)-m(k)&0&0\\ 0&0&\varepsilon(k)+m(k)&-A(k_{x}+ik_{y})\\ 0&0&-A(k_{x}-ik_{y})&\varepsilon(k)-m(k)\end{array}\right)\left(\begin{array}[]{c}\Psi_{\uparrow}\\ \Psi_{\downarrow}\end{array}\right)=E\left(\begin{array}[]{c}\Psi_{\uparrow}\\ \Psi_{\downarrow}\end{array}\right)$$

$${\left[F{\partial }_{\mathbf{n}}\Psi +G\Psi \right]}_{\Gamma }=0$$

$$\left(\begin{array}[]{cc}M(\bm{k})-E&v\bm{\sigma}\cdot\bm{k}\\ v\bm{\sigma}\cdot\bm{k}&-M(\bm{k})-E\end{array}\right)\left(\begin{array}[]{cc}\Psi_{c}\\ \Psi_{v}\end{array}\right)=0$$

$${\left[\frac{B}{M_0}{\partial }_{\mathbf{n}}\Psi +Q\Psi \right]}_S=0$$

$$\left(\begin{array}[]{cc}v\bm{\sigma}\cdot\left(\bm{k}-\bm{K}\right)&0\\ 0&v\bm{\sigma}\cdot\left(\bm{k}-\bm{K^{\prime}}\right)\end{array}\right)\left(\begin{array}[]{cc}\Psi_{K}\\ \Psi_{K^{\prime}}\end{array}\right)=E\left(\begin{array}[]{cc}\Psi_{K}\\ \Psi_{K^{\prime}}\end{array}\right)$$

$${\left[{\Psi }_K\sin \beta +i{e}^{i\phi }\left(e^{i\gamma \mathbf{\sigma }\cdot \mathbf{n}}+{\sigma }_z\cos \beta \right){\sigma }_z{\Psi }_{K^{\prime }}\right]}_{\Gamma }=0$$

$$v\mathbf{\sigma }\cdot \mathbf{p}{\Psi }_{K(K^{\prime })}=E{\Psi }_{K(K^{\prime })}$$

$${\left[{\psi }_{1K(K^{\prime })}+i{a}^{\tau }{\psi }_{2K(K^{\prime })}\right]}_{\Gamma }=0$$

$$E_{\tau }(k_{||})=\frac{2a}{1+a^2}\tau v{k}_{||},\tau k_{||}>0$$

$$E=\frac{2a\hslash v}{R_0}\left(|j|-\frac{\tau }{2}\right)-i\frac{2\pi \hslash v}{R_0}\frac{|a|{\left(|a|(|j|-1/2)\right)}^{2|j|-1}}{{\Gamma }^2(|j|-1/2)}.$$

Full text

I THE ENVELOPE-FUNCTION METHOD. INTRODUCTION TO THE HISTORY OF THE PROBLEM

By the middle of the 20-th century, the problem arose of explaining and predicting the electronic properties of semiconductors in external fields. Various versions of the single-band effective-mass method were developed. The Luttinger–Kohn envelope-function (EF) method Luttinger and Kohn, 1955 based on a generalization of the $\bm{kp}$-approach turned out to be very convenient. The formalism of EFs admitted a natural multiband generalization. Such a generalization was made by Keldysh in his theory of deep impurities Keldysh, 1964. In [Keldysh, 1964], Keldysh also noticed that, under certain conditions, EFs in a narrow-band semiconductor obey an effective Dirac equation (a system of four first-order differential equations). It is interesting that, in III–V semiconductors, this situation occurs under reversal of the sign of the strong spin–orbit splitting of the valence band. The idea of inversion of bands in crystals with strong spin–orbit interaction will be addressed below in this paper. In the same year as [Keldysh, 1964], it was shown [Wolff, 1964] that the spectrum of electrons and holes in bismuth near the $L$ points of the Brillouin zone should be described by an anisotropic Dirac Hamiltonian.

According to modern terminology, bismuth can be considered as a Dirac material, and the first Dirac material at that. These materials include graphene, bismuth– antimony alloys, lead chalcogenides, 2D and 3D topological insulators, Dirac semimetals, Weyl semimetals, and many other materials. One-particle excitations in these materials are called massless (for gapless materials) or massive Dirac fermions (DFs).

The energy spectrum $E(p)$ of free DFs in the three-dimensional isotropic case is analogous to the spectrum of a relativistic electron:

[TABLE]

Here $m$ is the electron effective mass, which is usually one or two orders of magnitude less than the mass of a free electron in vacuum, and $c$ is the effective velocity of light. This velocity is two orders of magnitude less than the velocity of light in vacuum; therefore, real Fermi excitations are, in fact, non-relativistic. However, relativistic corrections due to spin-orbit interaction in crystals with close bands may be rather large, which contributes to the formation of DFs. A finite mass in the spectrum of relativistic particles corresponds to a finite width of the forbidden band according to the well-known formula $E_{g}=2mc^{2}$.

The 1960s marked the emergence of the physics of 2D electron systems. Quantum size effect was started to be used to realize the 3D $\to$ 2D transition. The conventional theory of size quantization is based on the single-band effective mass approximation with zero boundary conditions for EFs. Physically, this is justified by the impenetrability of potential barriers on the surfaces (interfaces) that confine the motion of an electron. This approximation works well in crystals with parabolic band spectrum, for example, for electrons in silicon. However, the situation is qualitatively changed when one tries to correctly describe the size quantization in a system of DFs.

The possibility of implementation of the size quantization was predicted in 1962 [Sandomirskii, 1962] for films made of narrow-band semiconductors and related semimetals with low concentration of carriers. The experimental discovery of this phenomenon [Ogrin et al., 1966, Lutskii et al., 1966] and its explanation [Sandomirskii, 1967] occurred a few years later when studying bismuth films.

Strictly speaking, envelope wave functions strongly differ from real wave functions of an electron in a crystal: the former functions are the envelopes of the latter. However, this fact was not of fundamental importance for qualitative considerations at the initial stage. But the development of experiment raised the question of quantitative description of the results. Unexpectedly, the problem turned out to be rather difficult. A seemingly academic question of the simplest boundary conditions for effective wave functions on the surface of a bismuth-type narrow-band semiconductor inevitably entailed the problem of description of SSs. It turned out that the simplest and physically obvious boundary conditions—zero boundary conditions for all four EFs in the Dirac equation—lead to a trivial (everywhere zero) solution due to the overdeterminacy of the problem. This is how the theoretical problem of derivation of correct boundary conditions for the envelope-function method arose, which is especially important in its multiband version. It is important that these boundary conditions should also describe the spectrum of SSs in addition to volume states.

The point is that, on the real surface of any crystal, including a Dirac crystal, there always exist electron or hole SSs, which are extrinsic and/or intrinsic. The first are attributed to defects and contamination of the surface, while the second (we will consider only these SSs) exist on the ideal surface. Intrinsic SSs were theoretically considered at the early stage of band theory in different models [Tamm, 1932a, b; Maue, 1935; Shockley, 1939; Davison and Levine, 1970]. I.E. Tamm’s works were the first ones; therefore, these SSs are called Tamm states. W. Shockley considered another popular model, and the SSs obtained by him are often called Shockley states. However, there is no essential difference between the Tamm and Shockley SSs. Moreover, in the simplest models, these states turn into each other under the variation of model parameters [Davison and Levine, 1970], which may lead to confusion. Therefore, we will call intrinsic SSs Tamm–Shockley states, bearing in mind that both types of SSs are attributed to the sharp (at the atomic scale) discontinuity of the crystal potential at the surface.

The Tamm–Shockley states are not only situated (completely or partly) in the forbidden band of a bulk crystal. It is very important that these states should form a surface band of conducting states that are delocalized in the plane of the surface and are characterized by a 2D dispersion law. However, intrinsic SSs are very sensitive to surface roughness and contamination. Therefore, it is not surprising that a lot of time passed until the Tamm–Shockley states were reliably detected experimentally. As a rule, these states are investigated in ultra-high vacuum usually by local techniques (STM, ARPES, and so on). The problem of existence of an SS band under normal conditions requires special investigation. For instance, it is not quite clear how much stringent conditions should be imposed on the perfection of the surface or the interface, and what does this depend on.

In ordinary semiconductors (with wide forbidden gap and small spin–orbit interaction), a band of Tamm–Shockley states arises far from always, and the wave functions of these states are not usually mixed with the wave functions of size quantized states. The very existence conditions of these states strongly depend on the structure of the surface or the interface at the atomic scale; therefore, in spite of the 80-year history of the problem, the general principles of behaviour of these functions have been investigated rather fragmentarily, especially experimentally. The situation started to be changed in recent years, when topological materials described by a modified Dirac equation appeared. The theory of topological insulators predicts that ”topological” SSs protected from backscattering should exist in these materials for topological reasons. The very fact of existence of these states does not depend on the details of the structure of a surface region.

A more general question arises: are the requirements imposed on the perfection of the surface (interface) relaxed for the manifestation of SSs on the surface of any Dirac crystal, not necessarily a topological insulator? The theoretical answer in the simplest Dirac model [Volkov and Pinsker, 1981] is positive. Nevertheless, as expected for the Tamm–Shockley states, the spectrum of SSs for Dirac materials depends on the properties of the surface, although not in a quite obvious way.

To solve this problem, one had to derive boundary conditions of general form that are invariant under time reversal and describe DFs near an impenetrable wall. The solution of the genuine (rather than a modified, as in the theory of topological insulators later) Dirac equation in a half-space with these boundary conditions showed that SSs should appear on any surface (more precisely, for any parameters of the boundary conditions) and should be extremely strongly split with respect to spin as a result of spin–orbit interaction with the surface (this interaction is often called the Rashba interaction). These SSs are characterized by a conical dispersion law but they occurred not all the values of momenta. The ”strength” of surface spin–orbit interaction is characterized by a real phenomenological parameter $a_{0}$ that appears in the boundary conditions. Depending on the sign of $a_{0}$, there exist two classes of surfaces. For the surfaces of one of the classes, the conical (”Dirac”) point in the spectrum is located in the forbidden gap of the bulk material; i.e., SSs in this case have a 2D spectrum of massless spin-nondegenerate DFs, as in the case of a topological insulator. The physical meaning of the sign of $a_{0}$ remained unclear.

In 1985, in the famous work [Volkov and Pankratov, 1985], the authors considered another Dirac model: the model of inverse heterocontact. It is described by the Dirac equation in which the mass of a DF is varied in space smoothly (at the atomic scale) and symmetrically with respect to the electron–hole (e–h) transformation. The interface is determined by the position of the plane in which the sign of the mass is reversed (the conduction and valence bands are interchanged). It is this case where the band of nondegenerate heterointerface states possessing a conical 2D spectrum of massless DFs is formed. It is remarkable that this conclusion does not depend on the details of the interface potential; one just needs the inversion of the bands. The point is that, under these conditions, the equation for DFs in the inverse contact takes a form typical for supersymmetric quantum mechanics with zero mode corresponding to interface states. This nontrivial result is actively used in modern physics of topological insulators. The Dirac point in this symmetric model is located exactly at the center of the gap.

The reason for some disagreement between the results of [Volkov and Pankratov, 1985] and [Volkov and Pinsker, 1981] and the physical meaning of the sign of $a_{0}$ were elucidated in the survey [Volkov et al., 1995]. A significant e–h asymmetry was introduced into the model of an abrupt heterocontact (which, however, is still smooth at the atomic scale). As this asymmetry increases, the Dirac point was more and more shifted from the center of the gap and disappeared. For a greater discontinuity of bands at the heterojunction, the spectrum of interface states completely corresponds to the spectrum of SSs from [Volkov and Pinsker, 1981]. The comparison of the results yielded a model expression for the boundary parameter $a_{0}$. The sign of this parameter generally correlated to the reversal of the sign of the gap at the interface.

Unfortunately, the incipient state of the manufacturing technology of samples did not allow one to verify these conclusions experimentally 30–20 year ago. The works [Volkov and Pinsker, 1981; Volkov and Pankratov, 1985; Volkov et al., 1995] acquired importance in relation to the discovery of topological insulators, as well as in relation to the recently detected conducting edge states of non-topological type in nanoperfortated graphene [Latyshev et al., 2013, Latyshev et al., 2014].

In the present study, we discuss the problems touched upon above and, within the formalism of EFs, formulate a minimal model that describes SSs of various, topological and non-topological, types in a number of Dirac materials. The paper is organized as follows. In Section 2, we give a brief survey of works on Tamm–Shockley-type SSs within the formalism of EFs. In Section 3, we consider topological SSs and, in Section 4, Tamm–Shockley edge states in graphene. In Section 5, we formulate conclusions.

II TAMM–SHOCKLEY-TYPE SURFACE STATES

The method of effective wave-functions—envelope functions—is widely used to describe the behavior of electrons in multilayer semiconductor structures. The method of EFs can be applied to the description of smooth (at the atomic scale) fields and is inapplicable to the real case of atomically abrupt interfaces. Information on the microscopic structure of the interface can be taken into account by appropriate boundary conditions for EFs.

The problem of boundary conditions in bounded crystals has a long history. Theoretical studies on boundary conditions in semiconductor structures can be classified into two groups. The works of the first, largest, group are devoted to the derivation of ”two-sided” boundary conditions relating EFs to their left and right derivatives at the interface. They contain different approaches to the solution of mathematical problems related, in particular, to possible singular behaviour of EFs on a heterointerface [Ivchenko, 2005; Zawadzki and Pfeffer, 2004; Foreman, 2005; Takhtamirov and Volkov, 1999; Rodina et al., 2002; Takhtamirov and Volkov, 1997; Takhtamirov and Melnik, 2010].

We will focus only on the works of the second group, which are devoted to the derivation of one-sided boundary conditions at the crystal–high barrier interface (in particular, at the crystal–vacuum interface). The problems of this type arise when describing the Tamm–Shockley surface (interface) states.

Under the neglect of the spin–orbit interaction, a microscopic derivation of boundary conditions for EFs at a step-like boundary between semiconductor ($z>0$) and vacuum ($z<0$) was apparently first presented in [Volkov and Pinsker, 1976, Volkov and Pinsker, 1977]. The boundary conditions introduced there contain boundary parameters that are analytically (but in a complicated way) expressed in terms of a complete infinite-band structure of a semiconductor analytically continued to the domain of complex quasimomenta. The numerical determination of these parameters remains an unsolved problem.

III SHALLOW TAMM–SHOCKLEY STATES IN A SINGLE-BAND APPROXIMATION

In a single-band approximation, a boundary condition represents a linear relation between an EF and its normal derivative with a single boundary parameter with the dimensionality of length, which we denote below by $R$. This length characterizes the structure of the interface between semiconductor ($z>0$) and an impenetrable barrier ($z<0$) at the atomic scale and has the meaning of the localization depth of a shallow Tamm–Shockley state when it exists (to this end, the condition $R>0$ should be satisfied). Moreover, the length $R$ depends on the parameters of the bulk band structure.

In [Volkov and Pinsker, 1979], a much simpler derivation of the same boundary condition is presented from the Hermiticity of the effective Hamiltonian for EFs in a half-space bounded by an impenetrable barrier. Within such a phenomenological approach, the parameter $R$ should be determined from experiment. The high-barrier model is applicable when the interface length $R$ is large compared with the penetration depth into the barrier. The effect of the spin–orbit interaction on single-band boundary conditions and spin splitting of shallow Tamm–Shockley states in the conduction band of a semiconductor with inversion symmetry was considered in [Vasko, 1979]. In the model used in [Vasko, 1979], the spin splitting is controlled by the product of the length $R$ and the parameter of bulk spin–orbit interaction. A non-parabolic generalization of boundary conditions [Vasko, 1979] in an asymmetric quantum well with infinite barriers is presented in [Rodina and Alekseev, 2006].

We begin with the analysis of a single-band boundary condition as applied to heterostructures based on III–V semiconductors with the heterointerface orientation (001). In [Devizorova and Volkov, 2013,Devizorova et al., 2014], the authors considered the effect of an atomically abrupt heterointerface on the effective single-band Hamiltonian of a quantum well and the spin splitting in a conduction band of symmetry $\Gamma_{6c}$ in crystals with lack of inversion symmetry. The discontinuity of bands at the heterointerface is assumed to be large, and the heterobarrier, to be impenetrable. The latter is characterized by a certain boundary condition for EFs.

The dynamics of a conduction electron for $z>0$ within the multiband EF method is described by the system of Kohn–Luttinger $\bm{kp}$-equations

[TABLE]

where $n$ is the band number, $E_{n}(0)$ is the extremum energy of the $n$-th band, $\Phi_{n}$ is a set of EFs, $\bm{p}_{nn^{\prime}}$ is a matrix element of the momentum operator $\hat{\bm{p}}$ on Bloch functions at the center of the Brillouin zone, $m_{0}$ is the mass of a free electron, the last term in curly brackets is a matrix element of the operator of bulk spin–orbit interaction on Bloch functions, and $\sigma_{i}$, $i=x,y,z,$ are Pauli matrices.

The Hermiticity of Hamiltonian (III) in the half-space $z>0$ reduces, after integration by parts, to vanishing of the surface contribution. This is equivalent to setting the matrix element of the normal component of the current operator between any pair of states at the boundary to zero:

[TABLE]

where $\hat{v}_{z}$ is a non-diagonal matrix of velocity $\left(\hat{v}_{z}\right)_{nn^{\prime}}=\partial_{p_{z}}\left(H_{nn^{\prime}}\right)$.

Since only the spinor corresponding to the conduction band $\Gamma_{6c}$ is large in the multicomponent function $\Phi$, we apply the unitary transformation $\Phi=e^{S}\phi$ [Winkler, 2003] (with regard to $\bm{kp}$ terms up to the third order inclusively), which reduces the Hamiltonian to a single-band one with effective mass $m*$ and smooth potential $V(z)$. Now, the 3D Hamiltonian of the conduction band contains contributions $\hat{H}_{BIA}$ and $\hat{H}_{SIA}$ that describe the spin splitting due to the lack of inversion symmetry of the crystal and the asymmetry of the well:

[TABLE]

[TABLE]

[TABLE]

By the same transformation, condition (3) reduces to a certain constraint for a two-component EF. Next, we require the invariance of this constraint with respect to the operation of time reversal

[TABLE]

where $\hat{K}$ is the operator of complex conjugation. We obtain T-invariant boundary conditions that take into account the spin–orbit interaction in the bulk and on the symmetry interface $C_{2v}$, as well as the absence of the inversion center in the bulk crystal:

[TABLE]

Here $\sigma_{0}$ is the identity matrix. The real quantity $R$, which has the dimensionality of length, depends on the microscopic structure of the boundary. The constants $\gamma_{c}$ and $\chi$ are defined by the volume parameters; for GaAs, $|\gamma_{c}|=24.4$ eV$\cdot{\rm A}^{3}$ and $\chi=0.082$. The constants $\gamma_{c}^{int}$ and $\chi_{c}^{int}$ characterize the spin–orbit interaction with the interface crystal potential and are determined from a comparison with experiment [Devizorova et al., 2014].

The Tamm problem in the half-space $z>0$ in this case corresponds to the solution of a single-band effective mass equation with Hamiltonian (4) and boundary conditions (III). It is important that such an approach leads to the discontinuity of single-band EFs on the interface. This is attributed to the non-perturbative effect of the interface potential. Therefore, the approach allows one to describe, for example, shallow Tamm–Shockley states, if they exist, even in the single-band approximation [Volkov and Pinsker, 1976, 1977, 1979]. Unfortunately, this situation is not realized for electrons of the conduction band at the (001) GaAs/AlGaAs hererojunction because $R<0$. A quantitative comparison with precision measurements of electron spin resonance that are sensitive to the parameters of boundary conditions shows [Devizorova et al., 2014] that $R=-2.2$ nm.

On the surface or at the interface where $R>0$, a band of shallow SSs [Volkov and Pinsker, 1979] is formed that is described, with neglect of the spin–orbit interaction, by the first two terms in the boundary conditions (III) (Fig. 1a). These SSs are ”suspended” over the bottom of the conduction band, and their localization depth is equal to $R$. The volume and interface contributions to the spin–orbit interaction, which are described by the remaining terms in (III), anisotropically split the band of SSs with respect to spin (more precisely, with respect to chirality—the projection of spin onto the direction related to momentum) (Fig. 1b). This splitting was first introduced in [Vasko, 1979] in the simplest isotropic model. As the interface spin–orbit interaction increases, the splitting increases and becomes comparable to the bandgap $E_{g}$ (Fig. 1c), and the single-band approximation fails.

IV TAMM–SHOCKLEY STATES IN THE TWO-BAND APPROXIMATION

The two-band model is a minimal multi-band model. For Dirac materials such as bismuth, lead chalсogenides, and many other, the Tamm problem in the simplest approximations reduces to the solution of the Dirac equation in the half-space $z>0$,

[TABLE]

with the boundary conditions [Volkov and Pinsker, 1981]

[TABLE]

Here the two-component EFs $\Psi_{c}$ and $\Psi_{v}$ describe the conduction and valence bands, $a_{0}$ is a real phenomenological parameter describing the small-scale structure of the surface, $\bm{n}$ is a unit normal to the surface $S$ of the crystal, and $E_{g}=2mc^{2}$ is the bulk band gap. The boundary conditions (10) follow from the Hermiticity (self-adjointness) of the Dirac Hamiltonian in the half-space and the symmetry of the problem with respect to time reversal.

A perturbative account of remote bands gives rise to diagonal corrections to the Dirac Hamiltonian in (10) that are quadratic in momentum. The corresponding band gap model, sometimes called Dimmock’s model, is used in the theory of topological insulators (see below). In this section, the contribution of remote bands is neglected.

As a result of the strong spin–orbit interaction, the Tamm problem for DFs always gives rise to a non-degenerate spectrum of SSs that fills a part of the cone surface (Fig. 2). It is interesting to compare this result with the theory of topological insulators and semi-metals. The result qualitatively depends on the sign of the boundary parameter $a_{0}$.

For a surface with a positive value of $a_{0}$, a conical (”Dirac”) point in the spectrum lies in the forbidden gap of a bulk single-valley material; i.e., in this case, SSs are characterized by a 2D spectrum of massless DFs (Figs. 1c and 2a). The comparison with topological theory implies that this case is an analog of a topologically nontrivial case. Figure 2b corresponds to the topologically trivial phase, for which the parameter $a_{0}$ is negative. In the massless limit, we obtain SSs for Dirac’s single-valley 3D semimetal (Figs. 2c, 2d).

Analytically, the spectrum of SSs for the Dirac equation in the half-space is expressed as

[TABLE]

Here the quantum number $s=\pm 1$ describes the chirality, $\bm{k}_{||}=(k_{x},k_{y},0)$. The energy of the conical point is measured from the center of the gap and has the form

[TABLE]

The branches of the spectrum (11) with different chirality are implemented under the condition

[TABLE]

Another approach to the theory of SSs is related to the formulation of models whose bandgap parameters smoothly (at the atomic scale) vary near the interface (”heterocontact”). In this case, the multi-band method of EFs can be applied in the whole space, and no boundary condition problem arises. In [Volkov and Pankratov, 1985], the authors considered interface states of DFs in the model of a symmetric inverse heterocontact. These states are described by the Dirac equation in which the mass of DFs is inverted in space, the conduction and valence bands being interchanged (Fig. 3).

Near the inversion point, a 2D band of SSs arises. The equation for DFs at the inverse junction takes a form characteristic of supersymmetric (SUSY) quantum mechanics:

[TABLE]

The zero mode of this equation describes interface states,

[TABLE]

where $e^{i\theta}=(k_{x}+ik_{y})/k_{||}$. Solution (15) is stable against to perturbations of the function $E_{g}(z)$ in view of SUSY. The Dirac point in this symmetric model is located precisely at the center of the gap, and the spectrum of SSs is actually described by formulas (11) and (12) in which $a_{0}=1$.

The model of an asymmetric contact considered in the survey [Volkov et al., 1995] is more realistic (see Fig. 4). The asymmetry of bands is defined by the scalar potential $\varphi(z)$ that have the meaning of the bending of bands.

For a greater discontinuity of bands, the spectrum of interface states is given by

[TABLE]

where $\varphi=\varphi_{0}f(z)$, $E_{g}(z)=E_{g0}f(z)+E_{g1}$, and the model function $f(z)$ describes a correlated variation of the band gap and the bending of bands at the heterointerface.

It makes sense to compare the results of [Volkov and Pinsker, 1981] and [Volkov et al., 1995]. The heterointerface spectrum (16) turns into the spectrum of SSs (11) under the change

[TABLE]

This demonstrates the physical meaning of the boundary parameter $a_{0}$: its sign correlates with the sign of the gap, and the difference of its amplitude from 1 is the measure of the e–h asymmetry. An analog of the topologically non-trivial phase in Fig. 2a corresponds to the inversion of bands. An important consequence of (16) is as follows: in general, the inversion of bands is insufficient for the emergence of supersymmetric interface states. For higher e–h asymmetry, when $\varphi_{0}/E_{g0}>1$, these states cease to exist.

In principle, the conditions for the emergence of SUSY states can be controlled by external fields. For example, consider [Enaldiev and Volkov, 2015] a solution to the Dirac equation (9) with the boundary conditions (10) for a cylindrical nanowire of radius $R_{0}$ in a longitudinal magnetic field $B$. The spectrum of SSs is controlled by magnetic flux $\Phi=\pi BR_{0}^{2}$ through a section of the wire. For the positive sign of the boundary parameter $a_{0}$ in a strong field, the quasi-classical spectrum of SSs is given by

[TABLE]

where $k_{z}$ is a projection of quasi-momentum to the wire axis, $j$ are half-integer numbers, and $\Phi_{0}=hc/e$ is the flux quantum. When the magnetic flux is a half multiple of the flux quantum, the spectrum (18) becomes gapless. For any sign of $a_{0}$, the contribution of SSs to the total density of states oscillates as a function of the magnetic flux with period $\Phi_{0}$. This leads to Aharonov–Bohm-type oscillations of the wire resistance.

More interesting is the case of a wire with inhomogeneous cross section in a uniform longitudinal magnetic field. Magnetic flux through the section of a wire is inhomogeneous, whereby local SUSY states can be formed (Fig. 5).

In [Enaldiev and Volkov, 2015], it was found that such a SUSY solution with energy $E_{0}$ has the form

[TABLE]

where $v=2a_{0}c/(1+a_{0}^{2})$ is the velocity of surface states. The effective local gap $\Delta(z)=\hbar v\left(j+\Phi(z)/\Phi_{0}\right)/R(z)$ is inverted at the point $z=z_{0}$; therefore, solution (19) is localized precisely at this point.

V TOPOLOGICAL SURFACE STATES

The theory of band structure of semiconductors involving topological arguments, which were first introduced in the theory of the quantum Hall effect [Thouless et al., 1982], gave rise to a new physical concept—a topological insulator (TI) [Hasan and Kane, 2010,Qi and Zhang, 2011]. In topological insulators with a strong spin–orbit interaction, the specific features of the bulk band structure are such that a band of conducting surface (or edge) states of the type shown in Figs. 1c and 2a should exist on the surface of TIs for topological reasons. This is the most fundamental prediction of topological theory. Surface DFs are non-degenerate with respect to the spin quantum number (chirality) and have a conical dispersion law, which leads to the weakening of their scattering by non-magnetic impurities. It is proved that this property is a specific feature of precisely topological surface or edge states.

To implement a TI phase, one usually should impose three qualitative conditions on the parameters of the bulk band structure: symmetry with respect to time reversal (T symmetry), strong spin–orbit interaction, and inversion of bands. As a result, topological SSs arise in the bulk gap; according to topological theory, the very fact of emergence of these states does not depend on the details of the structure of the near-surface region. The possibility of implementation of 2D TI [Kane and Mele, 2005] and 3D TI [Fu et al., 2007] phases was predicted in 2005 and 2007, respectively.

First, consider a 1D band of edge states in a 2D crystal system. The topological invariant $Z_{2}$ is defined by the properties of the anti-unitary matrix $w_{nm}=\langle u_{n}(\bm{k})|\hat{T}|u_{m}(-\bm{k})\rangle$, which is composed of the matrix elements of the time-reversal operator on Bloch factors of all occupied bands [Hasan and Kane, 2010]. The T symmetry allows one to classify all 2D systems into two groups that differ by the value of the $Z_{2}$ invariant $\nu=\left\{0,1\right\}$ [Kane and Mele, 2005]. A crystal is said to be in the phase of a topological (trivial) insulator if the Hamiltonian characterizing it has $\nu=1$ ($\nu=0$). A TI phase is distinguished by the existence of an odd number of Kramers pairs of edge states at the Fermi level in the bulk gap. It is important that any T-invariant perturbations of the Hamiltonian of the crystal that do not close the gap in the spectrum cannot change the value of the $Z_{2}$ invariant of the system. Topological protection of SSs in a TI is understood precisely in this sense.

In the theoretical study [Bernevig et al., 2006], the authors predicted that a 2D TI phase should be implemented in a CdTe–HgTe double heterostructure for the width of the HgCdTe quantum well greater than 6.4 nm. It is these conditions under which the electron and hole size-quantized subbands are inverted.

A two-component EF describing the states in a quantum well with two different values of the spin quantum number satisfies a modified Dirac equation (the so-called BHZ model [Bernevig et al., 2006]):

[TABLE]

Here $\varepsilon(k)=C-D(k_{x}^{2}+k_{y}^{2})$, $m(k)=M-B(k_{x}^{2}+k_{y}^{2})$, $A,B,C,D$, and $M$are the parameters of the system, and $k_{x}$ and $k_{y}$ are the components of 2D momentum. The inversion of bands corresponds to a positive value of the product of parameters $MB>0$.

The Tamm problem in the half-plane $x>0$ in the model of zero boundary conditions gives rise to spin-split edge states with linear dispersion. Recently, 1D conducting edge channels associated with these states have been observed experimentally [Konig et al., 2007,Nowack et al., 2013]. However, the spectrum of edge states has not been measured.

In [Enaldiev et al., 2015], the authors considered the problem of sensitivity of edge states to the form of nonzero boundary conditions. They derived a general T-invariant boundary condition for the EF $\Psi=\left(\Psi_{\uparrow},\Psi_{\downarrow}\right)^{T}$ satisfying Eq. (20):

[TABLE]

where $G$, a matrix that mixes different components of the EF, is defined by six real phenomenological parameters and the matrix $F$ is defined by the bulk parameters $B,A$, and $D$; this matrix is introduced into the boundary condition (21) for convenience. The spatial symmetry reduces the number of independent phenomenological parameters in the matrix $G$ to four ones. The authors showed that the difference of boundary conditions from zero conditions significantly changes the spectrum of edge states and, generally speaking, can even lead to the vanishing of these states from the region of small quasi-momenta, where model (20) is applicable.

Note the following important feature of this system. Although non-zero values of the parameters $B$ and $D$ are assumed to be important for topological classification, they make a small contribution to the energy of size-quantized subbands for actual quantum wells of Hg(Cd)Te [Bernevig et al., 2006]. Therefore, one can neglect these parameters in the Tamm problem, and Eq. (20) reduces to a 2D version of the Dirac equation. For the latter equation, the spectrum of the Tamm problem is defined by formula (11) in which $\bm{k}_{||}$ is an conserving projection of quasimomentum to the 1D boundary of the 2D system considered. In this case, information on the inversion of bands is contained in the sign of the boundary parameter $a_{0}$ rather than in the bulk parameters of the BHZ equation (20).

The topological classification is complicated as the dimensionality increases. For a 3D TI, one usually defines four topological invariants ($\nu_{0};\nu_{1},\nu_{2},\nu_{3}$), each of which can take values 0 or 1. A 3D crystal characterized by $\nu_{0}=1$ is called a strong TI. In the simplest case, a gapless 2D band of non-degenerate massless DFs is formed on the surface of a strong 3D TI. 3D crystals with $\nu=0$ are called weak TIs for which the appearance of SSs in the forbidden band depends on both the values of $\nu_{1,2,3}$ and the translational symmetry of the surface [Fu and Kane, 2007]. In the general case, weak TIs have no gapless SSs.

It has been theoretically shown that a whole class of Bi2Se3-type compounds with inversion symmetry should belong to the TI phase. In this case, it is the strong spin–orbit interaction that leads to the inversion of two bands with different parity that are nearest to the Fermi level (Fig. 6). The existence of SSs in Bi2Se3, Bi2Te3, and Sb2Te3 compounds was experimentally established in [Hsieh et al., 2009a,Hsieh et al., 2009b]. A standard 3D TI for ARPES measurements is the crystalline compound Bi2Se3 in which the gap is about 0.3 eV. By this method, it was established [Hsieh et al., 2009a] that on the (111) surface there exist non-degenerate SSs with massless Dirac spectrum in the gap of the crystal.

In the approximation of EFs, these compounds are described near the center of the Brillouin zone by a modified Dirac equation, which, with neglect of the anisotropy and the e–h asymmetry, has the form

[TABLE]

where the significantly dispersive mass term $M(\bm{k})=M_{0}-B\bm{k}^{2}$ distinguishes this equation from the true Dirac equation. A topological phase is characterized by the positive value of the product of parameters $M_{0}B>0$, which is attributed to the inversion of bands. The break of the e–h symmetry is described by the introduction of an additional diagonal contribution to Eq. (22). We will neglect this asymmetry.

To describe SSs within Eq. (22), one should specify boundary conditions that take into account the atomic-scale structure of the surface. To this end, one often uses the fact that the above-mentioned modification of the Dirac equation doubles the order of the system of differential equations. Therefore, it is simplest to use zero boundary conditions for all components of EFs [Liu et al., 2010, Shan et al., 2010]: $\left[\Psi_{c}=\Psi_{v}\right]_{S}=0$. Such boundary conditions, which are often called ”open” do not take into account the details of the microscopic structure of the surface at the atomic scale, but nevertheless guarantee the existence of massless DFs on a plane surface.

However, it was demonstrated that non-zero boundary conditions may significantly affect the spectrum of SSs in the case of 3D TIs as well [Enaldiev et al., 2015, Men’shov et al., 2014]. General boundary conditions for a four-component EF mix the values of the EF with its normal derivative on the boundary $S$:

[TABLE]

where the $4\times 4$ matrix $Q$ should be determined from additional considerations. The T symmetry, along with the spatial symmetry, reduces the number of independent real phenomenological parameters in the matrix to three ones. For sufficiently strong surface spin–orbit interaction described by these parameters, SSs may disappear from the forbidden gap in the neighbourhood of the point $\Gamma$ at which the EF method is applicable [Enaldiev et al., 2015].

Experiments have shown that the properties of SSs are sensitive to the perturbations of the surface potential, which can be taken into account by an appropriate choice of the parameters of the boundary conditions. For example, the doping of the surface with non-magnetic and magnetic atoms shifts the position of the Dirac point and modifies the dispersion of SSs [Valla et al., 2012; Scholz et al., 2012; Bianchi et al., 2011; Roy et al., 2014].

Another Dirac system is given by a topological crystalline insulator (TCI). In a TCI, the protection of SSs is guaranteed by two symmetries at once: the T symmetry and the spatial symmetry of the crystal and its surface [Fu, 2011]. A topological classification is performed similar to the above-considered case of TIs, with the difference that the matrix $w_{nm}$ is constructed on the operator $\hat{U}\hat{T}$, where $\hat{U}$ is the operator of spatial symmetry.

Lead chalcogenides and Pb1-xSnxTe(Se) solid solutions pass to the TCI phase at certain values of the concentration of tin when the inversion of bands of different parities occurs at the $L$ points [Hsieh et al., 2012; Xu et al., 2012; Dziawa et al., 2012; Egorova et al., 2015]. The role of spatial symmetry in these compounds is played by mirror symmetry, which persists not for any orientations of the surface. In the Brillouin zone of this system, there are four inequivalent $L$ valleys in the forbidden band of each of which massless surface 2D DFs appear. On a (001) surface, pairs of Dirac cones of SSs are projected to a single point, leading to the degeneracy of the spectra of SSs. Atomically abrupt perturbations preserving mirror symmetry lift this degeneracy, moving apart the Dirac points in the reciprocal space without opening a gap in the spectrum of SSs [Hsieh et al., 2012; Xu et al., 2012; Dziawa et al., 2012]. However, perturbations breaking the mirror symmetry of the surface give rise to a gap in the spectrum of SSs [Zeljkovic et al., 2015], thus breaking the topological protection. The fundamental features of the spectrum of SSs in any phase of TCIs based on lead chalcogenides can be found when one neglects the intervalley interaction in the solution of an anisotropic Dirac equation with boundary conditions (10).

VI EDGE STATES IN GRAPHENE

The Tamm–Shockley states arise far from always. The following questions are fundamental: Under what conditions, except for those implemented in TIs, do Tamm–Shockley-type states arise? Do they exist on real interfaces under normal conditions (rather than in vacuum)? Does the conductivity through these states have a band character? The analysis carried out above shows that, in a number of Dirac materials, such states exist (in contrast to, for example, GaAs-type semiconductors).

Quite recently, edge states of this type have been experimentally observed by resonance methods in another 2D Dirac material, graphene. It was found [Latyshev et al., 2013] that the spectrum of Tamm–Shockley edge states in graphene is linear with respect to momentum, which is similar to a section of the conical spectrum of SSs in TIs. Moreover, the band character of the conductivity through these states in normal conditions [Latyshev et al., 2014] was proved by a direct transport experiment. Let us describe a theory necessary to explain these experiments.

2D DFs in graphene are described by a two-valley system of equations for two-component functions $\Psi_{K}$ and $\Psi_{K^{\prime}}$ describing EFs in two valleys, $K$ and $K^{\prime}$:

[TABLE]

where $\bm{\sigma}=(\sigma_{x},\sigma_{y})$ are the Pauli matrices in standard representation and $v\approx 10^{6}$ m/s.

General boundary conditions can be determined, as usual, from two requirements: the Hermiticity of Hamiltonian (24) and the invariance of the boundary conditions with respect to time reversal. They have a rather complicated form [Zagorodnev et al., 2011; McCann and Fal’ko, 2004; Zagorodnev et al., 2015; Akhmerov and Beenakker, 2008; Basko, 2009; Tkachov, 2009]

[TABLE]

Here $\beta,\varphi$, and $\gamma$ are real phenomenological parameters.

With the neglect of intervalley scattering, the system can be described by a pair of independent 2D Weyl equations with one-parameter boundary conditions [Volkov and Zagorodnev, 2009,van Ostaay et al., 2011] for the single-valley function $\Psi_{K(K^{\prime})}=\left(\psi_{1K(K^{\prime})},\psi_{2K(K^{\prime})}\right)^{T}$:

[TABLE]

where $a$ is a real phenomenological parameter that determines the properties of the edge of graphene at the atomic scale and the index $\tau=\pm 1$ numbers the valleys: $\tau=1$ ($\tau=-1$) in valley $K$ ($K^{\prime}$).

The Tamm problem (26) is easily solved. The wave function of edge states is exponentially localized near the boundary. A typical localization length determined from experiment [Latyshev et al., 2014] is 2 nm. The value of $a$ should be determined from comparison with experiment; however, from a comparison with model microscopic calculations [van Ostaay et al., 2011,Koskinen et al., 2008], one should expect that it is small: $|a|\ll 1$.

For a graphene sample in the form of a half-plane, the spectrum of these states in the absence of a magnetic field has the form [Basko, 2009; Tkachov, 2009; Volkov and Zagorodnev, 2009; van Ostaay et al., 2011] (Fig. 7):

[TABLE]

Here $k_{||}$ is the one-dimensional momentum of an electron along the edge, measured from the center of the valley. It is important that the sign of $k_{||}$ is determined by the number of the valley. For clarity, different valleys in the figure are shown in different colors. The localization length of an edge state is equal to $1/k_{||}$.

In the case of a perpendicular magnetic field, two types of edge states exist: Tamm–Shockley states and magnetic edge states due to skipping orbits. They interact with each other, and the result of their interference depends on the geometry of a sample. For the half-plane $x>0$, the resulting spectrum [Zagorodnev et al., 2011] is shown in Fig. 8.

Edge states also exist on the boundary of a circular hole (an antidot with radius $R_{0}$) in an infinite sheet of graphene. The finiteness of the perimeter of the antidot leads to the quantization of the momentum component parallel to the boundary. Now discrete edge states are characterized by the total angular momentum $j=\pm 1/2,\pm 3/2,\pm 5/2,\dots$, different signs of $j$ corresponding to different valleys. Within a quasiclassical approach, the spectrum of such edge states is obtained from (28) by replacing $k_{||}$ by a multiple of $1/R_{0}$. A more accurate asymptotes ($|a|\ll 1$) yields a decay-type complex spectrum:

[TABLE]

Here $\tau j>0$, $\Gamma(z)$ is a gamma function, and $j=\pm 3/2,\pm 5/2,\dots$. In a special case of $j=\pm 1/2$, the energy of edge states is determined by the root of the equation

[TABLE]

Although these states are quasistationary in the absence of a magnetic field, their lifetime with respect to the decay into bulk states is large in the actual case of small $a$. DFs trapped in edge states perform a cyclic motion along the perimeter of the antidot (clockwise in one valley, or counterclockwise in the other), acquiring an additional Aharonov–Bohm phase in a magnetic field.

The spectrum of edge states in the antidot in a perpendicular magnetic field is controlled by a magnetic flux $\Phi=\pi R_{0}^{2}B$ through the area of the antidot:

[TABLE]

The applicability conditions of this asymptotes is given by

[TABLE]

Asymptotes (31) is the more justified, the less the effect of the magnetic field on the orbital part of the wave function, i.e., when the localization length of an edge state is small compared with the magnetic length. In the general case, a numerically calculated spectrum is shown in Fig. 9 [Enaldiev and Volkov, 2011].

These states explain the recently revealed Aharonov–Bohm effect in the resistivity of nanoperforated graphene. From comparison with experiment, we can determine the value of the boundary parameter $a=-0.05$.

Thus, we can conclude that hole-type DFs slowly rotate around every antidot in nanoperforated graphene. The sign of rotation is determined by the number of the valley, and the speed of rotation is ten times less than the volume velocity of DFs in graphene.

VII CONCLUSIONS. MINIMAL MODEL OF SURFACE STATES FOR DIRAC ELECTRONS

One can draw the following qualitative conclusion from the survey of literature data on Dirac materials: the Dirac property favours the emergence of a Tamm–Shockley-type band of edge and surface states. However, this band is not always situated in the bulk gap where its spectrum is analogous to the spectrum of topological states. For some values of boundary parameters, this band is expelled from the gap and its energy overlaps with the energy of bulk bands.

The analysis carried out allows us to formulate a minimal model suitable for the analytic investigation of fundamental features of the spectrum of surface or edge states in a number of Dirac systems: a generally anisotropic version of the 3D Dirac equation (9) with boundary conditions (10) and their 2D generalizations.

Such a model can be applied with neglect of intervalley interaction to describe SSs in the following Dirac systems.

In 3D systems: lead chalcogenides and solid solutions based on lead chalcogenides Pb1-xSnxSe(Te)-type, bismuth- and Bi1-xSbx-type semimetals at the $L$ point, both in the inverted and ordinary, non-inverted, modes.

In 2D systems: graphene and 2D TIs in quantum wells based on a HgTe–CdTe heteropair.

The minimal model cannot be applied to describe SSs in Bi2Se3, Bi2Te3, and Sb2Te3 because of a significant contribution of remote bands that forms the dispersion of the mass term in a modified Dirac equation.

VIII ACKNOWLEDGMENTS

We are grateful to I.V. Zagorodnev for critical remarks.

This work was supported in part by the Russian Foundation for Basic Research.