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

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 [1–3], topological insulators and topological superconductors [4–6], Weyl semimetals [7], and many others. Traditional techniques to calculate the wave functions associated with the boundary modes include numerical techniques such as tight-binding exact diagonalization [8–13], solving the Schrödinger equation with the corresponding boundary conditions [14–21], using the bulk-boundary correspondence [22–26] and the extended Bloch theorem [27].

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 references [28–30], and subsequently developed in references [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 [33]. 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 [34–37], and analytical studies of the Schrödinger equation with specific boundary conditions in graphene [38–40] and Weyl semimetals [41, 42]. 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 [43–46] 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 section 2 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 section 3, and we leave the conclusions to section 4.

The simplest lattice model for graphene can be written as ${H}_{0}=\int \frac{\mathrm{d}\boldsymbol{k}}{{\left(2\pi \right)}^{2}}{{\Psi}}_{\boldsymbol{k}}^{{\dagger}}{\mathcal{H}}_{0}\left({k}_{x},{k}_{y}\right){{\Psi}}_{\boldsymbol{k}}$, with ${{\Psi}}_{\boldsymbol{k}}\equiv {\left\{{\psi }_{\boldsymbol{k}}^{A},\enspace {\psi }_{\boldsymbol{k}}^{B}\right\}}^{\text{T}}$, where the indices A and B refer to the corresponding sublattices, and

$$\begin{equation}{\mathcal{H}}_{0}\left({k}_{x},{k}_{y}\right)=t\left(\begin{matrix}\hfill 0\hfill & \hfill {\text{e}}^{-\text{i}{k}_{x}}\enspace +\enspace 2{\text{e}}^{\text{i}\frac{{k}_{x}}{2}}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\hfill \\ \hfill {\text{e}}^{\text{i}{k}_{x}}\enspace +\enspace 2{\text{e}}^{-\text{i}\frac{{k}_{x}}{2}}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation} \tag{ 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 $\left({k}_{x},{k}_{y}\right)=\left({\pm}\frac{2\pi }{3},{\pm}\frac{2\pi }{3\sqrt{3}}\right),\left(0,{\pm}\frac{4\pi }{3\sqrt{3}}\right).$ The bare Matsubara Green's function for the Hamiltonian in equation (1) can be calculated using the standard definition:

$$\begin{align}\hfill {G}_{0}\left({k}_{x},{k}_{y},i{\omega }_{\ell }\right)& \equiv {\left[i{\omega }_{\ell }-{\mathcal{H}}_{0}\left({k}_{x},{k}_{y}\right)\right]}^{-1}\hfill \\ \hfill & =-\frac{1}{3+2\enspace \mathrm{cos}\enspace {k}_{y}\sqrt{3}+4\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\enspace \mathrm{cos}\enspace \frac{3{k}_{x}}{2}-{\left(i{\omega }_{\ell }\right)}^{2}}\hfill \\ \hfill & \quad {\times}\left(\begin{matrix}\hfill i{\omega }_{\ell }\hfill & \hfill {\text{e}}^{-\text{i}{k}_{x}}+2{\text{e}}^{\text{i}\frac{{k}_{x}}{2}}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\hfill \\ \hfill {\text{e}}^{\text{i}{k}_{x}}+2{\text{e}}^{-\text{i}\frac{{k}_{x}}{2}}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\hfill & \hfill i{\omega }_{\ell }\hfill \end{matrix}\right).\hfill \end{align} \tag{ 2 }$$

where ${\omega }_{\ell }=\pi T\left(2\ell +1\right)$, 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 figure 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 figure 1. Such an impurity can be described by $V=U\sum _{n\in \mathbb{Z}}{\left[{\psi }_{\left(0,n\sqrt{3}\right)}^{A}\right]}^{{\dagger}}{\psi }_{\left(0,n\sqrt{3}\right)}^{A}$, with $\left(0,n\sqrt{3}\right)$ being the real-space lattice points. This choice of the impurity potential in the limit of U → ∞ 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 figure 1).

Zoom In Zoom Out Reset image size

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:

$$\begin{equation}{G}_{0}\left(x,{k}_{y},i{\omega }_{\ell }\right)\equiv \underset{-2\pi /3}{\overset{2\pi /3}{\int }}\frac{\mathrm{d}{k}_{x}}{4\pi /3}{G}_{0}\left({k}_{x},{k}_{y},i{\omega }_{\ell }\right){\text{e}}^{\text{i}{k}_{x}x}.\end{equation} \tag{ 3 }$$

The integration limits ±2π/3 and the numerical factor 4π/3 were derived in reference [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,

$$\begin{equation}{x}_{A}=\frac{3n}{2}\quad \text{and}\quad {x}_{B}=\frac{3n+1}{2},\quad \text{with}\enspace n\in \mathbb{Z}.\end{equation} \tag{ 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.

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

$$\begin{equation}{X}_{m}\left(x,{k}_{y},i{\omega }_{\ell }\right)\equiv -\enspace \underset{-2\pi /3}{\overset{2\pi /3}{\int }}\enspace \frac{\mathrm{d}{k}_{x}}{4\pi /3}\frac{{\text{e}}^{\text{i}{k}_{x}m/2}\enspace {\text{e}}^{\text{i}{k}_{x}x}}{D},\end{equation} \tag{ 5 }$$

where $D\equiv 3+2\enspace \mathrm{cos}\enspace {k}_{y}\sqrt{3}+4\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\enspace \mathrm{cos}\enspace \frac{3{k}_{x}}{2}-{\left(i{\omega }_{\ell }\right)}^{2}$, m ∈{0, ±1, ±2}. Due to the translational invariance in the y direction, the integrals above are periodic functions of k_y, with a period given by $\left[0,\frac{2\pi }{\sqrt{3}}\right]$. The Green's function in equation (3) can thus be written as

$$\begin{align}\hfill & {G}_{0}\left(x,{k}_{y},i{\omega }_{\ell }\right)=\left(\begin{matrix}\hfill i{\omega }_{\ell }{X}_{0}\hfill & \hfill 2{X}_{1}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}+{X}_{-2}\hfill \\ \hfill 2{X}_{-1}\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}+{X}_{2}\hfill & \hfill i{\omega }_{\ell }{X}_{0}\hfill \end{matrix}\right),\hfill \end{align} \tag{ 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 equation (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 {\left(\mathbb{I}-V{G}_{0}\right)}^{-1}V$, which in our case becomes:

$$\begin{align}\hfill T\left({k}_{y},i{\omega }_{\ell }\right)& \equiv {\left[\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)-\left(\begin{matrix}\hfill U\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\cdot {G}_{0}\left(x=0,{k}_{y},i{\omega }_{\ell }\right)\right]}^{-1}\cdot \left(\begin{matrix}\hfill U\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\hfill \\ \hfill & \underset{U\to \infty }{\overset{}{\rightarrow }}-\left(\begin{matrix}\hfill {\left[{G}_{0}^{11}\left(x=0,{k}_{y},i{\omega }_{\ell }\right)\right]}^{-1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\hfill \end{align} \tag{ 7 }$$

where ${G}_{0}^{11}$ denotes the 11 component of the corresponding Green's function. We perform the analytical continuation iω_ℓ → E + iδ, δ → +0, and we find that the imaginary part of the trace of the T-matrix,

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace \mathrm{t}\mathrm{r}\enspace T\left({k}_{y},E+i0\right)=\mathrm{I}\mathrm{m}\enspace \frac{\sqrt{{\left(3+2\enspace \mathrm{cos}\enspace {k}_{y}\sqrt{3}-{\left(E+i0\right)}^{2}\right)}^{2}-16\enspace {\mathrm{cos}}^{2}\enspace \frac{{k}_{y}\sqrt{3}}{2}}}{E+i0}\end{equation} \tag{ 8 }$$

has a pole at E = 0 for all ${k}_{y}\in \left[0,\frac{2\pi }{\sqrt{3}}\right]{\backslash}\left\{\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right\}$ (see figure 2). Thus, there exists a zero-energy edge state for all the possible values of k_y except for the special points $\left\{\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right\}$. 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 \left(\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right)$, whereas for the bearded edge we would have it for ${k}_{y}\in \left[0,\frac{2\pi }{\sqrt{3}}\right]{\backslash}\left[\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right]$. The special points $\frac{2\pi }{3\sqrt{3}}$ and $\frac{4\pi }{3\sqrt{3}}$ are the points where the zero-energy state cease 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation}\left[\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)\enspace -\enspace {G}_{0}\left(x=0,{k}_{y},E=0\right)\enspace \cdot \enspace \left(\begin{matrix}\hfill U\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\right]{\Psi}\left(x\enspace =\enspace 0,{k}_{y}\right)=0.\end{equation} \tag{ 9 }$$

We then use the propagation relation

$$\begin{equation}{\Psi}\left(x,{k}_{y}\right)={G}_{0}\left(x,{k}_{y},E=0\right)\cdot \left(\begin{matrix}\hfill U\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right){\Psi}\left(x=0,{k}_{y}\right)\end{equation} \tag{ 10 }$$

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

$$\begin{equation}{\Psi}\left(x,{k}_{y}\right)=\left(\begin{matrix}\hfill 0\hfill \\ \hfill {X}_{0}\left({x}_{B}+1\right)+2\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}{X}_{0}\left({x}_{B}-\frac{1}{2}\right)\hfill \end{matrix}\right),\end{equation} \tag{ 11 }$$

where we have set E = 0 and we have omitted the k_y dependence in X₀ for the sake of brevity. The explicit definitions of x_B and X₀ are given in equations (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 equation (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 equation (11) then simplifies to:

$$\begin{equation}{\Psi}\left(n,{k}_{y}\right)=\left(\begin{matrix}\hfill 0\hfill \\ \hfill -{\left(-2\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\right)}^{-n-1}\hfill \end{matrix}\right),\end{equation} \tag{ 12 }$$

which yields the zigzag-edge wave functions for n < 0, ${k}_{y}\in \left(\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right)$ and the bearded-edge wave functions for n > 0, ${k}_{y}\in \left[0,\frac{2\pi }{\sqrt{3}}\right]{\backslash}\left[\frac{2\pi }{3\sqrt{3}},\frac{4\pi }{3\sqrt{3}}\right]$. 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 figure 3 we plot the corresponding LDOS ρ(x, k_y) = |Ψ(x, k_y)|². 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 equation (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.

Zoom In Zoom Out Reset image size

In this section we turn to calculating the wave functions for the Fermi-arc states and the quasiparticle interference patterns on the surface of a Weyl semimetal described by the following model:

$$\begin{equation}{\mathcal{H}}_{0}\left(\boldsymbol{k}\right)=v\mathrm{sin}{k}_{x}{\sigma }_{x}+v\mathrm{sin}{k}_{y}{\sigma }_{y}+\left(m-t\mathrm{cos}{k}_{x}-t\mathrm{cos}{k}_{y}-t\mathrm{cos}{k}_{z}\right){\sigma }_{z}.\end{equation} \tag{ 13 }$$

The above Hamiltonian is defined on a cubic lattice with the lattice constant set to unity. Here k ≡ (k_x, k_y, k_z), the Pauli matrices are denoted by (σ_x, σ_y, σ_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 $\mathrm{cos}\enspace {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:

$$\begin{equation}{G}_{0}\left(\boldsymbol{k},i{\omega }_{\ell }\right)\equiv {\left[i{\omega }_{\ell }-{\mathcal{H}}_{0}\left(\boldsymbol{k}\right)\right]}^{-1}=-\frac{i{\omega }_{\ell }{\sigma }_{0}+\tilde {g}{\sigma }_{z}+\mathrm{sin}{k}_{x}{\sigma }_{x}+\mathrm{sin}{k}_{y}{\sigma }_{y}}{{\tilde {g}}^{2}+{\mathrm{sin}}^{2}{k}_{x}+{\mathrm{sin}}^{2}{k}_{y}-{\left(i{\omega }_{\ell }\right)}^{2}},\end{equation} \tag{ 14 }$$

where we denote $\tilde {g}\equiv \tilde {g}\left(\boldsymbol{k}\right)\equiv m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{y}-\mathrm{cos}\enspace {k}_{z}$.

3.1. 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 δ-function impurity with the following potential $V\left(y\right)=\left(\begin{matrix}\hfill U\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill U\hfill \end{matrix}\right)\delta \left(y\right)$. To solve the corresponding impurity problem we compute the Matsubara Green's function in the mixed coordinate space and momentum space representation, i.e.,

$$\begin{equation}{G}_{0}\left({k}_{x},y,{k}_{z},i{\omega }_{\ell }\right)\equiv \underset{-\pi }{\overset{\pi }{\int }}\frac{\mathrm{d}{k}_{y}}{2\pi }{G}_{0}\left(\boldsymbol{k},i{\omega }_{\ell }\right){\text{e}}^{\text{i}{k}_{y}y}.\end{equation} \tag{ 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 ⩾ 0 or y ⩽ 0), in other words, y = na, where a = 1 is the lattice constant, and $n\in \mathbb{Z}$. To perform the Fourier transform in equation (15) we define three integrals:

$$\begin{equation}{X}_{s}\equiv -\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{\left({\bullet}\right){\text{e}}^{\text{i}{k}_{y}y}}{{\tilde {g}}^{2}+{\mathrm{sin}}^{2}\enspace {k}_{x}+{\mathrm{sin}}^{2}\enspace {k}_{y}-{\left(i{\omega }_{\ell }\right)}^{2}},\end{equation} \tag{ 16 }$$

where (•) = 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

$$\begin{equation}{G}_{0}\left({k}_{x},y,{k}_{z},i{\omega }_{\ell }\right)=i{\omega }_{\ell }{X}_{0}{\sigma }_{0}+\left(g{X}_{0}-{X}_{1}\right){\sigma }_{z}+\mathrm{sin}\enspace {k}_{x}{X}_{0}{\sigma }_{x}+{X}_{2}{\sigma }_{y},\end{equation} \tag{ 17 }$$

where we denote g ≡ g(k_x, k_z) ≡ 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 equation (17), and assuming that the impurity potential amplitude U is the largest energy scale in the system, i.e., U ≫ 1. Thus, we get

$$\begin{align}\hfill T\left({k}_{x},{k}_{z},i{\omega }_{\ell }\right)& \equiv {\left[\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)-U{G}_{0}\left({k}_{x},y=0,{k}_{z},i{\omega }_{\ell }\right)\right]}^{-1}\hfill \\ \hfill & \underset{U\to \infty }{\overset{}{\rightarrow }}-{\left[{G}_{0}\left({k}_{x},y=0,{k}_{z},i{\omega }_{\ell }\right)\right]}^{-1}.\hfill \end{align} \tag{ 18 }$$

The poles of the imaginary part of the trace of the T-matrix found via analytical continuation iω_ℓ → E + iδ, δ → +0 define the energies of the impurity-bound states at finite values of U, whereas at U → ∞ they yield the surface modes, i.e., the Fermi-arc surface states. We perform the analytical continuation, and we find that

$$\begin{equation}\mathrm{t}\mathrm{r}\enspace T\left({k}_{x},{k}_{z},E+i0\right)=\frac{2\left(E+i0\right){X}_{0}}{{\left(g{X}_{0}+{X}_{1}\right)}^{2}+\left[{\mathrm{sin}}^{2}\enspace {k}_{x}-{\left(E+i0\right)}^{2}\right]{X}_{0}^{2}}\end{equation} \tag{ 19 }$$

whose imaginary part has two poles, at E = ± sin k_x, for k_z lying in the interval between the Weyl nodes, i.e., $\left(-{k}_{z}^{0},{k}_{z}^{0}\right)$, and k_x such that cos k_x ⩾ 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 ⩽ 0, and another one for that in the upper half-space, y ⩾ 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:

$$\begin{equation*}{G}_{0}\left({k}_{x},y=0,{k}_{z},E={\pm}\mathrm{sin}\enspace {k}_{x}\right){{\Psi}}_{{\pm}}\left({k}_{x},y=0,{k}_{z}\right)=0,\end{equation*}$$

whose non-trivial solutions can be found by imposing det G₀(k_x, y = 0, k_z, E = ±sin k_x) = 0, which in turn yields:

$$\begin{equation}E=+\mathrm{sin}\enspace {k}_{x}\to {{\Psi}}_{+}\left({k}_{x},y=0,{k}_{z}\right)\propto \left(\begin{matrix}\hfill 1\hfill \\ \hfill -1\hfill \end{matrix}\right),\end{equation} \tag{ 20 }$$

$$\begin{equation}E=-\mathrm{sin}\enspace {k}_{x}\to {{\Psi}}_{-}\left({k}_{x},y=0,{k}_{z}\right)\propto \left(\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right).\end{equation} \tag{ 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 Ψ(k_x, y, k_z) = G₀(k_x, y, k_z, E) ⋅ UΨ(k_x, y = 0, k_z):

$$\begin{equation}{{\Psi}}_{+}\left({k}_{x},y,{k}_{z}\right)=\left(g{X}_{0}-{X}_{1}+i{X}_{2}\right)\left(\begin{matrix}\hfill 1\hfill \\ \hfill 1\hfill \end{matrix}\right)\end{equation} \tag{ 22 }$$

$$\begin{equation}{{\Psi}}_{-}\left({k}_{x},y,{k}_{z}\right)=\left(g{X}_{0}-{X}_{1}-i{X}_{2}\right)\left(\begin{matrix}\hfill 1\hfill \\ \hfill -1\hfill \end{matrix}\right).\end{equation} \tag{ 23 }$$

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

$$\begin{equation}\rho \left({k}_{x}=0,y,{k}_{z},E=0\right)\equiv \vert {{\Psi}}_{+}\left({k}_{x}=0,y,{k}_{z}\right){\vert }^{2}\end{equation} \tag{ 24 }$$

as a function of y ⩾ 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$.

Zoom In Zoom Out Reset image size

3.2. 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}\left(x,z\right)=\left(\begin{matrix}\hfill {U}_{s}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {U}_{s}\hfill \end{matrix}\right)\delta \left(x,z\right)$ on the surface of a Weyl semimetal. Along the lines of reference [31], the surface Green's function for the Weyl semimetal in the upper half plane (y > 0) is given by:

$$\begin{equation}{G}_{s}\left({k}_{x},{k}_{z}\right)=\enspace \int \enspace \frac{\mathrm{d}{k}_{1y}\enspace \mathrm{d}{k}_{2y}}{{\left(2\pi \right)}^{2}}G\left({k}_{x},{k}_{1y},{k}_{2y},{k}_{z}\right){\left.{\text{e}}^{\text{i}\left({k}_{1y}-{k}_{2y}\right)y}\right\vert }_{y=1},\end{equation} \tag{ 25 }$$

with

$$\begin{equation}G\left({k}_{x},{k}_{1y},{k}_{2y},{k}_{z}\right)=2\pi {G}_{0}\left({k}_{x},{k}_{1y},{k}_{z}\right)\delta \left({k}_{1y}-{k}_{2y}\right)+{G}_{0}\left({k}_{x},{k}_{1y},{k}_{z}\right)T\left({k}_{x},{k}_{z}\right){G}_{0}\left({k}_{x},{k}_{2y},{k}_{z}\right),\end{equation} \tag{ 26 }$$

where the unperturbed bulk Green's function and the T-matrix are given by equations (14) and (18), respectively. We set y = 1 in the expression above because in the presence of an infinite-amplitude δ-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 equation (25) using equation (26), and we obtain the surface Green's function performing both integrals over momenta in equation (25):

$$\begin{equation}{G}_{s}\left({k}_{x},{k}_{z}\right)={G}_{0}\left({k}_{x},y=0,{k}_{z}\right)+{G}_{0}\left({k}_{x},y=1,{k}_{z}\right)T\left({k}_{x},{k}_{z}\right){G}_{0}\left({k}_{x},y=-1,{k}_{z}\right),\end{equation} \tag{ 27 }$$

where the mixed coordinate space-momentum space representation of the unperturbed Green's function is given by equation (17).

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

$$\begin{equation}A\left({k}_{x},{k}_{z}\right)\equiv -\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{t}\mathrm{r}\enspace {G}_{s}\left({k}_{x},{k}_{z}\right).\end{equation} \tag{ 28 }$$

In figure 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation}{\Delta}\rho \left({k}_{x},{k}_{z}\right)=-\frac{1}{2\pi i}\enspace \int \enspace \frac{\mathrm{d}{q}_{x}\enspace \mathrm{d}{q}_{z}}{{\left(2\pi \right)}^{2}}f\left({q}_{x},{q}_{z},{k}_{x},{k}_{z}\right)\end{equation} \tag{ 29 }$$

where

$$\begin{equation*}f\equiv \mathrm{t}\mathrm{r}\left[{G}_{s}\left({q}_{x},{q}_{z}\right){T}_{s}{G}_{s}\left({q}_{x}+{k}_{x},{q}_{z}+{k}_{z}\right)-{G}_{s}^{{\ast}}\left({q}_{x}+{k}_{x},{q}_{z}+{k}_{z}\right){T}_{s}^{{\ast}}{G}_{s}^{{\ast}}\left({q}_{x},{q}_{z}\right)\right].\end{equation*}$$

Also

$$\begin{equation*}{T}_{s}={\left[\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)-{U}_{s}\int \frac{\mathrm{d}{k}_{x}\enspace \mathrm{d}{k}_{z}}{{\left(2\pi \right)}^{2}}{G}_{s}\left({k}_{x},{k}_{z}\right)\right]}^{-1}{U}_{s},\end{equation*}$$

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

In figure 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 $\left[-\pi /2,\pi /2\right]$, 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 $-2{k}_{z}^{0}$ to $-2{k}_{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 $2{k}_{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 $\left[-\pi ,\pi \right]$.

Zoom In Zoom Out Reset image size

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

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 [34–36], or by solving the Schrödinger equation with specifically derived boundary conditions [17, 42], or other methods [25, 27].

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 calculations 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 [58]; 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 [59].

VK and JHB would like to acknowledge the ERC Starting Grant No. 679722. VK 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.

In this Appendix we calculate the integrals in equation (5) for m ∈{−2, −1, 0, 1, 2}:

$$\begin{equation}{X}_{m}\left(x,{k}_{y},i{\omega }_{\ell }\right)\equiv -\enspace \underset{-2\pi /3}{\overset{2\pi /3}{\int }}\enspace \frac{\mathrm{d}{k}_{x}}{4\pi /3}\frac{{\text{e}}^{\text{i}{k}_{x}m/2}\enspace {\text{e}}^{\text{i}{k}_{x}x}}{3+2\enspace \mathrm{cos}\enspace {k}_{y}\sqrt{3}+4\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}\enspace \mathrm{cos}\enspace \frac{3{k}_{x}}{2}-{\left(i{\omega }_{\ell }\right)}^{2}}.\end{equation} \tag{ 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 \left[0,\frac{2\pi }{\sqrt{3}}\right]$. 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 figure 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:

$$\begin{equation}{X}_{m}\left(x,{k}_{y},i{\omega }_{\ell }\right)=-\frac{1}{1-{\left(i{\omega }_{\ell }\right)}^{2}}\enspace \underset{-2\pi /3}{\overset{2\pi /3}{\int }}\enspace \frac{\mathrm{d}{k}_{x}}{4\pi /3}\enspace {\text{e}}^{\text{i}{k}_{x}m/2}\enspace {\text{e}}^{\text{i}{k}_{x}x}=-\frac{1}{1-{\left(i{\omega }_{\ell }\right)}^{2}}\frac{\mathrm{sin}\enspace \pi \left(\frac{2}{3}x+\frac{m}{3}\right)}{\pi \left(\frac{2}{3}x+\frac{m}{3}\right)}\end{equation} \tag{ A2 }$$

When additionally x = −m/2, we get: ${X}_{0}\left(0,\frac{\pi }{\sqrt{3}},i{\omega }_{\ell }\right)=-\frac{1}{1-{\left(i{\omega }_{\ell }\right)}^{2}}.$

In what follows we turn to the case of ${k}_{y}\ne \frac{\pi }{\sqrt{3}}$. In this case the integral in equation (A1) can be rewritten as a contour integral in the complex plane, using the substitution $z={\text{e}}^{\frac{3}{2}{k}_{x}}$:

$$\begin{equation}{X}_{m}\left(x,{k}_{y},i{\omega }_{\ell }\right)\equiv -\frac{1}{2\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}}\enspace \frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{\frac{2}{3}\left(x+\frac{m}{2}\right)}}{{z}^{2}+fz+1},\end{equation} \tag{ A3 }$$

where the circle |z| = 1 is oriented counter-clockwise, and we defined

$$\begin{equation}f\equiv f\left({k}_{y},i{\omega }_{\ell }\right)=\frac{3+2\mathrm{cos}{k}_{y}\sqrt{3}-{\left(i{\omega }_{\ell }\right)}^{2}}{2\mathrm{cos}\frac{{k}_{y}\sqrt{3}}{2}}.\end{equation} \tag{ A4 }$$

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

$$\begin{equation}\text{sublattice}\enspace A:\quad {X}_{m}^{A}\left(x=\frac{3n}{2},{k}_{y},i{\omega }_{\ell }\right)=-\frac{1}{2\mathrm{cos}\frac{{k}_{y}\sqrt{3}}{2}}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{{z}^{2}+fz+1},\end{equation} \tag{ A5 }$$

$$\begin{equation}\text{sublattice}\enspace B:\quad {X}_{m}^{B}\left(x=\frac{3n+1}{2},{k}_{y},i{\omega }_{\ell }\right)=-\frac{1}{2\mathrm{cos}\frac{{k}_{y}\sqrt{3}}{2}}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+\left(m+1\right)/3}}{{z}^{2}+fz+1},\end{equation} \tag{ 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}\left(x\right)={X}_{2}^{A}\left(-x\right),{X}_{-1}^{A}\left(x\right)={X}_{1}^{A}\left(-x\right),$ whereas for sublattice B ${X}_{-2}^{B}\left(x\right)={X}_{1}^{A}\left(\frac{1}{2}-x\right)$, ${X}_{-1}^{B}\left(x\right)={X}_{0}^{A}\left(\frac{1}{2}-x\right)$, ${X}_{0}^{B}\left(x\right)={X}_{1}^{A}\left(x-\frac{1}{2}\right)$, ${X}_{1}^{B}\left(x\right)={X}_{2}^{A}\left(x-\frac{1}{2}\right)$, and ${X}_{2}^{B}\left(x\right)={X}_{0}^{A}\left(x+1\right)$. Above we omitted the dependence on k_y and iω_ℓ. Therefore, it is sufficient to compute only the three integrals in equation (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 equation (A5) it is clear that the case of m = 0 and the case of m ≠ 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 equation (A5). A simple analysis of the roots given by
  $$\begin{equation}{z}_{{\pm}}=\frac{1}{2}\left[-f{\pm}\sqrt{{f}^{2}-4}\right]\end{equation} \tag{ 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,
  $$\begin{equation}\quad \text{if}\enspace 0{\leqslant}{k}_{y}{< }\frac{\pi }{\sqrt{3}},\quad \text{then}\quad {z}_{+}\in \left(-1,0\right),\enspace {z}_{-}{< }-1,\end{equation} \tag{ A8 }$$
  $$\begin{equation}\quad \text{if}\enspace \frac{\pi }{\sqrt{3}}{< }{k}_{y}{\leqslant}\frac{2\pi }{\sqrt{3}},\quad \text{then}\quad {z}_{-}\in \left(0,+1\right)\enspace {z}_{+}{ >}+1.\end{equation} \tag{ 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.

A.1. The sub-case of m = 0

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

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}={\ast}.\end{equation} \tag{ 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{\leqslant}{k}_{y}{< }\pi /\sqrt{3}$), and using the residue theorem, we get:

$$\begin{equation}{\ast}=\underset{z={z}_{+}}{\mathrm{r}\mathrm{e}\mathrm{s}}\enspace \frac{{z}^{n}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}=-\frac{{z}_{+}^{\vert n\vert }}{{z}_{-}-{z}_{+}}.\end{equation} \tag{ A11 }$$

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}{\leqslant}2\pi /\sqrt{3}$), we just need to exchange z₊ → 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:

$$\begin{equation}{X}_{0}^{A}\left(x=\frac{3n}{2},{k}_{y},i{\omega }_{\ell }\right)=-\frac{1}{2\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}}\begin{cases}-\frac{{z}_{+}^{\vert n\vert }}{{z}_{-}-{z}_{+}},\quad \hfill & 0{\leqslant}{k}_{y}{< }\frac{\pi }{\sqrt{3}}\hfill \\ +\frac{{z}_{-}^{\vert n\vert }}{{z}_{-}-{z}_{+}},\quad \hfill & \frac{\pi }{\sqrt{3}}{< }{k}_{y}{\leqslant}\frac{2\pi }{\sqrt{3}}\hfill \end{cases}.\end{equation} \tag{ A12 }$$

A.2. The sub-cases of m = 1 and m = 2

In this subsection we calculate the integral

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}={\ast}\end{equation} \tag{ 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 \left(-\infty ,0\right]$. In order to start the calculation, we first define auxiliary contours in the complex plane, Γ = C₁ ∪ γ₊ ∪ γ₋ (see left and right panels of figure 7). The unit circle is denoted C₁, and γ_± are the right and left banks of the branch cut, parametrized by t ± i0 with $t\in \left[-1,0\right]$, 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{\leqslant}{k}_{y}{< }\pi /\sqrt{3}$ and $\pi /\sqrt{3}{< }{k}_{y}{\leqslant}2\pi /\sqrt{3}$ we choose the left and the right panels of figure 7, 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 Γ we neglect the tiny line between γ₊ and γ₋, i.e., $\left[-i0,+i0\right]$, since it never contributes to the value of the integral figure 7.

Zoom In Zoom Out Reset image size

We start by considering the case of $0{\leqslant}{k}_{y}{< }\pi /\sqrt{3}$ and we write down the residue theorem for the contour Γ = C₁ ∪ γ₊ ∪ γ₋ defined in the left panel of figure 7. Since there are no poles inside the contour Γ we have:

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{{\Gamma}}{\oint }\mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}=0=\frac{1}{2\pi i}\left[\enspace \underset{{C}_{1}}{\int }+\underset{{\gamma }_{+}}{\int }+\underset{{\gamma }_{-}}{\int }\right]\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}.\end{equation} \tag{ A14 }$$

Therefore, we can express the sought-for line integral along C₁ ≡ {∀z: |z| = 1} as

$$\begin{equation*}{\ast}=-\frac{1}{2\pi i}\left[\enspace \underset{{\gamma }_{+}}{\int }+\underset{{\gamma }_{-}}{\int }\right]\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}={\ast}{\ast}.\end{equation*}$$

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

$$\begin{align*}\hfill \underset{{\gamma }_{+}}{\int }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}& =\underset{-1}{\overset{0}{\int }}\mathrm{d}t\frac{{\left(t+i0\right)}^{n+m/3}}{\left(t+i0-{z}_{+}\right)\left(t-{z}_{-}\right)}=\underset{-1}{\overset{0}{\int }}\mathrm{d}t\frac{{e}^{\left(n+m/3\right)\mathrm{L}\mathrm{n}\left(t+i0\right)}}{\left(t+i0-{z}_{+}\right)\left(t-{z}_{-}\right)}\hfill \\ \hfill & =\underset{-1}{\overset{0}{\int }}\mathrm{d}t\frac{{\mathrm{e}}^{\left(n+m/3\right)\left[\mathrm{ln}\vert t\vert +\pi i\right]}}{\left(t+i0-{z}_{+}\right)\left(t-{z}_{-}\right)}={\text{e}}^{\text{i}\pi \left(n+m/3\right)}\underset{-1}{\overset{0}{\int }}\mathrm{d}t\frac{\vert t{\vert }^{n+m/3}}{\left(t+i0-{z}_{+}\right)\left(t-{z}_{-}\right)}\hfill \\ \hfill & ={\text{e}}^{\text{i}\pi \left(n+m/3\right)}\underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{\left(w+{z}_{+}-i0\right)\left(w+{z}_{-}\right)}\hfill \\ \hfill & ={\mathrm{e}}^{\text{i}\pi \left(n+m/3\right)}\left[\mathcal{P}\enspace \underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{\left(w+{z}_{+}\right)\left(w+{z}_{-}\right)}+i\pi \underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{w+{z}_{-}}\delta \left(w+{z}_{+}\right)\right].\hfill \end{align*}$$

Similarly, we get:

$$\begin{equation*}\underset{{\gamma }_{-}}{\int }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}=-{\text{e}}^{-\text{i}\pi \left(n+m/3\right)}\left[\mathcal{P}\enspace \underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{\left(w+{z}_{+}\right)\left(w+{z}_{-}\right)}-i\pi \underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{w+{z}_{-}}\delta \left(w+{z}_{+}\right)\right].\end{equation*}$$

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

$$\begin{equation*}{\ast}{\ast}=\frac{{\left(-1\right)}^{n+1}}{\pi }\left({I}_{\mathcal{P}}\enspace \mathrm{sin}\enspace \frac{\pi m}{3}+\pi {I}_{\delta }\enspace \mathrm{cos}\enspace \frac{\pi m}{3}\right)={\ast}{\ast}{\ast},\end{equation*}$$

where we defined

$$\begin{align}\hfill {I}_{\mathcal{P}}& =\mathcal{P}\enspace \underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{\left(w+{z}_{+}\right)\left(w+{z}_{-}\right)}\hfill \\ \hfill & =\frac{{\left(-{z}_{-}\right)}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{-}},n+\frac{m}{3}+1,0\right]-{\left(-{z}_{+}\right)}^{n+\frac{m}{3}}\left(B\left[-{z}_{+},-n-\frac{m}{3},0\right]+\pi \enspace \mathrm{c}\mathrm{t}\mathrm{g}\enspace \frac{\pi m}{3}\right)}{{z}_{-}-{z}_{+}}\hfill \end{align} \tag{ A15 }$$

$$\begin{equation}{I}_{\delta }=\underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{w+{z}_{-}}\delta \left(w+{z}_{+}\right)=\frac{{\left(-{z}_{+}\right)}^{n+m/3}}{{z}_{-}-{z}_{+}}.\end{equation} \tag{ A16 }$$

Substituting ${I}_{\mathcal{P}}$ and I_δ into the equation above, we obtain after simplifications:

$$\begin{equation*}{\ast}{\ast}{\ast}=\frac{{\left(-1\right)}^{n+1}}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\frac{{\left(-{z}_{-}\right)}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{-}},n+\frac{m}{3}+1,0\right]-{\left(-{z}_{+}\right)}^{n+\frac{m}{3}}B\left[-{z}_{+},-n-\frac{m}{3},0\right]}{{z}_{-}-{z}_{+}}.\end{equation*}$$

Thus we get:

$$\begin{align*}\hfill & \frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}\hfill \\ \hfill & \quad =\frac{{\left(-1\right)}^{n+1}}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\frac{{\left(-{z}_{-}\right)}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{-}},n+\frac{m}{3}+1,0\right]-{\left(-{z}_{+}\right)}^{n+\frac{m}{3}}B\left[-{z}_{+},-n-\frac{m}{3},0\right]}{{z}_{-}-{z}_{+}},\hfill \end{align*}$$

where B[z, α, β] is the incomplete Beta-function defined as:

$$\begin{equation}B\left[z,\alpha ,\beta \right]\equiv \underset{0}{\overset{z}{\int }}{t}^{\alpha -1}{\left(1-t\right)}^{\beta -1}\mathrm{d}t.\end{equation} \tag{ A17 }$$

In what follows we consider the remaining case of $\frac{\pi }{\sqrt{3}}{< }{k}_{y}{\leqslant}\frac{2\pi }{\sqrt{3}}$. First, we rewrite the integral as follows:

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}=\frac{1}{{z}_{-}-{z}_{+}}\left[\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{-}}-\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}\right].\end{equation} \tag{ 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:

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{-}}=\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \frac{\mathrm{d}w}{{w}^{2}}\frac{{w}^{-n-m/3}}{1/w-{z}_{-}}=-\frac{1}{{z}_{-}}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}w\frac{{w}^{-n-m/3-1}}{w-1/{z}_{-}}=-\frac{1}{{z}_{-}}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}w\frac{{w}^{-n-m/3-1}}{w-{z}_{+}}.\end{equation} \tag{ A19 }$$

From the expression above we see that to obtain the first integral from the second one we need to replace n → −n − 2, m → 3 − m and multiply the result by −1/z₋. Next, we compute the second integral

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}\end{equation} \tag{ A20 }$$

by writing down the residue theorem for the contour Γ = C₁ ∪ γ₊ ∪ γ₋ defined in the panel panel of figure 7 for the second integral. Since there are no poles inside the contour Γ we have:

$$\begin{equation}\frac{1}{2\pi i}\enspace \underset{{\Gamma}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}=0=\frac{1}{2\pi i}\left[\enspace \underset{{C}_{1}}{\int }+\underset{{\gamma }_{+}}{\int }+\underset{{\gamma }_{-}}{\int }\right]\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}.\end{equation} \tag{ A21 }$$

Therefore, we can express the sought-for line integral along C₁ ≡ {∀z: |z| = 1} as

$$\begin{equation*}\frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}=-\frac{1}{2\pi i}\left[\enspace \underset{{\gamma }_{+}}{\int }+\underset{{\gamma }_{-}}{\int }\right]\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}={\ast}.\end{equation*}$$

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

$$\begin{equation}\underset{{\gamma }_{{\pm}}}{\int }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}={\pm}\underset{-1}{\overset{0}{\int }}\mathrm{d}t\frac{{\left(t{\pm}i0\right)}^{n+m/3}}{t-{z}_{+}}={\pm}{e}^{{\pm}i\pi \left(n+m/3\right)}\underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{w+{z}_{+}}.\end{equation} \tag{ A22 }$$

Thus, we get:

$$\begin{align*}\hfill {\ast}& =-\frac{1}{2\pi i}\left[\enspace \underset{{\gamma }_{+}}{\int }+\underset{{\gamma }_{-}}{\int }\right]\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{+}}=-\frac{{\left(-1\right)}^{n}}{2\pi i}\left[{\text{e}}^{\text{i}\frac{\pi m}{3}}-{\text{e}}^{-\text{i}\frac{\pi m}{3}}\right]\underset{0}{\overset{1}{\int }}\mathrm{d}w\frac{{w}^{n+m/3}}{w+{z}_{+}}\hfill \\ \hfill & =-\frac{1}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\enspace {\text{e}}^{-\text{i}\frac{\pi m}{3}}{z}_{+}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{+}},n+\frac{m}{3}+1,0\right].\hfill \end{align*}$$

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

$$\begin{align}\hfill \frac{1}{2\pi i}\enspace \underset{{C}_{1}}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{z-{z}_{-}}& =\frac{1}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi \left(3-m\right)}{3}\enspace {\text{e}}^{-\text{i}\frac{\pi \left(3-m\right)}{3}}{z}_{-}^{n+\frac{m}{3}}B\left[-{z}_{-},-n-\frac{m}{3},0\right]\hfill \\ \hfill & =-\frac{1}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\enspace {\text{e}}^{\text{i}\frac{\pi m}{3}}{z}_{-}^{n+\frac{m}{3}}B\left[-{z}_{-},-n-\frac{m}{3},0\right].\hfill \end{align} \tag{ A23 }$$

Finally:

$$\begin{align}\hfill & \frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{n+m/3}}{\left(z-{z}_{+}\right)\left(z-{z}_{-}\right)}\hfill \\ \hfill & \quad =\frac{1}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\enspace \frac{{\text{e}}^{-\text{i}\frac{\pi m}{3}}{z}_{+}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{+}},n+\frac{m}{3}+1,0\right]-{\text{e}}^{\text{i}\frac{\pi m}{3}}{z}_{-}^{n+\frac{m}{3}}B\left[-{z}_{-},-n-\frac{m}{3},0\right]}{{z}_{-}-{z}_{+}}.\hfill \end{align} \tag{ A24 }$$

Combining the results for different ranges of k_y, we present the final result:

$$\begin{equation*}{X}_{m}^{A}\left(x=\frac{3n}{2},{k}_{y},i{\omega }_{\ell }\right)=-\frac{1}{2\enspace \mathrm{cos}\enspace \frac{{k}_{y}\sqrt{3}}{2}}\begin{cases}\frac{{\left(-1\right)}^{n+1}}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\frac{{\left(-{z}_{-}\right)}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{-}},n+\frac{m}{3}+1,0\right]-{\left(-{z}_{+}\right)}^{n+\frac{m}{3}}B\left[-{z}_{+}-n-\frac{m}{3},0\right]}{{z}_{-}-{z}_{+}}\hfill & 0{\leqslant}{k}_{y}{< }\frac{\pi }{\sqrt{3}}\hfill \\ \frac{1}{\pi }\enspace \mathrm{sin}\enspace \frac{\pi m}{3}\enspace \frac{{\text{e}}^{-\text{i}\frac{\pi m}{3}}{z}_{+}^{n+\frac{m}{3}}B\left[-\frac{1}{{z}_{+}},n+\frac{m}{3}+1,0\right]-{\text{e}}^{\text{i}\frac{\pi m}{3}}{z}_{-}^{n+\frac{m}{3}}B\left[-{z}_{-},-n-\frac{m}{3},0\right]}{{z}_{-}-{z}_{+}}\hfill & \frac{\pi }{\sqrt{3}}{< }{k}_{y}{\leqslant}\frac{2\pi }{\sqrt{3}}\hfill \end{cases},\end{equation*}$$

where m = 1 or m = 2.

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

$$\begin{equation}{X}_{s}\equiv -\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{\left({\bullet}\right){\text{e}}^{\text{i}{k}_{y}y}}{{\tilde {g}}^{2}+{\mathrm{sin}}^{2}\enspace {k}_{x}+{\mathrm{sin}}^{2}\enspace {k}_{y}-{\left(i{\omega }_{\ell }\right)}^{2}},\end{equation} \tag{ B1 }$$

where (•) = 1, cos k_y and sin k_y for s = 0, 1 and 2, respectively. First, it is easy to see that integrals X₁ and X₂ can be expressed in terms of the integral X₀ in the following way:

$$\begin{align}\hfill {X}_{1}\left({k}_{x},y,{k}_{z},i{\omega }_{\ell }\right)& =-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{\mathrm{cos}\enspace {k}_{y}\enspace {\text{e}}^{\text{i}{k}_{y}y}}{D}=\frac{1}{2}\left[-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\enspace \frac{{\text{e}}^{\text{i}{k}_{y}\left(y+1\right)}}{D}-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\enspace \frac{{\text{e}}^{\text{i}{k}_{y}\left(y-1\right)}}{D}\right]\hfill \\ \hfill & =\frac{1}{2}\left[{X}_{0}\left({k}_{x},y+1,{k}_{z},i{\omega }_{\ell }\right)+{X}_{0}\left({k}_{x},y-1,{k}_{z},i{\omega }_{\ell }\right)\right]\hfill \end{align} \tag{ B2 }$$

$$\begin{align}\hfill {X}_{2}\left({k}_{x},y,{k}_{z},i{\omega }_{\ell }\right)& =-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{\mathrm{sin}\enspace {k}_{y}\enspace {\text{e}}^{\text{i}{k}_{y}y}}{D}=\frac{1}{2i}\left[-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{{\text{e}}^{\text{i}{k}_{y}\left(y+1\right)}}{D}+\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{{\text{e}}^{\text{i}{k}_{y}\left(y-1\right)}}{D}\right]\hfill \\ \hfill & =\frac{1}{2i}\left[{X}_{0}\left({k}_{x},y+1,{k}_{z},i{\omega }_{\ell }\right)-{X}_{0}\left({k}_{x},y-1,{k}_{z},i{\omega }_{\ell }\right)\right],\hfill \end{align} \tag{ B3 }$$

where for the sake of brevity we denoted $D\equiv {\tilde {g}}^{2}+{\mathrm{sin}}^{2}\enspace {k}_{x}+{\mathrm{sin}}^{2}\enspace {k}_{y}-{\left(i{\omega }_{\ell }\right)}^{2},\enspace \tilde {g}=m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{y}\enspace \mathrm{cos}\enspace {k}_{z}$. Therefore, we need to compute only the integral X₀:

$$\begin{align*}\hfill {X}_{0}\left({k}_{x},y,{k}_{z},i{\omega }_{\ell }\right)& =-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{{\text{e}}^{\text{i}{k}_{y}y}}{{m}^{2}+{\omega }_{\ell }^{2}+\left(\mathrm{cos}\enspace {k}_{x}+\mathrm{cos}\enspace {k}_{y}+\mathrm{cos}\enspace {k}_{z}\right)\left(-2m+\mathrm{cos}\enspace {k}_{x}+\mathrm{cos}\enspace {k}_{y}+\mathrm{cos}\enspace {k}_{z}\right)+{\mathrm{sin}}^{2}\enspace {k}_{x}+{\mathrm{sin}}^{2}\enspace {k}_{y}}\hfill \\ \hfill & =-\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{{\text{e}}^{\text{i}{k}_{y}y}}{2\left(-m+\mathrm{cos}\enspace {k}_{x}+\mathrm{cos}\enspace {k}_{z}\right)\mathrm{cos}\enspace {k}_{y}+1+{m}^{2}+{\omega }_{\ell }^{2}+\left(\mathrm{cos}\enspace {k}_{x}+\mathrm{cos}\enspace {k}_{z}\right)\left(-2m+\mathrm{cos}\enspace {k}_{x}+\mathrm{cos}\enspace {k}_{z}\right)+{\mathrm{sin}}^{2}{k}_{x}}={\ast}.\hfill \end{align*}$$

If −m + cos k_x + cos k_z = 0, then we have:

$$\begin{equation}{\ast}=-\frac{1}{1+{\mathrm{sin}}^{2}\enspace {k}_{x}-{\left(i{\omega }_{\ell }\right)}^{2}}\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\enspace {\text{e}}^{\text{i}{k}_{y}y}=-\frac{1}{1+{\mathrm{sin}}^{2}\enspace {k}_{x}-{\left(i{\omega }_{\ell }\right)}^{2}}{\delta }_{y,0}.\end{equation} \tag{ B4 }$$

If −m + cos k_x + cos k_z ≠ 0, then

$$\begin{align*}\hfill {\ast}& =\frac{1}{2\left(m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{z}\right)}\enspace \underset{-\pi }{\overset{\pi }{\int }}\enspace \frac{\mathrm{d}{k}_{y}}{2\pi }\frac{{\text{e}}^{\text{i}{k}_{y}y}}{\mathrm{cos}\enspace {k}_{y}-f\left({k}_{x},{k}_{z}\right)}\hfill \\ \hfill & =\frac{1}{2\left(m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{z}\right)}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \frac{\mathrm{d}z}{z}\frac{{z}^{y}}{\frac{1}{2}\left(z+{z}^{-1}\right)-f\left({k}_{x},{k}_{z}\right)}\hfill \\ \hfill & =\frac{1}{m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{z}}\frac{1}{2\pi i}\enspace \underset{\vert z\vert =1}{\oint }\enspace \mathrm{d}z\frac{{z}^{y}}{{z}^{2}-2f\left({k}_{x},{k}_{z}\right)z+1}={\ast}{\ast}.\hfill \end{align*}$$

Above we introduced

$$\begin{equation*}f\left({k}_{x},{k}_{z}\right)\equiv \frac{1+{g}^{2}+{\mathrm{sin}}^{2}\enspace {k}_{x}-{\left(i{\omega }_{\ell }\right)}^{2}}{2g},\quad g\left({k}_{x},{k}_{z}\right)\equiv m-\mathrm{cos}\enspace {k}_{x}-\mathrm{cos}\enspace {k}_{z}.\end{equation*}$$

Note, that since m ∈ (1, 3), g(k_x, k_z) ∈ (−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 ⩾ 0. The calculation for y < 0 is not needed, since the integral is symmetric with respect to changing y → −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:

$$\begin{equation}{z}_{{\pm}}=f{\pm}\sqrt{{f}^{2}-1}.\end{equation} \tag{ B5 }$$

It is easy to show that

$$\begin{align*}\hfill {z}_{+}\in \left(-1,0\right),{z}_{-}{< }-1\quad \text{when}\enspace g\in \left(-1,0\right),\\ \hfill {z}_{-}\in \left(0,+1\right),{z}_{+}{ >}+1\quad \text{when}\enspace g\in \left(0,+1\right).\end{align*}$$

Therefore, for the integral above we get:

$$\begin{equation}{\ast}{\ast}=\frac{1}{g}\begin{cases}\frac{{z}_{+}^{\vert y\vert }}{{z}_{+}-{z}_{-}}\quad \hfill & \quad \text{for}\enspace g\in \left(-1,0\right)\hfill \\ \frac{{z}_{-}^{\vert y\vert }}{{z}_{-}-{z}_{+}}\quad \hfill & \quad \text{for}\enspace g\in \left(0,+1\right)\hfill \end{cases}.\end{equation} \tag{ B6 }$$