Surface Green's functions and boundary modes using impurities: Weyl semimetals and topological insulators

TL;DR

This paper introduces a new analytical method to compute surface Green's functions for 3D topological systems, enabling precise analysis of boundary states like Fermi arcs without numerical approximations.

Contribution

The authors develop a direct, non-numerical technique based on T-matrix formalism to obtain surface Green's functions for topological insulators and Weyl semimetals.

Findings

Calculated surface Green's functions for Weyl semimetals.

Analyzed Fermi-arc surface states.

Extended method to Kane-Mele and Chern insulators.

Abstract

In this work we provide a new direct and non-numerical technique to obtain the surface Green's functions for three-dimensional systems. This technique is based on the ideas presented in Phys. Rev. B 100, 081106(R), in which we start with an infinite system and model the boundary using a plane-like infinite-amplitude potential. Such a configuration can be solved exactly using the T-matrix formalism. We apply our method to calculate the surface Green's function and the corresponding Fermi-arc states for Weyl semimetals. We also apply the technique to systems of lower dimensions, such as Kane-Mele and Chern insulator models, to provide a more efficient and non-numerical method to describe the formation of edge states.

Equations

$$G\left({\mathbf{k}}_1,{\mathbf{k}}_2,i{\omega }_n\right)$$

$$+G_0\left({\mathbf{k}}_1,i{\omega }_n\right)T\left({\mathbf{k}}_1,{\mathbf{k}}_2,i{\omega }_n\right)G_0\left({\mathbf{k}}_2,i{\omega }_n\right)$$

$$T\left(i{\omega }_n\right)={\left[\mathbb{I}-V\cdot \int \frac{d{k}_x}{L_k}{G}_0\left(k_x,i{\omega }_n\right)\right]}^{-1}\cdot V$$

$$T$$

$$={\delta }_{k_{1y},k_{2y}}{\left[\mathbb{I}-V\cdot \int \frac{d{k}_x}{L_k}{G}_0\left(k_x,k_{1y},i{\omega }_n\right)\right]}^{-1}\cdot V$$

$$T$$

$$\times {\left[\mathbb{I}-V\cdot \int \frac{d{k}_x}{L_k}{G}_0\left(k_x,k_{1y},k_{1z},i{\omega }_n\right)\right]}^{-1}\cdot V,$$

$${\mathcal{G}}_s(k_y,k_z)\equiv \mathcal{G}(x=x^{\prime }=\pm 1;k_y,k_y;k_z,k_z)$$

$$=\int \frac{d{k}_{1x}}{L_k}\int \frac{d{k}_{2x}}{L_k}\mathcal{G}(k_{1x},k_{2x};k_y,k_y;k_z,k_z)e^{i{k}_{1x}x}{e}^{-i{k}_{2x}{x}^{\prime }}.$$

$$A(k_y,k_z,E)=-\frac{1}{\pi }\mathrm{Im}\{\mathrm{tr}[{\mathcal{G}}_s\left(k_y,k_z,E\right)]\}.$$

$${\mathcal{G}}_l(k_y)=\mathcal{G}(x=x^{\prime }=\pm 1;k_y,k_y)$$

$$=\int \frac{d{k}_{1x}}{L_k}\int \frac{d{k}_{2x}}{L_k}\mathcal{G}(k_{1x},k_{2x};k_y,k_y)e^{i{k}_{1x}x}{e}^{-i{k}_{2x}{x}^{\prime }}.$$

$$A(k_y,E)=-\frac{1}{\pi }\mathrm{Im}\{\mathrm{tr}[{\mathcal{G}}_l\left(k_y,E\right)]\}.$$

$$\delta N(k_y,E)=\int \frac{d{k}_x}{L_k}\delta A(k_x,k_y,E).$$

$$\delta A(\mathbf{k},E)=-\frac{1}{\pi }\mathrm{Im}\{\mathrm{tr}[{\mathcal{G}}_0\left(\mathbf{k}\right)T\left(\mathbf{k},\mathbf{k}\right){\mathcal{G}}_0\left(\mathbf{k}\right)]\}.$$

$$H_{1,2}=\sum _{\mathbf{k}}{\psi }^{†}(\mathbf{k}){\mathcal{H}}_{1,2}(\mathbf{k})\psi (\mathbf{k}),$$

$${\mathcal{H}}_1(\mathbf{k})$$

$$g_0(\mathbf{k})$$

$$g_1(\mathbf{k})$$

$$g_2(\mathbf{k})$$

$$g_3(\mathbf{k})$$

$$S_A^x(k_x,k_z)=-\frac{1}{\pi }\mathrm{Im}\left[{\mathcal{G}}_{12}^s\left(k_x,k_z\right)+{\mathcal{G}}_{21}^s\left(k_x,k_z\right)\right]$$

$$S_B^x(k_x,k_z)=-\frac{1}{\pi }\mathrm{Im}\left[{\mathcal{G}}_{34}^s\left(k_x,k_z\right)+{\mathcal{G}}_{43}^s\left(k_x,k_z\right)\right]$$

$$S_A^y(k_x,k_z)=-\frac{1}{\pi }\mathrm{Re}\left[{\mathcal{G}}_{12}^s\left(k_x,k_z\right)-{\mathcal{G}}_{21}^s\left(k_x,k_z\right)\right]$$

$$S_B^y(k_x,k_z)=-\frac{1}{\pi }\mathrm{Re}\left[{\mathcal{G}}_{34}^s\left(k_x,k_z\right)-{\mathcal{G}}_{43}^s\left(k_x,k_z\right)\right]$$

$$S_A^z(k_x,k_z)=-\frac{1}{\pi }\mathrm{Im}\left[{\mathcal{G}}_{11}^s\left(k_x,k_z\right)-{\mathcal{G}}_{22}^s\left(k_x,k_z\right)\right]$$

$$S_B^z(k_x,k_z)=-\frac{1}{\pi }\mathrm{Im}\left[{\mathcal{G}}_{33}^s\left(k_x,k_z\right)-{\mathcal{G}}_{44}^s\left(k_x,k_z\right)\right],$$

$${\mathcal{H}}_2(\mathbf{k})$$

$${\mathcal{H}}_{\mathrm{TB}}=\sum _{⟨ij⟩,\alpha }t{c}_{i,\alpha }^{†}{c}_{j,\alpha }+\sum _{⟨⟨ij⟩⟩,\alpha ,\beta }i{t}_2{\nu }_{ij}{s}_{\alpha \beta }^z{c}_{i,\alpha }^{†}{c}_{j,\beta }$$

$${\mathcal{H}}_{\mathbf{k}}=\left(\begin{array}{cccc}{h}_{NNN}& 0& h_{NN}& 0\\ 0& -h_{NNN}& 0& h_{NN}\\ h_{NN}^{\ast }& 0& -h_{NNN}& 0\\ 0& h_{NN}^{\ast }& 0& h_{NNN}\end{array}\right),$$

$$h_{NN}=t\left[1+e^{i\sqrt{3}{k}_y}+e^{i\frac{\sqrt{3}}{2}{k}_y}{e}^{-i\frac{3}{2}{k}_x}\right],\mathrm{and}$$

$$h_{NNN}=2t_2\left[2\cos \left(\frac{3}{2}{k}_x\right)\sin \left(\frac{\sqrt{3}}{2}{k}_y\right)-\sin \left(\sqrt{3}{k}_y\right)\right]$$

$$V=U\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),$$

$$\delta N(k_y,E)=\int \frac{d{k}_x}{L_k}\delta A(k_x,k_y,E),$$

$${\mathcal{H}}_{\mathbf{k}}$$

$$+B(2-M-\cos k_x-\cos k_y){\sigma }_z$$

$$T(k_y,i{\omega }_n)$$

$$=-{\left[\underset{-\pi }{\overset{\pi }{\int }}\frac{d{k}_x}{2\pi }{G}_0(k_x,k_y,i{\omega }_n)\right]}^{-1}$$

Full text

I Introduction

Boundaries of certain condensed matter systems host unique phenomena. For instance, graphene exhibits zero-energy zigzag-edge modes Fujita et al. (1996), and topological insulators exhibit conducting edge or surface states Kane and Mele (2005); Andrei Bernevig and Hughes (2013); Hsieh et al. (2008). In order to describe boundary effects, several techniques were developed, including the exact diagonalization of tight-binding Hamiltonians Slater and Koster (1954); Busch and Penson (1987), iterative methods to compute boundary Green’s functions Sancho et al. (1984, 1985); Peng et al. (2017), solving the Schrödinger equation Duncan et al. (2018) and the bulk-boundary correspondence Rhim et al. (2018).

A new method describing the formation of boundary modes was recently introducedKaladzhyan and Bena (2019). This method can be generalized to any dimensions, and in certain situations it can yield fully analytical results, providing a deeper physical insight than numerical techniques. The general idea is as follows: instead of considering a finite system with a sharp boundary, we consider an infinite system with a strong delta-potential impurity emulating the shape of the boundary. For example, in order to recover end, edge or surface boundaries, the impurity potential should be chosen to be point-like, line-like and plane-like, respectively. In the limit of an infinite impurity potential such impurities divide a given system into two independent semi-infinite regions. Subsequently we use the $T$-matrix formalism Bena (2016); Balatsky et al. (2006) to study the impurity-induced states which transform into boundary states when the impurity strength is larger than any energy scale in the system.

Along the same lines, we present here a direct and non-numerical technique to calculate the surface Green’s functions of an arbitrary three-dimensional system. The boundary can once more be modeled as a plane impurity potential with an amplitude going to infinity. The corresponding full Green’s functions can be calculated exactly using the $T$-matrix formalism. The resulting Green’s function evaluated on the plane neighboring and parallel to the impurity plane becomes the surface Green’s function (see Fig. 1). We apply this technique to calculate the surface Green’s functions for Weyl semimetals described by two different modelsKourtis et al. (2016); Lau et al. (2017). We recover in each case the corresponding Fermi-arc states.

Moreover, in this work we apply the technique from Ref. [Kaladzhyan and Bena, 2019] to a new class of systems – topological insulators. In particular we consider a 2D honeycomb lattice described by the Kane–Mele modelKane and Mele (2005), as well as a 2D Chern insulator Andrei Bernevig and Hughes (2013). We show that impurity-induced states in these two models transform into helical or chiral edge modes, respectively, when the impurity potential is taken to infinity. While the Kane–Mele model requires performing a numerical integration, the 2D Chern insulator allows an exact closed-form solution and thus demonstrates the analytical power of the method.

The paper is organized as follows: in Sec. II we introduce the formalism and the notations. In Sec. III we present the calculation of the surface Green’s functions for a Weyl semimetal described by two different models, and the formation of the corresponding Fermi-arc states. In section IV and V we focus on two-dimensional topological insulators described by the Kane–Mele and Chern-insulator models, respectively, and we obtain the corresponding edge modes. We leave the conclusions to Sec. VI.

II $T$-matrix formalism for surface Green’s functions and edge states

Below we consider an infinite system described by a momentum-space Hamiltonian $\mathcal{H}_{\bm{k}}$. The unperturbed Matsubara Green’s function (GF) can be written as: $G_{0}\left(\bm{k},i\omega_{n}\right)=\left[i\omega_{n}-\mathcal{H}_{\bm{k}}\right]^{-1}$, where $\omega_{n}$ denote the Matsubara frequencies. In the presence of an impurity, the Green’s function is modified to:

[TABLE]

where the $T$-matrix $T\left(\bm{k}_{1},\bm{k}_{2},i\omega_{n}\right)$ embodies all-order impurity-scattering processes Balatsky et al. (2006); Mahan (2000). Note that due to the impurity-induced breaking of translational symmetry, and consequently of the momentum conservation, the generalized Green’s function depends no longer on one, but on two values of momentum. For a delta-function impurity $V_{\mathrm{imp}}\left(\bm{r}\right)\equiv V\delta\left(x\right)$, the form of the $T$-matrix in 1D is momentum independent and is given by Balatsky et al. (2006); Ziegler et al. (1996); Salkola et al. (1996); Bena (2008):

[TABLE]

while in 2D and 3D, for a line and plane impurity respectively, localized at $x=0$ and perpendicular to the $x$ axis, we have:

[TABLE]

and

[TABLE]

respectively, with $L_{k}$ being a normalization factor. The limits of integration are given by the boundaries of the first Brillouin zone, i.e., for a 1D system and for a 2D square lattice or 3D cubic lattice we integrate from $-\pi$ to $\pi$ (with $L_{k}=2\pi$), while for a honeycomb lattice with an impurity along $y$ we integrate from $-2\pi/3$ to $2\pi/3$ ($L_{k}=4\pi/3$) (for a full justification see Appendix A). Note that Eqs. (3) and (II) are independent of $k_{1x}$ and $k_{2x}$ due to the fact that the impurity potential is a delta-function centered at $x=0$. Reversely, note that the T-matrix contains the terms $\delta_{k_{1y},k_{2y}}$ (2D) and $\delta_{k_{1y},k_{2y}}\delta_{k_{1z},k_{2z}}$ (3D), since the impurity is independent of $y$ in 2D and of $y$ and $z$ in 3D, and therefore, the momenta in the corresponding directions are conserved in all scattering processes.

The exact same formalism can be applied for impurities perpendicular to the other axes of the systems.

In what follows we employ this formalism at zero temperature to calculate the retarded GF $\mathcal{G}(\bm{k}_{1},\bm{k}_{2},E)$ obtained by the analytical continuation of the Matsubara GF $G(\bm{k}_{1},\bm{k}_{2},i\omega_{n})$ (i.e., by setting $i\omega_{n}\rightarrow E+i\delta$, with $\delta\rightarrow 0^{+}$).

For a three-dimensional system the surface Green’s function can be extracted from the perturbed generalized Green’s function in Eq. (1). This can be related to the mixed Green’s function in which we keep the momentum coordinates in the two directions parallel to the impurity plane ($k_{y}$ and $k_{z}$), but we perform a Fourier transform to write down the Green’s function in real space coordinates in the $x$ direction. Note that for a plane impurity the generalized Green’s function depends on two different values of momentum only for the direction perpendicular to the impurity, in the other two directions we recover a simple dependence on momentum due to unbroken translational invariance. Thus we have:

[TABLE]

We fix $x=x^{\prime}=\pm 1$ since we are interested in describing the lattice planes one lattice constant away from the impurity (see Fig. 1). The boundaries of the two resulting semi-infinite systems correspond to the two planes at $x=\pm 1$, and the Green’s functions taken at these two positions are effectively the surface Green’s functions for the semi-infinite systems. The physics at $x=0$ (impurity position) is relatively trivial since the infinite-amplitude impurity potential pushes away all the wave function weight off the impurity plane. For simplicity we have omitted writing down explicitly the energy dependence of the Green’s functions. Note again that translational invariance holds within the planes parallel to the impurity plane, thus the surface Green’s functions depend only on one momentum in each of the in-plane directions. Furthermore, in order to obtain the surface physics the value of the impurity potential in Eq. (II) needs to be set to a value much larger than all the energy scales in the problem.

The surface Green’s function allows us to recover the formation of the surface states such as, for instance, the Fermi-arc states. Thus, we can study the surface spectral function

[TABLE]

The same analysis can be performed for a two-dimensional system with a line-impurity to find the line Green’s functions

[TABLE]

and the corresponding edge states given by:

[TABLE]

Alternatively, in order to visualize the impurity-induced states, as described in Appendix B, we may focus on the average correction to the spectral function:

[TABLE]

where

[TABLE]

Above $\mathcal{G}_{0}\left(\bm{k}\right)$ stands for $\mathcal{G}_{0}\left(\bm{k},E\right)$ and $T\left(\bm{k},\bm{k}\right)$ for $T\left(\bm{k},\bm{k},E\right)$. The integral over $k_{x}$ is performed along the same interval as the one defined in Eq. (3). This quantity corresponds to the average number of available electronic states with wavevector $(k_{y},k_{z})$, where the average is performed along the direction perpendicular to the impurity. A more detailed description of the significance of this quantity is provided in Appendix B.

III Weyl semimetals

In what follows we consider a Weyl semimetal: there is a large number of models of various degree of complexity describing such a system. Here we focus only on the tight-binding models described in Refs. [Kourtis et al., 2016] and [Lau et al., 2017], which we denote by $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, respectively. The Bloch Hamiltonians for these two systems are given by

[TABLE]

where $\psi(\bm{k})=(c_{\bm{k}A\uparrow},c_{\bm{k}A\downarrow},c_{\bm{k}B\uparrow},c_{\bm{k}B\downarrow})$ is a spinor with the index $A/B$ denoting a generic unspecified orbital component, and the $\uparrow/\downarrow$ the physical spin.

For the model in Ref. [Kourtis et al., 2016] written in the basis above we have

[TABLE]

where

[TABLE]

and $\alpha$, $\beta$ are real parameters. The $2\times 2$ identity matrices $\sigma_{0}/\tau_{0}$ and the Pauli matrices $\sigma_{i}/\tau_{i}$, $i=1,2,3$ act in the spin and the orbital spaces, correspondingly, and the multiplication of the $\sigma$ and $\tau$ matrices indicates a tensor product.

We consider the same values of parameters as those in Ref. [Kourtis et al., 2016], thus we take a) $a=b=1$, $t=-1$, $m=0.5$, $d=0.8$, $\alpha=\beta=0$ and b) $a=b=1$, $t=-1.5$, $d=m=0$, $\beta=0.9$, and $\alpha=0.3$. The former is characterized by two Weyl points, while the latter by four Weyl points, and thus we expect to have one and two Fermi arcs, respectively.

In order to obtain the Fermi-arc surface states we need to introduce a surface into the system in such a way that the vector connecting the Weyl nodes has a nonzero projection onto it. For example, for the above model we choose to have a plane-like impurity at $y=0$, hence perpendicular to the $y$ direction. The resulting surface Green’s functions are described by the formalism in Sec. II, where we consider an impurity $V=U\delta(y)\mathbb{I}_{4}$, with $U\rightarrow\infty$ (i.e, much larger than all energy scales in the problem). Note that in Sec. II we describe an impurity at $x=0$ and not $y=0$, however the $y=0$ formalism is obtained by simply interchanging $x$ and $y$ in the corresponding formulas. The spectral function for the surface states $A(k_{x},k_{z},E)=-\frac{1}{\pi}\operatorname{Im}\{\operatorname{tr}[\mathcal{G}_{s}\left(k_{x},k_{z},E\right)]\}$ is depicted in Fig. 2 for two chosen configurations of parameters.

We note that the Fermi-arc states calculated using our method agree exactly with those predicted in Ref. [Kourtis et al., 2016]. Moreover, the full surface Green’s function that we have obtained contains all the information required to describe these states, such as their spin and orbital distribution, the full energy dispersion, etc.. For instance, we have for the different spin and orbital components:

[TABLE]

where we omit the energy dependence for the sake of brevity. In Fig. 3 we plot the $x$ and $z$ spin components of the Fermi-arc states at zero energy, separately calculated for the A and B orbitals. The parameters chosen correspond to Fig. 2 b). We do not plot the $y$ component since it is zero for both orbitals.

The spins of opposite arcs are of opposite signs, as expected Kourtis et al. (2016).

We perform a similar analysis on a different Weyl semimetal model, introduced in Ref. [Lau et al., 2017]:

[TABLE]

We consider the same values of parameters as those in Ref. [Lau et al., 2017], thus we take a) $a=b=1$, $t=-1.5$, $\lambda=0.5$, $d=0.1$, $\alpha=0.3$ and $\beta=0.7$ and b) $a=b=1$, $t=-1.5$, $\lambda=0.5$, $d=0.1$, $\alpha=0.3$ and $\beta=0.4$. These configurations are characterized by four Weyl points, and thus one expects two Fermi arcs on the surface, however Ref. [Lau et al., 2017] indicates the possible existence of an electron pocket coming from the bulk bands. We apply the same techniques as above, and we show in Fig. 4 the resulting spectral function $A(k_{x},k_{z},E)$ for the surface states, for the two chosen configurations of parameters.

Our results agree exactly with those in Ref. [Lau et al., 2017]. Furthermore, we compute the spin and orbital properties for this model (see Fig. 5).

We have also checked that for the first set of parameters in Fig. 4 the spins of the two Fermi arcs are opposite, same as for the $\mathcal{H}_{1}$ model. The $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ models differ in this case mainly by a nonzero $y$ component in the $\mathcal{H}_{2}$ model and the asymmetry of the two Fermi arcs in $k_{z}$.

IV Kane–Mele model

We start with the Kane–Mele model of a topological insulator on a honeycomb lattice Kane and Mele (2005). Therefore, we employ the following tight-binding model:

[TABLE]

where $c_{i,\alpha}^{\dagger}$ creates an electron on the lattice site $i$, with spin $\alpha=\uparrow,\downarrow$. The first term in Eq. (16) is the standard nearest-neighbor hopping term corresponding to the tight-binding Hamiltonian of graphene, which yields a spectrum with bands touching at the Dirac points situated at the Brillouin zone corners $(\pm 4\pi/3\sqrt{3},0)$, $(\pm 2\pi/3\sqrt{3},\pm 2\pi/3)$. In order to turn graphene into an insulator we add a next-nearest-neighbor term with a spin-dependent amplitude $\nu_{ij}=-\nu_{ji}=\pm 1$, defined by the orientation of the hopping direction (see Fig. 6). The second term opens a bulk gap in the energy spectrum at the Dirac points.

First, we obtain the boundary modes numerically by diagonalizing the tight-binding Hamiltonian in Eq. (16) and considering periodic boundary conditions in the $y$ direction and open boundary conditions in the $x$ direction. This corresponds to a ribbon with zigzag edges. For convenience we set the lattice spacing $a$ to unity. The corresponding energy spectrum is shown in Fig. 7. Note the formation of two subgap states (we have verified that these are actually edge states). The momentum-space dispersions for the two edge states cross at zero energy for $k_{y}=\pi/\sqrt{3}$.

In what follows we reproduce the formation of these edge states by considering a line impurity in an infinite system and subsequently taking the impurity potential to infinity. We can rewrite the tight-binding Hamiltonian in Eq. (16) in momentum space. Thus in the basis $(c_{A\uparrow},c_{A\downarrow},c_{B\uparrow},c_{B\downarrow})$, where $c_{i\sigma}$ is an electron operator with spin $\sigma=\uparrow/\downarrow$ on the sublattice $i=A/B$, the Kane-Mele Hamiltonian is expressed as:

[TABLE]

with

[TABLE]

being the nearest-neighbor and the next-nearest-neighbor terms with amplitudes $t$ and $t_{2}$, respectively. Here ∗ simply denotes the complex conjugation.

To reproduce the zigzag edge states, we choose an impurity potential localized on two adjacent rows of atoms corresponding to two different sublattices.

[TABLE]

(see Appendix D for more details).

In order to visualize the impurity-induced states, in Fig. 8 we plot the average correction to the spectral function due to the impurity, as defined in Eq. (9),

[TABLE]

for the same range of $k_{y}$ as the one considered in Fig. 7.

For a weak impurity, i.e., $U=t=1$, the impurity states appear as a distinct band at energies concentrated mostly outside of the gap. We expect that the impurity bound states will evolve into edge states and acquire the same properties (i.e., the same momentum dispersion) as the edge states derived previously using numerical methods, when the impurity strength $U$ goes to infinity. Indeed, for a stronger impurity potential $U=100$, the impurity-induced spectral function exhibits subgap states with the same dispersion as the ones derived via exact diagonalization and depicted in Fig. 7. The agreement between the two methods is remarkable, confirming the validity of our analytical approach towards finding the edge states of a simple topological insulator system.

While here we consider only zigzag edges, in Appendix D we have also considered the case of an impurity localized only on one row of atoms, which splits the systems into two subsystems with different edges, one with a zigzag edge, and one with a bearded edge. We expect that we will recover two distinct sets of edge states, and in Appendix D we show that this is indeed the case.

V Chern insulator

Below we consider the simplest lattice model defining a Chern insulator

[TABLE]

where we set the lattice constant to unity and $t=1$. Here $\bm{k}\equiv(k_{x},\,k_{y})$ and $\bm{\sigma}=(\sigma_{x},\,\sigma_{y},\,\sigma_{z})$ are the Pauli matrices. The subspace in which they act may be very general and depends on the given model, for example in a lattice model with two orbitals per site, $s$ and $p$, the $\bm{\sigma}$ matrices act in the orbital subspace. The above model yields topologically nontrivial phases for $M\in(0,\,2)\cup(2,\,4)$ (see Ref. [Andrei Bernevig and Hughes, 2013]).

In what follows we introduce a line-like impurity at $x=0$ described by the potential $V_{\mathrm{imp}}(x)=V\delta(x)\mathbb{I}$, with $V\to\infty$ and $\mathbb{I}$ is the $2\times 2$ identity matrix. In this limit the $T$-matrix can be written as:

[TABLE]

We compute the integral in Eq. (21), setting $B=M=1$ for the sake of simplicity. We note that the calculation can be performed for arbitrary values of $B$ and $M$. Thus, we have

[TABLE]

Plugging Eq. (22) into Eq. (21) we obtain the $T$-matrix, which in turn defines the perturbed GF given by Eq. (1). Note that the poles of the $T$-matrix obtained by taking $i\omega_{n}\rightarrow E+i\delta$, with $\delta\to+0$, are given by $E=\pm\sin k_{y}$ for $k_{y}\in\left[-\pi/2,\,\pi/2\right]$, corresponding to two chiral edge modes (cf. Ref. [König et al., 2008]). For the sake of brevity, we leave the derivation of the poles of the $T$-matrix to Appendix C.

To verify our findings, in Fig. 9 we plot the average correction to the spectral function defined in Eq. (9) as a function of $E$ and $k_{y}$. As expected, on one hand we can see the bulk states, originating from the poles of the bare Green’s function i.e., the eigenvalues of the Hamiltonian in Eq. (20) which for the values considered here, $B=M=t=1$ and $k_{x}=0$ correspond to $E=\pm 1$. More importantly we can identify also the two counter-propagating chiral edge modes of the Chern insulator crossing at $k_{y}=0$, whose dependence on $k_{y}$ is consistent with the fully analytical form above. This demonstrates the strength of our approach to recover fully analytical results for the edge states of certain models for which the unperturbed Green’s function in the real space can be obtained in an analytical closed form.

VI Conclusions

We have generalized the technique to obtain boundary modes, introduced in Ref. [Kaladzhyan and Bena, 2019], to calculate the surface Green’s functions of an arbitrary three-dimensional system, and we have applied it to calculate the surface Green’s functions for Weyl semimetals and recovered the corresponding Fermi-arc states. We have also shown that the technique in Ref. [Kaladzhyan and Bena, 2019] can be applied to other topological systems, such as topological insulators. Furthermore, we have demonstrated that it functions also for systems with more than one sublattice, and that this method can be easily employed to study boundary modes in lower dimensions. In particular, using line-like impurities, we have applied the formalism to derive the helical edge states of the Kane–Mele model and the chiral edge states of a Chern insulator. For the latter we have shown that using this formalism a full analytical form can be obtained for the $T$-matrix, and analytical expressions for the energies of the edge states can be recovered.

Acknowledgements.

V.K. would like to acknowledge the ERC Starting Grant No. 679722 and the Roland Gustafsson foundation for theoretical physics.

Appendix A Derivation of the T-matrix momentum limits of integration for a honeycomb lattice

While for a square or cubic lattice there is no subtlety concerning the integration limits, for the honeycomb lattice this is much more subtle. We will start with writing down the form of the contribution of the impurity potential to the Hamiltonian, such as it has to be written in the continuum. We will start with a row of impurities localized on A atoms, but the conclusions do not depend on the type of impurity we consider:

[TABLE]

where the integral is performed over the entire space, with

[TABLE]

where the sum is performed over a row of lattice unit cells, each unit cell is denoted by the index $j$ (see Eq. (24) in Ref. [Bena and Montambaux, 2009]) and $R_{j}=ja\sqrt{3}$ (see Fig. 9), with the lattice constant $a$ having been set to 1. The $a^{\dagger}_{j}$ and $a_{j}$ operators describe the formation and annihilation of electrons at site $j$, they no longer live in the continuum but on the lattice and as such they are defined as (see Eq. (21) in Ref. [Bena and Montambaux, 2009]):

[TABLE]

where $\int_{{\bf k}\in BZ}\equiv\int_{BZ}\frac{d^{2}k}{S_{BZ}}$ with $S_{BZ}=\frac{8\pi^{2}}{3\sqrt{3}}$. In order to use the momentum space T-matrix formalism we need to write $\delta H_{\mathrm{imp}}$ in momentum space:

[TABLE]

Since both $k$ and $k^{\prime}$ are in the first BZ the only possibilities for $n$ are $0,1$ and $-1$. It appears that implementing this constraint is quite subtle, however things get much simpler if we deform the first BZ and instead we consider a rectangle with $-\frac{2\pi}{3\sqrt{3}}<k_{y}<\frac{4\pi}{3\sqrt{3}}$ and $-2\pi/3<k_{x}<2\pi/3$. This is allowed since $a_{k+{\bf Q}_{\mu\nu}}=a_{k}$, where ${\bf Q}_{\mu\nu}$ are all the reciprocal basis vectors (see Ref. [Bena, 2009]; in the tight-binding basis considered here the same relation is valid also for the $B$ atoms Bena and Montambaux (2009)). Under this construction it is clear that the only possible solution for $n$ is $n=0$, and thus we have

[TABLE]

where $L_{k}=4\pi/3$.

Appendix B Significance of a position-averaged spectral function

In 2D, in the presence of a line impurity proportional to $\delta(x)$, the correction to the number of available electronic states at position $x$ and having momentum $k_{y}$ is given by:

[TABLE]

Note that since the spatial translational invariance along $y$ is not broken $k_{y}$ is still a good quantum number.

Averaging this over the direction perpendicular to the impurity we obtain:

[TABLE]

where $L_{k}$ and the limits of integration are $-\pi$ to $\pi$ with $L_{k}=2\pi$ for a square lattice and $-2\pi/3$ to $2\pi/3$ with $L_{k}=4\pi/3$ for a honeycomb one. Also

[TABLE]

is the correction to the perturbed spectral function in momentum space.

Appendix C $T$-matrix poles for the Chern insulator

In order to calculate the energies of the bound states, here we calculate analytically the poles of the $T$-matrix defined by Eqs. (21-22). The latter can be found from the trace of the $T$-matrix given by

[TABLE]

where we replaced $i\omega_{n}\to E+i\delta$, with $\delta\to+0$. We obtain straightforwardly the zeros of the denominator, namely, $E=\pm\sin k_{y}$. However, to make sure that the latter are poles, we need to verify that they are not zeros of the numerator. The trivial zero of the numerator is $E=0$, we discard it below. Thus, we need to analyze the zeros of the expression in the square brackets:

[TABLE]

We represent the complex numbers under the square roots in the trigonometric form, and applying the limit $\delta\to+0$ we get:

[TABLE]

where we defined $\phi_{1}(E)=\frac{\pi}{2}\operatorname{sgn}E\,\Theta\left(1-E^{2}\right),\;\phi_{2}(k_{y},E)=\frac{\pi}{2}\operatorname{sgn}E\,\Theta\left(5-4\cos k_{y}-E^{2}\right).$ Since we are searching for subgap solutions, i.e., $|E|<1$, then $E^{2}<5-4\cos k_{y}\;\forall k_{y}\in\left[-\pi,\,\pi\right]$. Thus, we have $\phi_{1}(E)+\phi_{2}(k_{y},E)=\pi\operatorname{sgn}E$ and, therefore, $e^{i\left[\phi_{1}(E)+\phi_{2}(k_{y},E)\right]}=e^{i\pi\operatorname{sgn}E}=-1$. Eq. (33) then becomes

[TABLE]

The equation above is equivalent to solving the system:

[TABLE]

When $k_{y}=0$, we get $E=\pm 1$ at the edge of the gap, therefore, $E=\pm\sin k_{y}$ from the first equation. The second equation then yields:

[TABLE]

Thus, the numerator of Eq. (31) has zeros at $E=0$, and $E=\pm\sin k_{y}$ when $|k_{y}|\in(\pi/2,\,\pi]$, and therefore, we conclude that the trace of the $T$-matrix has poles at $E=\pm\sin k_{y}$ only when $|k_{y}|\leqslant\pi/2$. This means that the edge modes exist only for $k_{y}$ lying in the interval $|k_{y}|\leqslant\pi/2$, and their dispersion is given by $E=\pm\sin k_{y}$.

Appendix D Impurity potentials for the Kane-Mele model

While working with the Kane-Mele model, we considered three different delta-function impurities which we illustrate in Fig. 10 with the associated matrix representation. We see that in order to reproduce the zigzag edge states of the Kane-Mele model we must introduce an impurity on both $A$ and $B$ sites (right panel of Fig. 10). If the impurity is localized only on the $A$ or $B$ sites, it will create both a zigzag edge state and a bearded edge state.

The formation of both zigzag and bearded edges can be recovered using our method by applying either of the potentials given in the left and middle panels of Fig. 10. Fig. 11 shows the energy spectrum obtained by exact diagonalization of the Hamiltonian on a strip with one zigzag edge and one bearded edge, along with the correction to the averaged spectral function due to a line impurity localized on either one of the two sublattices. Here we have enlarged the horizontal axis to $k_{y}\in[-2\pi/\sqrt{3};2\pi/\sqrt{3}]$ to see the edge states more clearly. We recover the same dispersion as in the main text for the zigzag edge states and obtain in addition other subgap states which are localized on the bearded edge. We need however to keep in mind the fact that in this case the two semi-infinite systems are not fully decoupled since the bulk Hamiltonian of the Kane-Mele model contains spin-flip NNN terms: this leakage effect needs to be taken into account carefully especially when we consider the spin-properties of the systems, but will not affect the spectrum below.