Exciting Pseudospin Dependent Edge States in Plasmonic Metasurfaces

TL;DR

This paper investigates plasmonic metasurfaces supporting pseudospin dependent edge states, revealing their properties, propagation behavior, and optical responses, with potential implications for topological photonics.

Contribution

It introduces a detailed analysis of pseudospin dependent edge states in plasmonic metasurfaces considering full electrodynamic effects, expanding understanding beyond idealized models.

Findings

Edge states exhibit pseudospin dependence influenced by electromagnetic spin angular momentum.

Propagation direction is determined by local spin angular momentum sign.

Radiation and retardation significantly affect the optical response of the metasurface.

Abstract

We study a plasmonic metasurface that supports pseudospin dependent edge states confined at a subwavelength scale, considering full electrodynamic interactions including retardation and radiative effects. The spatial symmetry of the lattice of plasmonic nanoparticles gives rise to edge states with properties reminiscent of the quantum spin Hall effect in topological insulators. However, unlike the spin-momentum locking characteristic of topological insulators, these modes are not purely unidirectional and their propagation properties can be understood by analysing the spin angular momentum of the electromagnetic field, which is inhomogenous in the plane of the lattice. The local sign of the spin angular momentum determines the propagation direction of the mode under a near-field excitation source. We also study the optical response under far-field excitation and discuss in detail the effects of radiation and retardation.

Equations

$$\frac{1}{\alpha (\omega )}{\text{p}}_i={\text{E}}_i+\sum _{i\ne j}\hat{\text{G}}({\text{d}}_i-{\text{d}}_j,\omega ){\text{p}}_j,$$

$$\hat{\text{G}}({\text{d}}_i-{\text{d}}_j,\omega )=k^2\frac{e^{ikd}}{d}\left[\left(1+\frac{i}{kd}-\frac{1}{k^2{d}^2}\right)\hat{\text{I}}-\left(1+\frac{3i}{kd}-\frac{3}{k^2{d}^2}\right)\text{n}\otimes \text{n}\right].$$

$${\alpha }_s(\omega )=\frac{V}{4\pi }\frac{ϵ(\omega )-1}{1+L\left(ϵ(\omega )-1\right)},$$

$$ϵ(\omega )={ϵ}_{\infty }-\frac{{\omega }_p^2}{{\omega }^2+i\omega \gamma }.$$

$${\alpha }_{\text{MLWA}}(\omega )=\frac{{\alpha }_s(\omega )}{1-\frac{k^2}{l_E}D{\alpha }_s(\omega )-i\frac{2k^3}{3}{\alpha }_s(\omega )},$$

$$\left(\frac{1}{\alpha (\omega )}-\hat{\text{H}}(\text{k},\omega )\right)\cdot \text{p}=\text{E},$$

$$H_{ij}=\{\begin{array}{ll}\sum _{\text{R}}\hat{\text{G}}({\text{d}}_i-{\text{d}}_j+\text{R},\omega )e^{i\text{k}\cdot \text{R}}& i\ne j\\ \sum _{|\text{R}|\ne 0}\hat{\text{G}}(\text{R},\omega )e^{i\text{k}\cdot \text{R}}& i=j\end{array},$$

$$P_M=\frac{{\mu }_0}{4\pi }\frac{{\omega }^4|\text{m}{|}^2}{3c^3},$$

$$S_{\text{incl}}=\sum _{\text{R}}{e}^{i\text{k}\cdot \text{R}}\hat{\text{G}}(\text{r},\omega )=S_L+S_{SM},$$

$$S_{SM}\approx k^2\sum _{\text{R}=0}^{{\text{R}}_{min}}{e}^{i\text{k}\cdot \text{R}}{e}^{ikd}\left(\frac{ik}{d^2}-\frac{1}{d^3}\right)+k^2{\int }_{R_{min}}^{\infty }{e}^{i\text{k}\cdot \text{R}}{e}^{ikd}\left(\frac{ik}{d^2}-\frac{1}{d^3}\right).$$

$$S_{\text{excl}}=\sum _{\text{R}\ne 0}{e}^{i\text{k}\cdot \text{R}}\hat{\text{G}}(\text{R},\omega ),$$

$$\left(\hat{\text{H}}(\text{k},\omega )-\frac{1}{\alpha (\omega )}\hat{\text{I}}\right)\cdot \text{p}=0.$$

$$\left(\hat{\text{H}}(\text{k},{\omega }_{sp})-\frac{1}{\alpha (\omega )}\hat{\text{I}}\right)\cdot \text{p}=0,$$

$$\left(\hat{\text{H}}(\text{k},\omega )-\frac{1}{\alpha (\omega )}\hat{\text{I}}\right)\cdot \text{p}={\text{E}}_{\text{inc}}.$$

$$E_z=\frac{-k_{\parallel }}{k_z}=\frac{-k_{\parallel }}{k_0^2-k_{\parallel }^2}.$$

$${\sigma }_{\text{ext}}=\frac{4\pi k{\sum }_i\text{Im}({\text{p}}_i\cdot {\text{E}}_{\text{inc}}^{\ast })}{|\text{E}{|}^2}.$$

$${\sigma }_{\text{spectral}}=4\pi k\sum _i\text{Im}({\alpha }_{\text{eff}}^{(i)}),$$

Full text

Imperial College London]Department of Mathematics, Imperial College London, London, SW7 2AZ, UK \alsoaffiliationDepartment of Physics, Imperial College London, London, SW7 2AZ, UK

Imperial College London]Department of Mathematics, Imperial College London, London, SW7 2AZ, UK

Munich]Chair in Hybrid Nanosystems, Nanoinstitut München, Faculty of Physics, Ludwig-Maximilians-Universität München, 80539, München, Germany \alsoaffiliationDepartment of Physics, Imperial College London, London, SW7 2AZ, UK

Madrid]Instituto de Estructura de la Materia (IEM), Consejo Superior de Investigaciones Científicas (CSIC), Serrano 121, 28006, Madrid, Spain

Lisbon]Instituto de Telecomunicações, Insituto Superior Tecnico-University of Lisbon, Avenida Rovisco Pais 1,1049-001 Lisboa, Portugal \alsoaffiliationDepartment of Physics, Imperial College London, London, SW7 2AZ, UK

\abbreviations1D, 2D, SSH, LSP, QSA, CDA, PTI

Topological insulators are materials which are insulating in the bulk but which have conduction surface states protected against disorder (1). The remarkable properties of these states in electronic systems has inspired the search for photonic topological insulators (PTIs), which aim to guide and manipulate photons with the same level of control and efficiency (2, 3, 4, 5, 6). Systems which possess these effects whilst preserving time reversal symmetry are appealing as they do not require complicated experimental setups such as strong magnetic fields or bianisotropic coupling. Motivated by this, a proposal to emulate the quantum spin Hall (QSH) effect in photonic crystals was presented by Wu and Hu in Ref. 7. Effects reminiscent of the QSH phase such as a band inversion between dipolar and quadrupolar modes, and pseudospin dependent edge states are realised but, rather than relying on the time reversal symmetric pairs characteristic of electronic systems, they instead rely on the spatial symmetry of the lattice structure. As a result, the edge states have a reduction in backscattering over trivial ones (8, 9). The method has since been applied to a variety of bosonic systems (10, 11, 12, 13, 14), and has recently experimentally been demonstrated in the visible regime (15).

The combination of topological effects with plasmonics offers the possibility of precisely controlling light on the the nanoscale. The strong enhancement and localisation of electric fields due to localised surface plasmon (LSP) resonances (16) is a widely employed platform for light confinement on the nanoscale (17, 18). Plasmonic metasurfaces can be formed by arranging plasmonic nanoparticles in two-dimensional (2D) lattices, where the LSPs become delocalised across the whole metasurface as collective resonances. The optical properties of metasurfaces are then determined by the individual nanoparticle elements as well as the geometry of the lattice (19, 20). The tunable optical properties of metasurfaces makes them versatile tools for the manipulation of light on the nanoscale (21, 22). For instance, appropriately designed plasmonic metasurfaces can host spin dependent directional states which can couple to valley excitons when interfaced to 2D materials (23).

One-dimensional (1D) chains of dielectric and metallic nanoparticles were some of the first systems used for the investigation of topological phases in nanophotonics (24, 25, 26, 27, 28, 29, 30), in particular systems analogous to the Su-Schrieffer-Heeger (SSH) model, which hosts topologically protected edge states in 1D. In the plasmonic chain, initial studies into these topological states were limited to the quasistatic approximation (QSA) (31) despite the well known radiative effects which are of great importance for large enough nanoparticles and retardation which is important at large lattice periods (32, 33) or at very small periods, where higher order multipolar effects occur (34). Indeed, beyond the quasistatic limit, the plasmonic SSH chain becomes non-Hermitian and band structures are distorted compared to the quasistatic model, with effects such as polariton splitting at the light line (35). In addition, the ubiquitous bulk-edge correspondence of topological insulators has been shown to break down due to retardation (36). 2D plasmonic systems, including graphene and arrays of plasmonic nanoparticles have also been considered for hosting topological states. Time reversal broken effects reliant on magnetic fields have been proposed for graphene plasmons (37, 38). Honeycomb lattices of plasmonic nanoparticles have been investigated in the quasistatic limit, where a direct analogy between the tight binding model in graphene and nearest neighbour approximation in plasmonics can be made (39, 40, 41). More recently, theoretical and experimental investigations have shown how the long range, retarded interactions affect the physical behaviour of the system (42, 43). Here, we consider the lattice geometry proposed in Ref. 7, shown in Figure 1, to study pseudospin dependent edge states in plasmonic metasurfaces. This scheme has been considered in 2D arrays of plasmonic nanoparticles in the QSA (44). In this work, we show the importance of going beyond this approximation and study the realisation of electromagnetic modes resembling the QSH effect on a plasmonic metasurface, including full electrodynamic interactions in our description of the system.

We use semi-analytical techniques to investigate plasmonic metasurfaces. These consist of a triangular lattice with unit cells containing six nanoparticles arranged in a hexagon, Figure 1a. We begin by outlining the coupled dipole method used to model these arrays and then we investigate the behaviour of the modes supported by the metasurface. By calculating the eigenmodes of the infinite lattice, we characterise the band inversion process which occurs between a shrunken phase (where nanoparticles in the unit cell are displaced towards the centre) and an expanded phase (where they are displaced outwards from the centre). The response of the infinite system is also studied under far-field excitation.

We then consider the interface between two phases in a semi-infinite ribbon layout to elucidate the nature of the edge states that emerge due to the band inversion between the two phases. We characterise their pseudospin dependence, showing that these modes are more appropriately characterized by means of the spin angular momentum of the electromagnetic fields, and probe these states in the near-field with a magnetic dipole source. By investigating the effect of the source position, we unambiguously show that, unlike the unidirectional edge states characteristic of topological insulators, the directionality of these edge states depends on the source position. Finally, we highlight the radiative and retardation effects on the edge states.

Design and set up of the plasmonic metasurface

The metasurface we consider here consists of a 2D array of metal nanorods (modelled as spheroidal nanoparticles) arranged in the lattice shown in Figure 1a. The unit cell contains six nanoparticles of radius $r$ and height $h$ arranged in a hexagon separated by nearest neighbour spacing $R$ with lattice vectors $\textbf{a}_{1}$ and $\textbf{a}_{2}$, as shown in the figure. In the regime $R>3r$, the nanoparticles can be considered as point dipoles and higher order resonances can be neglected (45, 32). A nanoparticle at position $\textbf{d}_{i}$ with polarisability $\alpha(\omega)$ supports a dipole moment $\textbf{p}_{i}$ when excited by an external field $\textbf{E}_{i}$. For a lattice of nanoparticles, a coupled dipole equation can be written which describes the dipole moment of a nanoparticle induced by an external field plus the sum of all neighbouring dipole moments,

[TABLE]

where the dyadic Green’s function $\hat{\textbf{G}}$, describes the interactions between point dipoles, and is given by,

[TABLE]

Here, $k=\sqrt{\epsilon_{B}}\omega/c$ is the wavenumber of the surrounding medium, with $\epsilon_{B}$ the permittivity of the surrounding environment (which is assumed to be the vacuum throughout this work, $\epsilon_{B}=1$), and $d=|\textbf{d}_{i}-\textbf{d}_{j}|$ is the distance between nanoparticles, with $\textbf{n}=(\textbf{d}_{i}-\textbf{d}_{j})/|\textbf{d}_{i}-\textbf{d}_{j}|$ the unit vector in the direction along the line that joins two nanoparticles.

In equation 1, the polarisability, $\alpha(\omega)$, describes the optical response of an individual nanoparticle. The static polarisability for a spheroidal nanoparticle is written,

[TABLE]

where $\epsilon(\omega)$ is the dielectric function of the metal, $V$ is the spheroid volume, and $L$ is the static geometrical factor which is dependent on the radius and height of the nanoparticle; for a sphere $L=\frac{1}{3}$ (46). The dielectric function of the nanoparticles is given by the Drude model,

[TABLE]

Throughout the work we consider silver nanoparticles, with $\epsilon_{\infty}=5$, $\omega_{p}=8.9$ eV and $\gamma=1/17$ fs $\approx 0.039$ eV (47). The static polarisability neglects radiative effects which are essential for describing larger nanoparticles. We take this into account by means of the modified long wavelength approximation (MLWA),

[TABLE]

where $l_{E}$ is the spheroid major axis half length and $D$ is a dynamic geometrical factor; $D=1$ for a sphere (46). The importance of the radiative correction for spheroidal silver nanoparticles is exemplified in Figure 1c. The extinction cross section $\sigma_{\text{ext}}$ for a spheroidal nanoparticle with radius $r=5$ nm and height $h=20$ nm, as well as a spherical nanoparticle with radius $r=5$ nm is shown. We normalise to the cross sectional area $A$ perpendicular to the dipole moment. The nanoparticle supports two resonance modes: one where the dipole is aligned with the minor axis (in-plane) and one with the major axis (out-of-plane). The out-of-plane resonance becomes redshifted and well separated in frequency from the in-plane resonance, which allows us to investigate the out-of-plane and in-plane modes separately. The increased radiative effect on a single nanoparticle is demonstrated by the extinction cross section for larger radii, up to $r=20$ nm, where the resonance becomes broader and continues to be redshifted, Figure 1d. Whilst the static polarisabilty adequately describes the behaviour of the smallest nanoparticle sizes, the MLWA incorporates the effects of dynamic depolarisation and the radiative correction, and is necessary to correctly model particles of radius above $\sim 10$ nm.

Spectral response of the metasurface

We start by considering the optical response of the plasmonic metasurface in the expanded and shrunken phases, Figure 1b. We do so by setting up an infinite lattice of nanoparticles and applying periodic boundary conditions to Equation 1 by writing the external electric field $\mathbf{E}$, and dipole moments $\mathbf{p}$, as periodic Bloch functions. The following system of equations can then be written,

[TABLE]

where the interaction matrix $\hat{\textbf{H}}(\textbf{k},\omega)$ has elements,

[TABLE]

with q, the Bloch wavevector and $\mathbf{R}$, the lattice site positions. We note that the interaction matrix is a $6\times 6$ matrix since we restrict our study to the out-of-plane modes, such that there is a single degree of freedom for each particle in the unit cell. The in-plane modes of a honeycomb lattice of plasmonic nanoparticles have been studied elsewhere (41) and being well shifted in frequencies, they are completely decoupled from the out-of-plane modes. Finally, we note that sums in the interaction matrix are conditionally convergent due to the slowly decaying $1/r$ term and so additional manipulation is required to converge these expressions (see Methods).

Figure 2 shows the spectral response of the plasmonic metasurface under study. We consider nanoparticles with radius $r=5$ nm and height $h=20$ nm, nearest neighbour spacing $R=20$ nm and lattice constant $a_{0}=60$ nm. To calculate the eigenvalues and eigenvectors of the periodic lattice we solve Equation 6 without an incident field. Initially, in order to keep the eigenvalue problem linear we take the QSA and only consider the quickly decaying $1/d^{3}$ term in the Green’s function. However, we choose to go beyond the nearest neighbour approximation and include interactions between all particles in the lattice (48). The band structure of a metasurface with the nanoparticles arranged in a honeycomb lattice is shown in Figure 2a. Instead of using the conventional rhombic unit cell, we take the larger hexagonal cell such that the Brillouin zone (BZ) of the honeycomb lattice becomes folded and the Dirac points at $K$ and $K^{\prime}$ are mapped onto each other to create a doubly degenerate point at $\Gamma$ (7), as shown in Figure 2a. Whereas the original honeycomb lattice is formed of two triangular sublattices, the system is now formed of six sublattices corresponding to the six nanoparticles in the unit cell; meaning there are six bands present in the band structure. In the QSA, assuming nearest neighbour interactions, the Dirac points of the honeycomb plasmonic lattice occur at the surface plasmon resonance frequency $w_{sp}$ (41) and the band structure is symmetrical about this frequency (44). This is no longer the case due to the sublattice symmetry breaking interaction term between particles of the same sublattice in neighbouring unit cells. Figure 2c shows the band structures for metasurfaces with shrunken (left) and expanded (right) unit cells. The bands are labelled as $s,\,p,\,d,\,f$, indicating monopolar, dipolar, quadrupolar and hexapolar characters. This was determined by calculating the overlap with the eigenstates of an isolated hexagon of nanoparticles. The dipole and quadrupole modes of the lattice are shown in Figure 2b. In both the shrunken and expanded phases, we see how the the ordering of the modes is opposite to that obtained in photonic crystals and other bosonic analogues (7, 49, 13, 14), with the monopolar mode being the highest in energy and the hexapolar being the lowest in energy. This is due to the different electromagnetic properties of the two systems: while in the photonic crystal the permittivity is positive and constant, in the plasmonic system the permittivity is negative and dispersive. For the metasurface considered here, the bonding or antibonding character of the coupling between the plasmonic nanoparticles determines the mode ordering, as we explain in more detail later. Next, we show that near the centre of the BZ, there is a mode inversion between the shrunken and expanded phases. This is evident from the colour scale, which encodes the dipolar/quadrupolar character of the modes. In the shrunken phase ($s<1$), the band above the band gap is dipolar and the one below is quadrupolar (Figure 2c). On the other hand, for the expanded phase ($s>1$) they become inverted at $\Gamma$ (Figure 2d). The degeneracy of the bands above and below the gap at $\Gamma$ for the shrunken and expanded lattices suggests linear combinations of these modes can be taken, $p_{\pm}=(p_{x}\pm ip_{y})$ and $d_{\pm}=(d_{x^{2}-y^{2}}\pm id_{xy})$, which correspond to pseudospins; taking positive or negative combinations gives clockwise or anticlockwise rotations (7).

After characterising the eigenstates in the QSA, we now consider the full electrodynamic interaction between the nanoparticles by including all terms in the Green’s function, and explore retardation and radiative effects by calculating the extinction cross section of the system when excited by an external field. We show the response of the system in the shrunken and expanded phase over the BZ in Figure 2e-f (contour plots). The incident plane wave is defined as $\textbf{E}=(E_{x},E_{y},E_{z})$, where the field components satisfy Maxwell’s equations and $k_{0}^{2}=k_{\parallel}^{2}+k_{z}^{2}$. Above the light line, the wave is propagating in the $z$ direction, $k_{z}=\sqrt{k_{0}^{2}-k_{\parallel}^{2}}$ and below the light line it is evanescent, $k_{z}=i\sqrt{k_{\parallel}^{2}-k_{0}^{2}}$, which allows us to probe the eigenstates in this region of the spectrum. Starting at $\Gamma$, the incident field has no $z$-component meaning out-of-plane dipoles cannot be excited and the extinction cross section is zero. As soon as we move away from $\Gamma$, we begin to excite out-of-plane polarised modes. The highest energy mode corresponds to a monopolar mode, with all dipole moments in the unit cell pointing in the same direction. The energy ordering can be understood from the ordering of bonding and anti-bonding modes in a plasmonic dimer. The bonding mode lies at a higher energy compared to the anti-bonding mode, as has been shown theoretically (50) and experimentally (51). This is a consequence of the relative orientation of the dipoles: in the antibonding mode the dipoles are antiparallel which minimizes the radiated electric field and hence the energy, whereas in the bonding mode the dipoles are parallel and hence sustain a larger electric field and appear at higher energies. When plasmonic nanoparticles are arranged in a lattice with more than one particle per unit cell, this ordering is mantained, as has also been shown elsewhere (42, 52). Notably, this is not an effect of retardation and will arise in any near-field coupled ensemble of plasmonic nanoparticles, as shown by our quasistatic results (Figure 2c,d) and in Refs. (50, 42, 52).

In the infinite lattice, the bonding nature of the monopolar mode across the whole system causes it to be highly radiant and so dominates the response of the lattice for propagating waves. Nevertheless, we are still able to examine the peaks in extinction cross section at lower frequencies. By symmetry arguments, only the dipolar mode is excitable by a plane wave. The extinction cross section of the system agrees with the characterisation of the modes in the QSA. In the shrunken phase, the higher energy dipolar mode is visible in the extinction cross section whereas in the expanded phase a band inversion occurs. This inversion between dipolar and quadrupolar modes close to the BZ centre constitutes a signature to distinguish both metasurface phases, since it can be detected experimentally by far-field measurements (49). In the right panels of Figure 2e-f we show the extinction spectrum at a fixed incident momentum, k: the resonance peak corresponding to the dipolar mode is visible for the shrunken structure at higher energies than for the expanded structure. In Figure 2e-f, we plot the loci of peaks in the spectral function based on the effective polarisability (see Methods), to make all of the modes visible, even those which are dark in an external field. These peaks qualitatively agree with the QSA dispersion relation with only a slight redshift due to the radiative correction given that all the length scales in the system are very subwavelength. Retardation effects are most apparent in the highest energy band. The coupling of the plasmonic mode with free photons causes the strong polariton-like splitting at the light line (33, 48). Since the monopolar (bonding) mode is strongly radiant, it shows the largest interaction with the light line. We also observe how the radiative broadening of the mode grows larger going from $\Gamma$ towards the light line; this has been observed in 1D plasmonic chains and the effect is much greater in this 2D lattice due to constructive interference between dipole moments in the monopolar mode (42). The polariton splitting and radiative broadening importantly do not affect the band inversion at $\Gamma$. Finally, the lowest energy band is a hexapolar mode where all dipole moments are in anti-phase meaning such a mode cannot be excited by a plane wave, and it is not visible in the extinction cross section for most of the BZ.

Pseudospin edge states in finite systems

After characterising the modes in the metasurface, we now look at the edge states between regions of different phases. As discussed earlier, the band inversion between different phases is reminiscent of the QSH effect and thus pseudospin dependent, directional edge states are anticipated. In order to study the edge states, we consider an interface between the two different phases of the metasurface. In the calculations, a finite ribbon of a region in the expanded phase is cladded with regions in the shrunken phase so that we have a finite structure along the direction $\textbf{a}_{2}$, and we apply Bloch periodicity along $\textbf{a}_{1}$. Figure 3a shows the band structure of the ribbon, calculated with retardation (black lines), and also in the QSA (grey lines). In the case of full electrodynamics interactions, the Green’s function was linearised by letting $\omega=\omega_{sp}$ and the static polarisability was assumed (see Methods). Two edge states can be seen clearly within the band gap (we colour the retarded solution, as we will explain below). Importantly, the edge states do not join the bulk modes, rather they reconnect with each other at the edges of the BZ. This is a fundamental distinction with topological edge states since it means it is possible, by a continuous change of parameters, to remove these states from the band gap; which indicates that these edge states are not resistant to backscattering (53). In Figure 3b, we also plot the spectral function of the system with radiative corrections in the polarisability, which accounts for the frequency shift (in this plot we let the Drude losses $\gamma=0.01$ eV to improve the visibility of the edge states). Similarly to the bulk dipolar and quadrupolar bands, we see the edge states do not interact strongly with the light line although the radiative broadening and redshift is apparent. Since the pseudospin effect is reliant on the $C_{6v}$ symmetry of the bulk system which is necessarily broken at the interface between the shrunken and expanded region there will always be a ‘minigap’ between the edge states at $\Gamma$ which will never fully close. The amount of perturbation between shrunken and expanded phase determines the size of the bulk band gap and as well as the size of the minigap. To minimise the gap, the perturbation along the edge can be graded; the edge states are then excitable across the whole band gap (44, 11).

To investigate the nature of the edge modesand to elucidate their excitation under point sources we now look at the two eigenstates of the expanded/shrunken interface with opposite group velocity, and at a frequency $\omega$ in the upper band. These are shown as points $s_{1}$ and $s_{2}$ in Figure 3a. We first plot the time averaged Poynting vector $\textbf{S}=\frac{1}{2}\text{Re}(\textbf{E}\times\textbf{H}^{*})$ along with the normalised out-of-plane electric field $E_{z}$, in the plane of the metasurface $z=0$, in Figure 3c. This demonstrates that the edge modes are confined to the interface between the two regions. The Poynting vector (arrows) characterises the flow of electromagnetic energy, which is in opposite directions for wavevectors with opposite sign. In inhomogeneous, dispersive media it can be argued that the Poynting vector does not necessarily accurately characterise a spin or pseudospin (54). More appropriately, we also consider the spin angular momentum in the plane of the lattice, $\textbf{T}=\text{Im}(\textbf{E}^{*}\times\textbf{E}+\textbf{H}^{*}\times\textbf{H})$ (55). Since the dipole moments are aligned out of the plane, at $z=0$ the field components are $E_{z},H_{x},H_{y}$ and the spin angular momentum is then $T_{z}=\text{Im}(\textbf{H}^{*}\times\textbf{H})$. This spin angular momentum quantifies the degree of elliptical polarization of the magnetic field and its handedness, similar to the studies for photonic crystals in Refs. (56, 57). We plot the normalised spin over the same region as the Poynting vector in Figure 3c, where its inhomogenous character across the lattice is evident, with varying magnitude and sign across the interface. To determine the pseudospin dependence of the edge states, we integrate the spin angular momentum $T$ over the region shown in red in Figure 3d, for wavevectors $k_{\parallel}$ across the whole BZ, and we plot the edge states with this colour code in the dispersion relation in Figure 3a. At the edges of the BZ there is maximum mixing between pseudospins and the integral of $T=0$, but as we move towards the centre we see either edge state acquires an opposite pseudospin. This pattern of opposite pseudospins travelling in opposite directions is akin to the QSH effect. Again, we note that since in this system the two edge modes are linked through a pseudo-time reversal operator (7), which is only rigourously defined at $\Gamma$, there is not complete protection against backscattering.

We now consider the excitation of edge modes with a localised source, and its relation to the spin angular momentum. The sign of the spin angular momentum $T$ is positive in the region within the hexagon of nanoparticles immediately above the edge for the mode $s_{1}$, as shown by the purple area in Figure 3d, left panel. However, at a point directly on the edge, shown by the black line in Figure 3d, the sign of $T$ switches. Importantly, this means that the local handedness of the elliptical polarization of the magnetic field for each of the edge modes is inhomogeneous in the plane of the array, as expected for a confined mode in a complex environment. This position dependence of the spin angular momentum has implications on the excitation of unidirectional modes by point sources, which we show by modelling a finite array of nanoparticles with an interface between an expanded and shrunken region. To fully understand how the propagation along the edge is dependent on the source and its position, we consider an interface with zero material losses; the edge state is prevented from reflecting off the hard boundary with the vacuum by slowly increasing material losses at these boundaries. We place a right circularly polarised magnetic dipole source with magnetic field $\textbf{H}=H_{x}+iH_{y}$ at the various positions shown in Figure 4a. We choose the excitation frequency from Figure 3b as $\omega=2.83$ eV. The emission of the magnetic dipole source placed in the plasmonic lattice will be modified due to the surrounding environment. This is characterised by the Purcell factor $P_{F}$, the ratio of emitted power of a magnetic dipole to the emitted power in free space $P_{M}$ (58),

[TABLE]

where $\mu_{0}$ is the vacuum permeability and m is the magnetic dipole moment. The Purcell factor as a function of source position is shown in Figure 4b (dotted orange line). The enhancement is greatest when the source is close to the metallic nanoparticles in the lattice, but we note it is generally modest due to the distance between the source and nanoparticles.

To characterise the directionality of the energy flow along the edge, we integrate the Poynting vector through a plane perpendicular to the interface and metasurface, and we plot this as solid lines in Figure 4b; with purple and green corresponding to the fraction of flow to the right ($P_{R}$) and left ($P_{L}$), respectively. Regions in which the coupled dipole approximation does not hold are shaded. Starting at $\textbf{r}_{1}$, and looking at the power flow for the right circular polarisation, the flow of energy is predominantly to the right ($P_{R}$, purple line) as expected from the polarisation of the source, which couples to the right-propagating pseudospin mode. Importantly there is still a fraction of energy travelling to the left ($P_{L}$, green line) which demonstrates the existence of pseudospin mixing. As the source is moved towards $\textbf{r}_{2}$ the flow is completely to the right before quickly flipping to the opposite direction. From $\textbf{r}_{2}$ to $\textbf{r}_{3}$ the energy flow again changes sign as the source moves from an area with negative spin angular momentum to an area with positive spin angular momentum. The propagation direction of the excited edge state can then be predicted from the interplay between the source polarisation and the local handedness of the polarisation of the mode, given by the spin angular momentum in real space, rather than by the pseudospin. More details on the directionality of the modes excited by dipole sources are given in the Supporting Information, Figure S1). We emphasise that we have used a right hand polarised source throughout, which suggests one should always expect energy flow in the right direction.

In Figure 4c-d we show the stark difference in directionality along the edge depending on the position of the source (magenta star). We consider the same right-circularly polarized source placed at two different positions. In Figure 4c we choose the optimal placement for excitation in the expected direction, and in Figure 4d the least optimal placement. All bosonic systems with these lattice symmetries, such as the photonic crystal, will also possess these position dependent directional modes when excited with a circularly polarised source (56). In the left and right panels, we show the electric field intensity, $|\textbf{E}|^{2}$, at plane cut perpendicular to the metasurface in the left and right directions. We note that the edge state is not only confined to the edge in the plane of the metasurface but also out of the plane at subwavelength scales. Finally, we stress that in this discussion we have considered sources placed in the plane of the metasurface. For sources above the nanoparticles, the excitation of directional edge modes will be determined by, first, the value of the Purcell factor, and, second, the interplay between the source polarization and the distribution of spin angular momentum. (Plots of the spin angular momentum in planes above the metasurface are given in the Supporting Information, Figure S2.)

Retardation and Radiative Effects

We finally discuss in detail the effect of retardation and radiation in the pseudospin edge states. While so far we have considered a very subwavelength period and nanoparticle size, we have already seen the effect of retarded interactions which cause the band structure to be altered with respect to the QSA, in particular close to the light line, and the radiative broadening and redshifting of resonances. Radiative effects become more apparent for larger nanoparticles, with very large broadenings and shifts as in the single particle extinction cross section shown in Figure 1d; effects which are not captured in the QSA. On the other hand, it is important to note that retardation can have striking consequences not only on bulk band structures but also on the properties of edge states. It has been shown for example how in the 1D plasmonic SSH model, retardation can result in the breakdown of bulk boundary correspondence and the disappearance of edge states (36). To investigate the effects of retardation in the 2D system studied here, we excite an interface between expanded and shrunken regions in a ribbon with a plane wave at a finite wavevector above the light line, $k_{x}=0.15\pi/a_{0}$, and calculate the extinction cross section for increasing lattice constants $a_{0}$.

In Figure 5a, we show the extinction cross section for a metasurface with an interface between the two phases for silver nanoparticles with radius $5$ nm and height $20$ nm. We let the Drude losses $\gamma=0.01$ eV, as in Figure 3b, and highlight the edge states with white dots for visibility. As in the infinite lattice, in the quasistatic and nearest neighbour approximations, the edge states are expected to be symmetrical about the plasma frequency $\omega_{sp}$ (dotted purple line) (44). In contrast, when radiative effects are taken into account, there are substantial shifts in the edge state frequencies. Initially, for $a_{0}=60$ nm the edge states are well separated from the bulk but as the lattice constant increases up to $90$ nm the bulk modes close up and the edge states are lost. We calculate the full width at half maximum (FWHM) of the highest energy edge state at $a_{0}=60$ nm where FWHM $=0.0059$ eV. As an aside, this edge state has a larger cross section as along the interface the nanoparticles form a bonding-like state (Figure 5c, top) whereas the lower energy edge state has an anti-bonding distribution which leads to lower cross section (Figure 5c, bottom). In Figure 5b we show the cross section for nanoparticles with radius $10$ nm and height $40$ nm, with lattice constant varying from $120$ to $180$ nm. Here, the edge states become significantly redshifted far below $\omega_{sp}$. Again we measure the FWHM, at $a_{0}=120$ nm, FWHM $=0.0077$ eV demonstrating the broadening of the mode for larger nanoparticles. We note that these widths are smaller than the Drude losses since the lattice structure modifies the optical response by increasing the quality factor (17, 18).

Conclusion

In this work we have presented a study of spin dependent edge states in a plasmonic metasurface. These states rely on lattice symmetries and are similar to the quantum spin Hall effect in topological insulators. By going beyond the quasistatic approximation and including retardation and radiative effects, we model the plasmonic system appropriately and show how the long range interactions result in a more complex band structure than the QSA. The bands involved in the band inversion are not greatly affected by retardation, and they provide a signature for far-field measurements. We note that the ordering in energy of the modes of the plasmonic metasurface are opposite to what is seen in photonic crystals (7), with the dipolar modes lying at higher energies than the quadrupolar modes for the non-inverted bands.

Importantly, we determine the spin angular momentum of the edge modes, and show that this is the quantity that characterises the directionality of the modes. Remarkably, this leads to a more complex behaviour which can be opposite to what would be expected from a direct analogy with the QSH effect. We emphasise that these conclusions are valid not only for the plasmonic metasurface but for any bosonic system with this lattice (56). We probe the spin angular momentum by looking at the excitation of the edge states in the near-field by a magnetic dipole source. By varying the position of the source, we highlight how the location is essential in exciting a purely unidirectional state. In agreement with the spin angular momentum characterisation of the edge in the ribbon we show that in some positions a source will excite a mode in completely the opposite way to the expected direction.

Finally, we have also considered the optical response of the plasmonic metasurface and the edge modes under far-field excitation, showing the effect of radiative and retardation effects for the larger radius nanoparticles and larger lattice constants. Although the edge states persist for larger nanoparticles, we observe a radiative broadening and shift in frequency which is naturally not captured in a quasistatic approximation. For increasing lattice constants, the bulk modes close up and eventually the edge states are indistinguishable, which provides a parameter regime for which these edge states could be experimentally observed.

Methods

Lattice sums

The summations in the interaction matrix in Equation 7 are conditionally convergent. Firstly, we will consider sums including the origin term and split these into long range and short/medium range terms,

[TABLE]

where $S_{L}=k^{2}\sum_{\textbf{R}}e^{i\textbf{k}\cdot\textbf{R}}e^{ikd}\frac{1}{d}$ and $S_{SM}=k^{2}\sum_{\textbf{R}}e^{i\textbf{k}\cdot\textbf{R}}e^{ikd}\left(\frac{ik}{d^{2}}-\frac{1}{d^{3}}\right)$. The slowly converging $S_{L}$ term is handled by using Ewald’s method (59). This splits the real space sum into two and then takes the Fourier transform of one part using Poisson’s summation, resulting in a sum over the reciprocal lattice. The sum is optimised with a Ewald parameter to ensure the real space and reciprocal space sum converge within approximately the same number of lattice constants. The $S_{SM}$ term converges rapidly above the light line. Outside of this region the first few terms within a radius $R_{min}$ are added and then the rest of the sum is calculated numerically by approximating the summation as an integral (48),

[TABLE]

For sums excluding the origin,

[TABLE]

the integral method from Ref. 48 is used. For ribbons which are infinite in only one direction, the summations converge easier and so no techniques are used to speed up convergence (41).

In the quasistatic approximation, we only consider the short range $1/d^{3}$ term which converges quickly when including all neighbours, in both the infinite lattice and semi-infinite ribbon. We note that this method of including all neighbours results in a kink at $\Gamma$ in one of the modes of the infinite lattice which is due to the group velocity necessarily being zero at the centre of the BZ (48).

Linearised Green’s Function

In Figure 3, we linearise the Green’s function to investigate the spin angular momentum semi analytically. For the retarded Green’s function and radiative polarisabilities the eigenvalue problem is non-linear and non-Hermitian,

[TABLE]

To avoid the computational complexity of searching for complex $\omega$ solutions we linearise the Green’s function by making the approximation $\omega=\omega_{sp}$, the surface plasmon frequency. In spherical nanoparticles, $\omega_{sp}=\omega/\sqrt{\epsilon_{\infty}+2}$ and for spheroids, $\omega_{sp}=\omega/\sqrt{\epsilon_{\infty}-1+1/L}$. For out of plane modes in the system we consider, $\omega_{sp}=\omega/3.12$. This is valid for the size of nanoparticles considered here since $\omega$ varies faster in the polarisabilty term than in the Green’s function. As particle size increases the approximation becomes less valid close to the light line. We can then rewrite Equation 13 as,

[TABLE]

for which eigenvalues $\lambda$ and eigenvectors p are found at each point in the BZ, k, and band structures are calculated by rearranging $\lambda=1/\alpha(\omega)$ to find $\omega$, for the static polarisability.

Extinction Cross Section

To calculate the extinction cross section, $\sigma_{\text{ext}}$, we use the following system of equations for a non-zero external field,

[TABLE]

From Maxwell’s equations, we have $k_{\parallel}E_{\parallel}+k_{z}E_{z}=0$ and we assume $E_{\parallel}=1$. Rearranging, we can then write,

[TABLE]

The total incident field on each particle includes an additional phase due to the position within the unit cell d, $\textbf{E}_{inc}=\textbf{E}\exp{(i\textbf{k}\cdot\textbf{d})}$. After calculating dipole moments using Equation 14, the extinction cross section is given by the optical theorem,

[TABLE]

Spectral Function

The spectral function method relies on an effective polarisability formulation of the system (48, 41). We rewrite the system of equations for a non-zero external field as $\textbf{p}=\alpha_{\text{eff}}\textbf{E}$. The effective polarisability $\alpha_{\text{eff}}=1/\lambda$ for eigenvalues of M, $\lambda$. The spectral function is analogous to the extinction cross section but rather than describing the system when excited by a well defined external field instead characterises all modes in the system in the retarded, radiative regime, regardless of whether they are bright or dark modes. This corresponds to the forced oscillation of each mode of the lattice at some driving frequency $\omega$ and Bloch wavevector k. The spectral function is defined,

[TABLE]

where the sum is over the number of elements in the unit cell in the infinite lattice or the super cell in the semi-infinite ribbon. Peaks in the spectral function will correspond to the real part of the band structures from the linearised Green’s function.

{acknowledgement}

We acknowledge fruitful discussions with A. García-Etxarri and Mehul P. Makwana. M.P., R.V.C. and P.A.H. acknowledge funding from the Leverhulme Trust. P.A.H. also acknowledges funding from Fundação para a Ciência e a Tecnologia and Instituto de Telecomunicações under project CEECIND/03866/2017. S.A.M. and R.V.C. acknowledge funding from EPSRC Programme Grant “Mathematical Fundamentals of Metamaterials” (EP/L024926/1). S.A.M. additionally acknowledges the Lee-Lucas Chair in Physics. We thank Jude Marcella for the TOC figure.

{suppinfo}

Position dependence on directionality of the edge states, with left and right polarised magnetic dipole sources. Spin angular momentum calculated at planes above the metasurface. This material is available free of charge on the ACS Publications website.