{'fraction_non_alphanumeric': 0.045547396897256114, 'fraction_numerical': 0.0131765398700749, 'mean_word_length': 3.9584816132858838, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing edge states without hard truncation

We present a numerical method which accurately computes the discrete spectrum and associated bound states of Hamiltonians which model electronic"edge"states localized at boundaries of one and two-dimensional crystalline materials. The problem is non-trivial since arbitrarily large finite"hard"truncations of the Hamiltonian in the infinite bulk direction tend to produce spurious bound states partially supported at the truncation. Our method, which overcomes this difficulty, is to compute the Green's function of the Hamiltonian by imposing an appropriate boundary condition in the bulk direction; then, the spectral data is recovered via Riesz projection. We demonstrate our method's effectiveness by studies of edge states at a graphene zig-zag edge in the presence of defects modeled both by a discrete tight-binding model and a continuum PDE model under finite difference discretization. Our method may also be used to study states localized at domain wall-type edges in one and two-dimensional materials where the edge Hamiltonian is infinite in both directions; we demonstrate this for the case of a tight-binding model of distinct honeycomb structures joined along a zig-zag edge.

1. Introduction. When computing the electronic states of crystalline materials it is typical to first treat the material as infinite in all directions of periodicity, leading to an periodic bulk Hamiltonian which is naturally analyzed using the Bloch transform (see [3], for example). The eigenstates of this Hamiltonian, known as bulk states, are quasi-periodic with respect to the crystal lattice and hence extensive throughout the material. Materials may also host electronic states which decay rapidly away from the physical edge of the material known as edge states. Edge states are bound states of the semi-infinite Hamiltonian (known as the edge Hamiltonian) obtained by truncating the infinite bulk Hamiltonian in one direction and imposing a Dirichlet boundary condition at the truncation. Such states play an important role in the theory of the quantum Hall effect and of topological insulators, and are of independent interest for wave-guiding applications because of their robustness to certain classes of local perturbations [4,14,16,19,23].
In this work we propose a numerical method for computing the discrete spectrum and associated bound states of such semi-infinite edge Hamiltonians. The problem is non-trivial by the fact that arbitrarily large finite truncations of the Hamiltonian produce spurious bound states partially supported at the truncation which do not correspond to bound states of the original Hamiltonian. The key idea of our method, which overcomes this difficulty, is to instead compute the Green's function of the edge Hamiltonian along an appropriate contour by truncating the problem while imposing an appropriate boundary condition at the truncation. The spectrum and eigenstates of the original Hamiltonian may then be recovered using the Riesz projection formula. Our method relies on eventual periodicity of the Hamiltonian which implies that exact boundary conditions for computing the Green's function can be derived from a transfer matrix; similar ideas have been explored for numerical computation of local defects in electronic structures in [18], and for analytical calculations e.g., [9,24]. We provide self-contained proofs of the existence of exact boundary conditions for computing the Green's function for arbitrary finite range discrete Hamiltonians and that all spectral data can be recovered exactly from the Green's function. Hence our method is exact up to rounding errors for discrete models and exact up to discretization error of continuum models.
The rest of the paper is organized as follows. In Section 2 we motivate our method by demonstrating the failure of a naïve method using a one-dimensional SSHtype model which has the essential features of the general case. In Section 3 we introduce our method and show that it accurately computes spectral data of this model. In Section 4 we generalize our method to a class of models whose Hamiltonians are "eventually periodic", and prove the existence of exact boundary conditions for computing the Green's function: first under a mild simplifying assumption on the Hamiltonian (Section 4.1) and then in generality for finite range Hamiltonians (Appendices A and B). In Section 5 we present results of computations using our method of states localized at edges of two-dimensional honeycomb structures modeled both in the tight-binding limit and by a continuum Schrödinger equation under finite difference discretization.
Acknowledgements: This work is supported in part by National Science Foundation through grants DMS-1454939 and ACI-1450280 and Department of Energy through grant DE-SC0019449. The authors would also like to thank Jacob Shapiro and Michael I. Weinstein for stimulating discussions. Early stages of this project were carried out by Qinyi Zhu as an undergraduate researcher visiting Duke University during the Summer of 2017.
2. Spurious bound states in a one-dimensional model: the SSH model. In this section, in order to motivate and clearly lay out the ideas behind our method, we restrict attention to a discrete one-dimensional model which has the essential features of the general case: the Su-Schrieffer-Heeger (SSH) model with real, nearestneighbor hopping [27]. The structure of this section is as follows.
In Section 2.1, we review the definition and spectrum (known as the bulk spectrum) of the SSH bulk Hamiltonian. In Section 2.2, we introduce an SSH edge Hamiltonian with defects whose essential spectrum is equal to that of the bulk Hamiltonian but whose discrete spectrum and associated bound states are difficult to compute by hand. In Section 2.3 we then demonstrate the appearance of spurious bound states in naïve numerical computations.
The stage will then be set for us to present in Section 3 our main method which we refer to as the Green's function method. This method accurately computes the discrete spectrum and associated bound states of the SSH edge Hamiltonian with defects within the gap between essential spectrum.
2.1. SSH bulk Hamiltonian and bulk spectrum. In this section, we review the definition and spectrum of the SSH bulk Hamiltonian. We consider an electron (we ignore spin) "hopping" along a one-dimensional infinite periodic lattice (chain) with hopping amplitudes that alternate between sites along the chain. A fundamental cell of the chain consists of any two neighboring sites, which we label A and B. We introduce the notation to represent the restriction of the electron wavefunction ψ ∈ l 2 (Z; C 2 ) to the mth fundamental cell of the chain. Denoting the intra-and inter-cell hopping amplitudes by t 1 and t 2 (for simplicity, we ignore all hopping other than between nearest-neighbors and assume that t 1 and t 2 are real and non-zero), the bulk Hamiltonian H bulk is Periodicity of the bulk Hamiltonian (2.2) implies, via Bloch's theorem (see [3,25], for example), that all (generalized) eigenfunctions Φ of H bulk are k-quasi-periodic: . A standard calculation then shows that the associated eigenvalues of these eigenfunctions are given by The spectrum of (2.2), σ(H bulk ), is therefore the union of the real intervals (known as spectral bands) swept out by the functions (2.3) as k is varied over the interval [−π, π]. If |t2| |t1| = 1, then σ(H bulk ) consists of a single interval. If |t2| |t1| = 1, then σ(H bulk ) consists of two intervals, separated by a gap which is symmetric about 0.
The (generalized) eigenvectors of H bulk corresponding to the eigenvalues E ± (k) are known as bulk states. These states do not decay as m → ±∞.
2.2. SSH edge Hamiltonian and edge spectrum. We define a SSH left edge Hamiltonian which acts on l 2 (N; C 2 ) as follows Here and in what follows we adopt the convention that the natural numbers start at 1, i.e. N = {1, 2, ...}. Note that we allow the edge to have defects, modeled by (non-zero) hopping amplitudes whose values depend on the cell index m, and for possibly non-zero real onsite potentials V σ (m), σ ∈ {A, B}. We make the following assumption on the behavior of these functions for large m.
Assumption 2.1 (Periodicity of bulk medium). There exists a positive integer M ≥ 1 such that for all m ≥ M : where H bulk is defined by (2.2) with t 1 and t 2 replaced by t ∞ 1 and t ∞ 2 . Regarding the discrete spectrum of (2.4) we have the following standard result on existence of bound states with eigenvalue E = 0.
where c is an arbitrary complex constant. If instead The bound state (2.6) is known as an edge state because it decays away from the edge of the material. The existence of edge states can be related to topological invariants associated with the associated bulk Hamiltonian in a link known as the bulk-boundary correspondence (see [2,7,12,16], and references within, for example). We will adopt the convention that any bound state of the edge Hamiltonian is an edge state even if it is not associated with any bulk invariant.
An SSH right edge Hamiltonian H R edge can be defined on l 2 (−N; C 2 ) analogously to H L edge . Under an analogous assumption to Assumption 2.1 it is clear that σ ess (H R edge ) = σ(H bulk ) and that H R edge has a bound state with eigenvalue 0 whenever |t2| |t1| > 1 and V σ (m) = 0 for σ ∈ {A, B} and all m ∈ −N. While (2.6) is entirely supported on A sites, the bound state of H R edge with eigenvalue 0 is entirely supported on B sites. We emphasize that when V A (m) and V B (m) are non-zero and t 1 (m) and t 2 (m) are not constant, H L edge may have multiple eigenvalues in the gap between essential spectrum and 0 may or may not be one of them (see Figure 3.1). In the next section we consider a naïve approach to computing all such eigenvalues of (2.4) and their associated bound states.
2.3. Spurious bound states associated with finite "hard" truncation of the semi-infinite Hamiltonian. Since bound states of (2.4) must decay as m → ∞, it is tempting to try to approximate such bound states and their associated eigenvalues by imposing a Dirichlet boundary condition at some large value of the cell index m, i.e. by computing solutions of the truncated eigenequation for some choice of M M , where M is the integer appearing in Assumption 2.1 and where H L trunc is defined by where V (m) and A(m) are as in (2.4). We call this method the hard truncation method. This approach produces inaccurate results because H L trunc can have bound states partially supported at the truncation for arbitrarily large M; see Remark 2.1. The problem described above can be avoided when |t ∞ 2 | |t ∞ 1 | > 1 by truncating the chain in the middle of a cell (i.e. after an A site) rather than at the end . The situation is analogous to the case of two identical wellseparated potential wells in the semi-classical regime [15,26] where one observes splitting (width O(e −cγ ) where γ is the distance between wells and c > 0 is a constant) of the single potential well eigenvalues.
of a cell (i.e. after a B site). The right edge is then well described by H R edge , the right-edge Hamiltonian, but with t 1 replaced by t ∞ 2 and t 2 replaced by t ∞ 1 which has no edge state in this case. The method we describe below has the advantage that it may be applied without modification to all possible choices of t ∞ 1 and t ∞ 2 equally well and also generalizes naturally to the more general cases we consider below.
In the next section we introduce a method for computing spectral data of H L edge which overcomes this problem.
3. Green's function method. In this section we introduce our main method which we refer to as the Green's function method. Rather than trying to compute eigenvectors of H L edge directly we compute the associated spectral projection operator. This is done by computing the resolvent (Green's function) of the Hamiltonian and then recovering the projection via the Riesz projection formula.
is an orthogonal projection operator; in particular, it is the spectral projection onto the associated eigenspace of the eigenvalues enclosed in the contour C.
Note that our method only yields information about eigenvalues within the contour C. In practice this is not a problem: typically one is only interested in eigenvalues within the gap between essential spectrum (recall that for the SSH model σ ess (H L edge ) In spectral plots, essential spectrum (computed using Bloch theory) is shown as a blue line, while eigenvalues are labeled by ×s. The associated eigenvalue of each bound state which is plotted in the above figures is plotted with the same color in the figures below. In both cases, we drew the values of the hopping amplitudes t 1 (m) and t 2 (m) for m ≤ 40 from normal distributions with standard deviation 0.5 and means equal to t ∞ 1 and t ∞ 2 , with t ∞ 1 = 1 and t ∞ 2 = 1.6, respectively. To generate (a), we set all onsite potentials V σ (m) = 0. In this case, we find that (2.4) has an eigenvalue which is precisely 0 (confirming Proposition 2.1). In addition, we see bound states with non-zero eigenvalues whose associated eigenfunctions are localized away from the edge. To generate (b), we drew the values of the onsite potentials V A (m) and V B (m) for m < M from a normal distribution with standard deviation 0.5 and mean 0. In this case, we find again one bound state which decays rapidly away from the edge and multiple bound states whose associated eigenvectors are localized away from the edge. Note that the eigenvalue of the eigenvector which decays rapidly away from the edge is no longer 0 in this case. has a gap containing 0 whenever |t ∞ 2 | |t ∞ 1 | = 1). One therefore chooses the contour C in such a way as to encircle a large subinterval of the gap without intersecting with the essential spectrum.
In Section 3.1 we explain how to compute the upper left block of the Green's function of H L edge for any z / ∈ σ(H L edge ) and hence the upper left block of the spectral projection. In Section 3.2 we show how all relevant spectral information can be recovered using this method. In Figure 3.1 we present results computed using our method. It is natural to label the entries of G(z) similarly to those of the wavefunction ψ so that the mth 2×2 block of the part of G(z) which acts on the m th fundamental cell of ψ is written .
We now express (3.2) column by column as follows. For notational clarity, we suppress z-dependence of entries of G(z) and use the vector notation The defining equation of the Green's function (3.2) then implies that the σ , m th column of G(z) satisfies 5) and the condition that G σ m,m → 0 as m → ∞ sufficiently fast that {G σ m,m } m∈N ∈ l 2 (N, C 2 ). Here the matrices A(m) and V (m) are as in (2.4), and δ σ,σ m,m denotes the vector with 1 in its m , σ th entry and 0 in all its other entries.
For m ≤ M and m ≥ M + 2, the system (3.5) can be solved via a transfer matrix which depends on z Since the determinant of T ∞ (z) is 1 for all z, its spectrum is constrained as follows.
Note that whenever z is such that T ∞ (z) has an eigenvalue of norm 1 (i.e. either possibility (A2) or (A3) of Lemma 3.2 is realized) we can produce a Weyl sequence and hence z ∈ σ ess (H L edge ). Hence the following Lemma.
For the solution of (3.5) to decay as m → ∞ it must be that G σ M +1,m is proportional to the associated eigenvector of this eigenvalue. Hence the following.
,m , the M +1th entry of the σ , m th column of the Green's function is proportional to the associated eigenvector of the eigenvalue of T ∞ (z) with norm less than 1 whose existence is guaranteed by Lemma 3.3. From Lemma 3.4 we conclude that we can compute the first M + 1 entries of the solution of (3.5) (for m ≤ M ) by truncating the system and imposing the boundary condition where [ξ A (z), ξ B (z)] T denotes the eigenvector of T ∞ (z) whose associated eigenvalue has norm less than one (whose existence is guaranteed by Lemma 3.3) and c denotes an arbitrary constant. In this way we can compute the first M +1 entries of the (m , σ )th column of G(z) for each m ∈ {1, ..., M } and σ ∈ {A, B} for any z / ∈ σ(H L edge ). The choice to impose the boundary condition (3.7) in the M + 1th cell is to some extent arbitrary; the results of our method do not change if this boundary condition is imposed in any cell whose index is greater than M . However, choosing M + 1 clearly minimizes the size of the matrix problem to solve and is hence optimal. For a numerical verification that our results do not depend on this choice see Figure 3.2.
Remark 3.2. Under Assumption 2.1, the eigenequation H L edge ψ = Eψ can be solved directly using a transfer matrix similarly to the system (3.5). E is then an eigenvalue of H L edge whenever the solution of this system is such that the M th entry of ψ is proportional to an eigenvector of the transfer matrix whose associated eigenvalue has norm < 1. However this approach leads to a eigenvalue problem whose boundary condition (the analog of (3.7)) depends on E and is hence nonlinear. Such problems can be solved by iteration. We do not consider this idea further since it is inferior to the Green's function method.
It remains to recover the eigenvalues and eigenfunctions of H L edge from the upper left 2(M + 1) × 2M block of G(z); we discuss this in Section 3.2. In fact there is no loss in at this point discarding the 2(M + 1)th row of this block and in what follows we will consider only the upper left 2M × 2M block of G(z).

3.2.
Recovery of spectral data from the truncated Green's function. We now discuss how to recover the eigenvalues and eigenstates of H L edge from the truncated Green's function computed in the previous section. The upper left 2M × 2M block of the spectral projection Γ may be found by numerically evaluating the contour integral (3.1) using the upper-left 2M × 2M block of G(z) computed in the previous section (see [13] for methods which dramatically speed up the evaluation of this integral). For notational simplicity, let us define The truncated spectral projection, Γ trunc , is not sufficient for recovering all eigenvalues and eigenfunctions of H L edge . The missing information can be recovered by also using the upper-left 2M × 2M block of the operator H L edge Γ, which may be computed via a contour integral, similarly to Γ.
where the second equality is clear by analyticity of the operator (z − H L edge )G(z) = I. The upper left 2M × 2M block of H L edge Γ may now be computed similarly to that of Γ. We define Let ψ j , where j ∈ {1, ..., K}, denote the bound states of H L edge with associated eigenvalues E j within the contour γ and ψ j,trunc the truncations of these bound states to the first M cells. Then, using Dirac notation for simplicity, we have that Our goal now is to recover each E j and corresponding ψ j,trunc from the matrices Γ trunc and (H L edge Γ) trunc . We consider in detail the case K = 2; the result for general K is an easy generalization.
Let v 1 , v 2 denote any basis (not necessarily orthogonal) of the range of Γ trunc . Since the matrices Γ trunc and (H L edge Γ) trunc map this space to itself, we can represent their action on vectors in the space v = c 1 v 1 + c 2 v 2 by matrices as follows When we choose the basis v j = ψ j,trunc then using the definitions of Γ trunc and (H L edge Γ) trunc we have that Equation (3.13) implies that the eigenvalues E j are exactly the eigenvalues of the matrix M which relates the coefficients α ij and β ij , which satisfies where T denotes matrix transpose. The existence of such a matrix is guaranteed by the fact that Γ trunc is non-singular on its range, so that M can be defined by where −T denotes the inverse of the matrix transpose. With the eigenvalues E j in hand, the ψ j,trunc can then be found by computing each of the null spaces of the matrices (H L edge Γ) trunc − E j . 4. Green's function method for arbitrary eventually periodic finite range Hamiltonians. In this section we describe how the basic idea introduced in Section 3 generalizes in a natural way to a class of one-dimensional discrete Hamiltonians which are finite range and eventually periodic. This generalization is necessary for us to consider the models of honeycomb structure edges we consider in Section 5. We start with the following definition.
Definition 4.1. We call any operator H which is self-adjoint on l 2 (N; C N ) a semi-infinite one-dimensional discrete Hamiltonian. In addition: • We say H is finite range if there exists a non-negative integer R such that the action of H on elements ψ of l 2 (N; C N ) can be written in terms of the components ψ m ∈ C N of ψ as Here A * := A T denotes the conjugate transpose of the matrix A.
• We say H is eventually periodic if there exists a non-negative integer M such that for all m ∈ N such that m ≥ M . If H is both finite range and eventually periodic we will usually abbreviate and call such a Hamiltonian a semi-infinite eventually periodic Hamiltonian.
Note that the SSH edge Hamiltonian with nearest neighbor hopping which we considered in Sections 2 and 3 is a semi-infinite eventually periodic Hamiltonian with R = 1 and N = 2.
Remark 4.1. By enlarging the fundamental cell sufficiently we can always take R = 1 in (4.1). However, we prefer not to do this so that we keep 1 as the minimal period of the models we consider.
It is important to note that our method can be generalized further and applied to discrete Hamiltonians defined on Z under appropriate conditions. Definition 4.2. We call any operator H which is self-adjoint on l 2 (Z; C N ) an infinite one-dimensional discrete Hamiltonian. In addition: • We call H finite range if its action on elements ψ of l 2 (Z; C N ) can be written in terms of the components ψ m ∈ C N of ψ as for all m ∈ Z such that m ≤ M − . If H is both finite range and eventually periodic we will usually abbreviate and call such a Hamiltonian an infinite eventually periodic Hamiltonian.
In this case the method is almost identical except that when computing the Green's function the problem must be truncated in both directions with boundary conditions imposed at each end. We consider an example of this type in Section 5.3.
Since the method of recovering spectral data from the Green's function in general is identical to that given in Section 3.2 and the method is similar for infinite eventually periodic Hamiltonians we focus in this section on the problem of computing the upper left block of the Green's function of a general semi-infinite eventually periodic Hamiltonian.

Constructing boundary conditions for
Green's functions computation for semi-infinite eventually periodic Hamiltonian. For values of z ∈ C in the resolvent set, the Green's function G(z) of an arbitrary semi-infinite eventually periodic Hamiltonian H is defined as the bounded inverse of zI − H on l 2 (N; C N ), Just as in Section 3, we represent the mth N × N block of G(z) which acts on the m th fundamental cell as Equation (4.6) can then be written column-by-column as follows. We again suppress z-dependence of entries of G(z) and then assemble vectors We have then that the entries of the m , σ th column of the Green's function satisfy the system 9) and the condition that G σ m,m → 0 as m → ∞ sufficiently fast that {G σ m,m } m∈N ∈ l 2 (N; C 2 ). Here the matrices A * J (m) are those which appear in (4.1), and δ σ,σ m,m denotes the vector whose σ , m th entry is one and all others are zero. Since we require only the first M columns of the Green's function we need only solve (4.9) for m ≤ M .
For clarity we at this point restrict to the case R = 1 (see Remark 4.2 for the case where R > 1), and drop the subscript on A 1 → A so that (4.9) reduces to together with the condition that G σ m,m → 0 as m → ∞. Rather than study (4.10) directly, we will study the model difference equation for a general right hand side where our only assumption on the f m is that f m = 0 for m > M and we impose g m → 0 as m → ∞. Our goal is to find appropriate boundary conditions which will allow us to truncate (4.11) to a finite problem whose solution can be computed.
Under the eventual periodicity assumption (4.2), for m ≥ M +2 the system (4.11) can be written redundantly as Here V ∞ := V (M ) and A ∞ := A(M ). We will prove the following theorem, which generalizes the basic ideas developed in Section 3. what follows by V decaying , such that entries g 1 , g 2 , ..., g M +2 of the solution of the system (4.11) subject to {g m } m∈N ∈ l 2 (N; C N ) can be found by truncating the system after M + 2 cells and imposing the boundary condition that [g M +1 , g M +2 ] T ∈ V decaying .
The proof makes use of the following Lemma and hence whenever λ ∈ C \ 0 is an eigenvalue of T with multiplicity µ so is 1 λ . Proof of Lemma 4.4. The proof is by direct computation using standard identities satisfied by the determinant: We now prove Theorem 4.3.
Proof of Theorem 4.3. In the simplest case, S(z) is non-singular and we can obtain a transfer matrix relating g m−2 , g m−1 with g m , g m+1 By assumption, z / ∈ σ(H). If T ∞ (z) has an eigenvalue of magnitude 1 we can generate a Weyl sequence which is a contradiction. On the other hand, when z is real we may apply Lemma 4.4 to the matrix T ∞ (z) to see that T ∞ (z) must have the same number of (possibly generalized) eigenvectors corresponding to eigenvalues with norm less than 1 as corresponding to eigenvalues with norm greater than 1, and hence must have precisely N (possibly generalized) eigenvectors corresponding to eigenvalues with norm less than 1. The statement of the theorem now follows for such z by setting V decaying equal to the span of these vectors.
When z / ∈ σ(H) is complex, the conclusion still holds by the following argument. Consider a smooth path between any real z 0 / ∈ σ(H) and z which does not intersect with σ(H). Then the number of eigenvalues of T ∞ (z) with norm less than 1 cannot change along this path since if it did at the point of change an eigenvalue would have to lie on the unit circle in the complex plane. But we have a contradiction since this point must then be in the spectrum of H since we can produce a Weyl sequence.
Often, however, the matrix S(z) is singular. For the SSH model with nearestneighbor hopping for example, S(z) has the form which has an obvious null space spanned by [0, 0, 0, 1] T . In such a case the system can still be solved using a transfer matrix, but the derivation is more involved. We require the following assumption which is clearly satisfied by (4.17).
Assumption 4.1. We assume that if S(z) is not invertible, it has the form For simplicity we now suppress dependence on z and introduce the notation (4.22) A n := P Null S * g n g n+1 , B n := P ⊥ g n g n+1 , C n := P Null S g n g n+1 .
The system (4.12) can then be written row by row as  Assuming that B m−2 and C m−2 satisfy the solvability condition (4.24), we see that the system (4.23) uniquely specifies B m through (4.27), which in turn specifies A m through (4.25). In fact, C m is also uniquely specified by B m through the condition that the system (4.23) be solvable when m is replaced by m + 2, i.e. (4.29) It now follows that if B m−2 lies in the associated eigenspace of the eigenvalues of (4.28) with norm less than 1, then the system (4.23) has a solution which decays as m → ∞. In terms of the original notation, we have that as long as (1) g m−2 , g m−1 T satisfies the condition that and (2) P ⊥ g m−2 , g m−1 T lies in the image of the projection onto the space of associated eigenvectors of eigenvalues of T (4.28) with norm less than one, then the system (4.12) has a unique solution which lies in l 2 (N; C N ). We therefore define V decaying as the space of vectors satisfying both of these conditions. By Lemma 4.4 this space has dimension precisely N . The generalization to complex z / ∈ σ(H) follows from the same argument given for the case where S(z) is non-singular.
Remark 4.2. When R > 1, the discussion is the same except that the matrix S must be taken larger. For example when R = 2, we take which now acts on 4 cells of the column {g m } m∈N . Assumption 4.1 holds for all the models we consider in this paper; so, Theorem 4.3 gives a relatively simple and direct proof of our method for these cases. However, we note that Assumption 4.1 is not necessary. In Appendices A and B, we generalize our boundary condition method so that only two assumptions are required --the Hamiltonian must be (1) of finite range and (2) eventually periodic. The more general algorithm is slightly different than the one given above due to some special situations that can occur in the more general setting. See the appendices for the detailed general algorithm, as well as the full proof that the more general algorithm works in this general setting.
5. Application to edge states of honeycomb structures. In this section, we apply our method to Hamiltonians which model electronic states at edges of twodimensional honeycomb structures with defects. The structure of this section is as follows.
In Section 5.1 we define the tight-binding bulk Hamiltonian of a general honeycomb structure with real, nearest-neighbor hopping. In Section 5.2 we consider Hamiltonians describing states at a zig-zag edge of a graphene-like structure when the edge has defects. We choose to focus on the zig-zag edge because this edge (as opposed to, for example, an armchair edge) is known to support edge states in the absence of disorder. In Section 5.3 we consider an edge Hamiltonian describing states localized along the interface of two distinct such structures. We assume in this case again that the interface occurs along a zig-zag edge. In both cases we employ a supercell-type approximation to reduce a two-dimensional problem to a onedimensional problem. In Section 5.4 we show that our method can also be applied to a continuum Schrödinger equation model of the boundary of a honeycomb structure under finite difference discretization.
5.1. Tight-binding honeycomb structure with nearest-neighbor hopping bulk Hamiltonian. We consider a single electron without spin hopping on a two-dimensional honeycomb lattice. Each fundamental cell, labeled by integers m and n, of such a lattice hosts two atoms, which we label A and B. We introduce the notation to denote the restriction of the electron wavefunction ψ ∈ l 2 (Z 2 ; C 2 ) to the (m, n)th fundamental cell of the lattice. The bulk Hamiltonian of a general honeycomb structure with nearest-neighbor hopping is: [H bulk ψ] m,n = A * ψ m−1,n + B * ψ m,n−1 + V ψ m,n + Aψ m+1,n + Bψ m,n+1 m, n ∈ Z × Z (5.2) where the onsite potentials V A , V B are assumed real, and the hopping amplitudes t 0 , t 1 , t 2 are assumed real and non-zero. Bloch's theorem in this context implies that bounded eigenfunctions of H bulk satisfy: for k j ∈ [−π, π], j ∈ {1, 2}. It is then a standard calculation to find the Bloch bands E ± (k 1 , k 2 ) of H bulk . Depending on the values of the t j and V σ the Hamiltonian (5.2) may or may not have a spectral gap. When t 0 = t 1 = t 2 and V σ = 0, σ ∈ {A, B} the Hamiltonian (5.2) models graphene whose Bloch bands are non-degenerate other than at K := 2π 3 , − 2π 3 and −K, known as the "Dirac points", where the bands touch at E = 0 [22].

5.2.
Tight-binding honeycomb structure edge Hamiltonians with defects. We now model the edge of a honeycomb structure possibly subject to defects as follows. We allow for inter-atom hopping amplitudes t j (m, n), j ∈ {0, 1, 2} and for onsite potentials V σ (m, n), σ ∈ {A, B} which depend on the cell indices m, n, making the following assumptions.
Assumption 5.1 (Periodicity of bulk medium). There exists a non-negative integer M such that for all m ≥ M , V σ (m, n) = V σ for some real V σ and t j (m, n) = t j for some real t j , for all n ∈ Z, σ ∈ {A, B}, and j ∈ {0, 1, 2}.
Remark 5.1. Our method may be used to study the effect of more general disorder not satisfying Assumption 5.2 for any N . This leads to an approximate numerical method known as the supercell method (see for example [1]). Under Assumption 5.2 our method is exact.
In summary, we consider the Hamiltonian where the matrices A(m, n), B(m, n), V (m, n) are defined by  and is hence amenable to our method. An illustration of this notation is given in Figure 5.1, while results of computations using our method are shown in Figure 5.2.

5.3.
Honeycomb structure tight-binding domain wall edge Hamiltonian. We now consider a second example: a "domain wall" interface between distinct honeycomb structures with defects. For simplicity we again assume that the interface occurs along a "zig-zag" edge. We make analogous assumptions to Assumptions 5.1-5.2. Specifically, we allow for inter-atom hopping amplitudes t j (m, n), j ∈ {0, 1, 2} and onsite potentials V σ (m, n), σ ∈ {A, B} which are non-trivial functions of the cell indices m, n for |m| sufficiently small. The restrictions we make are summed up in the following assumptions.   , 2} for all m, n ∈ N × Z is plotted. We find that this bound state has eigenvalue precisely 0. The absolute value |ψ| and the argument arg ψ of the wave function at each lattice point is represented by the radius and color, respectively, of the circle at that point. In (b), an edge state of the zig-zag edge Hamiltonian (5.7) with multiple atomic vacancies, modelled by a large onsite potential at each vacancy, is plotted. We observe that the edge state survives the perturbation: the Hamiltonian retains a bound state with eigenvalue 0. Note that although the spectrum of the full two-dimensional bulk Hamiltonian is not gapped, for fixed values of k / ∈ 2π 3 , − 2π 3 the essential spectrum of H ZZ edge (k ) has a gap which is symmetric about 0.
data of the Hamiltonian (5.9) for the case (5.11) This case, studied for example in [5], is of particular interest to us because hard truncation yields inaccurate results in this case (see Lee-Thorp [17], in particular Figure 26.7). In Figure 5.3 we demonstrate this by comparing the edge dispersion curves computed with hard truncation and the Green's function method. In Figure  5.4 we study a bound state of (5.9) in the presence of defects.

Honeycomb structure continuum edge Hamiltonian.
Our method can also be used to compute spectral data of continuum Schrödinger edge models discretized by a finite difference scheme. Following [11], we consider the eigenvalue problem H cont edge ψ = Eψ for the Schrödinger operator   acting on L 2 (R 2 ), where the real potential function V (x) models a "half-plane" honeycomb structure. More precisely, introduce the lattice vectors Then we assume where Λ + is the "half-plane" honeycomb lattice Λ + = {mv 1 + nv 2 : m ∈ Z, n ∈ N}. For simplicity we take the potential at each lattice vertex V 0 (x) to be a Gaussian well (see Figure 5.5). Such a Hamiltonian models electronic states at a "zig-zag" edge of a graphene-like structure. We model atomic vacancy defects at the edge by omitting terms in the sum (5.14). Under assumptions analogous to Assumptions 5.3 and 5.4 we obtain a Hamiltonian which is periodic under translation by N v 1 (parallel to the edge) for sufficiently large N > 1 and eventually periodic (see Remark 5.3) with respect to translation by v 2 for x sufficiently far from the edge, both into the bulk of the material and into free space in the opposite direction. Imposing Bloch periodicity in the edge direction Φ(x + N v 1 ; k ) = e iN k Φ(x; k ), we derive a family of eigenvalue problems on L 2 ([0, N ] × R) parameterized by k . Discretizing any one of these eigenvalue problems by a five-point stencil we obtain an infinite eventually periodic Hamiltonian in the sense of Definition 4.2 whose bound states can be computed using the Green's function method. Our results are shown in Figures 5.5 and 5.6. Remark 5.3. In fact the Hamiltonian cannot be exactly eventually periodic even in the absence of defects because of the exponential tail of the potential at each vertex. Since the error accumulated by assuming that it is in fact eventually periodic is exponentially small (indeed can be made arbitrarily small by truncating further away from any defects) we ignore this issue.

Conclusion.
In this work we have introduced a numerical method for computing spectral data of a class of discrete Hamiltonians which model edge states: electronic states localized at edges of one and two-dimensional materials. Such states have generated considerable interest because of their potential as robust wave-guides. Our method is exact up to rounding error for discrete models and exact up to discretization error for continuum models and is hence a significant improvement on the naïve "hard truncation" method which may yield spurious states localized at the truncation rather than the edge of the structure. We have verified the effectiveness of our method through studies of the robustness of edge states to defects in various one and two-dimensional models.
Much remains to be explored related with edge states and their stability. For instance, it should be possible using the ideas presented in this paper to realize a scheme for studying the propagation in time of wave-packets propagating along edges which doesn't require "hard truncation" into the bulk. Another interesting direction would be to adapt our method to compute resonances resulting from defects (localized in both directions as well as along one-dimensional "edges") in periodic media. Locating such resonances (specifically, computing their imaginary part) would give information about the life-times of approximate bound states created by such defects. Finally, our methods can already be directly applied to the computation of "surface" states of three dimensional materials in the presence of defects under analogous assumptions to Assumptions 5.1 and 5.2.
A. Proof of existence of boundary conditions for arbitrary finite range eventually periodic Hamiltonians. In this section we discuss a more general framework for proving the existence of appropriate boundary conditions (i.e. those which ensure that the Green's function we compute maps l 2 (N; C N ) to itself) for computing the solution of (4.9).
The only assumptions we need on the Hamiltonian H ∈ C N×N are made explicit next. Note that we do not even require that H is Hermitian; our method works for non-Hermitian operators as well. However, if H is not Hermitian, then Theorem 3.1 is no longer guaranteed. Unfortunately, the conclusion of Theorem 3.1 does not hold for general operators, so one would have to be more careful in justifying that step in the non-Hermitian case. Let us define some notation. Under Assumptions A.1 and A.2, the matrix z − H can be split up as shown next.
We call (z − H) per the periodic-part of z − H; it consists of the matrix P (representing one period) periodically repeating off to infinity. We call (z − H) fin the finite-part; it consists of the defect-part (z − H) defect (which is where defects are allowed) and the last portion of a P matrix, which we have denoted P. We assume that P has no "extra" columns of zeros on either the left or right end: It has a nonzero entry in its leftmost and rightmost column. We now introduce some more notation by defining the sizes of these matrices: We finally define one more piece of notation: The first P overlaps with b − columns of (z − H) defect . That is, the leftmost nonzero entry of (z − H) per occurs in column N d − b − + 1. Now that we have defined notation, let's look at the problem we actually want to solve. We wish to compute the upper-left N d × N d block of G(z); that is, we wish to compute the first N d entries of the solution to (A.2) (z − H)g j = e j for j = 1, 2, ..., N d . First, we remark that this solution is unique in 2 (N).
Lemma A.1. If z is not in the spectrum of H, then the linear system (z − H)g j = e j has a unique solution in 2 (N).
Proof. Under Assumptions A.1 and A.2, it is straightforward to show that H is a bounded operator on 2 (N). The conclusion then follows.
Our goal is to compute the finite-part of g j . In order to do this, there are three conditions g j must to satisfy: 1. g j must satisfy all the finite-part equations: (z − H) fin g j = e j . 2. g j must satisfy all the periodic-part equations: (z − H) per g j = 0. 3. g j must be in 2 . By Lemma A.1, g j is the only vector which satisfies all three of these conditions. We can reduce these three conditions down to just two by combining Conditions 2 and 3: We see that g j satisfies both conditions if and only if g j is in the nullspace of (z − H) per , considered as an operator on 2 (N). This leaves us with just two requirements: 1'. g j must satisfy all the finite-part equations: (z − H) fin g j = e j . 2'. g j must be in the nullspace of (z − H) per , considered as an operator on 2 (N). Note that these requirements are again both necessary and sufficient: A vector satisfies these two conditions if and only if it is g j .
The rest of this section is organized as follows. First, we characterize the nullspace of (z − H) per , considered as an operator on 2 (N). Once we have this nullspace, we show how to construct the boundary conditions. Lastly, we prove that this procedure for incorporating the boundary conditions always results in a system whose unique solution is the finite-part of g j .
To compute the nullspace of (z − H) per , we proceed in a few steps. First, Lemma A.2 characterizes the eigenspaces of the operator that shifts a vector by one period, S p : [ψ 1 , ψ 2 , ...] T → [ψ p+1 , ψ p+2 , ...] T . Second, using the fact that (z − H) per and S p commute, we show in Lemma A.3 that the nullspace of (z − H) per is spanned by a finite set of (generalized) eigenvectors of S p , which fall into two particular types. (Note that the standard "commuting operators have a common set of eigenvectors" theorem does not apply here because neither (z − H) per nor S p is known a priori to be diagonalizable; neither operator is normal.) Using these particular forms for the basis vectors of Nul(z − H) per , we show how a polynomial eigenvalue problem can be used to actually compute such a basis.
Lemma A.2. (Characterization of the eigenspaces of S p , as an operator on 2 (N)) Let η ∈ (−1, 0) ∪ (0, 1) and ξ ∈ C p , ξ = 0. Then, is an eigenvector of S p corresponding to eigenvalue η. This characterizes the nonzero eigenspaces. The generalized eigenspace of the zero-eigenvalue is the space of vectors in C N with only a finite number of nonzero entries.
Proof. The nonzero claim is straightforward. For the zero-eigenspace, note that ψ is a generalized eigenvector if and only if S k p ψ = 0 for some k. The ψ which satisfy this are precisely the ones with only a finite number of nonzero entries. Proof. See Appendix B.1.
Next, we show how to find such a basis for Nul(z − H) per by solving a generalized eigenvalue problem. First, let us define two submatrices R and S of (z − H) per as follows. R and S are square matrices of identical size and have side lengths that are equal to an integer multiple of p. R and S are next to each other (with S immediately to the right of R), and both are aligned with the top of (z − H) per . Finally, we must choose R and S large enough so that they contain all the nonzero entries of (z − H) per in the rows they cover. Additionally, R must be chosen so that it contains all the all the nonzero entries of Type 1b vectors. With these definitions, the basis of Lemma A.3 can be found by solving the generalized eigenvalue problem R + λS. In particular, the eigenvectors of R + λS can be extended (in the sense of the following definition) to basis vectors of Nul(z − H) per .
Remark A.1. It is not always known a priori how large R should be chosen. But if a choice fails (i.e., you do not get enough eigenvectors of R + λS), then you simply try again with a larger choice.
Definition A.4. Let x be an eigenvector of R + λS; that is, for some λ 0 , (R + λ 0 S)x = 0. Then, we define the extension of x to 2 (N) as follows. If λ 0 = 0, then the extension is ξ = [ 0 N b x 0 ∞ ] T . If λ 0 = 0, then we define the N b + 1 through N b + N R entries to be x; i.e., in MATLAB notation, ξ N b +1:N b +N R = x. The next N R entries are defined to be λ 0 x, then λ 2 0 x, and so on. The first N b entries of ξ are defined by continuing this pattern backwards. Conversely, we define the restriction v restr of a vector v ∈ 2 (N) to be v restr To see that it is enough to solve the generalized eigenproblem R + λS, we first note that, since we chose the sides of R to be an integer multiple of p and chose them large enough, restrictions of the Type 1b and Type 2 basis vectors form solutions of (R + λS)ξ = 0, corresponding to |λ| < 1. Second, we show that these are the only such solutions. We note that there could also be generalized eigenvectors of R + λS. However, we do not want these --they correspond to generalized zero-eigenvectors of (z − H) per .
Remark A.2. Note that the restriction that we choose R large enough to encompass all the nonzero entries of the finitely-many nonzero entry basis vectors is important; it guarantees that the Type 1b vectors will be zero-eigenvectors of R + λS. An example of why we might have to choose R to be large is illustrated with the following example. Suppose This vector has nonzero values past the end of P. We note that this example can easily be generalized to make the last nonzero value go arbitrarily far past the end of P by choosing P to be for any positive integer k, where I k is the k × k identity matrix.
Using Lemma A.5, we can now detail how to construct the boundary conditions. The idea is simple: We restrict the solution to lie in the nullspace of (z − H) per . We do this by requiring that the solution is orthogonal to the orthogonal complement of the nullspace.
First, note that since {e j } N b j=1 ∈ Nul(z−H) per , there are no constraints on the first N b entries of the solution. Therefore, we only need to enforce the condition that the truncation of the solution is orthogonal to the eigenvectors of R + λS. This is easily done as follows. Let V ∈ C N R ×(N R −N d +N b ) be a basis for the orthogonal complement of the eigenvectors of R + λS. Then we can solve for g j using the system, Finally, we prove that g j is the unique solution of this system.
Lemma A.6. Using our method for constructing the boundary conditions, (A.7) has a unique solution. That unique solution is the first N R + N b entries of g j .
Proof. First, note that existence of a solution is clear: g j is a solution. Next, let v be a solution of (A.7). Then, v can be extended to 2 (N) and would satisfy Conditions 1' and 2'. Therefore, by Lemma A.1 and the discussion following it, v must equal g j . Lemma B.1. Let A be an invertible operator on 2 (N). In the standard basis for 2 (N), consider A as an infinite dimensional matrix. Let A n be the operator A with n rows removed. Then, dim(Nul(A n )) = n.
Proof. Let R be the n × ∞ matrix of the rows removed from A. We first show that R H and A H n span disjoint subspaces, Note that the set {P v i } n i=1 is linearly independent: If it was not then there would be a nonzero vector v in Ran(R H ) such that P v = 0; but this would imply that v ∈ Ran(A H n ), which is not possible. Finally, using a well-known property of linear operators, we obtain the desired conclusion, Then, the matrix V H S p V ∈ C N d ×N d has a Jordan normal form, Since S p is invariant on the range of V (i.e., on the nullspace of (z − H) per ), we can multiply on the left by V (and the right by M ) to obtain This is essentially a partial Jordan normal form of S p . Since the columns of V form a basis for Nul(z − H) per and M is nonsingular, (B.8) shows that there is a basis for Nul(z − H) per that is made up of (generalized) eigenvectors of S p . In particular, such a basis is given by the columns of V M . Lemma A.2 then implies that the nullspace of (z − H) per is spanned by vectors that are either of form (A.3) or that have a finite number of nonzero entries. The vectors with a finite number of nonzero entries can be split into Type 1a and Type 1b by taking proper linear combinations.

B.2. Matrix pencils.
In order to prove Lemma A.5, we will use results from the theory of matrix pencils. We review the definition and relevant properties of matrix pencils here. For more information, see [8].
Definition B.2. A matrix pencil is a set of matrices of the form {A − λB : λ ∈ C}. A matrix pencil is said to be regular if det(A − λB) is not identically zero. If det(A − λB) is identically zero, then the matrix pencil is said to be singular. Note that if a matrix pencil is singular, then the corresponding generalized eigenvalue problem has infinitely many generalized eigenvalues (any λ ∈ C is an eigenvalue). A concrete example of this phenomenon is given by the generalized eigenvalue problem, This phenomenon is very different from that encountered in the standard eigenvalue problem, (A − λI)x = 0, where there are always precisely n eigenvalues if A ∈ C n×n ; additionally, the eigenspaces corresponding to distinct eigenvalues are not all independent. This suggests that the theory of generalized eigenvalue problems is more complicated than that of the standard eigenvalue problem. However, the generalized eigenproblems corresponding to regular matrix pencils maintain many nice properties that hold for standard eigenvalue problems. In order to show this, we first introduce a concept integral to our proof of Lemma A.5 --the Kronecker canonical form. where P and Q are invertible square matrices, and A − λ B is block diagonal. Each block of A − λ B is in one of four forms, where J j (α) and N j are j × j matrices, L j is j × (j + 1) and L T j is (j + 1) × j. The value of j may change from block to block.
The J j (α) blocks correspond to finite generalized eigenvalue α. The N j blocks correspond to infinite generalized eigenvalues (i.e., correspond to generalized eigenvalue zero of the problem λA − B). L j and L T j are called singular blocks. Note that the Kronecker canonical form of a matrix pencil is very similar to the Jordan canonical form of a matrix. However, in Jordan canonical form, there are only J j (α) blocks. The similarity between the two canonical forms becomes even closer when one restricts to regular matrix pencils.
Lemma B.5. If A − λB is a regular matrix pencil, then its Kronecker canonical form does not contain any singular blocks (L j or L T j ). From this lemma, it is easy to see that generalized eigenproblems corresponding to regular matrix pencils satisfy two properties familiar from standard eigenproblem theory: (1) If A, B ∈ C n×n , then there are precisely n eigenvalues (including the infinite eigenvalues and counting algebraic multiplicity); and (2) The eigenspaces corresponding to distinct eigenvalues are all independent.
B.3. Proof of Lemma A.5. Before proving Lemma A.5, we first show that the matrix pencil considered in that lemma is regular.
Lemma B.6. The matrix pencil R + λS in Lemma A.5 is regular.
Proof. The idea is as follows: If R + λS is singular, then its Kronecker canonical form contains an L j (or L T j ) block; this fact can then be used to show that columns (or rows) of (z − H) must be linearly dependent, which contradicts the assumption that z / ∈ σ(H). The details now follow. First, we will only consider the case that the Kronecker canonical form (B.11) R + λS = P ( R + λ S)Q −1 has an L j block; the L T j case is similar. Assume without loss of generality that the top-left block of R + λ S is an L j block. This implies that Se j λ j+1 + ( Re j + Se j−1 )λ j + ... + ( Re 2 + Se 1 )λ + Re 1 = 0.
Since this must be true for all λ, the coefficient in front of each λ k must be zero. This then gives the set of equations,  Since Q has full rank, the vector is nonzero. The columns of the matrix on the left are identical to the columns of (z − H) as the matrix goes off to infinity (with the non-relevant zero-rows above and below removed). This equation then implies that the columns of (z − H) are linearly dependent.
Using this lemma and the nice properties that follow from it (as mentioned below Lemma B.5), we can now prove Lemma A.5.
Proof of Lemma A.5. As noted before the statement of the lemma, it is simple to check that restrictions of the Type 1b and Type 2 basis vectors satisfy (R + λS)ξ = 0 with |λ| < 1. Since R + λS is a regular pencil (Lemma B.6), the Kronecker canonical form implies that eigenspaces corresponding to different eigenvalues are independent. Now suppose that the set of restrictions of the Type 1b and Type 2 basis vectors of Nul(z − H) per is linearly dependent. This would imply that there is a dependent set would corresponds to the same eigenvalue. However, this then implies that the extensions of these vectors are linearly dependent --contradicting the assumption that they were not. Therefore, the restrictions of Type 1b and Type 2 vectors form a linearly independent set of N d − N b eigenvectors of R + λS.
Let ξ be an eigenvector of R + λS with |λ| < 1. Then, it is straightforward to show that the extension of ξ is in Nul(z − H) per . Therefore, the extension of ξ must be a linear combination of the basis vectors for Nul(z − H) per . Hence, ξ is a linear combination of the restrictions of these basis vectors. Therefore, there are precisely N d − N b independent eigenvectors of R + λS with |λ| < 1.
Finally, it is straightforward to show that the extensions of any N d − N b independent eigenvectors of R + λS (combined with {e j } N b j=1 ) gives a basis for Nul(z − H) per .  when the edge does not have defects. Note that we choose co-ordinates such that the edge is at y = 0. To generate figures (c) and (d), we have removed two atoms as can be seen in the plot of the external potential (c). An edge state with k = 3.8π 3 is also plotted. The edge state plotted in (b) has eigenvalue -1.382, while the edge state plotted in (d) has eigenvalue -1.411. Essential spectrum (computed via two-dimensional Bloch theory) is shown in blue while edge states computed via the Green's function method are marked with red ×s. Note that the curve of edge state eigenvalues is not flat, in contrast to the prediction of the analogous tight-binding model.