The effects of retardation on the topological plasmonic chain: plasmonic edge states beyond the quasistatic limit

We study a one-dimensional plasmonic system with non-trivial topology: a chain of metallic nanoparticles with alternating spacing, which is the plasmonic analogue to the Su-Schreiffer-Heeger model. We extend previous efforts by including long range hopping with retardation and radiative damping, which leads to a non-Hermitian Hamiltonian with frequency dependence. We calculate band structures numerically and show that topological features such as quantised Zak phase persist due to chiral symmetry. This predicts parameters leading to topologically protected edge modes, which allows for positioning of disorder-robust hotspots at topological interfaces, opening up novel nanophotonics applications.


Introduction
Plasmonic systems take advantage of subwavelength field confinement and the resulting enhancement to create hotspots, with applications in medical diagnostics, sensing and metamaterials [1,2]. Arrays of metallic nanoparticles support surface plasmons that delocalise over the structure, and whose properties can be manipulated by tuning the dimensions of the particles and their spacing [3][4][5][6]. In particular 1D and 2D arrays have significant uses in band-edge lasing [7,8], and can be made to strongly interact with emitters [9,10]. Configurations of nanoparticle dimers have been shown to exhibit interesting physical properties [11]; in the following we consider a nanoparticle dimer array in the context of topological photonics.
The rise of topological insulators (TIs), materials with an insulating bulk and surface states which are protected from disorder, has inspired the study of analogous photonic and plasmonic systems [12][13][14][15][16][17][18][19][20][21][22][23]. Topological photonics shows exciting potential for unidirectional plasmonic waveguides [24], lasing [25], and field enhancing hotspots with robust topological protection, which could prove useful for nanoparticle arrays on flexible substrates [26]. Plasmonic and photonic systems provide a powerful platform to examine TIs without the complication of interacting particles, and with interesting additional properties like non-Hermiticity [27][28][29][30][31][32]. The lack of Fermi level simplifies the excitation of states, and the tunability made available by the larger scale allows for the study of disorder and defects in greater depth than electronic systems [33][34][35]. They also simplify the study of topology in finite systems [36].
One of the simplest topologically non-trivial models is that of Su, Schrieffer and Heeger (SSH) [37,38], which features a chain of atoms with staggered hopping. The one-dimensional chain of metallic nanoparticles with alternating spacing (see figure 1(a)) has been studied in analogue to this by taking the quasistatic (QS) approximation, where the dimensions of the chain are much less than the wavelength kd ! 1 [39][40][41]. In fact, when damping is neglected and only nearest neighbours considered the two systems are physically equivalent. However, the QS limit has been shown to be insufficient for describing the band structures of larger particles and spacing in equally spaced chains [42,43], as it neglects the effects of retardation and radiative losses.
In the following work we present a more complete treatment of the staggered plasmonic chain which takes into account retardation and radiative effects over long range, providing a natural extension to the SSH model that complements those already in the literature [44][45][46]. The model breaks chiral symmetry only trivially by adding an identity term to the Bloch Hamiltonian, which is non-Hermitian with frequency dependence. We calculate band structures and compare to the QS approximation, showing that the system is indeed still topologically non-trivial because it shares eigenvectors with a chiral system. The transverse and longitudinal modes are shown to have notably different band structures, and in the transverse case it is shown that the Zak phase is not the always the same as predicted by the QS approximation. In addition, we compare two methods of calculating the Zak phase and further confirm that for inversion symmetric systems it is possible to apply Zak's original results even in the case of non-Hermiticity. We go on to study the effects of disorder on the topologically protected edge states, which we find to be extremely robust.

Topological plasmonic chain
The plasmonic analogue of the SSH model is a chain of metallic nanoparticles with alternating spacing, as in figure 1(a). Particles have radius a and unit cells length d, with the spacing between the A and B sublattices given by t " βd{2, where β acts as a tuning parameter. If the spacing of the particles is large enough compared to the radius of the spheres (t, d´t ě 3a) the nanospheres can be treated as dipoles with dipole mo-ments p n , and the system is described by the coupled dipole equations: where Gpr nj , ωq is the free space 3ˆ3 Green's dyadic, which depends on the separation of the dipoles, r nj " r n´rj and complex frequency ω.
The properties of the individual nanospheres are represented by the polarizability αpωq, which in the quasistatic approximation is given by, where pωq is the dielectric function of the metal, 0 the permittivity of free space and B the background dielectric. This neglects radiative damping, which is essential for the model to be consistent with the optical theorem [47]. This is addressed by the radiative correction, αpωq " α qs pωq 1´i k 3 6π 0 α qs pωq . (3) Throughout this work we consider gold nanospheres using an extended Drude model with 8 " 9.1, ω P " 1.38ˆ10 16 Hz and γ D " 1.08ˆ10 14 Hz [48], embedded in a material like glass, B " 2.25. Figure 1(b) shows the normalised extinction cross section of a single nanoparticle, where the solid lines are quasistatic and the dashed make use of the radiative correction. While for particles of 5 nm radius the QS approximation agrees with the radiative correction, for particles with radius 20 nm radiative losses strongly effect the extinction cross section by reducing and broadening the resonance over wavelength. Particles with radius as small as 5 nm are well described classically and do not require quantum size effects to be taken into account [49]. The chain is confined to the x-axis, so the longitudinal (x) and transverse (y, z) parts of equation 1 decouple, where ν " x, y, z labels the orientation of the dipoles, and the hopping between dipoles is scalar, G y,z pr nj , ωq "´e with ω dependence contained within the wavenumber k " x`k 2 y`k 2 z . Unlike this work, previous studies have ignored Drude damping (γ D " 0) and taken the QS approximation with nearest neighbour hopping to solve the system [39][40][41], which neglects retardation by assuming that the dimensions of the chain are very small compared to the wavelength kd ! 1. In this limit, where m ν " 2 for the longitudinal case and´1 for transverse. This removes all ω dependence from G ν and neglects the intermediate and long range dipolar interactions. The resulting nearest neighbour, real, staggered hopping provides a close analogue to the SSH model, apart from a transformation from the eigenvalues to the frequency ω.
For gold nanoparticles embedded in glass the non-radiative surface plasmon resonance ω sp " ω P { ? 8`2 B corresponds to the wavelength λ sp " 504 nm. Figure 1(c) shows the dramatic difference between the QS approximation and the retarded treatment for an evenly spaced chain (β " 1) when d " 200 nm, on the same order as λ sp . Green dashed lines show bands resulting from the QS approximation, which are therefore symmetric around ω sp . The blue solid lines show the result of including retardation and radiative effects; as has previously been shown for the evenly spaced chain, retardation leads to polariton splitting and discontinuities at the light lines k "˘k x [42,43], which are completely absent in the QS approximation. The difference is even greater for the transverse polarisation due to the extremely long range " exppikrq{r term in the full dipolar interactions. When including retardation and radiative losses hopping becomes complex and long range, giving rise to a non-Hermitian topologically non-trivial Hamiltonian.

Bulk Bloch Hamiltonian
Topologically non-trivial systems exhibit a bulkboundary correspondence, where properties of the bulk, here the Zak phase [50], predict the existence or absence of edge states in the finite case [51]. We study the bulk by way of an infinite chain, where we relabel the two particles in the unit cell A and B as in figure 1(a), apply Bloch's theorem, and arrive at the equations where G ν pk x , ωq acts as an ω-dependent non-Hermitian Bloch Hamiltonian which is, in matrix form, The eigenvalues of G ν are 1{αpωq, but the band structure is given by ω. Topological properties of the system are associated with the eigenvectors p ν .
As is the case for any two band Hamiltonian, it is possible to write G in terms of the Pauli matrices tσ i u, with g 0 and g " pg x , g y , g z q given by examining G ν . The QS nearest neighbour approximation has strict chiral symmetry, or sublattice symmetry, because there is no hopping from sites A to A or from B to B. This can be expressed by the equation σ zĤ σ z "´Ĥ, which is true when g 0 " 0 " g z . For the retarded treatment we still have g z " 0, but g 0 ‰ 0, which we will call 'trivial' chiral symmetry breaking. Strict chiral symmetry leads to eigenvalues that are symmetric around 0, but here they are symmetric around g 0 pk x , ωq. This system still has inversion symmetry in the x-direction expressed by σ x G ν pk x qσ x " G ν p´k x q, with inversion centres marked by crosses in figure 1(a). This guarantees that the band structure is symmetric in k x . Calculating the band structure is a matter of fixing real k x and finding corresponding complex ω numerically, to solve detˆG ν´1 αpωq I˙" 0.
This is complicated by the fact that the elements of G ν are infinite, slowly converging sums. Faster evaluation can be achieved by writing the sums in terms of polylogarithms and the Lerch transcendental function, detailed in the supporting information (SI). These analytical expressions show that when β " 1, the offdiagonals of G ν are zero at the edge of the Brilluoin zone (BZ), at k x d{2 " π{2. Therefore the eigenvalues are degenerate here, leading to a band crossing as in figure 1(c). This signifies a topological phase transition at β " 1. Figure 2 shows numerically calculated band structures for various choices of chain parameters d and a, displaying only the real part of ω, with some discussion of Impωq in the SI. Blue dashed lines show the β " 1 case, and red solid lines show the β " 0.9 (identical to β " 1.1) case. For β ‰ 1 a complex valued gap opens at the edge of the BZ, which increases in magnitude with increasing |β´1|.
For small chain geometry (d " 50 nm, a " 5 nm) in figure 2(a) the band structure is well approximated by the QS model (yellow line), which makes a reasonable prediction of the band gap but fails to predict the small deviations of the band structure at the light line in the transverse case. Already some asymmetry in Repωq exists due to the trivial breaking of chiral symmetry. Figures 2(b,c) demonstrate band structures with larger particles and spacing with d " λ sp , well away from the QS limit. Once again the polariton splitting and Repωq asymmetry are present, as well as discontinuities at the light line, as in the upper band for (b) longitudinal. The QS approximation, not shown for clarity, is poor here and completely fails to predict that a gap in Repωq does not always open, such as in (b) transverse and (c) longitudinal. It is important to note that in these cases there still exists an gap in Impωq, such that there is no band crossing and associated topological phase transition (see SI). This means these cases can still be topological or trivial.

The Zak phase
The relevent topological number is the Zak phase, which for an Hermitian system like the SSH model is given by However, the Hamiltonian G ν is non-Hermitian, so we must be more careful. The generalisation of the Berry phase for non-Hermitian systems [29], written in 1D for our system, is where p R and p L are normalised biorthogonal right and left eigenvectors, solving equation 8 and its Hermitian conjugate respectively. It has been shown that chiral symmetry quantises this non-Hermitian Zak phase [30]. As discussed previously, our Hamiltonian breaks chiral symmetry trivially due to an additional identity term g 0 . Since all vectors are eigenvectors of the identity, the system shares eigenvectors with a chirally symmetric counterpart (see SI) and the chiral symmetry result quantising the Zak phase applies for this system too.
In fact, the inversion symmetry of the system leads to quantisation of the Hermitian Zak phase as well [52], so that in this case both calculations have the same result, where φ is the relative phase difference between p A and p B . It follows that γ is either 0 or π modulo 2π. It is important to note that topological systems with chiral symmetry but no inversion symmetry do exist, so the topological nature here arises because of the (trivially broken) chiral symmetry. When the Zak phase is γ " π we expect topologically protected edge modes. Figure 3(a) shows how φpk x q changes across the Brillouin zone for the lower bands of figure 2(b). We examine β either side of the topological phase transition at β " 1. The longitudinal case has the same property as the SSH model and QS approximation, that γ " π when β ą 1 and γ " 0 when β ă 1 [39]. Surprisingly, the transverse case is in the opposite topological phase to the longitudinal case for the same choice of unit cell when it has the same β. Of the example band structures given in figure 2, all bands have the same topological properties as the SSH model and (b) longitudinal, except (b) transverse.
Due to the existence of this unusual case and the analytically difficult elements of G ν it is unclear how to predict the Zak phase for a given set of parameters of the chain without performing the full numerics. Despite this, any interface between chains with the same d and a, and with β either side of the topological phase transition β " 1, will still be a topologically non-trivial interface featuring a topological edge mode.
M. Xoai et al. showed that, in photonic systems, the Zak phase is also given by the behaviour of the electric field at the inversion centres x " t{2, pd`tq{2 of the chain, at the centre (k x " 0) and edge (k x " π{d) of the Brillouin zone (BZ) [53]. Considering the inversion centre at x " t{2, if |E t{2 pk x " 0q| and |E t{2 pk x " π{dq| are both either zero or non-zero, we have γ " 0. If |E t{2 pk x " 0q| and |E t{2 pk x " π{dq| are opposite (one is zero while the other is non-zero), the Zak phase is given by γ " π. The normalised magnitude of the electric field at x " t{2 is shown across the BZ in figure 3(b), which agrees with the calculations of figure 3(a). If we take the opposite inversion centre the Zak phases are switched. This method of determining the topological nature of the system relies on Zak's results [50] originally for Hermitian systems, but inversion symmetry assures that any results for the Hermitian Zak phase is identical to the non-Hermitian Zak phase.

Finite chains and disorder
We now consider the implications of the topological phases in finite chains. Figure 4(a) shows the eigenmodes of a finite chain with an example choice of parameters so that there is a gap in Repωq. The gap as defined by the bulk modes (blue) increases symmetrically away from β " 1, where there is a topological phase transition. As expected, edge modes (yellow) appear in the gap when the Zak phase γ " π.
Bulk modes can be identified by their mode profiles, which are typically similar to normal modes of a chain, as in figure 5(b). They can therefore be ordered by assigning a mode number n, the number of times the sign of Repp ν q changes plus 1, and a k x given by where N is the number of particles in the chain [42], and N " 60 in our calculations. These pk x , ωq pairs are plotted for an example set of parameters in 4(b), where the finite QS approximation (red circles) and retarded system (blue dots) are compared. Bulk modes of the finite chain approximate the Bloch bulk band structure in both cases, while the two edge modes exist instead in the gap. Figure 4(c) shows the dipole moments of a set of particles near the edge, with a comparison between the QS edge mode (yellow dashed) and the retarded edge mode (green solid). The QS edge mode is fully supported on only one (the A) sublattice due again to chiral symmetry, while the retarded edge mode exists on both sublattices. This is due to the long range nature of the hopping, forcing the retarded case to be further from the fully dimerised limit than the QS case. This also explains why in the QS case the edge modes have energies fixed to Repωq " ω sp but the retarded edge modes' energies are slightly different. The QS edge mode decays exponentially into the chain whereas the retarded edge mode decreases more slowly into the chain. The real part of p x has a minimum at particle 5 before increasing and then decreasing again, because of the longer range, out of phase, dipole-dipole interactions, but the absolute value |p x | still decays monotonically on each sublattice into the chain as illustrated by the inset log |E| field.
When the gap has no real part, edge mode frequencies have an imaginary part so that they still sit in the imaginary valued gap for γ " π (see SI), and have similar profiles to figure 4(c). For the transverse case the extremely long range interactions " exppikrq{r necessitate a very long chain in order to distinguish between an edge and a bulk, which is numerically challenging.
One of the most relevant properties that arises due to topology is the protection of the edge modes from disorder in the axis of the chain. In figure 5 we apply disorder in the form of a random positive or negative shift to each particle's position in the chain axis, and measure the root mean square of the disorder as a percentage of Λ " |β´1|. When the disorder is 50% the system is within one standard deviation of the topological phase transition, where the particles are equally spaced.
In figure 5(a) a random choice of disorder is scaled smoothly for different choices of β, causing the bulk modes (blue) to enter and eventually close the gap at around 50%. These bulk modes also become Anderson localised, as in any 1D system with disorder, with example mode profiles in figure 5(b). The edge modes (yellow) separate in energy until they are lost in the bulk, but survive for high levels (sometimes ą 20%) of disorder, especially for larger |β´1|. Figure 5(c) shows the mode profiles of the two edge modes for two joined chains with 20% disorder and opposite Zak phases, which illustrates the continued existence of the edge modes in disordered systems. These disorder-protected edge modes act as plasmonic hotspots, which can be positioned anywhere at the interface of two chains with opposite Zak phase.

Conclusion
Here we have presented a detailed study of the 1D topological plasmonic chain beyond the quasistatic limit. We have discussed how appropriately modelling the interaction between the plasmonic nanoparticles by including the effects of retardation and radiative damping, as well as losses in the metal, leads to fundamental differences with its original electronic analogue, the SSH model. In particular, the plasmonic chain has a non-Hermitian Hamiltonian with long range hopping which breaks chiral symmetry in a 'trivial' way. This implies that, because the system has the same eigenvectors as a chirally symmetric counterpart, it is still a topologically non-trivial system that supports edge modes at interfaces between topological phases.
While in the QS limit the behaviour of the SSH model is recovered, specifically band gaps opening both for transverse and longitudinal modes and edge modes confined to a single sublattice in the topological phase, we have shown that as the size of the particles and their separation increases a richer phenomenology appears. In particular, the bulk band structures deviate strongly from the QS prediction, and a band gap opening for the dimerized chain does not always appear in the real frequency axis. However, in these cases there is still a gap in the imaginary part of the frequency so that it is possible to calculate topological invariants and define topological phases. We have calculated the Zak phase and discussed a remarkable case where the topological phase of the chain is opposite to that predicted by the SSH model, suggesting that there is still greater understanding of the topological plasmonic chain to seek.
Finally, we showed that the edge states survive positional disorder in the axis of the chain to a great extent. This hints at potential uses for the 1D topological plasmonic chain in plasmonic systems, which could take advantage of this robustness against fabrication imperfections to design plasmonic hotspots.

Supporting Information Lerch transcendent form
One of the challenges associated with the solving of the band structure of the infinite system is how slowly the sums converge in the Bloch Hamiltonian G ν , given bŷ where the primed sums have n ‰ 0. Even for as many as 1000 elements of the sum calculations can be quite far off. To do this for the case of evenly spaced particles, corresponding here to the on-diagonal sums, Citrin [54] and Koenderink and Polman [43] used polylogarithms defined by the sum, Li s pzq :" These have the advtange of being implemented as standard in several scientific packages. It has been shown that we can write the on-diagonal terms in polylog form, and The off diagonal terms, with x n " nd˘t, no longer fit the form of a polylogarithm. We resort to a more general special function, the Lerch transcendent, Φpz, s, νq :" For our numerical calculations we make use of the fact that the Lerch transcendent can be expressed as an integral which converges much faster than the sum, Φpz, s, νq " 1 Γpsq where Repνq ą 0, and either |z| ď 1, z ‰ 1, Repsq ą 0 or z " 1, Repsq ą 1 [55]. This imposes the restriction that it's not possible to calculate the Lerch transcendent on the light lines for the transverse case, as can be seen from the following equations. The off diagonal elements in terms of the Lerch transcendent are and Substituting k x "˘π{d and t " d{2 we see from the above equations that which predicts the closing of the gap at k x "˘π{d, t " d{2 since g x " g y " 0 in equation 26.

Trivial chiral symmetry breaking
As is the case for any 2ˆ2 matrix, it is possible to write G in the form where σ is the vector of Pauli spin matrices, and for this system in particular Chiral symmetry or sublattice symmetry exists when there are no interactions between sites A to A or B to B, which can be expressed in the form σ zĤ σ z "´Ĥ [38]. This is the case when g 0 " 0 " g z in equation 25. In this model we break this chiral symmetry in a 'trivial' way.
Interactions between A to A and B to B are introduced, but the A and B sublattices remain indistinguishable, because the difference between the retarded plasmonic system and a chirally symmetric system is an identity term g 0 ‰ 0. This is important because it means that the system has the same eigenvectors as a chirally symmetry system with Bloch Hamiltonian G´g 0 I, Since the eigenvectors are the same we can apply Zak phase results about chirally symmetric systems to our trivially non-chiral system. We also have a pseudochiral equation for G from the original chiral equation, σ z pG´g 0 Iqσ z "´pG´g 0 Iq. This explains why the eigenvalues are not symmetric around 1{α " 0, and goes one step further: so that much like the chirally symmetric case every eigenmode p has a counterpart σ z p, with eigenvalues related by These equations relate the upper and lower bands of G, although an additional non-trivial mapping is required from these eigenvalues to ω. Notably G also has ω dependence, which is why the eigenvalues aren't perfectly symmetric around g 0 . This can also be used to explain the fact that the edge modes of our system do not appear to be supported on only one of the sublattices, which is the case for the SSH model. Here the sublattice projection operators are given byP A " pI`σ z q{2 andP A " pI´σ z q{2. In the case of the fully dimerised SSH model edge modes |ψ n y are zero energy E n " 0, so that HP A{B |ψ n y "Ĥ p|ψ n y˘σ z |ψ n yq " 0. (30) An equivalent equation for our model would be ClearlyP A{B p is only an eigenvector of G when p and σ z p have the same eigenvalues-which we would expect in a 'fully dimerised limit', except that such a thing doesn't exist due to the long range dipole-dipole interactions of the model. This suggests that the nature of the hopping makes it difficult to see a 'purely' chiral behaviour, and we can expect stronger confinement to one sublattice when |β´1| is larger. It is possible to break chiral symmetry in a non-trivial way by adding a term in σ z to the Bloch Hamilotnian, for example by changing the relative sizes of A and B particles, or making A and B different metals.
Complex valued frequencies Figure 6 shows the Bloch band structure (red) in the case where there is a crossing in the real frequency axis, and the gap still has an imaginary part. In this case it is still possible to calculate a Zak phase and, for a finite chain, find bulk modes (blue) which follow the bands and edge modes (yellow) which sit in the gap. It is important to note that all of the retarded and radiative band structures presented in this work have small imaginary components to the frequency. The presence of such imaginary part implies that it may be necessary to use an evanescent wave to excite the modes in practice, which could be acheived with, for example, the usual Otto or Kretschmann configurations [56,57]. Figure 6: Bulk and finite 60 particle chain band structure for d " 330 nm and a " 20 nm, β " 1.1. This is in the γ " π topological phase, leading to edge modes. (a) The band structure projected onto the Repωq axis, with finite bulk modes (blue) following the Bloch (red) bands, and edge modes with difference frequency. (b) The imaginary part of the band structure, which shows that there still exists a gap in which the edge modes sit (c) A 3D visualisation of the band structure. Discontinuities in the bulk band structure exist at the light line.