We study the topological properties of a Kagome plasmonic metasurface, modeled with a coupled dipole method that naturally includes retarded long range interactions. We demonstrate that the system supports an obstructed atomic limit phase through the calculation of Wilson loops. Then, we characterize the hierarchy of topological boundary modes hosted by the subwavelength array of plasmonic nanoparticles: both one-dimensional edge modes and zero-dimensional corner modes. We determine the properties of these modes, which robustly confine light at subwavelength scales, calculate the local density of photonic states at edge and corner modes frequencies, and demonstrate the selective excitation of delocalized corner modes in a topological cavity, through nonzero orbital angular momentum beam excitation.

Two-dimensional (2D) higher-order topological insulators (HOTIs) host protected zero-dimensional (0D) localized corner states in addition to conventional one-dimensional (1D) edge states. Following their theoretical proposal in condensed matter^{1} and subsequent experimental demonstration in electrical circuits,^{2} many of the concepts have been transferred to classical wave systems including photonics.^{3} Early proposals in photonics focused on realization of corner modes in evanescently coupled waveguides,^{4,5} which map directly to tight-binding models with nearest neighbor interactions due to the short-range nature of the coupling. More recent literature has turned to photonic crystals, where the interactions are necessarily long-range.^{6} A number of experimental investigations have shown how corner modes can arise in photonic crystals with different crystalline symmetries, such as *C*_{3},^{7} *C*_{4},^{8} and *C*_{6}.^{9}

In particular, the *C*_{3}-symmetric “breathing Kagome” lattice and its hierarchy of topological modes have attracted considerable attention.^{10–14} It hosts topological valley edge states due to the localization of the Berry curvature with opposite sign at opposite valleys.^{15,16} Notably, it also has a topologically distinct bulk as a result of displaced Wannier centers,^{10} which result in corner modes in finite lattices. These corner modes have been probed experimentally in coupled waveguides^{4} and recently in a photonic crystal.^{17} There it was shown that by adding *ad hoc* long-range interactions to a next-nearest neighbor tight-binding model, an additional set of trivial corner modes appears.

In this Letter, we investigate higher-order topology in a subwavelength dipolar array: the breathing Kagome plasmonic metasurface. To study the hierarchy of subwavelength topological modes emerging in this system and the effects of slowly decaying photonic interactions, we consider a coupled dipole model that naturally incorporates full electromagnetic interactions. We start by calculating Wilson loops for the system in order to characterize the bulk topology. To appropriately capture the physical behavior of the metasurface, we calculate the optical response including retarded interactions. Next, we reveal the spin angular momentum properties of the valley topological edge states and their directional excitation. Moving to a finite system, we show how topological corner states emerge and explore their robustness. Finally, we study methods of exciting these modes, first by probing the local density of states (LDOS) of the corner and edge modes and then by showing the selective excitation of corner modes.

We model the metasurface, a 2D lattice of metallic nanoparticles (NPs), with the coupled dipole method (CDM). We consider silver nanorods modeled with the analytical spheroid polarizability,^{18} which includes depolarization and radiative effects to appropriately describe the optical response. The dielectric permittivity is given by a Drude model, $\u03f5(\omega )=\u03f5\u221e\u2212\omega p2/(\omega +i\gamma )$, with $\omega p=8.9$ eV, *γ* = 38 meV, and $\u03f5\u221e=5$.^{19} We only consider the out-of-plane modes here and note that the in-plane modes, which occur at a different frequency due to the spheroidal NPs, have been studied previously in a dipolar breathing Kagome lattice,^{20} although only nearest and next-nearest neighbor interaction terms were included. Importantly, interactions between NPs are given by the full dyadic Green's function, which incorporates retarded interactions that are important in describing the physical behavior of a plasmonic metasurface as can be seen by comparing to the quasistatic (QS) approximation.^{21} The CDM is also a suitable model for any system of subwavelength dipolar scatterers, such as silicon carbide NPs,^{20} quarter-wavelength resonators,^{22} or cold atoms.^{23} More details can be found in the supplementary material.

We begin by calculating the spectrum in the QS approximation; however, we go beyond a nearest-neighbor (NN) approximation and include all neighbors in the lattice sum. The Kagome unit cell is shown in Fig. 1(a). In Fig. 1(b), we plot the band structure of the unperturbed Kagome lattice (*δ* = 0) in black. The cusp at Γ is unphysical, as we shall show, but is typical of including all neighbors in the QS lattice sum.^{21,24} At $K(K\u2032)$, we see a Dirac cone between the two highest energy bands, associated with the $C3v$ symmetry of the lattice.^{25} To gap the Dirac cone, we break inversion symmetry, *σ _{v}*, by perturbing the trimer in the unit cell inwards or outwards by a scale factor

*δ*, such that $R=(1+\delta )R0$. We will refer to these as the “contracted” and “expanded” lattices, respectively. We plot the bands for these lattices with $\delta =\xb10.1$ as dashed blue lines in Fig. 1(b), noting that both perturbed lattices have identical band structures.

^{12,15,16}

Even though these lattices share band structures, they can be distinguished by their Wilson loops.^{26} Eigenvalues of the Wilson loop operator describe the evolution of Wannier centers along closed paths in the Brillouin zone and provide information about the band's topological characters.^{27} We calculate Wilson loops in the QS approximation for the two connected bands below the topological bandgap, for the two lattice configurations. While retardation can have significant effects on topological properties, we emphasize here that we are in a subwavelength parameter regime.^{28,29} For the contracted lattice, Fig. 1(c), the Wilson loops are displaced about *W* = 0. However, for the expanded lattice, Fig. 1(d), they are displaced about *W* = 0 as well as away from zero. Lattices with this behavior are in an obstructed atomic limit phase.^{27} We note that unlike previous works, our calculations of Wilson loops go beyond a NN approximation.^{13}

To capture the optical response of this fully electrodynamical system, we use the retarded Green's function and the radiative polarizability to calculate the spectral function $\sigma sf$. We plot this for the unperturbed, Fig. 1(e), and perturbed lattices, Fig. 1(f). First, $\sigma sf$ shows the broadening and slight redshifting of the modes due to the radiative effects. Starting at Γ, the highest energy mode has a monopolar character and the two lowest energy modes are degenerate dipolar modes. At the light line, a strong polariton-like splitting occurs in the highest energy monopolar band due to the coupling to free photons, similar to effects in 1D^{30} and other 2D plasmonic lattices.^{31–33} Notably, the fully electrodynamical interactions remove the cusp at Γ observed in the QS approximation.

The system also has a band inversion at $K/K\u2032$ between the contracted and expanded phases.^{16} We plot the out-of-plane electric field *E _{z}* for the two bands involved in the Dirac cone gapping, for contracted, Fig. 1(g), and expanded, Fig. 1(h), phases. In the contracted phase, the highest energy band is monopolar and the next band is dipolar, while the ordering is flipped in the expanded phase. Finally, we note that the overall ordering here is opposite to the dielectric system,

^{15}and the flatband is the lower frequency band. This is due to the negative and dispersive permittivity of metallic NPs and is not a result of the dipolar nature of the model, as claimed elsewhere.

^{20}

We now move onto investigating the topological valley edge states, which arise at the interface between contracted and expanded regions. These are guaranteed to occur via the bulk boundary correspondence, provided the perturbation between the two regions *δ* is small enough and preserves the localization of the Berry curvature at the $K/K\u2032$ valleys.^{34,35} Propagation along these interfaces is protected against disorder, which does not mix the valleys due to intervalley separation.^{36} To model these edge states, we set up a ribbon with interfaces between the two phases as shown in Fig. 2(a), with $\delta =\xb10.1$—see the supplementary material. The ribbon spectral function, Fig. 2(b), shows how within the bandgap there are two edge states; the red band corresponds to the expanded over contracted (E/C) interface and the blue band to the contracted over expanded (C/E) one.

To calculate the edge mode profiles, we linearize the Green's function, see the supplementary material. In Fig. 2(c), we plot the out-of-plane electric field *E _{z}* for the E/C interface, which shows how the edge state is strongly confined to the interface. Then, we show the in-plane time averaged Poynting vector $\u27e8S\u27e9=12Re(E*\xd7H)$ (d), out-of-plane electric field phase $\varphi (Ez)$ (e), and spin angular momentum $T=12Im(E*\xd7E+H*\xd7H)$

^{37}(f), which encodes the degree of elliptical polarization of in-plane magnetic fields and is out of the plane $T=Tzz\u0302$ for

*z*= 0.

From the Poynting vector $\u27e8S\u27e9$, we see how the energy flow is confined to the interface, with direction given by the group velocity $vg<0$ of the mode at $kx>0$, while $\varphi (Ez)$ and *T _{z}* determine the propagation direction of modes upon excitation.

^{32}For example, beams with nonzero orbital angular momentum will excite a directional mode if the phase vortex of the beam matches the phase vortex in $\varphi (Ez)$ of the edge eigenmode: this occurs at the center of an expanded unit cell next to the interface and has been shown in a similar photonic crystal.

^{38,39}Additionally, if a circularly polarized in-plane magnetic dipole source is placed at the center of an expanded unit cell, it will excite a mode along a direction determined by the source polarization: right(left)-propagating for right(left) circular polarization. However, at the edge of the unit cell, the local handedness of polarization changes sign and sources will couple to the mode propagating in the opposite direction. We show examples of directional propagation in the supplementary material.

Before studying the emergence of corner modes, let us briefly review the concept of chiral symmetry. In the limit of short range NN interactions, lattices with an even number of elements in the unit cell are chirally symmetric.^{40} In condensed-matter systems, chiral symmetry pins corner modes to a specific frequency in HOTIs.^{41} This imposes a subtle distinction between HOTIs and photonic systems with long-range couplings, where chiral symmetry is necessarily broken due to interactions between sites in the same sublattice. By continuously reducing the interaction range, we have previously shown that retarded systems are deformable to ones with chiral symmetry,^{6} meaning they possess an approximate chiral symmetry. We refer to corner modes in these systems as “higher-order topological modes” (HOTMs).

The breathing Kagome lattice has an odd number of elements in the unit cell, meaning it is incompatible with chiral symmetry. On the other hand, a “generalized chiral symmetry” has been proposed for this system.^{12,13} This acts similarly to chiral symmetry in that, in the limit of short-range interactions, it pins the corner modes at a mid-gap energy (or the localized surface plasmon frequency, $\omega =\omega lsp=2.845$ eV for these plasmonic NPs) and forces the excitation on a single sublattice in real space. With long-range interactions, the breathing Kagome lattice still has an approximate generalized chiral symmetry^{17} as a result of critical interactions between the same sublattice being negligible compared to interactions between different sublattices. As such, long-range interactions only slightly shift the frequency of the corner modes away from mid-gap energy, $\omega =\omega lsp$.

To investigate corner modes, we setup a periodic supercell containing a region in the expanded phase (the “topological particle”) surrounded by the contracted phase, Fig. 3(a). The finite nature of the topological particle means there will be a discrete set of modes, which arise on the interface.^{42}

Since HOTMs are reliant on the displaced Wannier centers,^{10} rather than the localization of the Berry curvature as the edge modes are, we are free to choose the perturbation $\delta >0$. This controls the size of the bandgap as well as the localization of modes to the corners, see the supplementary material. Again, we linearize the Green's function to calculate the spectrum of the topological particle. Figure 3(b) shows the spectrum for $\delta =0.15$ with corner modes (red) and edge modes (cyan) appearing within gapped bulk modes (gray). In Fig. 3(c), we plot the real part of the out-of-plane dipole moments $Re(pz)$ of the three corner eigenmodes. For a system with short-range interactions, the eigenmodes are triply degenerate,^{4} whereas here the long-range interactions split the degeneracy of the eigenmodes into $(1+2)$ as dictated by symmetry of the topological particle. Mode **s** is a monopolar-like mode, with the NPs at the three corners oscillating in phase, while modes $p1$ and $p2$ are dipolar-like modes. Although the effects of long-range interactions on corner modes in the breathing Kagome lattice have been studied previously in photonic crystals,^{17} here we start with a model that includes fully retarded photonic interactions such that there is coupling between all sites in the lattice.

The HOTMs are protected by the approximate generalized chiral symmetry and crystalline symmetry of the particle. We demonstrate the robustness by plotting the mode spectrum for increasing random positional disorder in Fig. 4(a). To apply the disorder, we shift each NP randomly in the *xy*-plane up to a percentage $\kappa =|\kappa |$ of the lattice constant *a*_{0}, $r\u2192r+\kappa a0$. Since positional disorder breaks *C*_{3} and mirror symmetries across the topological particle, the $(1+2)$-degeneracy is broken. Despite this, the HOTMs are robust up to $\kappa \u22484%$ where they eventually merge with the bulk. In Fig. 4(b), we plot a HOTM for $\kappa =1%$, indicated by the vertical black dotted line in panel (a). The slowly decaying nature of retarded interactions allows the excitation of modes delocalized over the cavity even if the crystalline symmetry is broken, provided the size of the cavity is not too large. In comparison, in the QS approximation (not shown), the modes become localized to single corners. This highlights the limitations of tight-binding models when studying nanophotonic systems,^{4} and explains how cavity delocalized corner modes can be observed in experimental setups,^{7} despite disorder due to fabrication errors. Finally, we note that another indicator of higher-order topology is a fractional corner anomaly.^{41} The breathing Kagome plasmonic metasurface presents this behavior, as we show in the supplementary material, which justifies referring to localized corner states as examples of higher-order topology.

Following the understanding of how topological boundary modes emerge, we now investigate methods of exciting them. We consider the coupling of a near-field source to the edge and corner modes of the topological particle. We characterize this through calculations of the partial LDOS for a vertical source, $\rho z(\omega )$, which captures the available photonic states which the source can couple to.^{43,44} From this, we calculate the Purcell factor,

with the decay rate enhancement $\Gamma /\Gamma 0$ and LDOS enhancement $\rho z/\rho 0$ with respect to vacuum. We consider a vertical dipole source, with dipole moment $\mu =(0,0,1)T$, and calculate the electric field scattered by the metasurface at the source position, $Es(r0)$. We take *γ* = 1 meV to make the modes more visible. For higher losses, the modes will become broader but will still be visible, provided they are well separated from the bulk modes.

We plot *F _{P}* for a frequency within the edge states band in Fig. 5(a). As expected, the

*F*is the highest at the edges of the topological particle and we note that peaks in

_{P}*F*are not restricted to a single sublattice. On the other hand, in Fig. 5(b) we show how when exciting at the HOTM frequency

_{P}*F*peaks at the corners and on only one sublattice due to the generalized chiral symmetry respected by the corner modes. As we showed earlier, another consequence of this symmetry is the HOTM frequency being close to $\omega lsp$. Therefore, the increase in

_{P}*F*over the NP directly at the corner is a combination of coupling both to a single NP localized surface plasmon (LSP) and to a HOTM. Critically, the coupling to the HOTM allows for the excitation of delocalized modes across the topological particle: A point source at one corner couples to the other two corners through the HOTM (see the supplementary material). This phenomenon could be exploited to generate topological long range effective interactions between quantum emitters.

_{P}^{45}

We will now show how to selectively excite HOTMs from the far-field using Laguerre–Gaussian (*LG _{pl}*) beams. These are higher-order laser modes, which can have nonzero orbital angular momentum due to optical vortices in the phase of the electric field.

^{43}Laguerre–Gaussian beams are parameterized by an azimuthal index

*l*, which controls the optical vortices, and a radial index

*p*, which controls the nodes of the beam. In the following, we will use

*LG*

_{00}(corresponding to a conventional Gaussian beam) and

*LG*

_{01}. The electric field intensity and phase of these are shown in the insets in Figs. 5(c) and 5(e). We excite with a beam with waist $w0=600$ nm such that the excitation extends across the topological particle and vary the frequency over the HOTM frequency. Additionally, we excite the system slightly off-normal incidence, $\theta =1\xb0$, in order to couple to out-of-plane modes (with the electric field polarized in

*xz*).

To characterize the selective excitation, we calculate the overlap *η* of the excited dipole moments $qexc$ across the system with the monopolar and dipolar eigenmodes calculated earlier and normalized to the norm of excited dipole moments,

In Fig. 5(c), we see that the *LG*_{00} beam predominantly excites a monopole (green line), compared to the dipole (purple line). This is due to the beam having orbital angular momentum (OAM) *l* = 0 and exciting each corner of the particle in phase. Next, in Fig. 5(e) we show that the *LG*_{01} beam preferentially excites the dipolar mode, which is due to the OAM of the beam *l* = 1. The phase vortex of the beam at the origin of the system results in a rotating phase pattern, which is no longer azimuthally symmetric across the topological particle. To confirm the monopolar and dipolar selective excitation, we excite the system at fixed frequency and plot the dipole moments in real space, in Fig. 5(d) for *LG*_{00} and Fig. 5(f) for *LG*_{01}. We note that any beam with a phase pattern that breaks the *C*_{3} symmetry of the system will allow for dipolar modes to be excited. This means, for example, at large angles of incidence, one can excite a dipolar mode with an *LG*_{00} beam.

In summary, we have investigated crystalline-symmetry-dependent higher-order topological corner modes in a plasmonic metasurface of metal nanoparticles arranged in a Kagome lattice. The system hosts a hierarchy of topological modes, which we characterize appropriately with long-range interactions. With Wilson loops, we identify an obstructed atomic limit bulk phase. We additionally calculate the optical response of the bulk with retardation to appropriately capture the physical response of the system. The interfaces between topologically distinct regions of the lattice host topological valley edge modes and we reveal the spin angular momentum mechanism for exciting directional modes. We complete the hierarchy with zero-dimensional localized corner modes, which are robust to random spatial perturbations due to their topological origin from the bulk. Motivated by typical experimental setups, we probe the modes of the system in both the near and far-field. We calculate the Purcell factor for a point source above the metasurface to demonstrate optimal positions for exciting edge and corner modes. Finally, we show how HOTMs can be selectively excited given the orbital angular momentum of a beam. Both edge and corner modes in the plasmonic metasurface could be exploited for mediating interactions between quantum emitters, with the benefits of topological protection.^{45}

See the supplementary material for (1) details on the coupled dipole method, including form of the polarizability, and Green's function linearization; (2) Berry curvature and directional propagation of edge states; and (3) corner mode localization and fractional corner anomaly.

M.P. and P.A.H. acknowledge funding from the Leverhulme Trust. P.A.H. further acknowledges funding from Fundação para a Ciência e a Tecnologia and Instituto de Telecomunicações under Project No. UIDB/50008/2020 and the CEEC Individual program with Reference No. CEECIND/03866/2017. D.B. and A.G.E. acknowledge funding from the Spanish Ministerio de Ciencia, Innovación y Universidades through Projects Nos. FIS2017-82804-P and PID2019-109905GA-C22, respectively. D.B. acknowledges funding from the Transnational Common Laboratory. A.G.E. received funding from the Gipuzkoako Foru Aldundia No. OF23/2019 (ES) project and from Eusko Jaurlaritza Grant Nos. IT1164-19 and KK-2019/00101. Quantum–ChemPhys.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.