# Exact lattice-model calculation of boundary modes for Weyl semimetals and graphene

Vardan Kaladzhyan, Sarah Pinon, Jens H. Bardarson, Cristina Bena

► To cite this version:

Vardan Kaladzhyan, Sarah Pinon, Jens H. Bardarson, Cristina Bena. Exact lattice-model calculation of boundary modes for Weyl semimetals and graphene. New Journal of Physics, 2020, 22 (10), pp.103042. 10.1088/1367-2630/abe52. cea-03196128

HAL Id: cea-03196128

https://cea.hal.science/cea-03196128

Submitted on 12 Apr 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exact lattice-model calculation of boundary modes for Weyl semimetals and graphene

Vardan Kaladzhyan,¹,* Sarah Pinon,² Jens H. Bardarson,¹ and Cristina Bena²

¹Department of Physics, KTH Royal Institute of Technology, Stockholm, SE-106 91 Sweden

²Institut de Physique Théorique, Université Paris Saclay,
CEA CNRS, Orme des Merisiers, 91190 Gif-sur-Yvette Cedex, France

(Dated: October 21, 2020)

We provide an exact analytical technique to obtain within a lattice model the wave functions of the edge states in zigzag- and bearded-edge graphene, as well as of the Fermi-arc surface states in Weyl semimetals described by a minimal bulk model. We model the corresponding boundaries as an infinite scalar potential localized on a line, and respectively within a plane. We use the T-matrix formalism to obtain the dispersion and the spatial distribution of the corresponding boundary modes. Furthermore, to demonstrate the power of our approach, we write down the surface Green's function of the considered Weyl semimetal model, and we calculate the quasiparticle interference patterns originating from an impurity localized at the respective surface.

## I. INTRODUCTION

Systems which exhibit edge states, be they topological or not, one-, two- or three-dimensional, superconducting or normal, have come under intense scrutiny over the past years. Among such examples one can mention graphene,¹⁻³ topological insulators and topological superconductors,⁴⁻⁶ Weyl semimetals,⁷ and many others. Traditional techniques to calculate the wave functions associated with the boundary modes include numerical techniques such as tight-binding exact diagonalization,⁸⁻¹³ solving the Schrödinger equation with the corresponding boundary conditions,¹⁴⁻²¹ using the bulk-boundary correspondence²²⁻²⁶ and the extended Bloch theorem.²⁷

Here we focus on a qualitatively different approach to solving the problem of finding boundary modes, which was first introduced for Majorana bound states in Refs. [28–30], and subsequently developed in Refs. [31] and [32] for topological insulators and Andreev bound states, respectively. In this work we extend the technique to derive analytically the boundary modes of Dirac systems. The core of the technique consists of modeling the boundary as a point-, line-, or plane-like localized scalar impurity potential with an infinite amplitude. This model can be solved exactly using the T-matrix formalism and allows to obtain in a very elegant and straightforward analytical manner the energy dispersion and the wave function of the bound states, as well as the surface/edge Green's functions of the system.

In this paper we focus on deriving exact closed-form analytical expressions for the edge modes of zigzag- and bearded-edge graphene, as well as for their three-dimensional generalization — Fermi-arc surface states in Weyl semimetals described by a minimal tight-binding model.³³ Our results are consistent with previous findings obtained using different techniques, both in graphene and in Weyl semimetals, including, e.g., a recursive evaluation of the edge states in a tight-binding model,³⁴⁻³⁷ and analytical studies of the Schrödinger equation with specific boundary conditions in graphene³⁸⁻⁴⁰ and Weyl

semimetals.⁴¹,⁴² Additionally, using the analytical result for the surface Green's functions, we calculate the quasiparticle interference patterns originating from a localized impurity at the surface of the considered Weyl semimetal. We should stress that our method provides a number of advantages with respect to more traditional ones, in that it does not require any numerical algorithms, such as recursive Green's function calculation⁴³⁻⁴⁶ or exact diagonalization, and therefore, can significantly speed up other calculations that require finding the boundary modes (e.g., in transport simulations). Moreover, in certain cases, like the ones described in this work, it yields exact closed-form expressions for the surface Green's functions, the edge states and the surface states, which provides one with more insight into the physics of the problem. Also, our technique is general, being applicable to any tight-binding model on any type of lattice structure in any number of dimensions. Finally, we are not required to make any low-energy approximations, as we can employ the full lattice model for the system under consideration.

The paper is organized as follows: in Sec. II we obtain the energy dispersion and the wave functions of the boundary modes for both zigzag- and bearded-edge graphene. We derive the Fermi-arc surface states and compute quasiparticle interference patterns in a Weyl semimetal in Sec. III, and we leave the conclusions to Sec. IV.

## II. ZIGZAG- AND BEARDED-EDGE MODES IN GRAPHENE

The simplest lattice model for graphene can be written as $H_0 = \int \frac{dk}{(2\pi)^2} \Psi_k^\dagger \mathcal{H}_0(k_x, k_y) \Psi_k$, with $\Psi_k = \{\psi_k^A, \psi_k^B\}^T$, where the indices A and B refer to the corresponding sublattices, and

$$ \mathcal{H}_0(k_x, k_y) = t \begin{pmatrix} 0 & e^{-ik_x} + 2e^{i\frac{k_x}{2}} \cos \frac{k_y \sqrt{3}}{2} \\ e^{ik_x} + 2e^{-i\frac{k_x}{2}} \cos \frac{k_y \sqrt{3}}{2} & 0 \end{pmatrix}, \quad (1) $$
where we set the lattice constant to unity and $t$ is the hopping amplitude. In what follows we express all energies in units of the hopping amplitude, equivalently we set $t = 1$. The first Brillouin zone for this model is defined as a hexagon with corners located at $(k_x, k_y) = (\pm \frac{2\pi}{3}, \pm \frac{2\pi}{3\sqrt{3}})$, $(0, \pm \frac{4\pi}{3\sqrt{3}})$. The bare Matsubara Green's function for the Hamiltonian in Eq. (1) can be calculated using the standard definition:

$$ G_0(k_x, k_y, i\omega_\ell) \equiv [i\omega_\ell - \mathcal{H}_0(k_x, k_y)]^{-1} = \frac{1}{3 + 2\cos k_y\sqrt{3} + 4\cos\frac{k_y\sqrt{3}}{2}\cos\frac{3k_x}{2} - (i\omega_\ell)^2} \times \begin{pmatrix} i\omega_\ell & e^{-ik_x} + 2e^{i\frac{k_x}{2}}\cos\frac{k_y\sqrt{3}}{2} \\ e^{ik_x} + 2e^{-i\frac{k_x}{2}}\cos\frac{k_y\sqrt{3}}{2} & i\omega_\ell \end{pmatrix}. \quad (2) $$

where $\omega_\ell = \pi T (2\ell + 1)$, with $\ell \in \mathbb{Z}$, denote the fermionic Matsubara frequencies.

Our goal is to find the boundary modes for the two types of edges we are interested in: zigzag and bearded edges (shown in red in Fig. 1). In order to model an edge along the $y$ axis, we introduce an impurity potential localized solely on the atoms of sublattice A, as shown by large blue circles at $x = 0$ in Fig. 1. Such an impurity can be described by $V = U \sum_{n \in \mathbb{Z}} [\psi_{(0,n\sqrt{3})}^A]^\dagger \psi_{(0,n\sqrt{3})}^A$, with $(0, n\sqrt{3})$ being the real-space lattice points. This choice of the impurity potential in the limit of $U \to \infty$ divides the infinite graphene sheet into two halves, at $x < 0$ and $x > 0$, where the former has a zigzag edge and the latter a bearded one, both terminating on atoms of sublattice B (see the green lattice sites in Fig. 1).

To find the impurity-induced states, which evolve into boundary modes when $U$ is much larger than the other energy scales in the model, we calculate the bare Matsubara Green's function in the mixed coordinate space and momentum space representation, keeping in mind that the momentum along the $y$ direction $k_y$ is in this configuration a good quantum number:

$$ G_0(x, k_y, i\omega_\ell) \equiv \int_{-2\pi/3}^{2\pi/3} \frac{dk_x}{4\pi/3} G_0(k_x, k_y, i\omega_\ell) e^{ik_x x}. \quad (3) $$

The integration limits $\pm 2\pi/3$ and the numerical factor $4\pi/3$ were derived in Ref. [31] for a hexagonal lattice. Note also that in the expression above, $x$ can only take discrete values corresponding to the positions of the atoms in the honeycomb lattice. These values can be divided into two classes, corresponding to the atoms of the sublattices A and B,

$$ x_A = \frac{3n}{2} \text{ and } x_B = \frac{3n+1}{2}, \text{ with } n \in \mathbb{Z}. \quad (4) $$

We note that when $n < 0$ we are dealing with lattice sites on the left side of the impurity line, whereas for $n > 0$ with those on the right side.

<figure><img src="image_2.png" /><figcaption>Figure 1. Graphene lattice with an infinite-amplitude $\delta$-function impurity introduced at $x = 0$ (impurity sites are thus effectively disconnected from the lattice, as indicated by the dashed lines). Effectively we have a zigzag edge half-plane on the left side of the impurity, i.e., for $x < 0$, and a bearded edge half-plane on the right side of the impurity, i.e., for $x > 0$ (red lines).</figcaption></figure>

In order to perform the integration in Eq. (3), we first denote:

$$ X_m(x, k_y, i\omega_\ell) \equiv - \int_{-2\pi/3}^{2\pi/3} \frac{dk_x}{4\pi/3} \frac{e^{ik_x m/2} e^{ik_x x}}{D}, \quad (5) $$

where $D \equiv 3 + 2\cos k_y\sqrt{3} + 4\cos\frac{k_y\sqrt{3}}{2}\cos\frac{3k_x}{2} - (i\omega_\ell)^2$, $m \in \{0, \pm 1, \pm 2\}$. Due to the translational invariance in the $y$ direction, the integrals above are periodic functions of $k_y$, with a period given by $[0, \frac{2\pi}{\sqrt{3}}]$. The Green's function in Eq. (3) can thus be written as

$$ G_0(x, k_y, i\omega_\ell) = \begin{pmatrix} i\omega_\ell X_0 & 2X_1 \cos \frac{k_y \sqrt{3}}{2} + X_{-2} \\ 2X_{-1} \cos \frac{k_y \sqrt{3}}{2} + X_2 & i\omega_\ell X_0 \end{pmatrix}, \quad (6) $$

where for simplicity we have omitted the explicit arguments of the $X_m$ functions. The detailed calculation of the five $X_m$ integrals is presented in Appendix A.

Below, using the Green's function in Eq. (6) we compute the $T$-matrix in order to find its poles defining the energies of the edge states of graphene in the limit where the impurity potential amplitude $U$ is the largest energy scale in the system. The $T$-matrix can be found as follows$^{47,48}$ $T \equiv (I - V G_0)^{-1} V$, which in our case be-
comes:

$$ T(k_y, i\omega_\ell) \\ = \left[ \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - \begin{pmatrix} U & 0 \\ 0 & 0 \end{pmatrix} \cdot G_0(x=0, k_y, i\omega_\ell) \right]^{-1} \cdot \begin{pmatrix} U & 0 \\ 0 & 0 \end{pmatrix} \\ \xrightarrow{U \to \infty} - \begin{pmatrix} [G_0^{11}(x=0, k_y, i\omega_\ell)]^{-1} & 0 \\ 0 & 0 \end{pmatrix}, \qquad (7) $$

where $G_0^{11}$ denotes the 11 component of the corresponding Green's function. We perform the analytical continuation $i\omega_\ell \to E + i\delta$, $\delta \to +0$, and we find that the imaginary part of the trace of the T-matrix,

$$ \mathrm{Im} \operatorname{tr} T(k_y, E + i0) = \mathrm{Im} \frac{\sqrt{(3 + 2 \cos k_y \sqrt{3} - (E + i0)^2)^2 - 16 \cos^2 \frac{k_y \sqrt{3}}{2}}}{E + i0} \quad (8) $$

has a pole at $E = 0$ for all $k_y \in [0, \frac{2\pi}{\sqrt{3}}] \setminus \{\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}}\}$ (see Fig. 2). Thus, there exists a zero-energy edge state for all the possible values of $k_y$ except for the special points $\{\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}}\}$. It is known from literature$^{35-37}$ that if we had only one type of edge—bearded or zigzag—we would only recover edge modes for specific values of $k_y$. Namely, for the zigzag edge we would have a zero-energy edge state only for $k_y \in (\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}})$, whereas for the bearded edge we would have it for $k_y \in [0, \frac{2\pi}{\sqrt{3}}] \setminus [\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}}]$. The special points $\frac{2\pi}{3\sqrt{3}}$ and $\frac{4\pi}{3\sqrt{3}}$ are the points where the zero-energy state ceases to be edge states and merge with the bulk. For our particular configuration we have both zigzag and bearded edges, and thus it is not surprising that the existence of the edge states extends to the entire range of $k_y$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Imaginary part of the trace of the T-matrix, as computed in Eq. (8). At $E = 0$, we recover the edge states corresponding to both the zigzag and the bearded edges for $k_y \in (\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}})$ and $k_y \in [0, \frac{2\pi}{\sqrt{3}}] \setminus [\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}}]$, respectively.</figcaption></figure>

To find the wave functions corresponding to the edge states we use the same algorithm as for the Yu-Shiba-Rusinov states.$^{49-54}$ First, we obtain the wave function at the impurity position ($x=0$) using

$$ \left[ \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - G_0(x=0, k_y, E=0) \cdot \begin{pmatrix} U & 0 \\ 0 & 0 \end{pmatrix} \right] \Psi(x=0, k_y) = 0. \qquad (9) $$

We then use the propagation relation

$$ \Psi(x, k_y) = G_0(x, k_y, E=0) \cdot \begin{pmatrix} U & 0 \\ 0 & 0 \end{pmatrix} \Psi(x=0, k_y) \quad (10) $$

to recover the dependence of the wave function on $x$, and we get:

$$ \Psi(x, k_y) = \begin{pmatrix} 0 \\ X_0(x_B + 1) + 2 \cos \frac{k_y \sqrt{3}}{2} X_0(x_B - \frac{1}{2}) \end{pmatrix}, \qquad (11) $$

where we have set $E = 0$ and we have omitted the $k_y$ dependence in $X_0$ for the sake of brevity. The explicit definitions of $x_B$ and $X_0$ are given in Eqs. (4) and (5), and the exact expression for the latter can be found in Appendix A. Note that this result is in qualitative agreement with previous studies of zigzag edge states$^{34-38,55-57}$ in that the edge wave function is nonzero only on the atoms of the sublattice B. The wave function is thus defined only for $x = x_B \equiv \frac{3n+1}{2}$ (see Eq. (4)), and it is equal to zero for $x = x_A$.

Further simplifications can be made assuming that $k_y$ lies in one of the two previously mentioned intervals. The wave function in Eq. (11) then simplifies to:

$$ \Psi(n, k_y) = \begin{pmatrix} 0 \\ -(-2 \cos \frac{k_y \sqrt{3}}{2})^{-n-1} \end{pmatrix}, \qquad (12) $$

which yields the zigzag-edge wave functions for $n < 0$, $k_y \in (\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}})$ and the bearded-edge wave functions for $n > 0$, $k_y \in [0, \frac{2\pi}{\sqrt{3}}] \setminus [\frac{2\pi}{3\sqrt{3}}, \frac{4\pi}{3\sqrt{3}}]$. Other combinations of $n$ and $k_y$ yield non-normalizable solutions, and thus can be discarded. Note the exact equivalence to previous results obtained using recursive calculations$^{34-37}$ and specific boundary conditions.$^{17,38}$

In Fig. 3 we plot the corresponding lattice local density of states (LDOS) $\rho(x, k_y) = |\Psi(x, k_y)|^2$. For values of $k_y$ between $\frac{2\pi}{3\sqrt{3}}$ and $\frac{4\pi}{3\sqrt{3}}$ the LDOS calculated from the wave function in Eq. (11) is nonzero only for $x < 0$, corresponding to the zigzag edge state. Contrary to that, if the value of $k_y$ is chosen to be outside of the aforementioned range, the LDOS is nonzero only for $x > 0$, and therefore corresponds to the bearded edge state.

### III. WEYL SEMIMETAL: FERMI ARCS AND QUASIPARTICLE INTERFERENCE PATTERNS

In this section we turn to calculating the wave functions for the Fermi-arc states and the quasiparticle in-
<figure><img src="image_4.png" /><figcaption>Figure 3. The LDOS calculated for the wave function in Eq. (11) for four different values of $k_y$. In the left and right columns we show the zigzag and bearded edges, respectively. The values of $k_y$ on the upper two panels correspond to the cases in which the edge states are the most localized, i.e., the localization length as a function of $k_y$ reaches a minimum, whereas the lower two panels show the opposite case, where the edge states are the most delocalized. The size of the circles reflects the magnitude of the LDOS.</figcaption></figure>

terference patterns on the surface of a Weyl semimetal described by the following model:

$$ \mathcal{H}_0(\mathbf{k}) = v \sin k_x \sigma_x + v \sin k_y \sigma_y + (m - t \cos k_x - t \cos k_y - t \cos k_z) \sigma_z. \quad (13) $$

The above Hamiltonian is defined on a cubic lattice with the lattice constant set to unity. Here $\mathbf{k} = (k_x, k_y, k_z)$, the Pauli matrices are denoted by $(\sigma_x, \sigma_y, \sigma_z)$, $v$ characterizes the group velocity of the low-energy Weyl fermions, $m$ is the mass term, and $t$ is the hopping amplitude. In what follows we express all quantities with the dimensionality of energy in terms of the hopping amplitude, and thus we set $t = 1$. We assume also for simplicity that $v = 1$. For $1 < m < 3$ the Hamiltonian above describes a Weyl semimetal phase. The Weyl nodes appear along the $k_z$ axis, at the two points defined by $\cos k_z^0 = m - 2$. Note that this is a simple minimal model exhibiting only two Weyl nodes in the Brillouin zone. The bare Matsubara Green's function is defined as follows:

$$ G_0(\mathbf{k}, i\omega_\ell) = [i\omega_\ell - \mathcal{H}_0(\mathbf{k})]^{-1} = \frac{i\omega_\ell\sigma_0 + \tilde{g}\sigma_z + \sin k_x\sigma_x + \sin k_y\sigma_y}{\tilde{g}^2 + \sin^2 k_x + \sin^2 k_y - (i\omega_\ell)^2}, \quad (14) $$

where we denote $\tilde{g} \equiv \tilde{g}(\mathbf{k}) \equiv m - \cos k_x - \cos k_y - \cos k_z$.

### A. Fermi-arc surface states

In order to find the Fermi-arc states we introduce into the system a plane-like boundary at $y = 0$, emulated by a $\delta$-function impurity with the following potential $V(y) = \begin{pmatrix} U & 0 \\ 0 & U \end{pmatrix} \delta(y)$. To solve the corresponding impurity problem we compute the Matsubara Green's function in the mixed coordinate space and momentum space representation, i.e.,

$$ G_0(k_x, y, k_z, i\omega_\ell) = \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} G_0(\mathbf{k}, i\omega_\ell) e^{ik_y y}. \quad (15) $$

In the expression above $y$ is not a continuous variable since the system is defined on a lattice, therefore, it admits only integer values (both positive and negative, depending on whether the system lies at $y \ge 0$ or $y \le 0$), in other words, $y = na$, where $a = 1$ is the lattice constant, and $n \in \mathbb{Z}$. To perform the Fourier transform in Eq. (15) we define three integrals:

$$ X_s = \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{(\bullet) e^{ik_y y}}{\tilde{g}^2 + \sin^2 k_x + \sin^2 k_y - (i\omega_\ell)^2}, \quad (16) $$

where $(\bullet) = 1, \cos k_y$ and $\sin k_y$ for $s = 0, 1$ and $2$, respectively. We leave the step-by-step calculations of the integrals above to Appendix B. In terms of these integrals the Green's function can be written as

$$ G_0(k_x, y, k_z, i\omega_\ell) = i\omega_\ell X_0\sigma_0 + (gX_0 - X_1)\sigma_z + \sin k_x X_0\sigma_x + X_2\sigma_y, \quad (17) $$

where we denote $g \equiv g(k_x, k_z) \equiv m - \cos k_x - \cos k_y$. Here we have omitted the arguments of the $X_s$ functions for the sake of brevity.

In what follows we compute the $T$-matrix using the mixed real and reciprocal space representation of the bare Green's function given in Eq. (17), and assuming that the impurity potential amplitude $U$ is the largest energy scale in the system, i.e., $U \gg 1$. Thus, we get

$$ T(k_x, k_z, i\omega_\ell) = \left[ \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - U G_0(k_x, y = 0, k_z, i\omega_\ell) \right]^{-1} \xrightarrow{U \to \infty} -[G_0(k_x, y = 0, k_z, i\omega_\ell)]^{-1}. \quad (18) $$

The poles of the imaginary part of the trace of the $T$-matrix found via analytical continuation $i\omega_\ell \to E + i\delta, \delta \to +0$ define the energies of the impurity-bound states at finite values of $U$, whereas at $U \to \infty$ they yield the surface modes, i.e., the Fermi-arc surface states. We perform the analytical continuation, and we find that

$$ \mathrm{tr} T(k_x, k_z, E + i0) = \frac{2(E + i0)X_0}{(gX_0 + X_1)^2 + [\sin^2 k_x - (E + i0)^2] X_0^2} \quad (19) $$
whose imaginary part has two poles, at $E = \pm \sin k_x$, for $k_z$ lying in the interval between the Weyl nodes, i.e., $(-k_z^0, k_z^0)$, and $k_x$ such that $\cos k_x \ge m - 1 - \cos k_z$. Outside of these regions, the T-matrix does not have any poles, and therefore, the Fermi arcs cease to exist. It is worth noting that there are two solutions, $E = +\sin k_x$ and $E = -\sin k_x$, due to the fact that a plane-like impurity introduces both a 'top' surface and a 'bottom' surface, one for the Weyl semimetal in the lower half-space, $y \le 0$, and another one for that in the upper half-space, $y \ge 0$.

To study the spatial dependence of the Fermi-arc wavefunctions we use the same procedure as in the previous section.49–54 First, we find the wave function at $y = 0$ from the equation:

$$G_0(k_x, y = 0, k_z, E = \pm \sin k_x) \Psi_{\pm}(k_x, y = 0, k_z) = 0,$$

whose non-trivial solutions can be found by imposing $\det G_0(k_x, y = 0, k_z, E = \pm \sin k_x) = 0$, which in turn yields:

$$E = + \sin k_x \rightarrow \Psi_+(k_x, y = 0, k_z) \propto \begin{pmatrix} 1 \\ -1 \end{pmatrix}, \quad (20)$$

$$E = - \sin k_x \rightarrow \Psi_-(k_x, y = 0, k_z) \propto \begin{pmatrix} 1 \\ 1 \end{pmatrix}. \quad (21)$$

The solution for $E = +\sin k_x$ corresponds to the Weyl semimetal lying in the upper half-space, whereas the one for $E = +\sin k_x$ to that in the lower half-space. The $y$-dependence of the wave function can be found using $\Psi(k_x, y, k_z) = G_0(k_x, y, k_z, E) \cdot U \Psi(k_x, y = 0, k_z)$:

$$\Psi_+(k_x, y, k_z) = (gX_0 - X_1 + iX_2) \begin{pmatrix} 1 \\ 1 \end{pmatrix} \quad (22)$$

$$\Psi_-(k_x, y, k_z) = (gX_0 - X_1 - iX_2) \begin{pmatrix} 1 \\ -1 \end{pmatrix}. \quad (23)$$

The physics encoded in these two wave functions is qualitatively the same, and thus in Fig. (4) we plot the LDOS for one of the solutions at $E = 0$ (equivalently, $k_x = 0$), i.e.,

$$\rho(k_x = 0, y, k_z, E = 0) = |\Psi_+(k_x = 0, y, k_z)|^2 \quad (24)$$

as a function of $y \ge 0$, at different values of $k_z$. It is clear that the localization length of the Fermi-arc surface state changes with $k_z$, and as expected, the state becomes completely delocalized exactly at the Weyl nodes, i.e., at $k_z^0 = \pm 2\pi/3$.

## B. Quasiparticle interference pattern

To demonstrate the power of our approach for solving boundary problems, in what follows we calculate the quasiparticle interference patterns associated to the presence of a single localized impurity $V_s(x, z) =$

<figure><img src="image_5.png" /><figcaption>Figure 4. The LDOS for the solution at $E = 0$, as defined by Eq. (24), plotted as a function of $y$ and $k_z$. The mass term is set to $m = 1.5$, yielding Weyl nodes at $k_z^0 = \pm 2\pi/3$, in the vicinity of which the Fermi-arc surface state becomes more delocalized, and the corresponding LDOS becomes more bulk-like.</figcaption></figure>

On the lines of Ref. [31], the surface Green's function for the Weyl semimetal in the upper half plane ($y > 0$) is given by:

$$G_s(k_x, k_z) = \int \frac{dk_1 y dk_2 y}{(2\pi)^2} G(k_x, k_{1y}, k_{2y}, k_z) e^{i(k_{1y} - k_{2y})y} \Big|_{y=1}, \quad (25)$$

with

$$G(k_x, k_{1y}, k_{2y}, k_z) = 2\pi G_0(k_x, k_{1y}, k_z)\delta(k_{1y} - k_{2y}) + G_0(k_x, k_{1y}, k_z)T(k_x, k_z)G_0(k_x, k_{2y}, k_z), \quad (26)$$

where the unperturbed bulk Green's function and the T-matrix are given by Eqs. (14) and (18), respectively. We set $y = 1$ in the expression above because in the presence of an infinite-amplitude $\delta$-potential impurity all sites at $y = 0$ are cut out of the system, thus moving the surface of the system to $y = 1$. Above we omitted writing down the explicit energy dependence in all Green's functions. We rewrite the integrand in Eq. (25) using Eq. (26), and we obtain the surface Green's function performing both integrals over momenta in Eq. (25):

$$G_s(k_x, k_z) = G_0(k_x, y = 0, k_z) + G_0(k_x, y = 1, k_z)T(k_x, k_z)G_0(k_x, y = -1, k_z), \quad (27)$$

where the mixed coordinate space-momentum space representation of the unperturbed Green's function is given by Eq. (17).
<figure><img src="image_6.png" /><figcaption>Figure 5. Spectral function at $E = 0$, as defined by Eq. (28), plotted as a function of $k_x$ and $k_z$. We set $m = 2.5$, thus Weyl nodes appear at $k_z^0 = \pm \pi/3$. We can clearly see the Fermi arc (in this toy model, a line) connecting the two nodes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Quasiparticle interference pattern at $E = 0$, as defined by Eq. (29), plotted as a function of $k_x$ and $k_z$. The mass term is set to $m = 2.5$, yielding Weyl nodes at $k_z^0 = \pm \pi/3$. The impurity potential amplitude is set to $U_s = 1$. The strong line-like feature in the middle reflects the surface-surface scattering processes, while the small background ($\Delta\rho \approx 10^{-2}$) accounts for the surface-bulk and bulk-bulk ones.</figcaption></figure>

In what follows, we study the spectral function defined through the surface Green's function in Eq. (27) as follows:

$$A(k_x, k_z) = -\frac{1}{\pi} \operatorname{Im} \operatorname{tr} G_s(k_x, k_z). \quad (28)$$

In Fig. 5 we plot the spectral function taken at $E = 0$, and we observe a line in the momentum space, corresponding to the Fermi arc. Since we employ one of the simplest lattice models of a Weyl semimetal (i.e., obtained by stacking Chern insulators in reciprocal space), the Fermi-arc surface states at a fixed value of energy appear as lines in the momentum space. Therefore, a priori we expect a very simple quasiparticle interference pattern in this toy-model case, namely, the scattering at the surface occurs mostly between the surface states, along with some residual scattering into the bulk as well.

Below we define the quasiparticle interference patterns in the momentum space via

$$\Delta\rho(k_x, k_z) = -\frac{1}{2\pi i} \int \frac{dq_x dq_z}{(2\pi)^2} f(q_x, q_z, k_x, k_z) \quad (29)$$

where

$$f = \operatorname{tr} [G_s(q_x, q_z) T_s G_s(q_x + k_x, q_z + k_z) - G_s^*(q_x + k_x, q_z + k_z) T_s^* G_s^*(q_x, q_z)].$$

Also

$$T_s = \left[ \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - U_s \int \frac{dk_x dk_z}{(2\pi)^2} G_s(k_x, k_z) \right]^{-1} U_s,$$

$U_s$ is the amplitude of the impurity potential, and * denotes complex conjugation.

In Fig. 6 we plot the quasiparticle interference pattern for the case of a shorter Fermi arc, i.e., when the Weyl nodes lie within the interval $[-\pi/2, \pi/2]$, or in other words, when $k_z^0 < \pi/2$. Most of the weight in the plot is concentrated, as expected, in a narrow line lying from $-2k_z^0$ to $-2k_z^0$ at $k_x = 0$, which corresponds to surface-surface scattering processes. Indeed, most of the scattering occurs within the Fermi arc, and the maximum change in momentum along the $k_z$ axis is $2k_z^0$, while along the $k_x$ axis it is 0. The very small and undifferentiated background corresponds to surface-bulk and bulk-bulk scattering processes. In the case of a longer Fermi arc (e.g., $m = 1.5$), the overall change in $k_z$ can go beyond the first Brillouin zone, and thus the narrow line is expected to lie in the interval $[-\pi, \pi]$.

More complex and realistic models for Weyl semimetals exhibiting topological Fermi-arc states were treated in Ref. [58], with a stark focus on spin-resolved components of quasiparticle interference patterns and the surface Green's function derivation employing the technique from Refs. [28–32]. Most importantly, no closed-form analytical solutions for the wave functions of the Fermi-arc states are presented in Ref. [58], contrary to the present work which mostly focuses on analytical derivations of boundary-mode wave functions.

## IV. CONCLUSIONS

To summarize, we have proposed an analytical route to calculate the boundary modes of graphene and Weyl semimetals within a lattice model. For the Weyl
semimetals we have considered a minimal tight-binding model exhibiting two cones and a line-like Fermi arc, and we have also computed the quasiparticle interference patterns via a calculation of the surface Green's function. Our results are obtained by modeling the boundaries as localized infinite-amplitude impurity potentials and treating the problem exactly within the T-matrix formalism. More specifically, we have recovered an exact closed form for the wave functions of the zigzag and the bearded edge states for graphene, as well as for the Fermi-arc surface states and surface Green's function for a Weyl semimetal model considered. The results presented in this work are in agreement with previous calculations performed either by recursive tight-binding methods,³⁴⁻³⁶ or by solving the Schrödinger equation with specifically derived boundary conditions,¹⁷,⁴² or other methods.²⁵,²⁷

The technique we have employed is very general and can be applied to any lattice model no matter its complexity or its dimensionality, and it is in no way limited to only minimal tight-binding models, but can be used for *any* tight-binding model. It allows to recover energies of the boundary modes, as well as exact forms for their wave functions without requiring any numerical calcula-

tions such as, e.g., exact diagonalization, except at most a numerical integral of the Green's function to Fourier transform it from momentum space to real space; however, in the present work we have obtained all integrals in a closed analytical form.

Potential applications include, but are not limited to: using the technique in setups with pseudomagnetic fields;⁵⁹ studying disordered systems numerically by treating a single disorder realization via the generalized T-matrix or analytically by introducing a boundary-emulating impurity after averaging over disorder; considering more realistic models for studying quasiparticle interference patterns in Weyl and Dirac semimetals.⁵⁸

ACKNOWLEDGMENTS

V.K. and J.H.B. would like to acknowledge the ERC Starting Grant No. 679722. V.K. acknowledges also the Roland Gustafsson foundation for theoretical physics, and the Karl Engvers foundation, as well as Mark O. Goerbig, Clément Dutreix and Loïc Herviou for fruitful discussions.

* vardank@kth.se

1. A. K. Geim and K. S. Novoselov, *Nature Materials* **6**, 183 (2007).

2. W. A. de Heer, C. Berger, X. Wu, P. N. First, E. H. Conrad, X. Li, T. Li, M. Sprinkle, J. Hass, M. L. Sadowski, M. Potemski, and G. Martinez, *Solid State Communications* **143**, 92 (2007).

3. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, *Rev. Mod. Phys.* **81**, 109 (2009).

4. C. Kane and E. Mele, *Phys. Rev. Lett.* **95**, 146802 (2005).

5. X.-L. Qi and S.-C. Zhang, *Rev. Mod. Phys.* **83**, 1057 (2011).

6. B. Andrei Bernevig and T. Hughes, *Topological Insulators and Topological Superconductors* (Princeton Press, 2013).

7. N. P. Armitage, E. J. Mele, and A. Vishwanath, *Rev. Mod. Phys.* **90**, 015001 (2018).

8. J. C. Slater and G. F. Koster, *Phys. Rev.* **94**, 1498 (1954).

9. R. Haydock, V. Heine, and M. J. Kelly, *Journal of Physics C: Solid State Physics* **5**, 2845 (1972).

10. R. Haydock, V. Heine, and M. J. Kelly, *Journal of Physics C: Solid State Physics* **8**, 2591 (1975).

11. L. M. Falicov and F. Yndurain, *Journal of Physics C: Solid State Physics* **8**, 147 (1975).

12. U. Busch and K. A. Penson, *Phys. Rev. B* **36**, 9271 (1987).

13. K. W. Kim, R. S. K. Mong, M. Franz, and G. Refael, *Phys. Rev. B* **92**, 075110 (2015).

14. M. V. Berry and R. J. Mondragon, *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences* **412**, 53 (1987).

15. S. Davison and M. Steślicka, *Basic theory of surface states*, Monographs on the physics and chemistry of materials (Clarendon Press, 1992).

16. J. Li, A. F. Morpurgo, M. Büttiker, and I. Martin, *Physical Review B* **82** (2010), 10.1103/physrevb.82.245404.

17. P. Delplace, D. Ullmo, and G. Montambaux, *Phys. Rev. B* **84**, 195452 (2011).

18. F. Zhang, C. L. Kane, and E. J. Mele, *Phys. Rev. B* **86**, 081303 (2012).

19. E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O. Sukhachov, *Phys. Rev. B* **91**, 235138 (2015).

20. Z. A. Devizorova and V. A. Volkov, *Phys. Rev. B* **95**, 081302 (2017).

21. C. W. Duncan, P. Öhberg, and M. Valiente, *Phys. Rev. B* **97** (2018), 10.1103/PhysRevB.97.195439.

22. Y. Hatsugai, *Phys. Rev. Lett.* **71**, 3697 (1993).

23. R. S. K. Mong and V. Shivamoggi, *Phys. Rev. B* **83**, 125109 (2011).

24. T. Fukui, K. Shiozaki, T. Fujiwara, and S. Fujimoto, *Journal of the Physical Society of Japan* **81**, 114602 (2012).

25. F. K. Kunst, M. Trescher, and E. J. Bergholtz, *Phys. Rev. B* **96**, 085443 (2017).

26. J.-W. Rhim, J. H. Bardarson, and R.-J. Slager, *Phys. Rev. B* **97**, 115143 (2018).

27. F. K. Kunst, G. van Miert, and E. J. Bergholtz, *Phys. Rev. B* **99**, 085427 (2019).

28. A. Zazunov, R. Egger, and A. Levy Yeyati, *Phys. Rev. B* **94**, 014502 (2016).

29. A. Zazunov, R. Egger, M. Alvarado, and A. L. Yeyati, *Phys. Rev. B* **96**, 024516 (2017).

30. V. Kaladzhyan and C. Bena, *Phys. Rev. B* **100**, 081106 (2019).

31. S. Pinon, V. Kaladzhyan, and C. Bena, *Phys. Rev. B* **101**, 115405 (2020).

32. S. Pinon, V. Kaladzhyan, and C. Bena, *Phys. Rev. B* **101**, 205136 (2020).
33 A. M. Turner and A. Vishwanath, in *Contemporary Concepts of Condensed Matter Science* (Elsevier, 2013) pp. 293–324.

34 M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, *Journal of the Physical Society of Japan* **65**, 1920 (1996).

35 K. Wakabayashi, S. Okada, R. Tomita, S. Fujimoto, and Y. Natsume, *Journal of the Physical Society of Japan* **79**, 034706 (2010).

36 K. Wakabayashi, K. ichi Sasaki, T. Nakanishi, and T. Enoki, *Science and Technology of Advanced Materials* **11**, 054504 (2010).

37 M. Bellec, U. Kuhl, G. Montambaux, and F. Mortessagne, *New Journal of Physics* **16**, 113023 (2014).

38 A. R. Akhmerov and C. W. J. Beenakker, *Phys. Rev. B* **77**, 085423 (2008).

39 B. Messias de Resende, F. C. de Lima, R. H. Miwa, E. Vernek, and G. J. Ferreira, *Phys. Rev. B* **96**, 161113 (2017).

40 A. L. Araújo, R. P. Maciel, R. G. F. Dornelas, D. Varjas, and G. J. Ferreira, *Phys. Rev. B* **100**, 205111 (2019).

41 R. Okugawa and S. Murakami, *Phys. Rev. B* **89**, 235315 (2014).

42 K. Hashimoto, T. Kimura, and X. Wu, *Progress of Theoretical and Experimental Physics* **2017** (2017), 10.1093/ptep/ptx053.

43 M. P. Lopez Sancho, J. M. Lopez Sancho, and J. Rubio, *Journal of Physics F: Metal Physics* **15** (1985), 10.1088/0305-4608/15/4/009.

44 M. P. L. Sancho, J. M. L. Sancho, and J. Rubio, *Journal of Physics F: Metal Physics* **14**, 1205 (1984).

45 A. Umerski, *Phys. Rev. B* **55**, 5266 (1997).

46 Y. Peng, Y. Bao, and F. von Oppen, *Phys. Rev. B* **95**, 235143 (2017).

47 G. D. Mahan, *Many-Particle Physics* (Springer US, 2000).

48 A. V. Balatsky, I. Vekhter, and J.-X. Zhu, *Rev. Mod. Phys.* **78**, 373 (2006).

49 L. Yu, *Acta Physica Sinica* **21**, 75 (1965).

50 H. Shiba, *Progress of Theoretical Physics* **40**, 435 (1968).

51 A. I. Rusinov, *Soviet Journal of Experimental and Theoretical Physics Letters* **9**, 85 (1969).

52 A. Sakurai, *Progress of Theoretical Physics* **44**, 1472 (1970).

53 F. Pientka, L. I. Glazman, and F. von Oppen, *Phys. Rev. B* **88**, 155420 (2013).

54 V. Kaladzhyan, C. Bena, and P. Simon, *Journal of Physics: Condensed Matter* **28**, 485701 (2016).

55 K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dresselhaus, *Phys. Rev. B* **54**, 17954 (1996).

56 L. Brey and H. A. Fertig, *Phys. Rev. B* **73**, 235411 (2006).

57 W. Yao, S. A. Yang, and Q. Niu, *Phys. Rev. Lett.* **102**, 096801 (2009).

58 S. Pinon, V. Kaladzhyan, and C. Bena, *Phys. Rev. B* **102**, 165117 (2020).

59 V. Peri, T. Dubček, A. Valenti, R. Ilan, and S. D. Huber, *Phys. Rev. B* **101**, 235117 (2020).

## Appendix A: Derivation of the edge states for zigzag- and bearded-edge graphene

In this Appendix we calculate the integrals in Eq. (5) for $m \in \{-2, -1, 0, 1, 2\}$:

$$X_m(x, k_y, i\omega_\ell) = - \int_{-2\pi/3}^{2\pi/3} \frac{dk_x}{4\pi/3} \frac{e^{ik_x m/2} e^{ik_x x}}{3 + 2\cos k_y \sqrt{3} + 4\cos \frac{k_y \sqrt{3}}{2} \cos \frac{3k_x}{2} - (i\omega_\ell)^2}, \quad (A1)$$

Translational invariance in the $y$ direction implies that the integral above is a periodic function of $k_y$, and we choose one period to be $k_y \in [0, \frac{2\pi}{\sqrt{3}}]$. The case of $k_y = \frac{\pi}{\sqrt{3}}$ should be considered separately since the integral simplifies significantly due to the fact that the denominator does not depend on $k_x$ anymore. Note also that the integral above is defined solely on the graphene lattice shown in Fig. 1, and thus $x = x_A = \frac{3n}{2}$ or $x = x_B = \frac{3n+1}{2}$, with $n \in \mathbb{Z}$.

We start with the case of $k_y = \frac{\pi}{\sqrt{3}}$, where we have:

$$X_m(x, k_y, i\omega_\ell) = - \frac{1}{1 - (i\omega_\ell)^2} \int_{-2\pi/3}^{2\pi/3} \frac{dk_x}{4\pi/3} e^{ik_x m/2} e^{ik_x x} = - \frac{1}{1 - (i\omega_\ell)^2} \frac{\sin \pi \left( \frac{2}{3}x + \frac{m}{3} \right)}{\pi \left( \frac{2}{3}x + \frac{m}{3} \right)} \quad (A2)$$

When additionally $x = -m/2$, we get: $X_0(0, \frac{\pi}{\sqrt{3}}, i\omega_\ell) = -\frac{1}{1-(i\omega_\ell)^2}$.

In what follows we turn to the case of $k_y \neq \frac{\pi}{\sqrt{3}}$. In this case the integral in Eq. (A1) can be rewritten as a contour integral in the complex plane, using the substitution $z = e^{\frac{3}{2}k_x}$:

$$X_m(x, k_y, i\omega_\ell) = - \frac{1}{2\cos\frac{k_y\sqrt{3}}{2}} \oint_{|z|=1} dz \frac{z^{\frac{2}{3}(x+\frac{m}{2})}}{z^2 + fz + 1}, \quad (A3)$$
where the circle $|z| = 1$ is oriented counter-clockwise, and we defined

$$f \equiv f(k_y, i\omega_\ell) = \frac{3 + 2 \cos k_y \sqrt{3} - (i\omega_\ell)^2}{2 \cos \frac{k_y \sqrt{3}}{2}}. \qquad (A4)$$

For simplicity we rewrite this integral for $x$ lying on sublattices A and B separately:

$$\text{sublattice A: } X_m^A(x = \frac{3n}{2}, k_y, i\omega_\ell) = -\frac{1}{2\cos\frac{k_y\sqrt{3}}{2}} \frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^{n+m/3}}{z^2 + fz + 1}, \qquad (A5)$$

$$\text{sublattice B: } X_m^B(x = \frac{3n+1}{2}, k_y, i\omega_\ell) = -\frac{1}{2\cos\frac{k_y\sqrt{3}}{2}} \frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^{n+(m+1)/3}}{z^2 + fz + 1}, \qquad (A6)$$

Before we proceed with the calculation, several important simplifications are worth pointing out:

* It is easy to notice that for sublattice A we have $X_{-2}^A(x) = X_2^A(-x)$, $X_{-1}^A(x) = X_1^A(-x)$, whereas for sublattice B $X_{-2}^B(x) = X_1^A(\frac{1}{2} - x)$, $X_{-1}^B(x) = X_0^A(\frac{1}{2} - x)$, $X_0^B(x) = X_1^A(x - \frac{1}{2})$, $X_1^B(x) = X_2^A(x - \frac{1}{2})$, and $X_2^B(x) = X_0^A(x + 1)$. Above we omitted the dependence on $k_y$ and $i\omega_\ell$. Therefore, it is sufficient to compute only the three integrals in Eq. (A5), for $m = 0$, $m = 1$ and $m = 2$, respectively, and for $x = x_A$ only, since all integrals on sublattice B can be expressed in terms of those on sublattice A.

* From the definition in Eq. (A5) it is clear that the case of $m = 0$ and the case of $m \neq 0$ (i.e., $m = 1$ and $m = 2$) should be considered separately. The reason for it is that when $m = 0$ there are no functions in the integrand requiring to make a branch cut in the complex plane, thus significantly simplifying the integration.

* Another important remark can be made about the roots of the quadratic polynomial in the denominator of the integrand in Eq. (A5). A simple analysis of the roots given by

$$z_{\pm} = \frac{1}{2} \left[ -f \pm \sqrt{f^2 - 4} \right] \qquad (A7)$$

shows that, first, both roots are always real for the values of the parameters in consideration. Second, one of the roots always belongs to the interval $(-1, 1)$, namely,

$$\text{if } 0 \le k_y < \frac{\pi}{\sqrt{3}}, \quad \text{then } z_+ \in (-1, 0), z_- < -1, \qquad (A8)$$

$$\text{if } \frac{\pi}{\sqrt{3}} < k_y \le \frac{2\pi}{\sqrt{3}}, \quad \text{then } z_- \in (0, +1), z_+ > +1. \qquad (A9)$$

This means that regardless of the value of $k_y$ one of the poles of the integrand always lies inside the unit circle in the complex plane, and has to be taken into account while applying the residue theorem. Additionally, Vieta's formula dictates $z_+z_- = 1$.

## 1. The sub-case of $m = 0$

Using the roots in Eq. (A7) we rewrite the integral in Eq. (A5) in the following form (omitting the prefactor):

$$\frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^n}{(z - z_+)(z - z_-)} = * \qquad (A10)$$

It is easy to prove that this integral is unchanged if one changes $n$ to $-n$, therefore, here we compute it only for $n > 0$. In that case one of the poles lies inside the unit circle, whereas the other one – outside. We start with the case in which $z_+$ lies inside the circle ($0 \le k_y < \pi/\sqrt{3}$), and using the residue theorem, we get:

$$* = \operatorname{res}_{z=z_+} \frac{z^n}{(z - z_+)(z - z_-)} = - \frac{z_+^{|n|}}{z_- - z_+} \qquad (A11)$$
<figure><img src="image_8.png" /><figcaption>Figure 7. Auxiliary contour $\Gamma$ for complex plane integration for $0 \le k_y < \frac{\pi}{\sqrt{3}}$ (left panel) and $\frac{\pi}{\sqrt{3}} < |k_y| \le \frac{2\pi}{\sqrt{3}}$ (right panel).</figcaption></figure>

To get the expression for the case where $z = z_-$ is inside the unit circle, and $z = z_+$ is outside $(\pi/\sqrt{3} < k_y \le 2\pi/\sqrt{3})$, we just need to exchange $z_+ \to z_-$ and vice versa in expression above. Taking into account the symmetry with respect to flipping the sign of $n$, we get the final expression:

$$X_0^A \left( x = \frac{3n}{2}, k_y, i\omega_\ell \right) = - \frac{1}{2\cos \frac{k_y \sqrt{3}}{2}} \begin{cases} - \frac{z_+^{|n|}}{z_- - z_+}, & 0 \le k_y < \frac{\pi}{\sqrt{3}} \\ + \frac{z_-^{|n|}}{z_- - z_+}, & \frac{\pi}{\sqrt{3}} < k_y \le \frac{2\pi}{\sqrt{3}} \end{cases}. \quad (A12)$$

## 2. The sub-cases of $m = 1$ and $m = 2$

In this subsection we calculate the integral

$$\frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = * \quad (A13)$$

for $m = 1$ or $m = 2$. In this case we have to introduce a branch cut on the non-positive part of the real axis, i.e., for $z \in (-\infty, 0]$. In order to start the calculation, we first define auxiliary contours in the complex plane, $\Gamma = C_1 \cup \gamma_+ \cup \gamma_-$ (see left and right panels of Fig. A 2). The unit circle is denoted $C_1$, and $\gamma_\pm$ are the right and left banks of the branch cut, parametrized by $t \pm i0$ with $t \in [-1, 0]$, correspondingly. Note, that to be entirely rigorous we should have added also a circle of infinitesimal radius around the point $z = 0$, to ensure that we treat correctly the divergence at that point when $n < -1$, however, it appears that the method of calculation we apply takes care of that problem automatically. For $0 \le k_y < \pi/\sqrt{3}$ and $\pi/\sqrt{3} < k_y \le 2\pi/\sqrt{3}$ we choose the left and the right panels of Fig. A 2, respectively. The difference between these panels is the position of the pole inside the unit circle: on the left panel, it falls into the branch cut, whereas on the right panel it lies within $(0, 1)$. We note that in the definition of the auxiliary contour $\Gamma$ we neglect the tiny line between $\gamma_+$ and $\gamma_-$, i.e., $[-i0, +i0]$, since it never contributes to the value of the integral.

We start by considering the case of $0 \le k_y < \pi/\sqrt{3}$ and we write down the residue theorem for the contour $\Gamma = C_1 \cup \gamma_+ \cup \gamma_-$ defined in the left panel of Fig. A 2. Since there are no poles inside the contour $\Gamma$ we have:

$$\frac{1}{2\pi i} \oint_{\Gamma} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = 0 = \frac{1}{2\pi i} \left[ \int_{C_1} + \int_{\gamma_+} + \int_{\gamma_-} \right] dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} \quad (A14)$$
Therefore, we can express the sought-for line integral along $C_1 \equiv \{\forall z : |z| = 1\}$ as

$$ * = -\frac{1}{2\pi i} \left[ \int_{\gamma_+} + \int_{\gamma_-} \right] dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = ** $$

In order to compute the integrals along $\gamma_+$ and $\gamma_-$ we use a parametrisation $z = t + i0, t \in [-1, 0]$ and $z = t - i0, t \in [0, -1]$, respectively. We start with $\gamma_+$ and we get:

$$ \begin{aligned} \int_{\gamma_+} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} &= \int_{-1}^{0} dt \frac{(t+i0)^{n+m/3}}{(t+i0-z_+)(t-z_-)} = \int_{-1}^{0} dt \frac{e^{(n+m/3)\ln(t+i0)}}{(t+i0-z_+)(t-z_-)} = \int_{-1}^{0} dt \frac{e^{(n+m/3)[\ln|t|+\pi i]}}{(t+i0-z_+)(t-z_-)} = \\ &= e^{i\pi(n+m/3)} \int_{-1}^{0} dt \frac{|t|^{n+m/3}}{(t+i0-z_+)(t-z_-)} = e^{i\pi(n+m/3)} \int_{0}^{1} dw \frac{w^{n+m/3}}{(w+z_+-i0)(w+z_-)} = \\ &= e^{i\pi(n+m/3)} \left[ \mathcal{P} \int_{0}^{1} dw \frac{w^{n+m/3}}{(w+z_+)(w+z_-)} + i\pi \int_{0}^{1} dw \frac{w^{n+m/3}}{w+z_-} \delta(w+z_+) \right] \end{aligned} $$

Similarly, we get:

$$ \int_{\gamma_{-}} dz \frac{z^{n+m/3}}{(z-z_{+})(z-z_{-})} = -e^{-i\pi(n+m/3)} \left[ \mathcal{P} \int_{0}^{1} dw \frac{w^{n+m/3}}{(w+z_{+})(w+z_{-})} - i\pi \int_{0}^{1} dw \frac{w^{n+m/3}}{w+z_{-}} \delta(w+z_{+}) \right] $$

Therefore, combining two integrals from left and right banks, we obtain:

$$ ** = \frac{(-1)^{n+1}}{\pi} \left( I_P \sin \frac{\pi m}{3} + \pi I_\delta \cos \frac{\pi m}{3} \right) = ***, $$

where we defined

$$ I_P = \mathcal{P} \int_0^1 dw \frac{w^{n+m/3}}{(w+z_+)(w+z_-)} = \frac{(-z_-)^{n+\frac{m}{3}} B[-\frac{1}{z_-}, n+\frac{m}{3}+1, 0] - (-z_+)^{n+\frac{m}{3}} (B[-z_+, -n-\frac{m}{3}, 0] + \pi \operatorname{ctg} \frac{\pi m}{3})}{z_- - z_+} \quad (A15) $$

$$ I_{\delta} = \int_{0}^{1} dw \frac{w^{n+m/3}}{w+z_{-}} \delta(w+z_{+}) = \frac{(-z_{+})^{n+m/3}}{z_{-}-z_{+}} \quad (A16) $$

Substituting $I_P$ and $I_\delta$ into the equation above, we obtain after simplifications:

$$ *** = \frac{(-1)^{n+1}}{\pi} \sin \frac{\pi m}{3} \frac{(-z_{-})^{n+\frac{m}{3}} B[-\frac{1}{z_{-}}, n+\frac{m}{3}+1, 0] - (-z_{+})^{n+\frac{m}{3}} B[-z_{+}, -n-\frac{m}{3}, 0]}{z_{-} - z_{+}}. $$

Thus we get:

$$ \frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = \frac{(-1)^{n+1}}{\pi} \sin \frac{\pi m}{3} \frac{(-z_-)^{n+\frac{m}{3}} B[-\frac{1}{z_-}, n+\frac{m}{3}+1, 0] - (-z_+)^{n+\frac{m}{3}} B[-z_+, -n-\frac{m}{3}, 0]}{z_- - z_+}, $$

where $B[z, \alpha, \beta]$ is the incomplete Beta-function defined as:

$$ B[z, \alpha, \beta] \equiv \int_{0}^{z} t^{\alpha-1} (1-t)^{\beta-1} dt. \quad (A17) $$
In what follows we consider the remaining case of $\frac{\pi}{\sqrt{3}} < k_y \le \frac{2\pi}{\sqrt{3}}$. First, we rewrite the integral as follows:

$$ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = \frac{1}{z_- - z_+} \left[ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z-z_-} - \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z-z_+} \right] \quad (A18) $$

The first integral in the sum can be expressed in terms of the second integral (with a parameter change) by means of a variable change $z = 1/w$. Note that such a variable change inverts the orientation of the integration contour, thus multiplying the result by $-1$:

$$ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z-z_{-}} = \frac{1}{2\pi i} \oint_{C_1} dw \frac{w^{-n-m/3}}{1/w - z_{-}} = -\frac{1}{z_{-}} \frac{1}{2\pi i} \oint_{C_1} dw \frac{w^{-n-m/3-1}}{w - 1/z_{-}} = -\frac{1}{z_{-}} \frac{1}{2\pi i} \oint_{C_1} dw \frac{w^{-n-m/3-1}}{w - z_{+}}. \quad (A19) $$

From the expression above we see that to obtain the first integral from the second one we need to replace $n \to -n-2$, $m \to 3-m$ and multiply the result by $-1/z_{-}$. Next, we compute the second integral

$$ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z - z_{+}} \qquad (A20) $$

by writing down the residue theorem for the contour $\Gamma = C_1 \cup \gamma_+ \cup \gamma_-$ defined in the panel panel of Fig. A 2 for the second integral. Since there are no poles inside the contour $\Gamma$ we have:

$$ \frac{1}{2\pi i} \oint_{\Gamma} dz \frac{z^{n+m/3}}{z - z_{+}} = 0 = \frac{1}{2\pi i} \left[ \int_{C_1} + \int_{\gamma_{+}} + \int_{\gamma_{-}} \right] dz \frac{z^{n+m/3}}{z - z_{+}}. \quad (A21) $$

Therefore, we can express the sought-for line integral along $C_1 \equiv \{\forall z : |z| = 1\}$ as

$$ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z - z_{+}} = -\frac{1}{2\pi i} \left[ \int_{\gamma_{+}} + \int_{\gamma_{-}} \right] dz \frac{z^{n+m/3}}{z - z_{+}} = * $$

In order to compute the integrals along $\gamma_+$ and $\gamma_-$ we use a parametrisation $z = t + i0, t \in [-1, 0]$ and $z = t - i0, t \in [0, -1]$, respectively:

$$ \int_{\gamma_{\pm}} dz \frac{z^{n+m/3}}{z - z_{+}} = \pm \int_{-1}^{0} dt \frac{(t \pm i0)^{n+m/3}}{t - z_{+}} = \pm e^{\pm i\pi(n+m/3)} \int_{0}^{1} dw \frac{w^{n+m/3}}{w + z_{+}} \quad (A22) $$

Thus, we get:

$$ * = -\frac{1}{2\pi i} \left[ \int_{\gamma_+} + \int_{\gamma_-} \right] dz \frac{z^{n+m/3}}{z - z_+} = -\frac{(-1)^n}{2\pi i} \left[ e^{i \frac{\pi m}{3}} - e^{-i \frac{\pi m}{3}} \right] \int_0^1 dw \frac{w^{n+m/3}}{w + z_+} = -\frac{1}{\pi} \sin \frac{\pi m}{3} e^{-i \frac{\pi m}{3}} z_+^{n+\frac{m}{3}} B \left[ -\frac{1}{z_+}, n+\frac{m}{3}+1, 0 \right] $$

Using the parameter substitution introduced above, we get the first integral:

$$ \frac{1}{2\pi i} \oint_{C_1} dz \frac{z^{n+m/3}}{z - z_{-}} = \frac{1}{\pi} \sin \frac{\pi(3-m)}{3} e^{-i \frac{\pi(3-m)}{3}} z_{-}^{n+\frac{m}{3}} B \left[ -z_{-}, -n-\frac{m}{3}, 0 \right] = -\frac{1}{\pi} \sin \frac{\pi m}{3} e^{i \frac{\pi m}{3}} z_{-}^{n+\frac{m}{3}} B \left[ -z_{-}, -n-\frac{m}{3}, 0 \right] \quad (A23) $$

Finally:

$$ \frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^{n+m/3}}{(z-z_+)(z-z_-)} = \frac{1}{\pi} \sin \frac{\pi m}{3} \frac{e^{-i \frac{\pi m}{3}} z_+^{n+\frac{m}{3}} B \left[ -\frac{1}{z_+}, n+\frac{m}{3}+1, 0 \right] - e^{i \frac{\pi m}{3}} z_-^{n+\frac{m}{3}} B \left[ -z_-, -n-\frac{m}{3}, 0 \right]}{z_- - z_+} \quad (A24) $$
Combining the results for different ranges of $k_y$, we present the final result:

$$X_m^A(x = \frac{3n}{2}, k_y, i\omega_\ell) = -\frac{1}{2\cos\frac{k_y\sqrt{3}}{2}} \begin{cases} \frac{(-1)^{n+1}}{\pi} \sin\frac{\pi m}{3} \frac{(-z_-)^{n+\frac{m}{3}} B[-\frac{1}{z_-}, n+\frac{m}{3}+1, 0] - (-z_+)^{n+\frac{m}{3}} B[-z_+, -n-\frac{m}{3}, 0]}{z_- - z_+} & 0 \le k_y < \frac{\pi}{\sqrt{3}}, \\ \frac{1}{\pi} \sin\frac{\pi m}{3} \frac{e^{-i \frac{\pi m}{3}} z_+^{n+\frac{m}{3}} B[-\frac{1}{z_+}, n+\frac{m}{3}+1, 0] - e^{i \frac{\pi m}{3}} z_-^{n+\frac{m}{3}} B[-z_-, -n-\frac{m}{3}, 0]}{z_- - z_+} & \frac{\pi}{\sqrt{3}} < k_y \le \frac{2\pi}{\sqrt{3}} \end{cases}$$

where $m=1$ or $m=2$.

## Appendix B: Derivation of the Fermi-arc states for a Weyl semimetal

In this Appendix we calculate the three integrals defined in Eq. (16):

$$X_s = - \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{(\bullet) e^{ik_y y}}{\tilde{g}^2 + \sin^2 k_x + \sin^2 k_y - (i\omega_\ell)^2}, \quad (B1)$$

where $(\bullet) = 1, \cos k_y$ and $\sin k_y$ for $s = 0, 1$ and $2$, respectively. First, it is easy to see that integrals $X_1$ and $X_2$ can be expressed in terms of the integral $X_0$ in the following way:

$$X_1(k_x, y, k_z, i\omega_\ell) = -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{\cos k_y e^{ik_y y}}{D} = \frac{1}{2} \left[ -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y(y+1)}}{D} - \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y(y-1)}}{D} \right] = \\ = \frac{1}{2} [X_0(k_x, y+1, k_z, i\omega_\ell) + X_0(k_x, y-1, k_z, i\omega_\ell)] \quad (B2)$$

$$X_2(k_x, y, k_z, i\omega_\ell) = -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{\sin k_y e^{ik_y y}}{D} = \frac{1}{2i} \left[ -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y(y+1)}}{D} + \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y(y-1)}}{D} \right] = \\ = \frac{1}{2i} [X_0(k_x, y+1, k_z, i\omega_\ell) - X_0(k_x, y-1, k_z, i\omega_\ell)], \quad (B3)$$

where for the sake of brevity we denoted $D = \tilde{g}^2 + \sin^2 k_x + \sin^2 k_y - (i\omega_\ell)^2$, $\tilde{g} = m - \cos k_x - \cos k_y \cos k_z$. Therefore, we need to compute only the integral $X_0$:

$$X_0(k_x, y, k_z, i\omega_\ell) = -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y y}}{m^2 + \omega_\ell^2 + (\cos k_x + \cos k_y + \cos k_z)(-2m + \cos k_x + \cos k_y + \cos k_z) + \sin^2 k_x + \sin^2 k_y} = \\ = -\int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y y}}{2(-m + \cos k_x + \cos k_z)\cos k_y + 1 + m^2 + \omega_\ell^2 + (\cos k_x + \cos k_z)(-2m + \cos k_x + \cos k_z) + \sin^2 k_x} = *$$

If $-m + \cos k_x + \cos k_z = 0$, then we have:

$$* = - \frac{1}{1 + \sin^2 k_x - (i\omega_\ell)^2} \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} e^{ik_y y} = - \frac{1}{1 + \sin^2 k_x - (i\omega_\ell)^2} \delta_{y,0} \quad (B4)$$

If $-m + \cos k_x + \cos k_z \neq 0$, then

$$\begin{aligned} * &= \frac{1}{2(m - \cos k_x - \cos k_z)} \int_{-\pi}^{\pi} \frac{dk_y}{2\pi} \frac{e^{ik_y y}}{\cos k_y - f(k_x, k_z)} = \frac{1}{2(m - \cos k_x - \cos k_z)} \frac{1}{2\pi i} \oint_{|z|=1} \frac{dz}{z} \frac{z^y}{\frac{1}{2}(z + z^{-1}) - f(k_x, k_z)} = \\ &= \frac{1}{m - \cos k_x - \cos k_z} \frac{1}{2\pi i} \oint_{|z|=1} dz \frac{z^y}{z^2 - 2f(k_x, k_z)z + 1} = ** \end{aligned}$$
Above we introduced

$$f(k_x, k_z) \equiv \frac{1 + g^2 + \sin^2 k_x - (i\omega_\ell)^2}{2g}, \quad g(k_x, k_z) \equiv m - \cos k_x - \cos k_z.$$

Note, that since $m \in (1, 3)$, $g(k_x, k_z) \in (-1, 1)$ for all values of $k_x$ and $k_z$. In what follows we assume that we compute the Fermi-arc states for the half-space above the impurity plane, i.e., for $y \ge 0$. The calculation for $y < 0$ is not needed, since the integral is symmetric with respect to changing $y \to -y$. In order to perform the integration above, we analyze the roots of the denominator in the complex plane, as a function of $k_x, k_z$ and $m$:

$$z_{\pm} = f \pm \sqrt{f^2 - 1} \qquad (B5)$$

It is easy to show that

$$z_+ \in (-1, 0), z_- < -1 \quad \text{when } g \in (-1, 0),$$

$$z_- \in (0, +1), z_+ > +1 \quad \text{when } g \in (0, +1).$$

Therefore, for the integral above we get:

$$** = \frac{1}{g} \begin{cases} \frac{z_{+}^{\vert y \vert}}{z_{+} - z_{-}} & \text{for } g \in (-1, 0) \\ \frac{z_{-}^{\vert y \vert}}{z_{-} - z_{+}} & \text{for } g \in (0, +1) \end{cases} \qquad (B6)$$