title: Near- and Far-Field Excitation of Topological Plasmonic Metasurfaces
authors: Matthew Proctor Xiaofei Xiao Richard V. Craster Stefan A. Maier Vincenzo Giannini Paloma Arroyo Huidobro
text: 
  
  
  
  
      matthew.proctor12@imperial.ac.uk
      Department of Mathematics, Imperial College London, London SW7 2AZ, UK
  
  
      Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom
      
  
      Department of Mathematics, Imperial College London, London SW7 2AZ, UK
      
  
      Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom
      Chair in Hybrid Nanosystems, Nano-Institute Munich, Faculty of Physics, Ludwig-Maximilians-Universität München, 80539 München, Germany
   
  
      Instituto de Estructura de la Materia (IEM), Consejo Superior de Investigaciones Científicas (CSIC), Serrano 121, 28006 Madrid, Spain
   
  
      Instituto de Telecomunicações, Insituto Superior Tecnico-University of Lisbon, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal
  
  
  
  
  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.
  
  
  Near- and Far-Field Excitation of Topological Plasmonic Metasurfaces
      Paloma Arroyo Huidobro
      December 30, 2023
  ====================================================================
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Topological nanophotonics offers a path towards efficient and robust control over light on the nanoscale <cit.>. Concepts borrowed from topological insulators, materials which host protected surface states for electrons, can also be applied to photonic systems. Following theoretical proposals <cit.>, protected photonic modes reliant on a explicit time-reversal symmetry breaking component have been demonstrated experimentally <cit.>. 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 <cit.>, 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 <cit.>, and pseudospin-dependent effects, which rely on the spin angular momentum texture of the electromagnetic fields  <cit.>. 
  
  First proposed in Ref. <cit.>, 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 ℤ_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 <cit.>, 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 <cit.>, 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 <cit.>. This result is true for any bosonic breathing honeycomb lattice <cit.> since it is rooted in time-reversal symmetry and the absence of Kramer's degeneracy for bosons as opposed to fermions. There have been a range of experimental investigations on the breathing honeycomb photonic crystal <cit.>, including the observation of edge modes in the visible regime <cit.>. 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 far-field 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.
  
  
  
  
  § 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 <cit.>. 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  <cit.>. The MWLWA polarisability α(ω) for spheroids is,
  
      α(ω) = α_s(ω)/1 - k^2/l_EDα_s(ω) - i2k^3/3α_s(ω),
  
  with the static polarisability α_s(ω),
  
      α_s(ω) = V/4πϵ(ω) - 1/1 + L(ϵ(ω) - 1).
  
  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,
  
      ϵ(ω) = ϵ_∞ - ω_p^2/ω^2 + iωγ.
  
  We use silver NPs, with ϵ_∞ = 5, ω_p = 8.9 eV and γ = 1/17 fs ≈38 meV <cit.>. 
  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 𝐄_inc,
  
      1/α(ω)𝐩_i = 𝐄_inc + 𝐆̂(𝐝_ij, ω)·𝐩_j.
  
  The interactions between dipoles are characterised by the dyadic Green's function <cit.>,
  
      Ĝ(d_ij, ω)    = k^2e^ikd/d[
          (
          1 + i/kd - 1/k^2d^2)Î - 
          (
          1 + 3i/kd - 3/k^2d^2)n⊗n],
  
  where 𝐝_ij is the separation between dipole i and j, d = |𝐝_ij| and 𝐧 = 𝐝_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-of-plane modes <cit.>. 
  For a periodic array of NPs in a plasmonic metasurface we can apply Bloch's theorem and write the following system of equations,
  
      (𝐈̂1/α(ω) - Ĥ(k_∥, ω))·p = 𝐄_inc.
  
  The vector 𝐩 contains the out-of-plane dipole moments p_z of all NPs in the unit cell. The interaction matrix Ĥ(k_∥, ω) has elements,
  
      H_pq=
          ∑_RĜ(d_pq + R, ω)  e^ik_∥·R    p ≠ q
      ∑_|R|≠0Ĝ(R, ω)  e^ik_∥·R    p = q
          ,
  
  where p, q are unit cell indices. The summations run over the lattice sites 𝐑=n𝐚_1 + m𝐚_2, with lattice vectors 𝐚_1 and 𝐚_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 <cit.>. When modelling plasmonic metasurfaces, it is necessary to include the long range, retarded interactions in the Green's function in <ref>. Whilst the quasistatic approximation can be used to model very subwavelength plasmonic chains and arrays <cit.>, 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 <cit.>. 
  
  By solving <ref>, we can calculate the extinction cross section σ_ext from the dipole moments and incident electric field using the optical theorem,
  
      σ_ext= k Im(𝐩·𝐄^*_inc)/|𝐄^2|.
  
  The incident field satisfies Maxwell's equations, |𝐤_∥| 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 - |𝐤_∥|^2) and when it is evanescent, below the light line, k_z = i√(|𝐤_∥|^2 - k^2). We can also solve <ref> as an eigenvalue problem and calculate the spectral function, by letting 𝐄_inc = 0, which allows us to probe modes whether they are bright or dark. The spectral function σ_sf is,
  
      σ_sf = k Im(α_eff),
  
  with the effective polarisability α_eff = ∑_i 1/λ^(i) for eigenvalues λ^(i) <cit.>. Finally, the coupled dipole method can also be applied to finite systems to model electromagnetic scattering <cit.>. 
  
  
  
  
  
  
  
  
  § 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 <ref>(a) with a unit cell with six NPs <cit.>. 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 = 120nm. 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 <ref>, <ref>). We note that the Dirac cone lies below the localised surface plasmon frequency ω_lsp due to chiral symmetry breaking, long range interactions <cit.>. By perturbing the nearest neighbour separation with a scale factor s, such that R = sR_0, a band gap is opened at Γ. In <ref>(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. The scale factors are chosen to ensure that the size of the band gaps for the two lattices approximately equal. We plot the bands from the spectral function as blue dots on top of the extinction cross section as a guide. Below the light line (white dotted line), where modes are confined, the modes have high quality factors whereas above the light line the modes are very broad corresponding to larger radiative losses.
  
     
  
  In the shrunken lattice, for 𝐤_∥ 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 <cit.>.
  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 <cit.>. We plot σ_ext for a fixed wavevector 𝐤_∥ in  <ref>(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 ℤ_2-invariant. However, they are still topologically distinct and one way to distinguish them is by their Wilson loops <cit.>. The shrunken lattice is topologically trivial whilst the expanded lattice is in an photonic `obstructed atomic limit' phase <cit.>; it is the C_6 symmetry which provides the topological character.
  
  Despite having trivial ℤ_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 <cit.>. 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 <cit.>. We model the edge states of the system by setting up a ribbon with an interface between the two phases, as shown in <ref>(b). The ribbon supercell is N=20 unit cells along the 𝐚_2 direction, with 10 unit cells in the expanded and shrunken phases, respectively. In <ref>(f), we plot the spectral function σ_sf for the ribbon for 𝐤_∥ = 0 to π/a_0. (The spectral function from 𝐤_∥ = 0 to -π/a_0 is identical.) Edge states appear in the band gap from ω≈ 2.74 eV to 2.81eV 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 <cit.>, including the `zig-zag' interface <cit.>. 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 <ref>, <ref>). Due to doubly-periodic supercell, we see multiple narrow bands bending down and crossing the band gap in <ref>(f). This is due to diffraction orders above the light line (white dotted line) <cit.> and does not affect the investigation of edge states which follows. 
  
  
  
  § CIRCULARLY-POLARIZED POINT SOURCES
  
  
  We begin by extending the investigation of our previous work <cit.> by exciting the system with point sources with circularly-polarized magnetic fields, H = H_x ± i H_y. 
  The directionality of modes is well understood for sources in the plane of the metasurface, h=0: The inhomogeneous spin angular momentum of the edge eigenmode determines the direction of propagation, rather than the polarization of the source. This means it is possible to selectively couple to either of the counter-propagating edge states by moving the source in the plane of the metasurface; a result which holds for any bosonic system with this lattice, including photonic crystals <cit.> (We show this explicitly in Appendix <ref>).
  
  
  It is important to understand whether it is still possible to excite edge modes for sources placed out of the plane, including whether directionality is still observed. This is relevant to experimental setups where quantum dots and emitters, or 2D materials are placed on top of a spacing layer above the NPs in the metasurface. 
  
  
  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 Zeeman-induced circularly-polarized excited states of quantum dots to the directional edge modes in a photonic crystal <cit.>. 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 <cit.>. 
  On the other hand, magnetic transitions in the optical regime can be realized with rare earth ions <cit.>, 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 <cit.>, 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 <cit.> 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 <cit.> or Mie resonators <cit.>. 
  
  
  The scattering setup is shown in <ref>(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 <ref>(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,
  
      β_i = P_i/P_0 P_F,
  
  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 <cit.>. For a magnetic dipole the emitted power in free space,
  
      P_0 = μ_0ω^4/3c^3|𝐦|^2,
  
  with vacuum permeability μ_0 and magnetic dipole moment 𝐦. 
  
  
  
  The source is placed
  
  initially at h=0 nm and then moved upwards to h=100 nm. We use a left circularly-polarized dipole throughout (H = H_x - iH_y); which, in combination with the xy-position at the centre of the expanded unit cell, means we expect it to couple to the left propagating edge mode. Additionally, we scan over the frequency range of the band gap to investigate whether exciting in the upper or lower band has an effect on the directionality. We measure the normalised power through yz-planes perpendicular to both the metasurface and the interface. First, in <ref>(b) we plot the Purcell factor. (A horizontal white dashed line serves as a guide to the height of the NPs in the metasurface and the band gap is highlighted with vertical dashed lines). We see that the Purcell factor is greatest near the metasurface since the array of NPs enhances the emission of the source by increasing the available local density of states (LDOS). For h > 50nm, there is almost no coupling to the metasurface and P_F ≈ 1. Scanning across frequencies at h close to zero, we can clearly see the minigap near 2.77 eV, where no NPs are excited which causes P_F to decrease. Lastly, P_F is largest at ω≈ 2.81 eV where the source begins to excite bulk modes.
  
  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 <ref>(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 > 50nm 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 <ref> (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 <cit.>, i.e. β_L≫β_R. This pattern of directionality holds as the source is moved upwards. 
  
  
  
     
  
  Following the understanding of directionality for sources at different heights, we now investigate the position dependence in a plane parallel to the metasurface.
  The sources are placed in a range of positions in the xy-plane at a fixed height, h=30 nm, close to the interface as shown in <ref>(a). We choose this height to ensure a 10nm gap between the top of the NPs and the excitation layer, and to ensure enough power couples into the edge. 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.795eV and in the lower band at ω = 2.76eV.
  
  In <ref>(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 <ref>(f) (and Appendix <ref>, <ref>), we showed that the upper and lower bands have a bonding and anti-bonding character. In <ref>(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 <ref>(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.
  
  
  
  
  
  § 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 <cit.>. 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 profile, 𝐄_inc = (0, 0, E_z) with,
  
      E_z ∝ e^-r^2/(2w^2) e^iϕ,
  
  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 <ref>(a), we plot the phase of the z-component of the electric field of the edge eigenmodes: for the lower band, at k_x > 0 and for the upper band, at k_x < 0. A phase vortex structure is present across the interface, with clearly distinct regions at the centres of unit cells and at the edge of unit cells.
  Mode matching and maximum directionality occurs when the phase vortex of the beam rotates in the the same direction as the phase vortex of the edge eigenmode <cit.>, these positions are highlighted as clockwise arrows. Mode mismatch occurs when the beam and eigenmode have vortices rotating in opposite directions and these positions are shown are anti-clockwise arrows. When there is a complete mode mismatch, the beam will excite a mode travelling in the opposite direction.
  The phase vortex behaviour of directional edge states has previously been studied in the context of valley modes <cit.>, where valleys at K and K' have vortices rotating in opposite directions. 
  
     
  
  
  
  To confirm this behaviour we will probe the finite system with the beam in <ref>. We let w = 50nm, meaning the beam covers approximately one unit cell and the centre of the beam is moved over the red and blue paths shown in <ref>(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 <ref>). 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≈-50nm) 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 <ref>, <ref>). 
  
  The vortex map in <ref>(a) shows that it is possible to choose a path which maximises directionality, by avoiding traversing over phase singularities.
  In <ref>(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 <cit.>. 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 <cit.>. Along with the results from <ref>, the directionality observed for far-field excitations suggests that a similar method could be employed to excite directional modes in the breathing honeycomb lattice. 
  
  
  
  
  
  
  
  § 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 non-zero 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.
  
  
  M.P. and P.A.H. acknowledge funding from the Leverhulme Trust. P.A.H. acknowledge financial support from Fundação para a Ciência e a Tecnologia and Instituto de Telecomunicações under projects UID/EEA/50008/2020 and the CEEC Individual program with reference CEECIND/03866/2017.
  
  
  
  
  
  
  
  
  § HONEYCOMB LATTICE: OPTICAL RESPONSE OF BULK MODES
  
  
   
  
  
  
  § EDGE STATE EIGENMODES
  
  
  The system of equations in <ref> for the ribbon interface is solved as an eigenvalue problem by letting 𝐄_inc = 0. Additionally, the Green's function is linearised by letting ω = ω_lsp, the localized surface plasmon frequency, in order to calculate dipole moments 𝐩. The dipole moments for the armchair and zig-zag interfaces are shown in <ref>.
  
     
  
  
  
  § DIRECTIONALITY OF MODES IN A PHOTONIC CRYSTAL
  
  
  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 <cit.> 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 <cit.>. 
  
  A left circularly-polarized magnetic dipole is used to excite modes at the three positions shown in <ref>(a). The electric field intensity |𝐄|^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.
  
  
  
  
  
  § FAR FIELD EXCITATIONS: POWER THROUGH LEFT AND RIGHT CHANNELS
  
  
  Power through the left and right channels is calculated in the same way as in <ref>, 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 <ref>(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. 
  
  
     
  
  
  
  