Near- and Far-Field Excitation of Topological Plasmonic Metasurfaces

The breathing honeycomb lattice hosts a topologically non-trivial bulk phase due to the crystalline-symmetry of the system. Pseudospin-dependent edge states which emerge at the interface between trivial and non-trivial regions can be used for directional propagation of energy. Using the plasmonic metasurface as an example system, we probe these states in the near and far-field using a semi-analytical model. We give the conditions under which directionality is observed and show that it is source position dependent. By probing with circularly-polarized magnetic dipoles out of the plane, we first characterize modes along the interface in terms of the enhancement of source emission due to the metasurface. We then excite from the far-field with non-zero orbital angular momentum beams. The position dependent directionality holds true for all classical wave systems with a breathing honeycomb lattice. Our results show that a metasurface in combination with a chiral two-dimensional material could be used to guide light effectively on the nanoscale.

Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal (Dated: August 28, 2020) The breathing honeycomb lattice hosts a topologically non-trivial bulk phase due to the crystalline-symmetry of the system. Pseudospin-dependent edge states which emerge at the interface between trivial and non-trivial regions can be used for directional propagation of energy. Using the plasmonic metasurface as an example system, we probe these states in the near and far-field using a semi-analytical model. We give the conditions under which directionality is observed and show that it is source position dependent. By probing with circularly-polarised magnetic dipoles out of the plane, we first characterize modes along the interface in terms of the enhancement of source emission due to the metasurface. We then excite from the far-field with non-zero orbital angular momentum beams. The position dependent directionality holds true for all classical wave systems with a breathing honeycomb lattice. Our results show that a metasurface in combination with a chiral two-dimensional material could be used to guide light effectively on the nanoscale.

I. INTRODUCTION
Topological nanophotonics offers a path towards efficient and robust control over light on the nanoscale [1]. Concepts borrowed from topological insulators, materials which host protected surface states for electrons, can also be applied to photonic systems. Following theoretical proposals [2,3], protected photonic modes reliant on a explicit time-reversal symmetry breaking component have been demonstrated experimentally [4]. However, these require strong magnetic fields or complex materials with a large magneto-optical response, which makes such systems difficult to miniaturise. More recently these effects have been proposed using graphene [5,6], which naturally has a large magneto-optical response. This is limited to the infrared regime however, which restricts potential applications. Methods of achieving topological protection through the crystalline symmetry of a system whilst preserving time-reversal symmetry are therefore appealing, since they do not require complex setups and are not restricted to a specific frequency regime. These effects fall into two main categories: valley effects, which rely on extrema in the band structures of materials [7][8][9], and pseudospin-dependent effects, which rely on the spin angular momentum texture of the electromagnetic fields [10,11].
First proposed in Ref. [12], the pseudospin-dependent effect relies on a triangular lattice of hexamers to produce states reminiscent of the quantum spin Hall effect (QSHE) in topological insulators. This breathing honeycomb lattice has two gapped phases: the shrunken phase, where the hexamers are perturbed inwards and the expanded phase where they are perturbed outwards. Despite both having a trivial Z 2 index, the phases of the breathing honeycomb are topologically distinct: while the shrunken phase is a trivial insulator, the expanded one is an instance of an 'obstructed atomic limit' phase [13,14], and edge states will appear between regions in either phase. A direct analogy of the QSHE would produce purely unidirectional edge states, in the absence of spin-mixing impurities. However, as we showed in a previous work [15], the edge mode directionality for near-field sources is source-position dependent and is determined by the spin angular momentum of the modes: it is the local handedness of the elliptical field polarisation which determines the propagation direction of the edge modes, rather than an absolute pseudospin quantity as in the QSHE [15]. This result is true for any bosonic breathing honeycomb lattice [16] since it is rooted in time-reversal symmetry and the absence of Kramer's degeneracy for bosons as op-posed to fermions. There have been a range of experimental investigations on the breathing honeycomb photonic crystal [17][18][19][20], including the observation of edge modes in the visible regime [21,22]. Despite this, there is no comprehensive theoretical study of the methods for exciting pseudospin edge modes and in particular, the necessary conditions for exciting unidirectional modes.
In this article, we consider a plasmonic metasurface consisting of a two-dimensional array of metallic nanoparticles with a breathing honeycomb lattice. We first characterise the optical response of the bulk modes and then investigate the propagation properties of edge states which emerge at the interface between trivial and and non-trivial regions. We then extend the understanding of excitation by near-field sources by probing edge states with sources out of the plane. Finally, we show that the propagation of edge states with farfield beams is determined by the electric field phase of the edge eigenmodes . Whilst the results we present are in the plasmonic metasurface, we emphasise that the properties and behaviour of the edge states are valid for any classical wave system.

II. METHODS
We model the system of subwavelength, metallic nanoparticles (NPs) using the coupled dipole method. When the nearest neighbour spacing R and NP radius r satisfy R > 3r, each NP can be treated as a point dipole [23,24]. To model nanorods, we use spheroidal NPs with radius r = 10 nm and height h = 40 nm. For NPs of this size and shape it is necessary to include depolarization and radiative effects and so we use the Meier Wokaun long wavelength approximation (MWLWA) to describe the the dipolar optical response of an individual NP [25,26]. The MWLWA polarisability α(ω) for spheroids is, with the static polarisability α s (ω), k = m ω/c is the wave number and the environment is a homogeneous vacuum with m = 1.
V is the NP volume and D and L are dynamic and static geometrical factors, and l E is the spheroid half-axis; D = 1, L = 1/3 and l E = r for a sphere. The dielectric function (ω) is given by the Drude model, We use silver NPs, with ∞ = 5, ω p = 8.9 eV and γ = 1/17 fs ≈ 38 meV [27]. For a system of multiple NPs, we can write a self-consistent coupled dipole equation which describes the dipole moment of a NP due to neighbouring NPs as well as an incident electric field E inc , The interactions between dipoles are characterised by the dyadic Green's function [28], where d ij is the separation between dipole i and j, d = |d ij | and n = d ij /d. (The· represents a dyadic operator.) We only take the zz-component of the Green's function, which corresponds to interactions between dipole moments perpendicular to the separation between the NPs (i.e. out-of-plane). This is a valid approximation since we use spheroidal NPs, causing the in-plane modes to become well separated in frequency from the out-ofplane modes [29]. For a periodic array of NPs in a plasmonic metasurface we can apply Bloch's theorem and write the following system of equations, The vector p contains the out-of-plane dipole moments p z of all NPs in the unit cell. The interaction matrixĤ(k , ω) has elements, where p, q are unit cell indices. The summations run over the lattice sites R = na 1 + ma 2 , with lattice vectors a 1 and a 2 . The interaction matrix has dimension N × N where N is the number of NPs in the unit cell. Since the lattice sums are slowly converging, we use Ewald's method to calculate them [30,31]. When modelling plasmonic metasurfaces, it is necessary to include the long range, retarded interactions in the Green's function in Equation 5. Whilst the quasistatic approximation can be used to model very subwavelength plasmonic chains and arrays [32,33], it fails to accurately capture the behaviour of modes at the light line and the radiative broadening and redshifting of modes, which becomes apparent in NPs with r > 10 nm. Additionally, retarded interactions can affect the topological properties of some plasmonic systems [34,35].
By solving Equation 6, we can calculate the extinction cross section σ ext from the dipole moments and incident electric field using the optical theorem, The incident field satisfies Maxwell's equations, |k |E + k z E z = 0 (we assume E = 1 and harmonic time dependence e −iωt ). When the field is propagating, above the light line, k z = k 2 − |k | 2 and when it is evanescent, below the light line, k z = i |k | 2 − k 2 . We can also solve Equation 6 as an eigenvalue problem and calculate the spectral function, by letting E inc = 0, which allows us to probe modes whether they are bright or dark. The spectral function σ sf is, with the effective polarisability α eff = i 1/λ (i) for eigenvalues λ (i) [33]. Finally, the coupled dipole method can also be applied to finite systems to model electromagnetic scattering [36].

III. OPTICAL RESPONSE OF BULK AND EDGE STATES
We begin by characterising the optical response of the periodic system, in order to later excite the finite system at the correct frequencies. The breathing honeycomb lattice setup is shown in Figure 1(a) with a unit cell with six NPs [12]. In a honeycomb lattice, R = R 0 = a 0 /3 where R is the nearest neighbour spacing and a 0 is the lattice constant. We let R 0 = 40 nm, so the lattice constant a 0 = 120 nm. The larger unit cell (compared to the rhombic Wigner-Seitz unit cell) folds the Brillouin zone and results in a double Dirac cone at Γ (as shown in Appendix A, Figure 5). We note that the Dirac cone lies below the localised surface plasmon frequency ω lsp due to chiral symmetry breaking, long range interactions [15].
By perturbing the nearest neighbour separation with a scale factor s, such that R = sR 0 , a band gap is opened at Γ. In Figure 1(c, d) we plot the extinction cross section σ ext across the Brillouin zone for the shrunken (R = 0.9R 0 ) and expanded (R = 1.065R 0 ) lattices. In the shrunken lattice, for k close to Γ, the NPs in the unit cell hybridise to form hexapoles, quadrupoles, dipoles and monopoles (from lowest energy to highest energy). The monopolar mode has dipole moments in phase, forming a bonding mode in the out-of-plane dipole moments across the whole lattice. As a result it couples very strongly with the light line and exhibits a polariton-like splitting. We note that the ordering in this plasmonic metasurface is opposite to the photonic crystal due to the metallic nature of the NPs [15].
A band inversion occurs around Γ between the shrunken and expanded phases, causing the dipolar and quadrupolar bands to flip. This is similar to the band inversion process in the QSHE. Although the bands become much broader above the light line, there is still a clear signature between shrunken and expanded phases which is observable in the far-field, as has been shown previously [15,37]. We plot σ ext for a fixed wavevector k in Figure 1(e), showing how the dipolar peak shifts between the two phases. The broad peak at the top corresponding to the monopolar mode does not shift. Unlike the QSHE, neither the shrunken or expanded lattices have a non-trivial Z 2 -invariant. However, they are still topologically distinct and one way to distinguish them is by their Wilson loops [38]. The shrunken lattice is topologically trivial whilst the expanded lattice is in an photonic 'obstructed atomic limit' phase [13]; it is the C 6 symmetry which provides the topological character.
Despite having trivial Z 2 -invariants, when regions in the shrunken phase are placed next to regions in the expanded phase, edge states will appear the band gap due to the topological origin of the band inversion [15]. These edge states are not topologically protected but do have a pseudospin character, which allows directional modes to be excited through a similar mechanism to the chiral light matter interactions in photonic crystals [39]. We model the edge states of the system by setting up a ribbon with an interface between the two phases, as shown in Figure 1(b). The ribbon supercell is N = 20 unit cells along the a 2 direction, with 10 unit cells in the expanded and shrunken phases, respectively. In Figure 1(f), we plot the spectral function σ sf for the ribbon for k = 0 to π/a 0 . (The spectral function from k = 0 to −π/a 0 is identical.) Edge states appear in the band gap from ω ≈ 2.74 eV to 2.81 eV and a small minigap appears between ω ≈ 2.77 eV to 2.78 eV. The mini-gap appears due to the interface breaking C 6 symmetry, which is the symmetry that protects the topological phase. We note that we use an 'armchair' interface here but it is also possible to define other terminations for this lattice [40], including the 'zig-zag' interface [11]. In the latter case, the C 6 symmetry breaking across the interface is slightly smaller and results in a smaller mini-gap, but otherwise the behaviour of the edge states is qualitatively similar. This is evident in the hybridisation of NPs close to the interface: the lower band has an anti-bonding character whilst the upper band has a bonding character (we show these in Appendix B, Figure 6). Due to doubly-periodic supercell, we see multiple narrow bands bending down and crossing the band gap in Figure 1(f). This is due to diffraction orders above the light line (white dotted line) [41] and does not affect the investigation of edge states which follows.

IV. CIRCULARLY-POLARIZED POINT SOURCES
We begin by extending the investigation of our previous work [15]  Experimentally, it is challenging to realise circularly-polarized magnetic dipoles at optical frequencies and in nanoscale setups, which would be required to excite directional modes in our metasurface by near-field sources. Previous experiments showed the coupling of Zeemaninduced circularly-polarized excited states of quantum dots to the directional edge modes in a photonic crystal [19]. However, the directionality in our metasurface is related to the circular polarisation of the in-plane magnetic field, meaning a magnetic rather than electric point source is required. At low frequencies these sources have been realized by means of coaxial cables [20,42]. On the other hand, magnetic transitions in the optical regime can be realized with rare earth ions [43], although these are usually linearly-polarized transitions.
A recent proposal makes use of two anti-parallel atomic dipoles to generate a magnetic dipole at optical frequencies [44], which could be extended to circular polarization, and possibly realized with quantum dots. Alternatively, the magnetic resonance of split-ring resonators could be engineered to provide the required source. In the following, we assume a simple model and use a dielectric NP as source, given that they support a magnetic dipole mode [45] which could be excited with a circularly polarised wave. Provided the system can be modelled in the dipole approximation with subwavelength resonators and point dipole sources, the following results are valid regardless of the physical nature of the source or frequency regime. Therefore our main conclusions also apply to metasurfaces made of microwave [20,42] or Mie resonators [46].
The scattering setup is shown in Figure 2(a). The source position is moved to various heights h in the z-direction and the xy-position is fixed at the centre of the expanded unit cell, as shown in Figure 2(a). The power through the left (P L ) and right (P R ) channels is calculated, as well as the power radiated by the source. The material losses along the interface are set to zero to test the directionality behaviour and they are gradually increased at the edge of the sample to prevent backscattering.
We define a β-factor to characterise the amount of power coupling into edge modes compared to the total power radiated by the source, where i = (L, R) corresponding to power through the left or right channel, respectively. P 0 is the power radiated by the source and P F is the Purcell factor. The total power coupling into the edge is then β edge = β L + β R . If β edge = 0, none of the energy radiated from the source couples into the edge mode, whereas if β edge = 1 all of the energy couples into the edge. β accounts for enhancement of the emission of the source due to the environment through the Purcell factor [43]. For a magnetic dipole the emitted power in free space, with vacuum permeability µ 0 and magnetic dipole moment m.
The source is placed initially at h = 0 nm and then moved upwards to h = 100 nm. We  Next, in panel (c) we plot β edge where it is clear that it is greatest for the edge state frequencies and smallest in the minigap and for bulk frequencies; in agreement with the spectral function of the edge states in Figure 1(f). A maximum of approximately 60% of the power emitted by the source couples into the metasurface at h = 0, and the rest radiates into free space as the source emits in all directions and the bulk of the metasurface is gapped.
As the source is moved upwards, β edge decreases and follows the same relationship as P F such that for h > 50 nm very little energy couples into the edge state. This confirms that the increase of P F is due to the source coupling to the edge. Finally, in Figure 2 (d, e) we plot the beta factors for the left β L and right β R edge channels. As explained previously, since we use a left circularly-polarized source which is initially placed at the centre of the unit cell, we expect and observe power flow predominantly through the left channel [15], i.e.
β L β R . This pattern of directionality holds as the source is moved upwards. In an experimental setup, this gap will be dependent on the type of source used. We will excite the system in the upper band at ω = 2.795 eV and in the lower band at ω = 2.76 eV.
In Figure 3(c) we plot the Purcell factor for the lower (left) and upper (right) bands: The pattern of each band is very different. From the ribbon eigenmodes in Figure 1(f) (and Appendix B, Figure 6), we showed that the upper and lower bands have a bonding and anti-bonding character. In Figure 3(b) we sketch the out-of-plane dipole moments for the two NPs closest to the edge. This shows how the electric field maxima and minima will be in different positions: For the anti-bonding mode there is a maximum (star) at some position above the NPs and for the bonding mode, the maximum is between the NPs. This is reflected in the peaks of the Purcell factor. Furthermore, for both bands the Purcell factor is greatest at the interface which again is due to the LDOS being greatest here. Away from the interface, it approaches unity but there is still some emission enhancement due to the proximity to the metasurface. However, this means that sources placed far from the interface will not excite edge modes.
Finally, in Figure 3(d, e) we plot the beta factors for the chosen frequencies in the lower and upper bands. β edge follows the same pattern as the Purcell factor. However, it is interesting to note that despite the upper and lower bands having different β edge patterns, the position dependent directionality can still be seen in β L (middle panels) and β R (right panels). Specifically, for both bands the coupling to the left propagating channel for a left circularly-polarized source is maximum when the source is above the centre of an expanded unit cell (lime star). As soon as the source is moved outside the hexamer of NPs (magenta star), it couples to the mode travelling in the opposite direction. For larger h this pattern will eventually be lost as the source couples less and less to the edge state.

V. FAR-FIELD CIRCULARLY-POLARIZED EXCITATIONS
A number of experimental works have investigated the directionality of these edge states under a far-field excitation, in a range of frequency regimes [17,21,22]. Despite this, there is no detailed theoretical investigation of directionality of edge modes with a far-field excitation and specifically whether the position-dependent directionality found for point sources also holds here. We investigate this by mimicking a circularly-polarized far-field Gaussian beam excitation with non-zero orbital angular momentum (OAM) incident on the metasurface. The NPs are excited with an incident electric field with a Gaussian intensity r is centre of the beam, w is the full width at half maximum (FWHM), φ is the angle in radians from the centre. The phase vortex in E z corresponds to a left or right circular polarization.
In Figure 4 To confirm this behaviour we will probe the finite system with the beam in Equation 12.
We let w = 50 nm, meaning the beam covers approximately one unit cell and the centre of the beam is moved over the red and blue paths shown in Figure 4(b). As in the previous section, we excite the system at frequencies in the lower and upper bands at ω = 2.76 eV and ω = 2.795 eV. In panels (c) and (d) we plot the fraction of power through the left (P L /P tot ) and right (P R /P tot ) channels for the two bands (plots of the total power coupling to the edge modes in each case are shown in Appendix D). Importantly, when the centre of the beam is at the centre of the expanded unit cell (y = 0 nm), the majority of power is through the left channel P L P R , as was the case with point sources. Similarly, when the beam moves across the interface (y ≈ −50 nm) the majority of power switches to the right channel P L P R . Although we have only demonstrated this for two specific frequencies here, the position dependent directionality holds for frequencies across the band gap (as we show in Appendix D, Figure 8). The vortex map in Figure 4(a) shows that it is possible to choose a path which maximises directionality, by avoiding traversing over phase singularities.
In Figure 4(d) we plot the power for the blue path shown in (a). This path is similar to scanning perpendicularly across a zig-zag interface, as in [21]. Compared to the red path, there is wider range of y-values that yield the expected directionality. Finally, we emphasise that it is the relationship between the FWHM and the lattice constant, rather than the FWHM value itself, which is important. For beams with a FWHM which is much larger than the lattice constant, directionality is mostly lost since the beam expands to cover and excite the majority of the interface.
In an experimental setup, edge modes in the metasurface could be excited with a combination of a 2D layer and far-field beam. For example, valley-selective modes can be excited using a transition metal dichalcogenide on top of plasmonic metasurface [52][53][54]. Along with the results from section IV, the directionality observed for far-field excitations suggests that a similar method could be employed to excite directional modes in the breathing honeycomb lattice.

VI. CONCLUSIONS
In this article we have provided a comprehensive study of the excitation of edge modes in a plasmonic metasurface with a breathing honeycomb lattice arrangement. The 2D lattice of metallic NPs hosts subwavelength pseudospin edge modes which arise due to the topology of the bulk. The plasmonic metasurface is a versatile system for testing the directionality of these modes through the coupled dipole method. Motivated by the excitation of chiral and valley selective modes in plasmonic metasurfaces with 2D layers, we probe the edge states of our system in the near-and far-field. With circularly-polarised magnetic dipole sources, we map the directionality of modes for sources out of the plane and show how, provided the source still couples to the interface, the pattern of directionality predicted by spin angular momentum is preserved. Additionally, we probe edge modes with far-field beams with nonzero orbital angular momentum. Here, the direction of propagation is predicted by the phase of the E z field of the edge eigenmodes. Importantly, although we have particularized to the plasmonic metasurface, the directionality behaviour of the edge modes is applicable to any classical wave system possessing the same breathing honeycomb lattice. Appendix B: Edge state eigenmodes The system of equations in Equation 6 for the ribbon interface is solved as an eigenvalue problem by letting E inc = 0. Additionally, the Green's function is linearised by letting ω = ω lsp , the localized surface plasmon frequency, in order to calculate dipole moments p.
The dipole moments for the armchair and zig-zag interfaces are shown in Figure 6. In the main text, we particularise to the subwavelength, plasmonic metasurface. We will show that the near-field position dependent directionality also holds in a photonic crystal, with constant permittivity, dielectric elements. We perform finite element simulations in COMSOL [55] on the same breathing honeycomb interface. The contracted region has s = 0.967 and the expanded s = 1.041, and the dielectric pillars are silicon with = 11.7 [12].
A left circularly-polarized magnetic dipole is used to excite modes at the three positions shown in Figure 7(a). The electric field intensity |E| 2 of the interface for the three source positions is shown in (b). When the source is at the centre of the expanded unit cell (magenta star), it couples to a left propagating mode. Whereas when the source is directly at the interface (green and cyan stars), it couples predominantly to a right propagating mode. Power through the left and right channels is calculated in the same way as in section V, however here we plot the total power through the edge (normalised to the maximum) rather than the fraction left/right. As the beam moves downwards, we see how peaks in P tot follow a qualitatively similar pattern to the beta factors calculated in Figure 3(b, d), which reflects the incident beam coupling to the edge modes. A maximum in P tot for the lower band occurs as the centre of the beam passes through the centre of two particles immediately at the interface, whereas for the upper band P tot peaks either side of this.