The Existence of Topological Edge States in Honeycomb Plasmonic Lattices

In this paper, we investigate the band properties of 2D honeycomb plasmonic lattices consisting of metallic nanoparticles. By means of the coupled dipole method and quasi-static approximation, we theoretically analyze the band structures stemming from near-field interaction of localized surface plasmon polaritons for both the infinite lattice and ribbons. Naturally, the interaction of point dipoles decouples into independent out-of-plane and in-plane polarizations. For the out-of-plane modes, both the bulk spectrum and the range of the momentum $k_{\parallel}$ where edge states exist in ribbons are similar to the electronic bands in graphene. Nevertheless, the in-plane polarized modes show significant differences, which do not only possess additional non-flat edge states in ribbons, but also have different distributions of the flat edge states in reciprocal space. For in-plane polarized modes, we derived the bulk-edge correspondence, namely, the relation between the number of flat edge states at a fixed $k_\parallel$, Zak phases of the bulk bands and the winding number associated with the bulk hamiltonian, and verified it through four typical ribbon boundaries, i.e. zigzag, bearded zigzag, armchair, and bearded armchair. Our approach gives a new topological understanding of edge states in such plasmonic systems, and may also apply to other 2D"vector wave"systems.


Introduction
Recent progress in topological phases [1][2][3][4][5][6][7] of condensed matters stimulates the studies of non-trivial topological properties in various physical systems [8][9][10][11]. In 2005, Haldane and Raghu introduced the topological concept in relation to the quantum Hall effect into photonic systems [12,13]. And then many physical phenomena in condensed matter physics which are difficult to be observed in the past have found their classical analogies [14][15][16][17][18][19][20]. Artificial graphene [21] is one of the typical examples, which offers a platform to investigate interesting phenomena related to the massless Dirac fermion dynamics near the Dirac cone [22]. As shown by recent studies, various physical systems possess the spectra similar to the π bands of electrons in graphene [23,24], such as 2D electron gas in quantum well [25], ultra cold atoms in optical lattices [26], photonic crystals [27][28][29][30][31][32][33], and plasmonic lattices composed of metal particles [34,35]. Although π bands in graphene are merely composed of p z orbital with only one freedom, it is also possible to construct band structures with multiple degrees of freedom and go beyond the "scalar wave" nature in some systems with honeycomb lattice. A simple idea to achieve this goal is to build "vector wave" bands by including the p x,y orbital of the lattice. Such schemes have been discussed for ultra cold atoms [26], polaritons of coupled micropillars [36], the TE (or in-plane polarized) modes in photonic systems [27,34], and the in-plane vibration modes in classical spring-mass systems [37].
More interestingly, there exist symmetry-protected edge states between the upper and lower π bands in graphene-like ribbon systems [38][39][40], which have been confirmed experimentally both in graphene [41] and artificial graphene systems [31,33,42]. According to previous theoretical studies [43][44][45][46], the emergence of edge states is associated with the chiral symmetry of the bulk Hamiltonian and can be determined by the Zak phase of the π bands, which is well-known as bulk-edge correspondence [7,[45][46][47][48][49]. However, previous analyses of the "vector wave" (in-plane polarized) bands of honeycomb lattices showed neither the existence ranges nor the number of the midgap edge states is identical with those of graphene [27,34], although they both possess the Dirac-like dispersion at the vertices of the first Brillouin zone. This is a strong hint that scalar wave bands and vector wave bands possess different topological properties, and the principle of bulk-edge correspondence for "scalar waves" should be modified in order to be used for the multifreedom systems. Therefore, a proper relationship of edge states and topological properties of the bulk bands of "vector wave" system is really needed.
In this paper, we mainly investigate the in-plane polarized modes in honeycomb plasmonic lattices (HPLs), as a representative example of vector wave systems. The interaction is simplified to form a concise Hermitian eigenvalue problem by using the coupled dipole method (CDM) and neglecting the dissipation and the retardation of the fields. Band structures of the infinite lattice and ribbons with four typical boundaries for out-of-plane and in-plane polarizations are calculated respectively in the tight-binding approximation. And the results taken the retardation and long range interactions into consideration also verify the effectiveness of our simplified model. It is shown that the ribbons with different boundary conditions possess different distributions of the edge states. Moreover, we develop a theorem of bulk-edge correspondence for the in-plane polarized modes which reveals the relations between the Zak phase, the winding number associated with bulk Hamiltonian and the existence of the zero-eigenvalue flat edge states in ribbons with various boundaries. Our result could be regarded as an extension of the theorem for graphene [44]. The existences of edge states predicted by our theorem are perfectly consistent with the calculated band structures in the four typical ribbons respectively.

Honeycomb plasmonic lattice and different boundaries
Metal nanoparticle lattice is a particular type of photonic systems which supports the collective resonance of localized surface plasmons of individual particles and have distinctive band dispersions [50][51][52][53][54][55][56]. Several studies have reported the interesting topological phenomena emerging in this type of systems [34,[57][58][59][60]. Specially, 2D honeycomb plasmonic lattices are made up of one-layer-thick metal particles and arranged in honeycomb structure just like graphene, which can be regarded as a compound structure nested with two triangular sublattices. The Bravais lattice is shown in Fig. 1(a), where a 1 , a 2 are two basis vectors, and two neighboring "atoms" consist of a unit cell (dimer). And a ribbon is a lattice with finite width in one direction with two parallel edges and is infinite along the direction of edges. Actually, there are myriad different types of ribbon edges, and a large class of ribbon edges can be characterized by the periodic vector T(m, n) = ma 1 + na 2 , where (m, n) is a pair of coprime integers [44]. According to Delplace's work [44], to obtain the boundary of a ribbon with the periodic vector T(m, n), we should firstly translate a dimer with the displacement a 1 (−a 1 if m < 0) |m| times (each time, we obtain a dimer on the edge), then translate with a 2 (−a 2 if n < 0) |n| times, and so on and so forth. There are three different orientations for a unit cell, the one adopted in our paper is shown in Fig. 1(a), the other two can be obtained from the former by a ±2π/3 rotation. For infinite systems, there are no differences between these three choices except an additional phase factor of the bulk Hamiltonian. However, they give rise to different results in finite systems even constructing the edge of a ribbon with the same T(m, n), in addition, the results are also different when calculating the Zak phases and winding numbers of bulk bands. Specifically, Fig. 1(b) shows the most studied four types of ribbon boundaries, i.e. zigzag, armchair, bearded zigzag, and bearded armchair. One choice (but not unique) of the periodic vectors of these four boundaries in Fig. 1(b) are, respectively, Zigzag: Bearded zigzag: Bearded armchair: T ba = T(2, −1).

Quasi-static and dipole approximation
In this paper, we use CDM to simplify the interactions between metal particles with treating metal particles as point dipoles, which means the results are effective only for r s /a 0 ≤ 1/3 [51][52][53], where r s and a 0 are the particle radius and the center-to-center distance between two adjacent particles respectively. The electric field generated by a harmonically radiating point dipole is characterized by the equation where r is the displacement vector from the source point to the field point, n is the unit vector along r, P is the dipole moment, k = ω/c is the wave vector (here ω is the angular frequency), and the surrounding medium is assumed to be vacuum. When the wavelength λ a 0 , the quasi-static approximation (QSA) can be applied, then Eq. (1) reduces to with the quasi-static Green tensor where ↔ I is identity matrix. For simplicity, we assume the metal particles are all spherical. The induced dipole moment of the particle on the lattice point R (Coupled Dipole  Equation) is where E R ,R denotes the electric field at R generated by the dipole on the lattice point R ,α(ω) = α(ω)/(4πε 0 ) = ε(ω)−1 ε(ω)+2 r 3 s is the polarizability of the spherical particles, ε(ω) is the frequency-dependent relative permittivity of the metal particles, and is the two-point correlation function. Indeed, the coupled dipole equation can be transformed to an eigenvalue problem whereα(ω) −1 serves as the eigenvalue.

Eigenvalue problem and Hamiltonian
Since the honeycomb lattice is a compound lattice with two inequivalent sublattices (see Fig. 1), we can combine the dipole moments of two sublattice points in one unit cell to a single polarization vector P R = (P A R , P B R ) T , then Eq. (5) can be written as an Hermitian eigenvalue equation where R(R ) denotes the index of a unit cell, the Hamiltonian element H R,R = H(R − R ) correlating the two cells R and R has the block form with ∆R = R − R . The hopping terms H AA and H BB represent the interaction from the same sublattices A and B respectively, between different unit cells, and they take the forms H AA Eq. (6) and Eq. (7) give the eigenequation and the component form of the Hamiltonian in Wannier basis |R . To consider the common eigenstates of both the Hamiltonian and the lattice translation operator subjected to the periodic boundary condition, we can obtain the Bloch state |ψ k = P R (k)|R with the component P R (k) = e ik·R ψ k , where k is the Bloch wave vector corresponding to the eigenvalue of the lattice translation operator. On the basis of Bloch states, the eigenequation of the Hamiltonian becomes H k · ψ k =α(ω) −1 ψ k with the component of the Hamiltonian in Bloch representation H k = ∆R H(∆R)e −ik·∆R . According to the tight-binding model (TBM), the Hamiltonian with only the nearest-neighbour interaction is where , t i (i = 1, 2, 3) denotes three nearest hopping vectors respectively (see Fig. 1). It is easy to know that the Hamiltonian possesses chiral symmetry: This property indicates that there may exist zero-eigenvalue edge states in corresponding ribbon systems [43]. Actually, the tight-binding Hamiltonian given in Eq. (9) is quite similar to the π electron Hamiltonian in graphene except the nearest neighbor hopping constant is replaced by the Green tensors ↔ G(t i ). Nevertheless, this difference has remarkable influences on band structures and the topological properties of the edge states.
Being the eigenvalue of an Hermitian matrix,α(ω) −1 is required to be real. However, if the permittivity ε(ω) of metal particles has imaginary part, a realα(ω) would lead to a complex frequency ω which means that the plasmonic resonance of the particles would decay over time. Without changing the topological properties, the imaginary part of the permittivity is omitted and the Drude model without damping ε(ω) = 1 − ω 2 p /ω 2 is used to characterize the metal paricles, where ω p is the plasmonic resonant frequency. Then we can obtain the relation betweenα and ω where ω 0 = ω p / √ 3 is the dipole resonant frequency of spheres. For ω > 0, α(ω) is a monotonic function, thus the correspondence between the eigenvalues and eigenfrequencies in Eq. (6) is a one-to-one mapping [59].

Bulk band structures of HPLs
Under the previous simplification, the band structure of honeycomb plasmonic lattice can be easily obtained from Eq. (6). Because H k is a 6 × 6 matrix, the system has 6 bands in total. On the other hand, the in-plane polarization P xy of the dipoles is orthogonal and decoupled with the out-of-plane polarization P z , since the Green tensor can be decomposed as the direct sum of the in-plane and out-of-plane parts xy )/r 3 and G z (r) = −1/r 3 for all r in xy plane. Therefore, it makes physics more clear if we investigate the out-of-plane and in-plane polarized eigenspectra of the two directly summed parts of the Hamiltonian The out-of-plane part of the Hamiltonian has a 2D massless Dirac-like form: where p(k) = − 1 is the triplet of Pauli matrices, and here p(k) = (Re p(k), −Im p(k), 0). It has two eigenvalues and the corresponding eigenstates are The Hamiltonian of the in-plane modes is a 4 dimensional matrix where where ψ and u are both normalized, and e −iγ is an additional phase, moreover, u satisfies the corresponding con-eigenvalue problem In addition, ψ is also the eigenvector of (H in k ) 2 = diag(AA † , A † A), thus u is the solution of the following equation where the effective Hamiltonian H eff = AA † which can also be expressed in the Dirac-like form Naturally, we get two eigenvalues of H eff and the corresponding eigenvectors where θ, ϕ are the components of g = (|g| , θ, ϕ) in spherical coordinates. Note that ϕ is not well-defined when g is along z axis, hence the phase e iϕ (e −iϕ ) in Eq. (20) is chosen to ensure the continuity of u 1 (u 2 ) at g x = g y = 0. According to Eq. (19), there are four eigenvalues of the ordinary Hamiltonian Eq. (14): which represent the four bulk bands of the in-plane modes. And the eigenvectors of these four bands are, respectively, where γ j can be determined by the argument ofλ j (λ j = u T j A † u j ). The bulk band structures of out-of-plane and in-plane polarizations are shown in Figs. 2(a) and 2(b) respectively. Here and in the following calculations, the parameters are chosen to be ω p = 6.18 eV, r s = 10 nm, and r s /a 0 = √ 3/6. Both of the in-plane and out-ofplane bands possess Dirac cones located at the K and K points (the corners of first Brillouin zone), and all of the Dirac points have the same frequency ω 0 corresponding to the zero-eigenvalue (α(ω 0 ) −1 = 0). For the out-of-plane case, the eigenvector on each "atom" has a fixed polarization along z direction just like the scalar wave function composed of p z orbital of electrons in graphene, therefore the out-of-plane polarized band of honeycomb plasmonic lattices pretty resembles the π bands of electron in graphene and acts as a classical photonic analog. Meanwhile, the in-plane polarized modes have two polarization degrees of freedom, this vector nature of eigenfield on all lattice points leads to its four-band peculiarity as shown in Fig. 2(b). However, in graphene, the p xy orbital hybridizes with the s orbital, composing the σ band, and its anisotropic feature is not so conspicuous because of the large s-orbital component. Moreover, the σ band is fully filled and inert. Therefore, the similar "vector wave" excitation is not observed in graphene. Another property of the in-plane bands is that the upper two bands and the lower two bands joint with a quadratic form at Γ point. In addition, it is obvious that both the out-of-plane and the in-plane band structures are symmetric about ω 0 plane due to chiral symmetry under the nearest neighbour tight-binding approximation.

Band structures and edge states of ribbons
In this section, we investigate the band structures of plasmonic ribbons with four types of boundary conditions, i.e. zigzag, bearded zigzag, armchair, and bearded armchair. For a finite ribbon with width N , the Hamiltonians are 2N order and 4N order square matrices for out-of-plane and in-plane polarizations respectively. With accuracy to the nearest neighbor hopping, we can obtain the Hamiltonians for zigzag and armchair ribbons respectively where k is the Bloch wave vector alone the periodic direction, and G i should be substituted by −1/a 3 0 for out-of-plane case, and by ↔ G xy (t i ) for in-plane case. And the Hamiltonians of bearded zigzag and bearded armchair ribbons are almost the same as their non-bearded counterparts except one more hopping term from the beard atoms of each edge. Immediately, we can obtain the band structures by solving the eigenvalue equations of these Hamiltonian matrices. Fig. 3 shows the band spectra of out-of-plane and in-plane polarized modes with four different types of ribbon boundaries. The ribbon band structures have some generic properties for both out-of-plane and in-plane cases. To be more specific, most of the blue lines lie in the projected bulk band area (gray region) with some additional bands outside the projected band which attribute to edge states. The blue lines tend to fill up the whole gray region as the widths of the ribbons increase towards infinity. As the boundaries at two sides of a ribbon are kept to be the same, the edge bands always appear in pairs. In each pair, there are two nearly overlapped bands that are excited on different sides  Table I. By checking the eigenvectors, the additional flat bands in the whole Brillouin zone for bearded zigzag and bearded armchair arise from the localized resonances at the beards. Another notable effect is that both the out-of-plane and in-plane modes of armchair ribbons are metallic, namely the gap at ω 0 is strictly closed at the Dirac points, only when ribbon width N = 3M − 1 (M ∈ N), and this is same to the electronic graphene system. When long range couplings are taken into consideration, the chiral symmetry is broken, consequently, both bulk and edge bands are distorted, as illustrated by the band structures (white curves) of QSA Hamiltonian including all long range hopping in Fig. 4. Nevertheless, the topology of the bands does not change, namely, the projected bulk bands are still degenerate at Dirac points which are connected by edge bands in corresponding regions. Also, we have broken the inversion symmetry by altering the particle size ratio of two sublattices for these four ribbons, then the gapless surface states split off naturally (see details in Appendix B).
Moreover, since the retardation effect usually influences the band dispersions significantly in plasmonic lattices [52][53][54][55][56], we have also studies this influence by calculating the band structures with the intact dipole radiation given in Eq. 1. The band dispersions with the influence of retardation can be characterized by the loci of the extinction cross section peaks which can be mapped out through the effective polarizability under an external driven field (see Appendix C for detailed derivations) [55]. On the other hand, the dissipation of nanoparticles also affects the dynamic dispersion and the Q factors of the resonant peaks. Here, we choose the Drude damping term with a small value 0.01 eV to obtain a rather clear band dispersion and also relative high Q factors. As shown in Fig. 4, the band structures are in good agreement with QSA results especially far away from the light cone (here only half of the Brillouin Zone is shown due to the symmetry, and the truncation number of hopping is chosen to be 200). The lower out-of-plane bulk band and the two upper in-plane bulk bands are affected little by retardation. However, the upper out-of-plane bulk band and the two lower in-plane bulk bands are dragged down as closing to the light cone and significantly blurred due to coupling with free photons. These observations are also revealed by the multiple scattering based calculations [34]. For the mid-gap edge bands corresponding to the flat edge bands in the nearest neighbor quasi-static approximation, outside the light cone, they are also satisfactorily consistent with the QSA results plotted by the red curves in Fig. 4. Even inside the light cone, these edge bands are also agree with the results of QSA on the whole, apart from two exceptions. For the in-plane modes in bearded zigzag and bearded armchair ribbons (Fig. 4(f) and (h)), a pair of edge bands is dragged down and submerged by the bulk bands in each case, while, in the bearded armchair case, another pair of edge bands is scarcely disturbed by retardation and still overlaps with the QSA bands inside the light cone. It has to be noted that the light cone will move to the right and occupy larger space in the Brillouin zone, when the lattice constant are increased relative to the wavelength. And the occupations of the light cone in the ribbon Brillouin zone are different for ribbons with different boundaries due to their different lengths of periods. As lattice constant increasing, the width of these dipole bands will be compressed, as a results, the narrowed bandwidth compels the resonant peaks to overlap with each other (supposing the damping is unchanged), and eventually the band gaps as well as the edge states are submerged by this compression. Therefore, a much smaller damping is needed to preserve the Q factors of the bands, when increasing the lattice size.
In terms of the above results, we have showed that the essential characters of the HPLs band structures are accurately caught by the simplified nearest neighbor QSA method. According to the simplified QSA results, it is natural to ask why the flat edge states of in-plane polarized modes appear in different ranges than out-of-plane cases, and what determines the existence of these flat edge states. In the previous work [44], the existence of edge states with a particular k in graphene ribbons for arbitrary boundaries has already been related with the Zak phase of π bands. This principle of bulk-edge correspondence can be used directly in our system for the out-of-plane cases, but it is invalid to the in-plane cases which have a 4 dimensional bulk Hamiltonian due to the vector nature of the eigenfields. For this purpose, we will extend this principle into 4 dimensional vector wave systems with a general form of the Hamiltonian in the next section.

Relation between number of edge states, Zak phase and Winding number
For out-of-plane polarized modes, the bulk Hamiltonian possesses a 2D massless Diraclike form. According to the bulk-edge correspondence in this kind of systems [44,45], the number of pairs of flat edge states at the mid-gap for a fixed k , in a ribbon with an arbitrary boundary, can be predicted by the winding number W (p), as well as the Zak phase of either one of the bulk bands divided by π, Z/π. And the Zak phase can be calculated by where the integration is along the path perpendicular to the edge with a fixed k , Γ N is the reciprocal vector perpendicular to the boundary.
In the following, we will derive the existence conditions of the in-plane polarized flat edge states for zigzag boundary as an example, since other boundaries can be analyzed in a similar way. We start from the Bloch states, the eigenstates of Hamiltonian for the infinite system: ψ k = (ψ A k , ψ B k ) T . On each site of a unit cell, the dipole moment of this Bloch state is (P T which consists of a common phase factor and the relative amplitudes on two sublattices. Further, we assume a semi-infinite lattice with a unique edge possesses such edge states that hold the form of Bloch waves but decay exponentially away from the edge to interior. In this case, the wave vector should be extended to complex κ = (k , κ ⊥ ) with analytical continuation. The parallel component k is real and still a good quantum number, while the normal component κ ⊥ is complex corresponding to the exponential decay factor. Namely, the phase factor associated with κ ⊥ is which gives rise to a decay field P n (κ) ∝ ξ n as long as |ξ| < 1, where n is the number of layers counted from the edge (n = 0) to interior (see Appendix A for detailed derivation). For zigzag boundary, ξ = exp (iκ ⊥ 3a 0 /2). Since the edge states ψ κ maintain the Blochlike form, they should also obey the eigenequation of bulk hamiltonian H in κ with the complex wave vector κ. However, the eigenequation of the edge atoms need to be modified due to the lack of hoppings from sites outside the edge. For zigzag boundary, the eigenequations of the the edge states with Bloch-like form read B : where ↔ G i is short for ↔ G xy (t i ) (i = 1, 2, 3), and it is assumed that the edge is terminated at the A sublattice. The missing term in Eq. (28), in comparison with Eq. (27), is corresponds to the hoppings of the missing B atoms cut off by the edge. Here, we are only interested in the flat zero-eigenvalue edge states, i.e., λ(κ) = 0. Substituting λ(κ) = 0 into Eqs. (27)-(29), then we could obtain the sufficient conditions of the flat edge states whereÃ(k , ξ) ≡ A(κ) denotes a complex-valued function of ξ for a special k . For Therefore, every zero point of detÃ(k , ξ) † inside the unit circle |ξ| = 1 corresponds to a zero-eigenvalue edge state of the semi-infinite ribbons. For a ribbon with two edges (ribbon width should be large enough so that the coupling between two edges can be neglected), the edge states appear in pairs localized at two sides of the ribbon, i.e. a solution of Eq. (33) corresponds to an edge state on each edge respectively. As ξ travels along the path |ξ| = 1, the trajectory of detÃ(k , ξ) † is also a closed curve on the complex plane. According to Cauchy's argument principle, the winding number of the trajectory of detÃ(k , ξ) † equals to the difference between the number of zeros and poles of detÃ(k , ξ) † inside the unit circle. Moreover, detÃ(k , ξ) † is analytic and thus possesses no pole inside the closed path. Therefore, with a fixed k , the winding number of detÃ(k , ξ) † is, gives the number of pairs of flat edge states in a finite ribbon system. Since the circle |ξ| = 1 corresponds to k ⊥ ∈ [0, |Γ N |) with a fixed k which denotes a closed path, marked as Γ ⊥ , perpendicular to the edge in Brillouin zone (see Appendix A), the winding number of det A † also can be integrated with respect to k ⊥ alone this path: =number of pairs of flat edge states.
This relation reveals the correspondence between the existence of the flat edge states and the bulk Hamiltonian. Moreover, the existence of the edge states is also related to the Zak phase of each band. Since A is a symmetric matrix, it can be regarded as a symmetric bilinear mapping, or a tensor of (2, 0) type. Then, the representation of A in the basis {u 1 , u 2 } given by Eq. (20) can be written as Thus the bulk Hamiltonian becomes where Λ j = (Reλ j , Imλ j , 0) T = λ j (cos γ j , sin γ j , 0) T , and {Σ j , Σ 0 j } is a set of "generalized Pauli matrices": which satisfies the Pauli algebra and iΣ j /2 gives a 4-dimensional representation of su(2) Lie algebra. In consideration of Eq. (21), H in therefore has been separated into two parts which correspond to bands 1, 4 and bands 2, 3 respectively. By analogy to the out-of-plane case, we can introduce the winding number of Λ j (k) (or the complex functionλ j (k)), W (Λ j ), via tracing the trajectory of Λ j as k varying along the closed path with a fixed k , k ⊥ ∈ [0, |Γ N |) in the Brillouin zone. On the other hand, the Zak phase of each band reads Obviously, Z j equals to π times the winding number of Λ j . In addition, according to Eq. (36), we have where U = (u 1 , u 2 ) is a unitary matrix satisfying U = U † , and det U = −1.
Consequently, the relation between the bulk hamiltonian and the Zak phase (winding number) of each band is: Note that W (Λ j ) and Z j can either be positive or negative, and the equality of Eq. (42) is established for arbitrary boundary conditions. By combining Eqs. (35) and (42), then the main theorem of the bulk-edge correspondence for our system could be obtained.
Namely, the following quantities are equivalent: 1. The number of pairs of flat edge states at the mid-gap; 2. The winding number of det A † , W (det A † ); 3. The sum of the winding numbers W (Λ j ) of bands 1 and 2; 4. The sum of the Zak phases of bands 1 and 2 divided by π, (Z 1 + Z 2 )/π. All of the winding numbers and Zak phases mentioned above are integrated along the path Γ ⊥ which is perpendicular to the edge and goes through the first Brillouin zone.
This theorem declares how does the topological properties of the bulk bands determine the existence of the flat edge states and gives their relation quantitatively. Though we only provide the derivation of the theorem for zigzag ribbon in the text, it is not difficult to give a similar proof to other types of ribbons, such as the bearded zigzag ribbon. In the following section, we will study four common boundaries to show the validity of the theorem.

Application of the theorem to four common boundaries
To verify our conclusion, the Zak phases and winding numbers associated with the in-plane polarized modes are calculated independently for the four typical ribbon boundaries. And we will not show the reduplicative results of out-of-plane polarization as they are similar to graphene ribbon systems [44]. The calculated results of these four ribbons are shown in Fig. 5-6, the directed loops are the trajectories of det A † (k , k ⊥ ) as k ⊥ travels alone the path Γ ⊥ in the Brillouin zone with several fixedk shown on the top of each illustration (k = k a/2π andk = k a /2π in Fig. 5 and 6 correspondingly). The number of times that a loop travels counterclockwise around the origin denotes the winding number of det A † for the correspondingk . As illustrated in Fig. 5(a), two loops ofk = 1/3 and 2/3 pass through the origin, which correspond to the two phase transition points of the winding number. In Fig. 5(b), the sum of the Zak phases of bands 1, 2 is non-zero and equals to π only in the range of k ∈ [0, 1/3) ∪ (2/3, 1], and it agrees with the winding number of det A † . This indicates the existence of one pair of flat in-plane polarized edge bands in this range for zigzag ribbon, which is consistent with the band structure illustrated in Fig. 3(e). And for the bearded zigzag ribbon, the sum of Zak phases of bands 1, 2 equals to π in the range ofk ∈ [0, 1/3) ∪ (2/3, 1], and equals to 2π in the range ofk ∈ (1/3, 2/3), meanwhile, the winding number of det A † equals to 1 and 2 in these two ranges respectively. This indicates the existence of one pair of flat edge bands in the range ofk ∈ [0, 1/3)∪(2/3, 1], and two pairs ink ∈ (1/3, 2/3) for in-plane polarization.  We have to be noted that the value of Zak phase (winding number) of each band may have a difference of ±2π from the present result for the different choices of the bulk unit cell, but these differences of bands 1 and 2 always have opposite sign, as a result, the sum of the Zak phases of the two bands is independent of the unit choice. In conclusion, the existence of edge states predicted by winding number or Zak phase are perfectly consistent with the statistics given in Table I for the four common boundary conditions. And these results confirm the validity of Theorem 1.

Distribution of dipole moments
Apart from the band structures, we can also obtain the dipole moment distributions of eigenfields on the lattices. Here the zigzag ribbon of in-plane polarization is shown as an example. As we see in Fig. 7, for bulk modes, the dipole moment on each lattice point is comparable for all k , however, the flat edge modes exponentially decay into the bulk, and only one sublattice, which the outmost atoms belong to, is excited on each side of the ribbon. Besides, the decay rate |ξ| tends to 1, when k approaches Dirac point. In addition, both the bulk modes and the edge modes are elliptically polarized in general, and the eccentricity of polarization ellipse varies in respect to k . At the boundary of Brillouin zone (k = 0), the edge mode is linearly polarized with the direction perpendicular to the edge (Fig. 7(e)), whereas the edge mode tends to circular polarization when k goes toward the projected Dirac point (Fig. 7(h)). And this effect reflects the vector nature of in-plane edge modes, and thus is a distinctive property different from out-of-plane cases. Hitherto, with our method, both the characters of dipole moment distribution and topological properties of the in-plane polarized midgap edge modes are clearly depicted.

Conclusion
In this paper, we have investigated the band structures and edge states in honeycomb plasmonic lattices with different boundaries. We have proposed a definite correspondence between the Zak phases of bulk bands, the winding number associated with the bulk Hamiltonian, and the existence of flat edge states with in-plane polarization for arbitrary ribbon boundaries, and further verified it through four typical ribbons. We emphasis that the vector nature of the in-plane polarized modes in our system makes the existence conditions of edge states differ from the out-of-plane polarized modes, the π electronic excitations in graphene, and all of the systems possessing a 2 dimensional Dirac-like Hamiltonian. Therefore, the theorem of bulk-edge correspondence obtained in this paper is a novel extension to previous works.
Besides, we also give some necessary remarks to our results. Firstly, we have adopted several assumptions and approximations, i.e. we suppose a proper scale of the system parameters so that the CDM and QSA are applicable, and we also omit the dissipation in the system. Nevertheless, according to our comparison between the QSA and the retardation-affected results, the simplified QSA model catches the main characters of the real dipole bands accurately. Secondly, this bulk-edge correspondence is derived in tight-binding approximation with only nearest neighbor interaction. When long range interactions are taken into consideration, both the bulk bands and edge bands will be distorted. In spite of some changes in band structures, the topological properties do not change (see Appendix B). Thirdly, we only give a proof of the theorem for zigzag in the text. Although similar proofs are easily obtained case-by-case for some other boundary conditions, a general proof for arbitrary boundaries is still hard, therefore, Eq. (43) is still a hypothesis to some extent. Fourthly, our theorem is only applicable to predict the flat edge bands at the mid-gap, but is helpless to predict the existence of the bent edge bands between bands 1, 2 and bands 2, 3 of in-plane polarized modes in these ribbons. It is still an open question whether there is a topological invariant corresponding to these bent edge states. Finally, the model proposed here is not limited to photonic systems, and we believe it can also be applied to other 2D wave systems which have an arbitrary even order Hamiltonian with chiral symmetry. On the other hand, the honeycomb plasmonic lattice not only is a beautiful analogy of classical wave to graphene but also owns additional interesting properties which may be easier to be realized experimentally. We hope our work could inspire further theoretical and experimental researches of similar systems.

Appendix A: Basis vector for different boundaries
In this Appendix, we briefly review how to select a appropriate pair of bases of the honeycomb Bravais lattice for different ribbon boundaries [44]. We choose a pair of vectors (T, N) to be the set of bases corresponding to the ribbon with the periodic vector T, where T = ma 1 + na 2 and N = m a 1 + n a 2 . The bases should satisfy a 1 × a 2 = T × N which leads to mn − nm = 1, and (n , m ) can be thus determined (the choice of (n , m ) is not unique). For zigzag ribbon, T z , N z are just the original bases (a 1 , a 2 ). The reciprocal vector bases (Γ T , Γ N ) corresponding to (T, N) satisfy Then, we have where (b 1 , b 2 ) are the reciprocal bases corresponding to (a 1 , a 2 ). Γ T and Γ N form a parallelogram that can be regarded as an alternative selection of Brillouin zone in the reciprocal space. Nevertheless, Γ N and an arbitrary vector Γ can also forms a bulk Brillouin zone as long as Γ × Γ N = Γ T × Γ N . We demand Γ is parallel to the edge, and obtain Γ = (Γ T ·t)t. And |Γ | = 2π/|T| is precisely the length of the Brillouin zone of the ribbon. The relations of these bases for zigzag ribbon is shown in Fig. 8. For an arbitrary wave vector k, it can be decomposed with respect to the directions of the two reciprocal bases, k = k Γ T + k Γ N , or equivalently, decomposed into the two components parallel and perpendicular to the edge of the ribbon, k = k t + k ⊥n , wheret andn are the unit vectors parallel and perpendicular to the edge. The relations between the two kinds of components are In all eigenequations above, the wave vector appears only in the fixed term k · R of the phase, where R = l 1 T+l 2 N is an arbitrary lattice vector. This term can be decomposed as follows k · R = k [l 1 |T| + l 2 (N ·t)] + l 2 k ⊥ 2π |Γ N | .
So in Wannier picture, the component of Bloch eigenfunction on lattice point R can be written as For the edge states in ribbon systems, k ⊥ should be replaced by the complex number κ ⊥ , and the first phase term in Eq. (48) corresponds to the decay of the fields from the edge to the interior, l 2 is the layer number at R. Thus, the decay rate of the edge state is ξ = exp(iκ ⊥ 2π |Γ N | ). For zigzag ribbon, Γ N = b 2 , so we have |Γ N | = 4π/(3a 0 ), and the decay rate ξ = exp(iκ ⊥ 3a 0 /2). The unit circle |ξ| = 1 on the complex plane corresponds to the path k ⊥ ∈ [0, |Γ N |) with a fixed k . Since the periodicity of the lattice makes the Brillouin zone composed by (Γ , Γ N ) be a torus, the path, k ⊥ ∈ [0, |Γ N |) with a fixed k , forms a closed path, marked as Γ ⊥ , on the torus, as shown in Fig. 8(c). For the in-plane polarization, the band structure also loses its mirror symmetry with respect to the mid-gap plane. Nevertheless, the topological properties of the bands do not change, concretely, the upper and lower bulk bands still degenerate at the Dirac points, despite their small shifts, and the mid-gap edge bands still appear in the ranges predicted by our theorem and connect the projected Dirac points, so the present result further illustrates our method works.
As the ribbon is made up of two sublattices, we can break the inversion symmetry by changing the size ratio of A type particles to B type particles. Here we define a ratio factor s = r 3 A /r 3 B , where r A and r B are the radii of A and B types of particles respectively. Apparently, s can be assigned arbitrary value, but we have to care about that the ratio must be reasonable so that the coupled dipole approximations are appropriate. In this sense, the Hamiltonian becomes Then we can map out the band structures in the similar way as before, the result is shown in Fig. 9. It is known that if either inversion symmetry or time-reversal symmetry is broken, the degeneracy at the Dirac points will immediately be lifted in 2D band structures according to the Von Neumann-Wigner theorem [12,13]. In contrast to graphene, the degeneracy will not open if we only consider nearest interactions, unless the second hopping influence is also included into the Hamiltonian, because there are no on-site energy terms in our eigen equation. For ribbon band structures, the original degenerate pairs of edge bands will completely separate into two gaped bands contacting the upper and lower bulk bands respectively, however, their topology is trivial since the Chern numbers of these bulk bands are zero when merely breaking the inversion symmetry.
As a result,M(k, ω) −1 serves as the effective polarizability ↔ α eff (k, ω) of the ribbon lattice in response to the driven field. And each eigenvalue α i eff (k, ω) of ↔ α eff (k, ω) indicates one forced oscillation mode of the lattice at the driven frequency ω and Bloch vector k. In terms of optical theorem, the total extinction cross section of all modes takes the following form Since the loci of the total extinction cross section peaks just correspond to the collective resonances of the plasmonic lattice, the band dispersion including the whole dipole interactions and retardation can accordingly be exhibited by the loci of the peaks of C ext . Besides,M(k, ω) is also composed by the direct sum of the in-plane polarization subspace and the out-of-plane polarization subspace, therefore, the band structures of in-plane modes and out-of-plane modes still can be investigated independently by this method.