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

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


I. INTRODUCTION
Boundaries of certain condensed matter systems host unique phenomena. For instance, graphene exhibits zeroenergy zigzag-edge modes, 1 and topological insulators exhibit conducting edge or surface states. [2][3][4] In order to describe boundary effects, several techniques were developed, including the exact diagonalization of tight-binding Hamiltonians, 5,6 iterative methods to compute boundary Green's functions, [7][8][9] solving the Schrödinger equation 10 and the bulk-boundary correspondence. 11 A new method describing the formation of boundary modes was recently introduced. 12 This method can be generalized to any dimensions, and in certain situations it can yield fully analytical results, providing a deeper physical insight than numerical techniques. The general idea is as follows: instead of considering a finite system with a sharp boundary, we consider an infinite system with a strong delta-potential impurity emulating the shape of the boundary. For example, in order to recover end, edge or surface boundaries, the impurity potential should be chosen to be point-like, line-like and plane-like, respectively. In the limit of an infinite impurity potential such impurities divide a given system into two independent semi-infinite regions. Subsequently we use the T -matrix formalism 13,14 to study the impurity-induced states which transform into boundary states when the impurity strength is larger than any energy scale in the system.
Along the same lines, we present here a direct and non-numerical technique to calculate the surface Green's functions of an arbitrary three-dimensional system. The boundary can once more be modeled as a plane impurity potential with an amplitude going to infinity. The corresponding full Green's functions can be calculated exactly using the T -matrix formalism. The resulting Green's function evaluated on the plane neighboring and parallel to the impurity plane becomes the surface Green's function (see Fig. 1). We apply this technique to calculate the surface Green's functions for Weyl semimetals described by two different models. 15, 16 We recover in each case the corresponding Fermi-arc states. Moreover, in this work we apply the technique from Ref. [12] to a new class of systems -topological insulators. In particular we consider a 2D honeycomb lattice described by the Kane-Mele model, 2 as well as a 2D Chern insulator. 3 We show that impurity-induced states in these two models transform into helical or chiral edge modes, respectively, when the impurity potential is taken to infinity. While the Kane-Mele model requires performing a numerical integration, the 2D Chern insulator allows an exact closed-form solution and thus demonstrates the analytical power of the method.
The paper is organized as follows: in Sec. II we introduce the formalism and the notations. In Sec. III we present the calculation of the surface Green's functions for a Weyl semimetal described by two different models, and the formation of the corresponding Fermi-arc states. In section IV and V we focus on two-dimensional topological insulators described by the Kane-Mele and Chern-insulator models, respectively, and we obtain the corresponding edge modes. We leave the conclusions to Sec. VI.

II. T -MATRIX FORMALISM FOR SURFACE GREEN'S FUNCTIONS AND EDGE STATES
Below we consider an infinite system described by a momentum-space Hamiltonian H k . The unperturbed Matsubara Green's function (GF) can be written as: where ω n denote the Matsubara frequencies. In the presence of an impurity, the Green's function is modified to: where the T -matrix T (k 1 , k 2 , iω n ) embodies all-order impurity-scattering processes. 14,17 Note that due to the impurity-induced breaking of translational symmetry, and consequently of the momentum conservation, the generalized Green's function depends no longer on one, but on two values of momentum. For a delta-function impurity V imp (r) ≡ V δ (x), the form of the T -matrix in 1D is momentum independent and is given by while in 2D and 3D, for a line and plane impurity respectively, localized at x = 0 and perpendicular to the x axis, we have: respectively, with L k being a normalization factor. The limits of integration are given by the boundaries of the first Brillouin zone, i.e., for a 1D system and for a 2D square lattice or 3D cubic lattice we integrate from −π to π (with L k = 2π), while for a honeycomb lattice with an impurity along y we integrate from −2π/3 to 2π/3 (L k = 4π/3) (for a full justification see Appendix A). Note that Eqs. (3) and (4) are independent of k 1x and k 2x due to the fact that the impurity potential is a delta-function centered at x = 0. Reversely, note that the T-matrix contains the terms δ k1y,k2y (2D) and δ k1y,k2y δ k1z,k2z (3D), since the impurity is independent of y in 2D and of y and z in 3D, and therefore, the momenta in the corresponding directions are conserved in all scattering processes. The exact same formalism can be applied for impurities perpendicular to the other axes of the systems.
In what follows we employ this formalism at zero temperature to calculate the retarded GF G(k 1 , k 2 , E) obtained by the analytical continuation of the Matsubara GF G(k 1 , k 2 , iω n ) (i.e., by setting iω n → E + iδ, with δ → 0 + ).
For a three-dimensional system the surface Green's function can be extracted from the perturbed generalized Green's function in Eq. (1). This can be related to the mixed Green's function in which we keep the momentum coordinates in the two directions parallel to the impurity plane (k y and k z ), but we perform a Fourier transform to write down the Green's function in real space coordinates in the x direction. Note that for a plane impurity the generalized Green's function depends on two different values of momentum only for the direction perpendicular to the impurity, in the other two directions we recover a simple dependence on momentum due to unbroken translational invariance. Thus we have: We fix x = x = ±1 since we are interested in describing the lattice planes one lattice constant away from the impurity (see Fig. 1). The boundaries of the two resulting semi-infinite systems correspond to the two planes at x = ±1, and the Green's functions taken at these two positions are effectively the surface Green's functions for the semi-infinite systems. The physics at x = 0 (impurity position) is relatively trivial since the infinite-amplitude impurity potential pushes away all the wave function weight off the impurity plane. For simplicity we have omitted writing down explicitly the energy dependence of the Green's functions. Note again that translational invariance holds within the planes parallel to the impurity plane, thus the surface Green's functions depend only on one momentum in each of the in-plane directions. Furthermore, in order to obtain the surface physics the value of the impurity potential in Eq. (4) needs to be set to a value much larger than all the energy scales in the problem. The surface Green's function allows us to recover the formation of the surface states such as, for instance, the Fermi-arc states. Thus, we can study the surface spectral function The same analysis can be performed for a twodimensional system with a line-impurity to find the line Green's functions and the corresponding edge states given by: Alternatively, in order to visualize the impurityinduced states, as described in Appendix B, we may focus on the average correction to the spectral function: where Above G 0 (k) stands for G 0 (k, E) and T (k, k) for T (k, k, E). The integral over k x is performed along the same interval as the one defined in Eq. (3). This quantity corresponds to the average number of available electronic states with wavevector (k y , k z ), where the average is performed along the direction perpendicular to the impurity. A more detailed description of the significance of this quantity is provided in Appendix B.

III. WEYL SEMIMETALS
In what follows we consider a Weyl semimetal: there is a large number of models of various degree of complexity describing such a system. Here we focus only on the tight-binding models described in Refs. [15] and [16], which we denote by H 1 and H 2 , respectively. The Bloch Hamiltonians for these two systems are given by where ψ(k) = (c kA↑ , c kA↓ , c kB↑ , c kB↓ ) is a spinor with the index A/B denoting a generic unspecified orbital component, and the ↑ / ↓ the physical spin. For the model in Ref. [15] written in the basis above we have where g 0 (k) = 2d(2 − cos k x − cos k y ) g 1 (k) = a sin k x g 2 (k) = a sin k y g 3 (k) = m + t cos k z + 2b(2 − cos k x − cos k y ). (13) and α, β are real parameters. The 2 × 2 identity matrices σ 0 /τ 0 and the Pauli matrices σ i /τ i , i = 1, 2, 3 act in the spin and the orbital spaces, correspondingly, and the multiplication of the σ and τ matrices indicates a tensor product. We consider the same values of parameters as those in Ref. [15], thus we take a) a = b = 1, t = −1, m = 0.5, d = 0.8, α = β = 0 and b) a = b = 1, t = −1.5, d = m = 0, β = 0.9, and α = 0.3. The former is characterized by two Weyl points, while the latter by four Weyl 3 Alternatively, in order to visualize the impurityinduced states, as described in Appendix B, we may focus on the average correction to the spectral function: where Above G 0 (k) stands for G 0 (k, E) and T (k, k) for T (k, k, E). The integral over k x is performed along the same interval as the one defined in Eq. (3). This quantity corresponds to the average number of available electronic states with wavevector (k y , k z ), where the average is performed along the direction perpendicular to the impurity. A more detailed description of the significance of this quantity is provided in Appendix B.

III. WEYL SEMIMETALS
In what follows we consider a Weyl semimetal: there is a large number of models of various degree of complexity describing such a system. Here we focus only on the tight-binding models described in Refs. [14] and [15], which we denote by H 1 and H 2 , respectively. The Bloch Hamiltonians for these two systems are given by where (k) = (c kA" , c kA# , c kB" , c kB# ) is a spinor with the index A/B denoting a generic unspecified orbital component, and the " / # the physical spin. For the model in Ref. [14] written in the basis above we have where g 0 (k) = 2d(2 cos k x cos k y ) g 1 (k) = a sin k x g 2 (k) = a sin k y g 3 (k) = m + t cos k z + 2b(2 cos k x cos k y ). (13) and ↵, are real parameters. The 2 ⇥ 2 identity matrices 0 /⌧ 0 and the Pauli matrices i /⌧ i , i = 1, 2, 3 act in the spin and the orbital spaces, correspondingly, and the multiplication of the and ⌧ matrices indicates a tensor product.
We consider the same values of parameters as those in Ref. [14], thus we take a) a = b = 1, t = 1, m = 0.5, d = 0.8, ↵ = = 0 and b) a = b = 1, t = 1.5, d = m = 0, = 0.9, and ↵ = 0.3. The former is characterized by two Weyl points, while the latter by four Weyl To make an exact correspondence with the spinless results in Ref. [14], in a) the trace is taken only over the spinup (first and third) components of the Green's function. We clearly see that there is a single Fermi arc emerging in a), whereas there are two Fermi arcs in b). We set U = 100.
points, and thus we expect to have one and two Fermi arcs, respectively.
In order to obtain the Fermi-arc surface states we need to introduce a surface into the system in such a way that the vector connecting the Weyl nodes has a nonzero projection onto it. For example, for the above model we choose to have a plane-like impurity at y = 0, hence perpendicular to the y direction. The resulting surface Green's functions are described by the formalism in Sec. II, where we consider an impurity V = U (y)I 4 , with U ! 1 (i.e, much larger than all energy scales in the problem). Note that in Sec. II we describe an impurity at x = 0 and not y = 0, however the y = 0 formalism is obtained by simply interchanging x and y in the corresponding formulas. The spectral function for the surface states Fig. 2 for two chosen configurations of parameters.
We note that the Fermi-arc states calculated using our method agree exactly with those predicted in Ref. [14]. Moreover, the full surface Green's function that we have obtained contains all the information required to describe these states, such as their spin and orbital distribution, the full energy dispersion, etc.. For instance, we have for 3 Alternatively, in order to visualize the impurityinduced states, as described in Appendix B, we may focus on the average correction to the spectral function: where Above G 0 (k) stands for G 0 (k, E) and T (k, k) for T (k, k, E). The integral over k x is performed along the same interval as the one defined in Eq. (3). This quantity corresponds to the average number of available electronic states with wavevector (k y , k z ), where the average is performed along the direction perpendicular to the impurity. A more detailed description of the significance of this quantity is provided in Appendix B.

III. WEYL SEMIMETALS
In what follows we consider a Weyl semimetal: there is a large number of models of various degree of complexity describing such a system. Here we focus only on the tight-binding models described in Refs. [14] and [15], which we denote by H 1 and H 2 , respectively. The Bloch Hamiltonians for these two systems are given by where (k) = (c kA" , c kA# , c kB" , c kB# ) is a spinor with the index A/B denoting a generic unspecified orbital component, and the " / # the physical spin. For the model in Ref. [14] written in the basis above we have where g 0 (k) = 2d(2 cos k x cos k y ) g 1 (k) = a sin k x g 2 (k) = a sin k y g 3 (k) = m + t cos k z + 2b(2 cos k x cos k y ). (13) and ↵, are real parameters. The 2 ⇥ 2 identity matrices 0 /⌧ 0 and the Pauli matrices i /⌧ i , i = 1, 2, 3 act in the spin and the orbital spaces, correspondingly, and the multiplication of the and ⌧ matrices indicates a tensor product.
We consider the same values of parameters as those in Ref. [14], thus we take a) a = b = 1, t = 1, m = 0.5, d = 0.8, ↵ = = 0 and b) a = b = 1, t = 1.5, d = m = 0, = 0.9, and ↵ = 0.3. The former is characterized by two Weyl points, while the latter by four Weyl To make an exact correspondence with the spinless results in Ref. [14], in a) the trace is taken only over the spinup (first and third) components of the Green's function. We clearly see that there is a single Fermi arc emerging in a), whereas there are two Fermi arcs in b). We set U = 100.
points, and thus we expect to have one and two Fermi arcs, respectively.
In order to obtain the Fermi-arc surface states we need to introduce a surface into the system in such a way that the vector connecting the Weyl nodes has a nonzero projection onto it. For example, for the above model we choose to have a plane-like impurity at y = 0, hence perpendicular to the y direction. The resulting surface Green's functions are described by the formalism in Sec. II, where we consider an impurity V = U (y)I 4 , with U ! 1 (i.e, much larger than all energy scales in the problem). Note that in Sec. II we describe an impurity at x = 0 and not y = 0, however the y = 0 formalism is obtained by simply interchanging x and y in the corresponding formulas. The spectral function for the Fig. 2 for two chosen configurations of parameters.
We note that the Fermi-arc states calculated using our method agree exactly with those predicted in Ref. [14]. Moreover, the full surface Green's function that we have obtained contains all the information required to describe these states, such as their spin and orbital distribution, the full energy dispersion, etc.. For instance, we have for To make an exact correspondence with the spinless results in Ref. [15], in a) the trace is taken only over the spinup (first and third) components of the Green's function. We clearly see that there is a single Fermi arc emerging in a), whereas there are two Fermi arcs in b). We set U = 100.
points, and thus we expect to have one and two Fermi arcs, respectively.
In order to obtain the Fermi-arc surface states we need to introduce a surface into the system in such a way that the vector connecting the Weyl nodes has a nonzero projection onto it. For example, for the above model we choose to have a plane-like impurity at y = 0, hence perpendicular to the y direction. The resulting surface Green's functions are described by the formalism in Sec. II, where we consider an impurity V = U δ(y)I 4 , with U → ∞ (i.e, much larger than all energy scales in the problem). Note that in Sec. II we describe an impurity at x = 0 and not y = 0, however the y = 0 formalism is obtained by simply interchanging x and y in the corresponding formulas. The spectral function for the surface states Fig. 2 for two chosen configurations of parameters.
We note that the Fermi-arc states calculated using our method agree exactly with those predicted in Ref. [15]. Moreover, the full surface Green's function that we have obtained contains all the information required to describe these states, such as their spin and orbital distribution, the full energy dispersion, etc.. For instance, we have for the different spin and orbital components: where we omit the energy dependence for the sake of brevity. In Fig. 3 we plot the x and z spin components of the Fermi-arc states at zero energy, separately calculated for the A and B orbitals. The parameters chosen correspond to Fig. 2 b). We do not plot the y component since it is zero for both orbitals. The spins of opposite arcs are of opposite signs, as expected. 15 We perform a similar analysis on a different Weyl semimetal model, introduced in Ref. [16]: We consider the same values of parameters as those in Ref. [16], thus we take a) a = b = 1, t = −1.5, λ = the di↵erent spin and orbital components: where we omit the energy dependence for the sake of brevity. In Fig. 3 we plot the x and z spin components of the Fermi-arc states at zero energy, separately calculated for the A and B orbitals. The parameters chosen correspond to Fig. 2 b). We do not plot the y component since it is zero for both orbitals. The spins of opposite arcs are of opposite signs, as expected. 14 We perform a similar analysis on a di↵erent Weyl semimetal model, introduced in Ref. [15]: Ref. [15] indicates the possible existence of an electron pocket coming from the bulk bands. We apply the same techniques as above, and we show in Fig. 4 the resulting spectral function A(k x , k z , E) for the surface states, for the two chosen configurations of parameters.
Our results agree exactly with those in Ref. [15]. Furthermore, we compute the spin and orbital properties for this model (see Fig. 5).
We have also checked that for the first set of parameters in Fig. 4 the spins of the two Fermi arcs are opposite, same as for the H 1 model. The H 1 and H 2 models di↵er in this case mainly by a nonzero y component in the H 2 model and the asymmetry of the two Fermi arcs in k z .

IV. KANE-MELE MODEL
We start with the Kane-Mele model of a topological insulator on a honeycomb lattice. 2 Therefore, we employ the following tight-binding model: where c † i,↵ creates an electron on the lattice site i, with spin ↵ =", #. The first term in Eq. (16) is the standard nearest-neighbor hopping term corresponding to the tight-binding Hamiltonian of graphene, which yields a spectrum with bands touching at the Dirac points situated at the Brillouin zone corners (±4⇡/3 p 3, 0), (±2⇡/3 p 3, ±2⇡/3). In order to turn graphene into an insulator we add a next-nearest-neighbor term with a the di↵erent spin and orbital components: where we omit the energy dependence for the sake of brevity. In Fig. 3 we plot the x and z spin components of the Fermi-arc states at zero energy, separately calculated for the A and B orbitals. The parameters chosen correspond to Fig. 2 b). We do not plot the y component since it is zero for both orbitals. The spins of opposite arcs are of opposite signs, as expected. 14 We perform a similar analysis on a di↵erent Weyl semimetal model, introduced in Ref. [15]: Ref. [15] indicates the possible existence of an electron pocket coming from the bulk bands. We apply the same techniques as above, and we show in Fig. 4 the resulting spectral function A(k x , k z , E) for the surface states, for the two chosen configurations of parameters.
Our results agree exactly with those in Ref. [15]. Furthermore, we compute the spin and orbital properties for this model (see Fig. 5).
We have also checked that for the first set of parameters in Fig. 4 the spins of the two Fermi arcs are opposite, same as for the H 1 model. The H 1 and H 2 models di↵er in this case mainly by a nonzero y component in the H 2 model and the asymmetry of the two Fermi arcs in k z .

IV. KANE-MELE MODEL
We start with the Kane-Mele model of a topological insulator on a honeycomb lattice. 2 Therefore, we employ the following tight-binding model: where c † i,↵ creates an electron on the lattice site i, with spin ↵ =", #. The first term in Eq. (16) is the standard nearest-neighbor hopping term corresponding to the tight-binding Hamiltonian of graphene, which yields a spectrum with bands touching at the Dirac points situated at the Brillouin zone corners (±4⇡/3 p 3, 0), (±2⇡/3 p 3, ±2⇡/3). In order to turn graphene into an insulator we add a next-nearest-neighbor term with a Ref. [16] indicates the possible existence of an electron pocket coming from the bulk bands. We apply the same techniques as above, and we show in Fig. 4 the resulting spectral function A(k x , k z , E) for the surface states, for the two chosen configurations of parameters.
Our results agree exactly with those in Ref. [16]. Furthermore, we compute the spin and orbital properties for this model (see Fig. III).
We have also checked that for the first set of parameters in Fig. 4 the spins of the two Fermi arcs are opposite, same as for the H 1 model. The H 1 and H 2 models differ in this case mainly by a nonzero y component in the H 2 model and the asymmetry of the two Fermi arcs in k z .

IV. KANE-MELE MODEL
We start with the Kane-Mele model of a topological insulator on a honeycomb lattice. 2 Therefore, we employ the following tight-binding model: where c † i,α creates an electron on the lattice site i, with spin α =↑, ↓. The first term in Eq. (16) is the standard nearest-neighbor hopping term corresponding to the tight-binding Hamiltonian of graphene, which yields a spectrum with bands touching at the Dirac points situated at the Brillouin zone corners (±4π/3 √ 3, 0), (±2π/3 √ 3, ±2π/3). In order to turn graphene into an insulator we add a next-nearest-neighbor term with a spin-dependent amplitude ν ij = −ν ji = ±1, defined by the orientation of the hopping direction (see Fig. 6). The second term opens a bulk gap in the energy spectrum at the Dirac points.
First, we obtain the boundary modes numerically by diagonalizing the tight-binding Hamiltonian in Eq. (16) and considering periodic boundary conditions in the y direction and open boundary conditions in the x direction. This corresponds to a ribbon with zigzag edges. For convenience we set the lattice spacing a to unity. The corresponding energy spectrum is shown in Fig. 7. Note the formation of two subgap states (we have verified that these are actually edge states). The momentum-space dispersions for the two edge states cross at zero energy for k y = π/ √ 3. In what follows we reproduce the formation of these edge states by considering a line impurity in an infinite system and subsequently taking the impurity potential with 2 ky e −i 3 2 kx , and being the nearest-neighbor and the next-nearest-neighbor terms with amplitudes t and t 2 , respectively. Here * simply denotes the complex conjugation.
To reproduce the zigzag edge states, we choose an impurity potential localized on two adjacent rows of atoms corresponding to two different sublattices.
(see Appendix D for more details). In order to visualize the impurity-induced states, in Fig. 8 we plot the average correction to the spectral function due to the impurity, as defined in Eq. (9), for the same range of k y as the one considered in Fig. 7. For a weak impurity, i.e., U = t = 1, the impurity states appear as a distinct band at energies concentrated mostly outside of the gap. We expect that the impurity bound states will evolve into edge states and acquire the same properties (i.e., the same momentum dispersion) as the edge states derived previously using numerical methods, when the impurity strength U goes to infinity. Indeed, for a stronger impurity potential U = 100, the impurity-induced spectral function exhibits subgap states with the same dispersion as the ones derived via exact diagonalization and depicted in Fig. 7. The agreement between the two methods is remarkable, confirming the validity of our analytical approach towards finding the edge states of a simple topological insulator system.
While here we consider only zigzag edges, in Appendix D we have also considered the case of an impurity localized only on one row of atoms, which splits the systems into two subsystems with different edges, one with a zigzag edge, and one with a bearded edge. We expect that we will recover two distinct sets of edge states, and in Appendix D we show that this is indeed the case.

V. CHERN INSULATOR
Below we consider the simplest lattice model defining a Chern insulator where we set the lattice constant to unity and t = 1.
Here k ≡ (k x , k y ) and σ = (σ x , σ y , σ z ) are the Pauli matrices. The subspace in which they act may be very general and depends on the given model, for example in a lattice model with two orbitals per site, s and p, the σ matrices act in the orbital subspace. The above model yields topologically nontrivial phases for M ∈ (0, 2) ∪ (2, 4) (see Ref. [3]).
In what follows we introduce a line-like impurity at x = 0 described by the potential V imp (x) = V δ(x)I, with V → ∞ and I is the 2 × 2 identity matrix. In this limit the T -matrix can be written as: We compute the integral in Eq. (21), setting B = M = 1 for the sake of simplicity. We note that the calculation can be performed for arbitrary values of B and M . Thus, we have π −π dk x 2π G 0 (k x , k y , iω n ) = 1 4 sin 2 ky 2 1 (iω n ) 2 − 1 (iω n ) 2 − 5 + 4 cos k y × × (iω n ) 2 − 2 cos k y + cos 2k y − (iω n ) 2 − 1 (iω n ) 2 − 5 + 4 cos k y σ z + 4 sin 2 k y 2 (iω n σ 0 + sin k y σ y ) . Note that the poles of the T -matrix obtained by taking iω n → E + iδ, with δ → +0, are given by E = ± sin k y for k y ∈ [−π/2, π/2], corresponding to two chiral edge modes (cf. Ref. [21]). For the sake of brevity, we leave the derivation of the poles of the T -matrix to Appendix C.
To verify our findings, in Fig. 9 we plot the average correction to the spectral function defined in Eq. (9) as a function of E and k y . As expected, on one hand we can see the bulk states, originating from the poles of the bare Green's function i.e., the eigenvalues of the Hamiltonian in Eq. (20) which for the values considered here, B = M = t = 1 and k x = 0 correspond to E = ±1. More importantly we can identify also the two counterpropagating chiral edge modes of the Chern insulator crossing at k y = 0, whose dependence on k y is consistent with the fully analytical form above. This demonstrates the strength of our approach to recover fully analytical results for the edge states of certain models for which the unperturbed Green's function in the real space can be obtained in an analytical closed form.

VI. CONCLUSIONS
We have generalized the technique to obtain boundary modes, introduced in Ref. [12], to calculate the surface Green's functions of an arbitrary three-dimensional system, and we have applied it to calculate the surface Green's functions for Weyl semimetals and recovered the corresponding Fermi-arc states. We have also shown that the technique in Ref. [12] can be applied to other topological systems, such as topological insulators. Furthermore, we have demonstrated that it functions also for systems with more than one sublattice, and that this method can be easily employed to study boundary modes in lower dimensions. In particular, using line-like impurities, we have applied the formalism to derive the helical edge states of the Kane-Mele model and the chiral edge states of a Chern insulator. For the latter we have shown that using this formalism a full analytical form can be obtained for the T -matrix, and analytical expressions for the energies of the edge states can be recovered.
Appendix B: Significance of a position-averaged spectral function In 2D, in the presence of a line impurity proportional to δ(x), the correction to the number of available electronic states at position x and having momentum k y is given by: δρ(x, k y , E) = − 1 π Im tr δG(x, x; k y ; E) (B1) Note that since the spatial translational invariance along y is not broken k y is still a good quantum number. Averaging this over the direction perpendicular to the impurity we obtain: where L k and the limits of integration are −π to π with L k = 2π for a square lattice and −2π/3 to 2π/3 with L k = 4π/3 for a honeycomb one. Also is the correction to the perturbed spectral function in momentum space.

Appendix C: T -matrix poles for the Chern insulator
In order to calculate the energies of the bound states, here we calculate analytically the poles of the T -matrix defined by Eqs. (21)(22). The latter can be found from the trace of the T -matrix given by tr T (k y , E + iδ) = − (E + iδ) (E + iδ) 2 − 2 cos k y + cos 2k y + (E + iδ) 2 − 1 (E + iδ) 2 − 5 + 4 cos k y where we replaced iω n → E + iδ, with δ → +0. We obtain straightforwardly the zeros of the denominator, namely, E = ± sin k y . However, to make sure that the latter are poles, we need to verify that they are not zeros of the numerator. The trivial zero of the numerator is E = 0, we discard it below. Thus, we need to analyze the zeros of the expression in the square brackets: (E + iδ) 2 − 2 cos k y + cos 2k y + (E + iδ) 2 − 1 (E + iδ) 2 − 5 + 4 cos k y = 0 (C2) We represent the complex numbers under the square roots in the trigonometric form, and applying the limit δ → +0 we get: where we defined φ 1 (E) = π 2 sgn E Θ 1 − E 2 , φ 2 (k y , E) = π 2 sgn E Θ 5 − 4 cos k y − E 2 . Since we are searching for subgap solutions, i.e., |E| < 1, then E 2 < 5 − 4 cos k y ∀k y ∈ [−π, π]. Thus, we have φ 1 (E) + φ 2 (k y , E) = π sgn E and, therefore, e i[φ1(E)+φ2(ky,E)] = e iπ sgn E = −1. Eq. (C3) then becomes E 2 − 2 cos k y + cos 2k y − 1 − E 2 5 − 4 cos k y − E 2 = 0. (C4) The equation above is equivalent to solving the system: E 2 − sin 2 k y sin 4 k y 2 = 0 (C5) E 2 − 2 cos k y + cos 2k y > 0 (C6) When k y = 0, we get E = ±1 at the edge of the gap, therefore, E = ± sin k y from the first equation. The second equation then yields: sin 2 k y − cos k y + cos 2k y > 0 ⇒ |k y | ∈ π 2 , π . (C7) Thus, the numerator of Eq. (C1) has zeros at E = 0, and E = ± sin k y when |k y | ∈ (π/2, π], and therefore, we conclude that the trace of the T -matrix has poles at E = ± sin k y only when |k y | π/2. This means that the edge modes exist only for k y lying in the interval |k y | π/2, and their dispersion is given by E = ± sin k y .
Appendix D: Impurity potentials for the Kane-Mele model While working with the Kane-Mele model, we considered three different delta-function impurities which we illustrate in Fig. 10 with the associated matrix representation. We see that in order to reproduce the zigzag edge states of the Kane-Mele model we must introduce an impurity on both A and B sites (right panel of Fig. 10). If the impurity is localized only on the A or B sites, it will create both a zigzag edge state and a bearded edge state. In each case, the impurity creates a "wall" in the system (dashed lines or shaded area) along with two boundaries. If the impurity is localized on a single sublattice (A or B), it creates one zigzag edge and one bearded edge. If it is localized on entire unit cells it will create two zigzag edges. The matrix representation in the insets is given in the basis (c A↑ , c A↓ , c B↑ , c B↓ ).
The formation of both zigzag and bearded edges can be recovered using our method by applying either of the potentials given in the left and middle panels of Fig. 10. Fig. 11 shows the energy spectrum obtained by exact diagonalization of the Hamiltonian on a strip with one zigzag edge and one bearded edge, along with the correction to the averaged spectral function due to a line impurity localized on either one of the two sublattices. Here we have enlarged the horizontal axis to k y ∈ [−2π/ √ 3; 2π/ √ 3] to see the edge states more clearly. We recover the same dispersion as in the main text for the zigzag edge states and obtain in addition other subgap states which are localized on the bearded edge. We need however to keep in mind the fact that in this case the two semi-infinite systems are not fully decoupled since the bulk Hamiltonian of the Kane-Mele model contains spin-flip NNN terms: this leakage effect needs to be taken into account carefully especially when we consider the spin-properties of the systems, but will not affect the spectrum below.