Theory of edge states based on the hermiticity of tight-binding Hamiltonian operators

We develop a theory of edge states based on the hermiticity of Hamiltonian operators for tight-binding models defined on lattices with boundaries. We describe Hamiltonians using shift operators which serve as differential operators in continuum theories. It turns out that such Hamiltonian operators are not necessarily hermitian on lattices with boundaries, which is due to the boundary terms associated with the summation by parts.The hermiticity of Hamiltonian operators leads to natural boundary conditions, and for models with nearest-neighbor (NN) hoppings only, there are reference states that satisfy the hermiticity and boundary conditions simultaneously. Based on such reference states, we develop a Bloch-like theory for edge states of NN models on a half plane. This enables us to extract Hamiltonians describing edge states at one end, which are separated from the bulk contributions. It follows that we can describe edge states at the left and right ends separately by distinct Hamiltonians for systems of cylindrical geometry. We show various examples of such edge state Hamiltonians (ESHs), including Hofstadter model, graphene model, and higher-order topological insulators (HOTIs), etc.


I. INTRODUCTION
In the studies of topological properties of matter [1][2][3][4][5], the bulk-edge correspondence [6,7] plays a central role. Although topological invariants are well-defined for the Bloch wave functions, such bulk topological invariants are not directly related with the physical observables, except for the integer quantum Hall effect [1,8]. Therefore, one usually judges topological phases of mater from observation of gapless edge states embedded in the bulk insulating states [4,5,9].
Usually, edge states mean localized states on d − 1 dimensional boundaries for d dimensional bulk systems. However, recent discovery of higher-order topological insulators (HOTIs) [10][11][12][13][14][15][16] have drawn our attention to unconventional d − 2, · · · dimensional higher-order edge states such as corner states, hinge states, etc. Therefore, the concept of edge states, including higher-order edge states, would become more and more important when one investigates topological phases for various kinds of systems.
Bulk states are described by the Bloch states. Let us consider a tight-binding model defined on a lattice which is characterized by a unit cell with n species as well as d dimensional translation vectors. Then, each unit cell is assigned by set of integers j = (j 1 , j 2 , · · · ) as its position. After the d-dimensional Fourier transformation for j, the Hamiltonian becomes a n × n matrix with wave numbers k = (k 1 , k 2 , · · · ) which gives n Bloch bands describing the bulk system. Their wave functions would define bulk topological invariants.
On the other hand, when one introduces boundaries to this model, one usually calculates the edge states and bulk states together. Let us consider d − 1 dimensional boundaries perpendicular to 1-direction, (1, 0, · · · , 0), e.g., at j 1 = 1 and N . Fourier-transforming all j except for j 1 , one obtains a one-dimensional Hamiltonian with N unit cells along 1-direction specified by j 1 . The resultant Hamiltonian is a nN × nN matrix, which includes not only the edge states but also the bulk states. In order to obtain the well-defined edge states, one should choose N as N n, so that among nN states of the Hamiltonian, the number of edge states is only of order n, and most of them are bulk states. Although edge states and bulk states are coupled together for finite N system, it is very desirable to separate the edge states from the bulk states and to derive an effective theory of the edge states only.
Contrary to such tight-binding Hamiltonians, the continuum Dirac models allow us to calculate edge states analytically without considering bulk states [15][16][17][18]. In particular, as pointed out by Witten [17] as well as Hashimoto et. al. [15,16], boundary conditions have intimate relationship with the hermiticity of continuum Hamiltonians. Along this line, higher-order topological insulators based on continuum models have been developed, introducing the idea of "edge of edge" states [16]. In Ref. [18], we investigated the edge states of the continuum Dirac model which is derived as a continuum limit of the tight-binding model defined on the square lattice. Then, it was pointed out that the boundary conditions of the lattice model also serve as those of the continuum model. This suggests an intimate relationship between the continuum models and tight-binding models.
The role of differential operators in continuum models are played by difference (or shift) operators in lattice models. It is then expected that Hamiltonians expressed by such difference (or shift) operators on lattices are not necessarily hermitian, as in the case of the continuum Dirac Hamiltonians, if systems have boundaries. Then, even for lattice models, the hermiticity of Hamiltonians would enable us to calculate edge states without considering the bulk states. This implies that the hermiticity gives us edge state Hamiltonians (ESHs) describing edge states only.
In this paper, we derive first-quantized Hamiltonian operators on lattices including shift operators. When we solve the Schrödinger equations using the Hamiltonian operators thus defined, it turns out that Hamiltonians are not necessarily hermitian when systems have boundaries. We show that like the continuum Dirac theories, the hermiticity of Hamiltonians would basically determine edge states. To be concrete, when we can choose boundary wave functions that satisfy the hermiticiy of Hamiltonians, we can develop a Bloch-like theory for edge states based on the boundary wave functions, which enables us to extract ESHs separated from bulk states. We show that this can be carried out for models with nearestneighbor (NN) hoppings only. The ESHs thus obtained are defined for a single end on a half line or a half plane. Therefore, for systems of cylindrical geometry with two ends, the edge states at the left end and right ends can be described separately by two distinct ESHs. Even if models include next-nearest-neighbor (NNN) hoppings, ESHs can be formally defined. However, in these models, simple Bloch-like wave functions ensuring the hermiticity of Hamiltonians are not enough: Their linear combinations are needed to ensure the boundary conditions, implying that for an edge state to exist, the ESH must allow degenerate two or more Bloch-like wave functions.
There have already been several papers dealing with extended Bloch theories for edge states defined on semiinfinite systems. Kunst et. al. [19][20][21][22] proposed systematic approach toward exact solutions of the edge states including higher-order edge states defined on semiinfinite lattice spaces. Pletyukhov et. al. [23,24] developed topological arguments based on exact edge states of a variant of the Hofstadter model solved also on semi-infinite lattice spaces. The method proposed in the present paper is a simple and unified theory based on the hermiticity of first-quantized Hamiltonian operators. In our approach, all information is included in the first-quantized Hamiltonian operators in any dimensions, and moreover, their relationship with continuum theories such as Dirac fermions is clear. This paper is organized as follows. In Sec. II, we formulate our method using a generalized version of the Su-Schrieffer-Heeger (SSH) model [21,25] as an example. Almost all concepts can be developed by this simple one-dimensional (1D) model, including the hermiticity of Hamiltonian, its relationship with the boundary condition, a Bloch-like edge state [21], and classes of models categorized according to how the present formulation can be applied, etc.
We then proceed to two-dimensional (2D) models [21,22].. We first discuss models with NN hoppings only. In Sec. III A, we apply our method to the Hofstadter model [26] and give explicit ESHs for the left as well as the right ends [23,24]. Sec. IV is devoted to edge states of the graphene in the absence/presence of a magnetic field [27,28]. As is well-known, the graphene allows various edge states depending on the shapes of edges. We present the ESHs describing zigzag and bearded edges. For an armchair edge in the presence of a magnetic field, it is too difficult to obtain the edge state by use of our method, since this is a NNN model. Instead, we present our attempt to obtain these edge states, by introducing some defects into the models.
Next, we apply out method to NNN models. In Sec. V, we consider the Wilson-Dirac model, as an example of simpler models including NNN hoppings. In Sec. VI, we consider the Haldane model [29] as a most difficult case for our formulation. From a practical point of view, our method is not useful for this class of models. Nevertheless, we would like to point out that the hermiticity of Hamiltonians basically determines edge states even in this class.
Finally, we switch to 2D HOTI models [20]. In Sec. VII, we study corner states as well as edge states for typical HOTI models defined on the square lattice and on the breathing kagome lattice [11,12,30,31]. Although edge states are governed by the same SSH Hamiltonian, different hermiticity conditions lead to different edge states. In Sec. VIII, we give summary and discussion.

II. HAMILTONIAN OPERATORS FOR LATTICE MODELS
For continuum models with boundaries, the hermiticity of their Hamiltonians is nontrivial, since momentum operators yield boundary terms associated with the integration by parts [16,17]. In tight-binding models defined on lattices, the differential operators are replaced by difference (or shift) operators, implying that the summation by parts plays a similar role. Therefore, even for lattice models, the hermiticity of Hamiltonians is expected also nontrivial when they are in the lattice space representation.
In II A, we introduce the shift operators and related summation by parts, and in Sec. II B, we demonstrate, using the generalized SSH model, that with a boundary the Hamiltonian becomes hermitian only if a suitable condition is imposed.

A. Summation by parts
In order to study conventional tight-binding models in condensed matter physics, it may be convenient to introduce the shift operators, instead of the difference operators. Let f j be a function of discrete sequence of integers denoted by j. Then, the forward and backward shift operators are defined by One can show the following summation by parts Thus, we have δ † = δ * for a bulk system, where bulk means a system defined over j = −∞ to +∞. If the system has a boundary, summation by parts yields a boundary term, These properties can be regarded as a lattice version of the integration by parts.
B. Typical example: 1D SSH model

First-quantized Hamiltonian operator
Let c j ≡ (c 1j , c 2j ) be the annihilation operators of the fermions on the j-th unit cell, as shown in Fig. 1. Then, the second-quantized Hamiltonian is given by where ← − δ acts on the left. The operator ← → H defined on the lattice will be referred to as the (first-quantized) Hamiltonian (operator).
For the bulk system, the Hamiltonian operator (4) can be rewritten as whereĤ defined bŷ operates on the right only, and hence, it is appropriate for the Schrödinger eigenvalue equation. The symbol hat on H emphasizes that H is not a simple matrix but a matrix-valued operator.Ĥ will be also referred to as the (first-quantized) Hamiltonian (operator). On the other hand, for the system with a boundary, there appears a boundary term. Let us introduce a boundary between j = 0 and j = 1 as in Fig. 1, and consider the system on the semi-infinite line, j 1 ≥ 1. Then, where the operator c 0 out of the boundary has been taken into account, and the subtracted term in right-hand-side is due to the summation by parts in Eq. (3). Thus, for a system with a boundary, ← → H andĤ differ in the boundary term. Generically, the boundary term can be written by nonhermitian matrix K associated with the operator In what follows, we will discuss various aspect of the Schrödinger equation forĤ, whereĤ is defined by Eq. (6), δ and δ * operate on j of ψ jn , and n denotes the energy quantum number.

Hermiticity ofĤ
First of all, let us discuss the hermiticity of this operator. For any wave functions ψ j and φ j including those at j = 0, ψ 0 and φ 0 , we have Therefore, only by imposing the condition will the Hamiltonian become hermitian. Now, assume that we have a set of eigenfunctions of H, ψ jn , with the hermiticity (10) imposed.
for all m, n. Then, they form a complete orthonormal set of functions ψ j on the semi-infinite line j ≥ 1. Hence, we have Define new operators d n and d † n such that and substitute these into the second quantized Hamiltonian (7), Note that the boundary term in the Hamiltonian vanishes due to the hermiticity Eq. (11) imposed on the wave functions, and we have a desired second-quantized Hamiltonian.

Boundary conditions
Let us now discuss the eigenvalue equation (8). For the bulk system, the Fourier transformation can be readily carried out just by replacing δ → e ik and δ * → e −ik . Then,Ĥ becomes a simple 2×2 matrix whose eigenstates are nothing but the Bloch states.
Let us next consider the system with a boundary. Assume that the system is defined on the semi-infinite line, j ≥ 1. Then, we need to specify the boundary condition at j = 1. To this end, let us write down the Schrödinger equation (8) explicitly, where the energy quantum number n has been suppressed. Here, the boundary condition associated with j = 1 above reads where we have supplementarily included ψ 0 which is out of the boundary. If one does not take ψ 0 into account or set ψ 0 = 0 from the beginning, the boundary condition at j = 1 is given by with a suitable initial value ψ 1 .
In this paper, we adopt Eq. (16) as a boundary condition. This may be much simpler than (17), if we can find ψ 0 . Based on such a ψ 0 , we can develop Bloch-like techniques, as we show below. Even if such a ψ 0 cannot exist, we can take the following route: Without considering the boundary condition (16), we solve the Schrödinger equation (8) only imposing the hermiticity (10). Then, using the wave functions thus obtained, we construct suitable eigenstates satisfying the condition (16) by taking their linear combination.
In this sense, we regard the hermiticity of Hamiltonians (11) as a guiding principle to choose ψ 0 . In what follows, we develop a Bloch-like theory for the edge states based on ψ 0 satisfying the hermiticity (10). In the case of NN hopping models, such a ψ 0 naturally satisfies the boundary condition (16), whereas for other models including NNN hoppings, we have to take linear combinations of degenerate Bloch-like wave functions to construct the wave functions satisfying the boundary condition (16), as will be argued below.

Edge states
For the system defined on the semi-infinite line j 1 ≥ 1, let us solve the eigenvalue equation (8) assuming wave functions decaying exponentially. To this end, we introduce a Bloch-like wave function with a complex wave number K, where we have utilized the supplemental wave function ψ 0 as a reference state, and by definition, κ > 0 is required. Since this model allows at the most one edge state, the quantum number n has been suppressed. Then, substituting Eq. (18) into Eq. (8), the eigenvalue equation becomes Since ψ 1 = ψ 0 e iK , the hermiticity condition (11) can be written solely by ψ 0 , This equation requires that The above vectors can also be interpreted as follows: Since {K, σ 3 } = 0, the eigenstates of σ 3 , ψ ↑,↓ , satisfy the hermiticity condition, as discussed in Refs. [16,17].
• NN model (t = 0) Let us first consider the case with nearest neighbor hopping only, setting t = 0. Here, consider the boundary in Fig. 1. The system defined on the infinite line can be separated at the dotted-line if the second component of the wave function at j = 0 is set 0. Thus, from the point of view of such a lattice termination, the appropriate wave function at j = 0 is the former ψ ↑ in Eq. (21). Namely, among the functions satisfying the hermiticity (20), we have to choose the reference function ψ 0 = ψ ↑ . What is important here is that holds. This guarantees that the boundary condition (16) is automatically satisfied. For NN models, the hermiticity, the lattice termination, and the boundary condition are consistently satisfied by a suitable reference state ψ 0 , as will be seen in various examples. Substituting ψ 0 = ψ ↑ into the eigenvalue equation (19), we have The latter equation gives Thus, an edge state exists when 0 < e −κ < 1. Explicitly, it is given by valid for |γ/λ| < 1.
• NNN model (t = 0). In this case, Kψ ↑,↓ = 0, implying that the hermiticity does not ensure the boundary condition. Even in this case, if Eq. (19) has degenerate two solutions K 1 = K 2 , we have a wave function of the form which becomesψ 0 = 0, and therefore, the boundary condition (16) holds. Let us choose ψ 0 = ψ ↑ . Then, similarly to Eq. (23), we have If we choose ψ 0 = ψ ↓ , we have When these equations have two solutions |e iK | = e −κ < 1, the model has an edge state. We show in Fig. 2 the phase diagram of the generalized SSH model.

Bulk states
For the bulk state, we know that the conditions associated with the boundaries are not very important, although they should be satisfied. Therefore, we firstly solve the eigenvalue equation for the Bloch states with a real wave vector k, where at the moment, we do not consider the hermiticity and the boundary condition. Then, from the Schrödinger equation, we have Note that ε ± (k) = ε ± (−k). Therefore, for each state with ε ± (k) we have left-going and right-going waves, and their superposition allows the wave function satisfying the boundary condition (16): where α n ψ k,n + β n ψ −k,n ∝ ψ ↑ or ψ ↓ , depending on the phases in Fig. 2. This is the bulk state satisfying the hermiticity. In what follows, we will not consider the bulk states.

C. Classes of models
In the following sections, we will apply the method developed above to various 2D models. For this purpose, it is useful to classify models into following three types.
A: Reference state ψ 0 which satisfies Kψ 0 = 0 can be chosen. Models with a boundary whose Hamiltonian can be represented only by NN hoppings belong to this class. Here, by NN hopping, we mean that if one of bonds connected by finite hoppings between unit cells is cut, the system is divided into two peaces. The generalized SSH model with t = 0 is one of examples, and various other models such as the Hofstadter model in Sec. III A, the graphene with the zigzag or bearded edges in the absence/presence of a magnetic field in Sec. IV, HOTIs in Sec. VII belong to A. In this class, it is very easy to extract Bloch-like ESHs free from the bulk contributions.
B: Models including NNN hoppings, but there exist a matrix A with A 2 = 1 which anti-commutes with hermiticity matrix K, {A, K} = 0. The generalized SSH model with t = 0 is one of typical examples, and the Wilson-Dirac model in Sec. V is another example in this paper. In this class, we choose the reference state ψ 0 as either eigenstates of A, Aψ 0 = ±ψ 0 , which guarantees the hermiticity of the Hamlitonians. Then, we can develop a Bloch-like theory based on ψ 0 . However, we need two independent Bloch-like states to yield edge states satisfying the boundary condition.
C: If one cannot find any matrix A that anti-commute with K, one has to solve the Bloch-like eigenvalue equation based on the reference state ψ 0 which satisfies the hermiticity (10). Like the class B, two independent solutions are needed. As an example of this class, we will discuss the Haldane model in Sec. VI. From a practical point of view, our method is not useful for these models. Nevertheless, the idea that the edge states are determined by the hermeticity seems to be important.
Without loss of generality, we consider the system defined on the half plane j 1 ≥ 1, and discuss the edge states along the boundary j 1 = 1. Then, we can make the Fourier transformation in the 2-direction, The last equality is the abbreviation for notational simplicity: The operators δ ≡ δ 1 and δ * ≡ δ * 1 act on j ≡ j 1 , and k 2 dependence is suppressed. Note that H k2 is nothing but a 1D Hamiltonian, so that we can obtain the edge states using the technique in Sec. II B.
In the following sections, we will show that the method in Sec. II B can apply various 2D models. In particular, for NN models (class A in Sec. II C) ESHs can be explicitly obtained, which is exemplified in detail by the Hofstadter model in the following Sec. III A.
3: The Hofstadter model on the square lattice with a uniform magnetic flux φ = 2πp/q per plaquette. The arrow with nφ denotes the link with the phase e inφ . Large boxes denote a magnetic unit cell.

A. Hofstadter model
As a typical example of 2D models, let us consider the Hofstadter model [7,26]. Despite its simplicity, the model has revealed various aspects and properties of the topological phase, including TKNN integers and associated Diophantine equation [1,8], the bulk-edge correspondence [6,7], the Streda formula [32,33], etc. In particular, using the transfer matrix method, Hatsugai showed that edge states give winding numbers around a complex energy surface [6,7]. On the other hand, since this model is a NN model, the present formulation enables us to treat the edge states at the left and right ends separately and to extract the ESH at each end, as will be shown momentarily. Based on these, we can discuss the bulk-edge correspondence, which will be published elsewhere [34].
We consider the square lattice with φ/(2π) = p/q flux per plaquette, where p and q are coprime integers. For the Landau gauge illustrated in Fig. 3, the firstquantized Hamiltonian operatorĤ(δ 1 , δ 2 ) in Eq. (32) is given bŷ After the Fourier transformation with respect to the 2-direction as in Eq. (33), the Hamiltonian operator for the 1D chain toward the 1-direction is given bŷ where the operators δ ≡ δ 1 and δ * ≡ δ * 1 act on j ≡ j 1 , as was already mentioned.

Edge states at the left end
Paying attention to the backward operator δ * of the above Hamiltonian, we find the hermiticity matrix K, The hermiticity condition Eqs. (11) or (20), i.e., ψ † m0 Kψ 0n = 0, leads to where χ n is a vector with q − 1 components. From the point of view of the lattice termination, we see that the former wave function matches the boundary condition for the system defined on the half plane j 1 ≥ 1. To see this, it may be more convenient not to use the magnetic unit cell, but to label each site as j. Then, the Hamiltonian can be expressed by , where c j=1 , c j=2 , · · · corresponds to c 1,j1=1 , c 2,j1=1 , · · · and v j = 2 cos(k 2 − φj). Using such a new labelling of sites and assuming the wave function |ψ = j c † j ψ j |0 , the eigenvalue equation, H|ψ = ε|ψ is explicitly given by tψ j−1 + v j ψ j + tψ j+1 = εψ j (j = 1, 2, · · · ). Now, let us consider the system defined on the half plane j ≥ 1. Since the above eigenvalue equation at j = 1 becomes tψ 0 + v 1 ψ 1 + tψ 2 = εψ 1 , it is natural to require ψ 0 = 0 as a boundary condition, where ψ 0 in the new notation corresponds to ψ q,j1=0 . In passing, we would like to mention that the eigenvalue equation is basically given by recurrence relation of order 2. Therefore, it is natural to use 2 × 2 transfer matrix, as was used in Ref. [6,7].
On the other hand, in the present formulation, we utilize the magnetic unit cell representation, and develop a Bloch-like theory. Here, it should be noted here that wave functions in (37) satisfy the boundary condition (16), implying that a single Bloch-like state ψ jn = ψ 0n e iKnj can be an eigenstate of the Hamiltonian, where K n = k n + iκ n with κ n > 0. Substituting this into the eigenvalue equation, we have This equation tells that the upper diagonal (q − 1) × (q − 1) part of the above equation determines the eigenvalues ε n (k 2 ) and eigenstates of the edge states, and the last equation yields the constraint e iKn χ 1,n + χ q−1,n = 0 which determines whether the edge states obtained above are localized near j 1 = 1. Namely, the edge states are eigenstates of the following ESH, As mentioned, although H e looks independent of K n , the last equation of Eq. (39) yields the following constraint ensuring the localization of the edge states, κ n > 0, As we have suppressed the k 2 -dependence in the above discussion, κ n (k 2 ) depends on k 2 through χ n (k 2 ). Therefore, although the energies ε n (k 2 ) and wave functions χ n (k 2 ) in Eq. (41) form q − 1 continuous functions of k 2 in the Brillouin zone, they do not necessarily describe the edge states: It is only those satisfying the condition Eq. (42) that can be the edge states localized near the In passing, we have to mention that the normalization of the wave functions can be determined by where the sum over j converges due to the condition e −κn < 1. With thin in mind, we will not pay attention to the normalization of the wave functions below. In Fig. 4, the spectra of the edge states for the system defined on the half plane 1 ≤ j 1 ≤ +∞ are shown on the background of the total spectra for the system of cylindrical geometry, 1 ≤ j 1 ≤ L. The red and yellow curves are the eigenvalues of the ESH Eq. (40) with e −κn(k2) < 1 and e −κn(k2) > 1, respectively. One can see that red curves describe the spectra of the edge states localized at the left boundary, whereas there are no states corresponding to the yellow curves, since they cannot be physical states. Thus, for the theory of the edge states, not only the Hamiltonian (40) but also the condition (42) are needed to describe the edge states at one end. In passing, we add a comment that if we consider a cylindrical system without the rightmost q-th site in the rightmost unit cell j 1 = L, the yellow curves are edge states localized at the right end. However, we will consider the edge states at the right end in the way given in the next Sec. III A 2.
2. Edge states at the right end As demonstrated above, when we numerically calculate systems of cylindrical geometry, 1 ≤ j 1 ≤ L, it is inevitable to have both edge states at the left and right ends. The Hamiltonian (40) is only for the edge states at the left boundary, j 1 = 1. This implies that there exists another Hamiltonian which would describe the edge states at the right edge, j 1 ∼ L. Below, we derive the ESH at the right end. It is indeed one of the merit of our method to be able to treat the left end and the right end separately.
So far we have studied the system defined on a half plane, j 1 ≥ 1. In order to find the right ESH, let us consider the same system defined on the opposite half plane, j 1 ≤ −1, for the edge states at the right end, j 1 = −1. Such edge states correspond to those localized at the right end, j 1 = L, mentioned above. In this case, the boundary condition should be instead of Eq. (16). Let ψ jn = ψ 0n e iKnj with K n = k n + iκ n be the eigenstate near j 1 = −1. Then, κ n < 0 is required, and from the point of view of the lattice termination, the latter type of ψ 0n in Eq. (37) is suitable in this case, which naturally satisfied the boundary condi-tion (44). Thus, the eigenvalue equation is decomposed into This equation tells that the lower (q − 1) × (q − 1) part of the same Hamiltonian (35) can be the Hamiltonian of the right edge states as long as the condition, holds, which is obtained in the first line of the above equation. Namely, the edge states at the right edge are determined by the following Hamiltonian, with the condition Eq. (46). This is the right ESH, actually different from the left ESH (40). In Fig. 5, we show, by the green curves, the eigenvalues of the Hamiltonian (47) with the condition (46) satisfied. One can find that the red curves and green curves fully describe the edge states for the system of cylindrical geometry. In passing, we mention that the Hamiltonian (47) differs from (40) in the constant sift of k 2 → k 2 + φ. Therefore, the yellow curves in (40) can be the edge states if one chooses a suitable boundary condition.

1/3-flux
In the case of φ/(2π) = 1 3 , we can obtain analytic solutions of the edge states. It follows from Eq. (40) that the left ESH is given by the following 2 × 2 matrix, which gives the eigenvalues The parameter κ ± is found to be We show in Fig. 6, the edge spectrum obtained in Eq. (49) as well as the plot of e −κ± as functions of k 2 . One sees that the edge state in the lower (higher) gap is the physical state in −π < k 2 < 0 (0 < k 2 < π).

IV. GRAPHENE
Another example in 2D systems treated in this section is (spinless) graphene. In Sec. IV A, we consider the model in the absence of a magnetic field and derive the famous edge states along the zigzag edge, etc [27], based on our Hamiltonian operator formalism. The zero energy edge states are protected by reflection symmetry characterized by quantized Berry phase [35]. In Sec. IV B, we discuss the edge states of graphene in the presence of a magnetic field [28]. In the case of the zigzag edge or bearded edge, we derive the ESHs, as discussed in Sec. IV B 1. Since our method treats edge states at the left end and right end separately, we can reproduce the edge states of the cylindrical system with the zigzag edge at the left end and bearded edge at the right end by use of left and right ESHs. In spite of the same graphene, the system with the armchair edge belongs to class B or C in Sec. II C, implying that simple Hamiltonian formalism is impossible. Nevertheless, if one deforms a model with the armchair edge kept unchanged, one can derive the ESH which reproduces well the edge states of the undeformed model, as demonstrated in Sec. IV B 3. This is due to localization property of edge states along the boundary. A. In the absence of a magnetic field

Zigzag edge
Let us start with the Hamiltonian in the lattice representation. It follows from Fig. 7 (a) that the Hamiltonian operator appropriate for the zigzag edge reads Then, this operator becomes 1D Hamiltonian operator after the Fourier transformion in the 2-direction by replacing δ 2 → e ik2 : where δ ≡ δ 1 . Form this Hamiltonian operator, the hermiticity matrix K reads which is the same matrix as the NN SSH model (t = 0), implying that this model belongs to class A in Sec. II C. Thus, as in the case of the SSH model in Sec. II B, we assume a Bloch-type wave function (18) based on ψ 0 = ψ ↑ in Eq. (21) which satisfies (22). Then, we have This equation separates into two components such that The last equation leads to In Fig. 8 (a), we show the edge states at zero energy by red lines which exists in the range Eq. (56). This reproduces the results in Ref. [27].

Bearded edge
Similarly to the zigzag edge, it follows from Fig.  7 (b) that the 2D Hamiltonian operator and Fouriertransformed 1D Hamiltonian operator appropriate for the bearded edge is given bŷ Then, the hermiticity matrix K reads Although K depends on k 2 , this is basically the same as (53), and this model is also a NN model in class A. Therefore, we can choose ψ 0 = ψ ↑ satisfying Eq. (22). Form the eigenvalue equation for the same Bloch-like wave function as in Sec. IV A 1, we have Hence, it follows from that the range of k 2 for the bearded edge is the complement of that for the zigzag edge (56). Lastly, in the case of the armchair edge, Fig.7 (c) leads us toĤ

Armchair edge
Then, the hermiticity matrix K reads which is the same as that for the generalized SSH model with t = 0 in Sec. II B 4, so that this case belongs to the class B in Sec. II C. Therefore, it is not enough to choose ψ 0 = ψ ↑ or ψ ↓ ensuring the hermiticity of the Hamiltonian: We need to find two independent solutions to obtain the wave function satisfying the boundary condition Eq. (16). When one chooses ψ 0 = ψ ↑ , the eigenvalue equation leads to and when one chooses ψ 0 = ψ ↓ , These equations for ψ 0 = ψ ↑,↓ give the same solutions for real part, |e iK | = e −κ . As shown in Fig. 9, we find one solution in e −κ < 1 and the other in e −κ > 1 at any k 2 , implying that it is impossible to get the wave function satisfying the boundary condition (16). Therefore, we conclude that the armchair edge allows no edge states. This also reproduces the results in Ref. [27]. It should be stressed that all the difference of the edge structure are incorporated in the difference of the Hamiltonian operators in Eqs. (52), (57), and (61).

B. In the presence of a uniform magnetic field
The edge states and the bulk-edge correspondence of graphene in a uniform magnetic field were studied in detail in Ref. [28], in which a model with hybrid edges, i.e., the zigzag edge at one end and the bearded edge at the other end was examined. The merit of our method developed in this paper is that we can examine these edge states separately. Let us first consider the zigzag edges.

Zigzag edge
For the choice of the magnetic unit cell and the labeling of sites in it in Fig. 10, the Hamiltonian operator with the zigzag edge at the left end under uniform magnetic flux φ/(2π) = p/q per hexagon is given bŷ wherê Here,ĥ φ = e iφ δ * 2 denotes the hopping towards the 2direction, and the hermiticity matrix K, which is given by the coefficient of δ * inĤ, is the same as that of the Hofstadter model in Eq. (36): Therefore, the reference wave function ψ 0 is given by Eq. (37), and assuming a Bloch-like edge state based on the reference state ψ jn = ψ 0n e iKnj with K n = k n + iκ n , the eigenvalue equation becomes where where D e is q × (q − 1) matrix defined by The equation of the last component in Eq. (68) yields the constraint e iKn χ 1,n + g * −k2 χ q,n = 0, from which it follows that the localization condition is given by Thus, we have successfully derived the ESH describing the edge states along the zigzag edge at the left end, which are completely separated from the bulk states. It should be stressed again that not only the Hamiltonian given by Eqs. (69) and (70) but also the localization condition (71) describe the edge states. We show in Fig. 11, two examples of the QHE of graphene. We plot, by red curves, the eigenvalues of the Hamiltonian (69) satisfying (71) on the background of the total spectra for the cylindrical systems with the zigzag edges at both ends. Half of the edge state spectra are well reproduced, from which we find that they are edge states localized at the left end. One of characteristic properties of the edge states of graphene is the existence of the zero energy flat bands [28]. Actually, in Fig. 11, at several intervals of k 2 , we find such bands. The zero energy edge states can be readily understood from the ESH in Eq. (69): Because of the reference state ψ 0 satisfying the hemiticity, the q × q square matrixD in Eq. (66) reduces to q × (q − 1) imbalanced matrix D e in Eq. (70). This reveals the existence of one zero-energy edge state at each k 2 , as long as the condition Eq. (71) is satisfied. To be more specific, let us write down the eigenvalue equation for the ESH, Note that D e D † e and D † e D e are q × q and (q − 1) × (q − 1) matrices and they have the same energy eigenvalues except for zero energy, from which it follows that the zero energy edge states are included in the upper part, D e D † e . Other nonzero energy edge states are obtained by diagonalizing the lower part, D † e D e η n = ε 2 n η n . Using Eq. (72), we find that the eigenfunction χ n for nonzero energy states is given by It then follows that the localization condition Eq. (71) can be rewritten as Therefore, if D † e D e has no zero energy eigenvalues, D e D † e has one zero energy state. Such a zero energy state should obey D † e ζ 0 = 0 and η 0 = 0, and the localization condition is given by Eq. (71), i.e., For small q systems, we can obtain exact spectra as well as exact wave functions of the edge states. These are so didactic that we describe separately below. • π-flux Let us first consider case with π flux per unit hexagon (p/q = 1/2). Although we have already derived the ESH for generic q, it may still be convenient to write down the total Schrödinger equation (68), since this includes the condition to determine e iK : where g π±k2 = 1 − e ±ik2 and g ±k2 = 1 + e ±ik2 , and the energy quantum number has been suppressed. The upper 3 × 3 matrix determine the spectrum of the edge states, from which we find 2 × 1 D e D e = g π−k2 1 .
The zero energy edge state lives in ζ sector satisfying D † e ζ 0 = 0, so that we set η 0 = 0. Then, we have the following edge state at zero energy, and e iK0 in this case in Eq. (76) or directly from the last component in Eq. (77) becomes which yields 0 < |k 2 | < π/6, 5π/6 < |k 2 | < π. In Fig. 12 (a), the spectra Eqs. (81) and (82) with (83) obtained from the ESH are shown as red curves on the background of the total spectrum of the cylindrical system with the zigzag edges at both ends. This figure suggests that the left and right ends yield the (non)zero energy edge states at the same (different) range of k 2 .

•2π/3-flux
In the case of φ/(2π) = 1/3, the eigenvalue equation for the edge states is where ξ n and η n denote three-and two-component vectors, respectively. It follows that D e matrix reads The nonzero energies of the edge states are thus determined by which gives two eigenenergies On the other hand, the zero energy state is determined by D † e ζ 0 = 0, from which one finds ζ 0 ∝ (1, −g * φ−k2 , g * φ−k2 g * −φ−k2 ) T . Therefore, the localization condition (76) is It follows that around three points k 2 ∼ π, ±π/3, the zero energy state appears. We show in Fig. 12 (b) the energies (87) and zero energies by red curves, whose localization conditions are numerically computed. One sees that half of the edge states are well reproduced, which are localized at the left end. Let us next consider the bearded edge in the presence of a uniform magnetic flux. When we choose the magnetic unit cell as illustrated in Fig. 13, we have the Hamiltonian Eq. (65) by use of the following new oper-atorD Thus, the K matrix is basically the same as the zigzag edge, although the single nonzero component depends on k 2 . Therefore, from a similar eigenvalue equation to Eq. (68), we find the Hamiltonian of the edge states at the left end is basically the same structure as Eq. (69), where q × (q − 1) matrix D L e is in this case defined by This is the ESH at the left end. Here, in Eq. (90), we have indicated the numbers of row and columns explicitly to compare the ESH at the right end below. The localization condition is Next, let us consider the right edge states. In the same way as in Sec. III A 2, we choose the reference states as in Eq. (45). Assuming a Bloch-like wave function for the edge states ψ jn = ψ 0 e iKnj , we find the following Schrödinger equation, This equation, from its lower part, leads to the following ESH at the right end, The localization condition, which is the equation of the first component in Eq. (93), becomes |e −iKn | = e −κn = χ q,n g k2 χ 2q−1,n < 1.
So far we have derived the two ESHs with the bearded edges, one is at the left end and the other is at the right end. In the previous Sec. IV B 1, we have also derived two ESHs with the zigzag edges. These are independent, since we have considered a half plane which has only one end to derive each ESH. In order to check this explicitly, we introduce a model with the zigzag edge at the left end and the bearded edge at the right end, which was studied in Ref. [28]. Figure 14 shows various spectra of such a model. The background dots are total energy eigenvalues of the cylindrical system with the zigzag and bearded edges. The red curves are eigenvalues of the left ESH for the zigzag edge defined by Eqs.   (69) and (94), respectively, whereas the background gray dots represent the total eigenvalues of the system of the cylindrical geometry with the zigzag and bearded edges.

Armchair edge
Irrespective of the presence or absence of a magnetic field, the choices of the unit cell for the zigzag and bearded edges lead to the Hamiltonian operators with NN hoppings only. This enables us to extract the exact ESHs at the left and right ends separately. In the case of the armchair edge, on the other hand, the model belongs to class B or C in Sec. II C, since the hermiticity matrix K has two nonzero matrix elements. In the absence of a magnetic field, the model is in class B. Indeed, we have a unique matrix A that anti-commutes with K. As a consequence, we have been able to solve the eigenvalue equation with the hermiticity condition imposed. However, if a magnetic field is introduced, the problem is much more difficult. One of the reasons is that the matrix A is not unique. Even if we treat the model as the one in class C, it is quite difficult to carry out numerical search for solutions in a huge parameter space.
In this section, we provide a different perspective on the edge states at the armchair edge. Let us try to deform the model so that the Hamiltonian belongs to class A without changing the edge. If such a deformation is far from the edge, edge states are not affected so much. Of course, such a deformation would modify the topological properties of the original model. Nevertheless, it could be useful to reproduce the edge state bands located at one end using a simple ESH. In this section, we propose such an attempt.
The Hamiltonian operator with the armchair edge at the left end is the same chiral structure as Eq. (65), whereD operator in the present choice of the magnetic unit cell in Fig. 15 is given bŷ with r = 1 for the uniform honeycomb lattice. As in the case of zero field, this model cannot be a simple NN model classified as A, since the K matrix is It follows that if one sets r = 0, above K matrix reduces to the one for the zigzag edge (67), and the edge states can be obtained in the same way as the zigzag and bearded edges. In this case, r = 0, it is easy to see that the ESH becomes Eq. (69) with In Fig. 16, we show spectra of graphene with the armchair edges. Fig. 16 (a) is the full spectrum of the uniform system with r = 1 of cylindrical geometry, which is compared with (b) of the system with r = 0. These two systems are different topological band structures, as can be seen from Figs. 16 (c) and (d), which are the Chern numbers as functions of the fermi energy [36,37] computed according to Ref. [37]. Nevertheless, the behavior of the edge states are rather similar. This is because for large q system, the left end and the defects at r = 0 bonds are well separated, so that edge states localized at the end and the impurities states localized at defects are almost decoupled. Although these edge states cannot be used for the discussion of the topological arguments such as the bulk-edge correspondence, they may be helpful, e.g., to determine at which end the edge states of the uniform system are localized, etc.  Fig. 16 (a) overwritten by edge state spectra (red curves) in Fig. 16 (b). (b) The same as (a) in the case of φ/(2π) = 1/5.
In fact, it may be interesting to overlay the edge states (the red curves) of the defect model in (b) and the full spectrum of the uniform system in (a), which is carried out in Fig, 17 (a). Of course, the red curves include impurity states, for example, those at zero energy near k 2 ∼ ±π are localized states at the defects of the r = 0 bonds. Nevertheless, overall spectrum of the edge states are well reproduced. On the other hand, for small q systems, one expects that the coupling between the edge states and impurity states becomes stronger. In Fig. 17  (b), we show the model with flux φ/(2π) = 1/5. One can indeed observe slight deviation of the red curves from the background spectrum, as expected.

V. WILSON-DIRAC MODEL
So far we have investigated 2D models mainly in class A defined in Sec. II C. In this section, we study a simple example, Wilson-Dirac model, in class B which allows some matrix A anti-commuting with the hermiticity matrix K.
The Wilson-Dirac model is a typical model whose edge states have been obtained analytically [9,38]. The Hamiltonian operator for the Wilson-Dirac model is given byĤ where in the second line j 2 has been Fourier-transformed by δ 2 → e ik2 , and δ ≡ δ 1 acts on j ≡ j 1 . Without loss of generality, we assume t > 0. Then, from the coefficients of δ * , we find For this K matrix, one finds {K, σ 2 } = 0, so that this model is classified into B. Actually, the K matrix (101) is unitary equivalent to that of the SSH model in the case of t = 0 in Eq. (20). Thus, one can choose the reference state as ψ 0 = ψ ± , where ψ ± is the eigenstates of σ 2 , σ 2 ψ ± = ±σ ± . Using these states, let us introduce a Bloch-like wave function for the edge states, ψ j = ψ ± e iKj , where K = k + iκ. Here, it should be noted that for a given model parameters, either one of ψ ± can be the wave function, provided that it allows two independent solutions, K 1 and K 2 , for the boundary condition (16) to be satisfied. If both ψ ± do not allow two solutions, one concludes that the model has no edge states. Substituting a Bloch-like wave function into the eigenvalue equation,Ĥψ j = εψ j , one has where M (k 2 ) = cos k 2 − 2 + m b , and the energy quantum number n has been suppressed. One sees that ψ ± is indeed the eigenstate ofĤ with energy eigenvalue ε = ±t sin k 2 , provided that The imaginary parts of the above equation yields In what follows, we calculate the case ψ 0 = ψ + (the case of upper sign in the above equations). Let us substitute the solutions (104) into the real part of Eq. (103). Then, one has where ± in the upper equation corresponds to k = 0, π. It then turns out that there exist two solutions when |M (k 2 )| < 1 for 0 < b and when |M (k 2 )| < 1 − (t/b) 2 for b < −t, and otherwise, two solutions are not allowed. Likewise, in the case of ψ 0 = ψ − , an edge state exist when |M (k 2 )| < 1 for b < 0 and when |M (k 2 )| < 1 − (t/b) 2 for t < b. As a consequence, edge states exist when |M (k 2 )| < 1, and the condition that there exists k 2 that satisfies |M (k 2 )| < 1 becomes This determine the topological phase of the Wilson-Dirac model. Moreover, from the types of K, one can get more information on the wave function of the edge state. In addition to |M (k 2 )| < 1, further conditions below classify the types of K as follows: For |b| < t, the two solutions K 1,2 are of the types and for |b| > t, and for k 2 that satisfies M (k 2 ) < 0, the wave function of the edge state is whereas k 2 that satisfies M (k 2 ) > 0, In Fig. 18, we show some examples that have different types of edge states dependent on k 2 in one model.

VI. HALDANE MODEL
In this section, we study the Haldane model as an example in class C in Sec.II C. For practical purpose, the present method may not be very useful in order to obtain the edge states in class C. Nevertheless, we would like to claim that the hermiticity of Hamiltonians would determine the edge states in principle.
From the unit cell defined in Fig. 19, the Hamiltonian operator readŝ whereĤ 1,2 denote the NNN hopping operators defined byĤ In order to find the edge states localized near the left end j 1 = 1, the 2-direction is Fourier-transformed to obtain where δ ≡ δ 1 operates to j ≡ j 1 and m 1 = m + 2t cos(k 2 + φ), m 2 = −m + 2t cos(k 2 − φ), The coefficient of δ * leads to the following K matrix, This matrix never anti-commutes with any other nontrivial matrix, so that the model belongs to class C. Thus, there is no simple way to choose ψ 0 ; rather, it should be determined using the hermiticity condition. To be more specific, we first set the wave function of the edge state ψ j = ψ 0n e iKnj (K n = k n + iκ n ) and require the hermiticity where we have set ψ 0n = (1, χ n ) T for simplicity. Under this condition, we solve the eigenvalue equation Thus, we have three equations (116) and (117) for three unknown parameters K, χ, and ε. When these equations have degenerate two solutions with |e iKn | = e −κn < 1, we can obtain the wave functions satisfying the boundary condition. Since the equations to be solved are complexvalued nonlinear equations, it is not so easy to determine whether they allow two independent solutions or not. In Fig. 20, we show the spectra of the edge states by red dots, which are obtained by numerical calculations. Figure 20 (a) and (b) are in the nontrivial and trivial phases, respectively. Even in the trivial phase, there appear edge states localized at boundaries, which can actually be captured by the present method based on the hermiticity of the Hamiltonian (116).

VII. HIGHER-ORDER TOPOLOGICAL INSULATORS
This section is devoted to the application of our formulation to the HOTIs. The HOTIs have been attracting much current interest, providing us a new concept of the bulk-edge correspondence. We study two typical models on the square lattice [11,12,30] and on the breathing kagome lattice [31]. Let us start with the square lattice with bondalternating hoppings towards the 1-and 2-directions, as illustrated in Fig. 21. On this lattice, we introduce zero flux and π-flux, which defines two different models. These models have nontrivial gapped edge states with nontrivial polarizations, and in particular, in the so-called BBH model with π-flux, these edge states yield corner states within a bulk gap [11,12,30].
From Fig. 21, it follows that the Hamiltonian operator is defined byĤ where∆ j = γ j + λ j δ j and∆ * j = γ j + λ j δ * j , and − signs representing π-flux are for the BBH model. From the coefficients of δ * j , we derive the K matrices for the 1-and 2-directions, Since each K matrix has two nonzero matrix elements, the models seem not to belong to class A in Sec. II C. However, due to the particularity of the models, we can apply our formulation to these models as class A, as shown below.

Edge states
To begin with, let us consider the system defined on the half plane, j 1 ≥ 1 and discuss the edge states along the left end. The Fourier transformation toward the 2direction can be carried out by replacing δ 2 → e ik2 , Then, when we choose where χ n = (χ 1,n , χ 2,n ) T is a two-component vector, the hermiticity of the Hamiltonian and the boundary condi-tion at j 1 = 1 are ensured by K 1 ψ 0n = 0. Generically, such a reference state yields two constraints on the eigenfunction of the ESH, which allows no solutions. However, for the present models, two constraints reduce to only one.
Generically, two conditions on the eigenfunctions χ n are incompatible with the eigenvalue equation of the Hamiltonian (122), H e χ n = ε n χ n . However, in the present case, the two conditions (123) reduce to one, which imposes moreover no constraint on the eigenfunction, i.e., γ 1 + λ 1 e iK1,n = 0. Thus, it turns out that when the edge states localized near j 1 = 1 exists, which is governed by the SSH Hamiltonian toward the 2-direction (122), regardless of whether the model includes π-flux or not. The difference of the two models lies in the bulk spectrum. Since there are no constraints on the wave functions χ n , all of the eigenstates of the ESH (122) are edge states of the present models, whose eigenvalues are ε ± (k 2 ) = ±|γ 2 + λ 2 e ik2 |.
In Fig. 22, we show by red curves the spectra ε ± (k 2 ) of the edge states (122). In both figures (a) and (b), they are of course the same and the difference is just the bulk spectra shown by gray dots in the background. The edge states for the cylindrical system, i.e., red curves in Fig.  22, are doubly-degenerate: One is at the left end, and the other is at the right end. In the 0-flux model, these edge states are partially embedded in the bulk bands for larger γ j -parameters. The background gray dots denote the total spectrum obtained numerically for cylindrical system.

Corner state
In order to study the corner state, we consider the system on the quarter plane defined by j 1 ≥ 1 and j 2 ≥ 1. For the hermiticity of the Hamiltonian, we have to impose two conditions associated with K 1 and K 2 in Eq. (119) on the reference wave function ψ 0 , K 1 ψ 0 = K 2 ψ 0 = 0. Then, we have where χ is just one component vector. Generically, these conditions are too strong to yield some meaningful solutions. However, the present models actually allow one solution, which is nothing but the corner state, as shown below. Let ψ j = ψ 0 e iK1j1+iK2j2 be the wave function of the corner state. For such a Bloch-like wave function, the eigenvalue equation of the Hamiltonian (118) becomes Thus, it turns out that the Hamiltonian for the corner state is the middle diagonal element denoted explicitly by 0 in the above equation, provided that the third and fourth components are satisfied for χ = 0: which reduces to Namely, a single corner state appears at ε = 0 for the system with the parameters (128). It is indeed a in-gap state for the BBH model, whereas it is embedded in the bulk band for the zero flux model. As another example of corner states, we consider the model reported in Ref. [31]. The unit cell illustrated in Fig. 23 leads to the Hamiltonian operator given bŷ from which the K matrices for 1-and 2-directions read We have written them as operator forms, but when we apply them to edge or corner states, δ 2 inK 1 is replaced by e ik2 , and so on. As already mentioned, two or more nonzero components of K matrices reduce models into class B or C generically. A major feature of Eq. (130) is that although each K matrix includes two nonzero elements, they are alined vertically. This enables us to treat this model as the one in class A, as shown below.

Edge states
Let us consider the system defined on the half plane j 1 ≥ 1, and study the edge state localized around the left end j 1 = 1. Then, we can replace the shift operator δ 2 by e ik2 . It follows from Eq. (130) that the K 1 matrix allows the following reference state which satisfies the hermiticity Eq. (10) or (20) where χ n is a two-component vector. Note that this state satisfies K 1 ψ 0n = 0, implying that the model with the present edge, i.e., a straight line toward the 1-direction, belongs to class A. Assume a Bloch-like wave function of the edge state as ψ j1n = ψ 0n e iK1,nj1 using Eq. (131). Then, the eigenvalue equation becomes Thus, the ESH is the upper 2 × 2 matrix defined by which gives the eigenvalues ε ± = ±|γ + λe ik2 |. This is the same SSH Hamiltonian as the ESH of the model on the square lattice in Eq. (122). However, the condition on the wave function given by the last component in (132) yields nontrivial constraint, (γ + λe iK1,± )χ 1,± + (γ + λe iK1,± e −ik2 )χ 2,± = 0, from which it follows |e iK1,± | = e −κ1,± = γ λ ε ± − (γ + λe ik2 ) ε ± − (γe −ik2 + λ) < 1. (134) In Fig. 24, we show the eigenvalues of Eq. (133) with the  The background spectra of the cylindrical system in Fig. 24 show the edge states at the right end as well. For example, in Fig. 24 (a), there appear a band in between the edge state at the left end and the flat band. This is the edge state at the right end. Since we consider the unit cell illustrated in Fig. 23, the right end has a zigzaglike shape. Unfortunately, the problem to determine the right edge states of such a shape is of class C type in Sec. II C, and hence, it may be too difficult for the present approach.

Acute angle corner state
Suppose that the system is defined on j 1 ≥ 1 and j 2 ≥ 1. This corresponds to the corner with angle π/3 on the kagome lattice. In this case, the reference state ψ 0 can be chosen as where χ is a one component vector. This reference state satisfies both K 1 ψ 0 = K 2 ψ 0 = 0. Assuming a Blochlike state ψ j1,j2 = ψ 0 e i(K1j1+K2j2) on the reference state (135), the eigenvalue equation becomes This equation tells that the upper diagonal component in the Hamiltonian matrix denoted explicitly by 0 is the corner state Hamiltonian. Therefore, a single corner state localized around the corner j = (1, 1) exist at zero energy, provided that the following conditions are satisfied, which are due to the second and third components in (136). While it appears within a bulk gap for −1 < γ/λ < 1/2, it is embedded in the bulk states for 1/2 < γ/λ < 1 [31].

Obtuse angle corner state
In order to discuss whether the model allows a corner state at the corner with obtuse angle 2π/3, we replace the unit cell by the dashed box in Fig. 24. The Hamiltonian operator is then from which the K matrices read Edge states: In exactly the same way as in the previous Sec. VII B 1, a single edge state exists, governed by the same Hamiltonian H e in Eq. (133) but with γ ↔ λ and k 2 → −k 2 , and by the same localization condition (134), which of course gives the same energies, as it should be.
Corner state: Let us consider the system defined on j 1 ≥ 1 and j 2 ≤ −1. The reference state should satisfy K 1 ψ 0 = K † 2 ψ 0 = 0, from which it follows that the same state as in Sec. VII B 2 (135) can be the reference state. However, due to the different localization condition, the corner with the obtuse angle has no corner state for any parameters.

VIII. SUMMARY AND DISCUSSION
We studied the edge state of various 2D tight-binding systems with particular emphasis on the hermiticity of the first-quantized Hamiltonian operators. We showed that even for tight-binding models defined on latices, Hamiltonians described by the shift operators, i.e., discrete versions of the differential operators are quite useful. Then, it turned out that with boundaries, Hamiltonians on lattices are also not necessarily hermitian, as is often the case with the Dirac Hamiltonians in the continuum space. We pointed out that the hermiticity of Hamiltonians has close relationship with boundary conditions. In particular, for models with NN hoppings only, both conditions can be satisfied simultaneously by specific wave functions. It then follows that we can develop a Bloch-like theory for edge states based on such wave functions. This enables us to extract Hamiltonians describing edge states at one end from the total states including the bulk contributions. Namely, we can study edge states appearing at one end for semi-infinite system, or in other words, we can treat edge states at the left and the right ends separately.
On the other hand, for models including NNN hoppings, we found it still difficult to apply our formulation. This may be due to the fact that while the boundaries of NN models are realized by simple lattice terminations, for those including NNN hoppings, how to realized the boundary conditions is not so obvious. Nevertheless, we showed, using Haldane model, that even if the boundary conditions are not so clear, the hermiticity of Hamiltonians serves as a guiding principle to determine edge states of tight-binding models defined on lattices.
Our formulation would have many potential applications. For example, the ESHs derived by our method enable us to investigate the bulk-edge correspondence for various kind of tight-binding models, as carried out by Pletyukhov et. al. for a variant of the Hofstadter model [23]. In the present paper, we focused our attention to technical details of deriving ESHs for well-known models. On the other hand, in the discussion of the bulk-edge correspondence, it is important to reveal how ESHs are embedded in the bulk Hamiltonians. In the subsequent paper, we would like to clarify this point [34].
As demonstrated by Kunst et. al. [20][21][22], it may also be interesting to construct various models which yield ESHs by exploring the first quantized Hamiltonian oper-atorsĤ and corresponding K matrices. Even if models do not look NN models, some specific forms of K matrices allow ESHs, as shown in Sec. VII. Based onĤ, one can construct a lot of tight-binding models to give ESHs.