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

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, 1-3 topological insulators and topological superconductors, [4][5][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][9][10][11][12][13] solving the Schrödinger equation with the corresponding boundary conditions, [14][15][16][17][18] using the bulk-boundary correspondence [19][20][21][22][23] and the extended Bloch theorem. 24 Here we focus on a qualitatively different approach to solving the problem of finding boundary modes, which was introduced in Ref. [25]. 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 zigzagand bearded-edge graphene, as well as for their threedimensional generalization -Fermi-arc surface states in Weyl semimetals described by a minimal tight-binding model. 26 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, [27][28][29][30] and analytical studies of the Schrödinger equation with specific boundary conditions in graphene [31][32][33] and Weyl semimetals. 34,35 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 [36][37][38][39] 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 tightbinding 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 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 arXiv:2005.01719v2 [cond-mat.mes-hall] 14 May 2020 The bare Matsubara Green's function for the Hamiltonian in Eq. (1) can be calculated using the standard definition: where ω = πT (2 + 1), with ∈ 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 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 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: The integration limits ±2π/3 and the numerical factor 4π/3 were derived in Ref. [40] 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, 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 Eq. (3), we first denote: x y A B Figure 1. Graphene lattice with an infinite-amplitude δ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).
where D ≡ 3 + 2 cos k y √ 3 + 4 cos Due to the translational invariance in the y direction, the integrals above are periodic functions of k y , with a period given by 0, 2π 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 41,42 −1 V , which in our case becomes: where G 11 0 denotes the 11 component of the corresponding Green's function. We perform the analytical contin- to recover the dependence of the wave function on x, and we get: 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 [27][28][29][30][31][49][50][51] 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 ≡ 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: which yields the zigzag-edge wave functions for n < 0, and the bearded-edge wave functions 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 [27][28][29][30] and specific boundary conditions. 17,31 In Fig. 3 we plot the corresponding lattice local density of states (LDOS) ρ(x, k y ) = |Ψ(x, k y )| 2 . For values of k y between 2π 3 √ 3 and 4π 3 √ 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 interference patterns on the surface of a Weyl semimetal described by the following model: 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 cos k 0 z = 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: 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 δ-function impurity with the following potential To solve the corresponding impurity problem we compute the Matsubara Green's function in the mixed coordinate space and momentum space representation, i.e., 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 ∈ Z. To perform the Fourier transform in Eq. (15) we define three integrals: 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 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 1. Thus, we get The poles of the imaginary part of the trace of the Tmatrix 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 whose imaginary part has two poles, at E = ± sin k x , for k z lying in the interval between the Weyl nodes, i.e., −k 0 z , k 0 z , 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. [43][44][45][46][47][48] First, we find the wave function at y = 0 from the equation: whose non-trivial solutions can be found by imposing det G 0 (k x , y = 0, k z , E = ± sin k x ) = 0, which in turn yields: 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 0 (k x , y, k z , E) · U Ψ(k x , y = 0, k z ): 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., 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 0 z = ±2π/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) = U s 0 0 U s δ(x, z) on the surface of a Weyl semimetal.
Along the lines of Ref. [40], the surface Green's function for the Weyl semimetal in the upper half plane (y > 0) is given by: with G(k x , k 1y , k 2y , k z ) = 2πG 0 (k x , k 1y , k z )δ(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 ), where the unperturbed bulk Green's function and the Tmatrix are given by Eqs. (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 Eq. (25) using Eq. (26), and we obtain the surface Green's function performing both integrals over momenta in Eq. (25): where the mixed coordinate space-momentum space representation of the unperturbed Green's function is given by Eq. (17).
In what follows, we study the spectral function defined through the surface Green's function in Eq. (27) as follows: 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 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 [−π/2, π/2], or in other words, when k 0 z < π/2. Most of the weight in the plot is concentrated, as expected, in a narrow line lying from −2k 0 z to −2k 0 z at k x = 0, which corresponds to surfacesurface scattering processes. Indeed, most of the scattering occurs within the Fermi arc, and the maximum  29), plotted as a function of kx and kz. The mass term is set to m = 2.5, yielding Weyl nodes at k 0 z = ±π/3. The impurity potential amplitude is set to Us = 1. The strong line-like feature in the middle reflects the surface-surface scattering processes, while the small background (∆ρ ≈ 10 −2 ) accounts for the surface-bulk and bulk-bulk ones.
change in momentum along the k z axis is 2k 0 z , 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 [−π, π].

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, [27][28][29] or by solving the Schrödinger equation with specifically derived boundary conditions, 17,35 or other methods. 22,24 The technique we have employed is very general and can be applied to any lattice model no matter its complexity or its dimensionality. 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 52 or considering more realistic models for studying quasiparticle interference patterns in Weyl and Dirac semimetals. 53 • It is easy to notice that for sublattice A we have , and X B 2 (x) = X A 0 (x + 1). Above we omitted the dependence on k y and iω . 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 = 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 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, 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): 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 k y < π/ √ 3), and using the residue theorem, we get: * = res To get the expression for the case where z = z − is inside the unit circle, and z = z + is outside (π/ √ 3 < k y 2π/ √ 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: 2. The sub-cases of m = 1 and m = 2 In this subsection we calculate the integral 1 2πi 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 ∈ (−∞, 0]. In order to start the calculation, we first define auxiliary contours in the complex plane, Γ = C 1 ∪γ + ∪γ − (see left and right panels of Fig. A 2). The unit circle is denoted C 1 , and γ ± are the right and left banks of the branch cut, parametrized by t ± i0 with t ∈ [−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 k y < π/ √ 3 and π/ √ 3 < k y 2π/ √ 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 Γ we neglect the tiny line between γ + and γ − , i.e., [−i0, +i0], since it never contributes to the value of the integral.