Topological Surface-plasmon-polaritons on Corrugated Metal-dielectric Surfaces

We study topological surface-plasmon-polaritons at optical frequencies in diffraction gratings formed by bipartite corrugated metal-dielectric gratings. To do so we implement the theory as developed in Della Valle and Longhi to describe the amplitude of the field by an emergent Schr\"odinger-like equation. The tri-harmonic grating generates a bipartite Kronig-Penney model. Topologically protected localised modes are then predicted to occur at the edges of the grating and at defects formed by the combination of two mirror antisymmetric corrugations, whose bulk invariant is a step-wise varying Zak phase in both cases.

The standard approach in all of these systems is to assume the tight-binding formalism.Then the constructed low-energy effective Hamiltonian, which generates the Schrödinger equation, may possess certain discrete symmetries.When one such symmetry, the chiral one, is present in one-dimension then topologically protected edge modes will be present.The appeal of such modes is both theoretical, for their fundamental interest due to the phenomenon of bulk-boundary correspondence, and experimental, for their robustness against bulk defects and disorder due to the topological invariance of the Hamiltonian through chiral-symmetry-preserving adiabatic deformations.
In the context of topological plasmonics [30][31][32][33][34][35][36] , the tightbinding approximation may be made if the lattice in question is built of nanoparticles that are capable of hosting localised-plasmon-polaritons (LPPs).These modes form the basis 'atomic wavefunctions' of the tight-binding model and the hopping parameters are simply the tunnelling amplitudes that arise from the dipole interactions of these LPPs between neighbouring nanoparticles.
Surface-plasmon-polaritons (SPPs), much like LPPs, are quasi-particles formed from the electromagnetic interaction of photons with plasmons.As such, they may only exist at a dielectric-metal interface whereat the sign of the dielectric function changes; a fact which has been recently shown to be of topological origin 37 .Work has been conducted on the appearance of SPPs at the surface of 3D topological insulators [38][39][40][41] whilst study into their own potential topological characteristics has begun in ernest [42][43][44] .Herein, we consider an extended corrugated surface that is made to be bipartite in nature through a periodic variation of the radii of curvature between neighbouring peaks and troughs.The 'bipartite-ness' here is constructed to effectively mimic the SSH model by keeping all peaks the same and varying the troughs between them.Such systems may be fabricated easily and with great accuracy using modern experimental techniques 45,46 .
By taking the Fourier transform of a typical bipartite grating, as shown in Fig. 2 dominant contributions to the periodic structure are the lowest three Fourier frequencies: k, 2k, 3k, with weights A 1,2,3 and the constant A 0 , which is not crucial.Therefore, constructing the profile g(x) = 3 n=0 A n cos(nkx) and plotting it as in Fig. 2(b) reveals that the bipartite structure is born from a triharmonic pattern.Note that panel (b) is not to scale, the widths of the corrugation are in fact much larger than the heights.
Within the first section, the emergent Kronig-Penney model is introduced.Its derivation is given in the supplementary material.The construction of the effective Schrödinger equation satisfied by the amplitude of the electromagnetic field follows Ref. 1 .The bulk system is then solved in the following section.There, the bulk invariant is identified.Then the finite systems are solved and degenerate edge modes are observed in concordance with the bulk invariant.Finally, conclusions are drawn with respect to the results as presented.

II. THE EMERGENT KRONIG-PENNEY MODEL
As detailed within the supplementary material and in Ref. 1 , the amplitude, F (σ), of the SPP fields, {E, B} ∝ F (σ)e [i(ne+∆n)z−γη]/λ , may be described through the following Schrödinger-like equation: in the presence of a corrugated metal-dielectric interface, where (σ, η) are curvilinear coordinates that are local to the surface as in Fig. 1(a).In addition, λ = λ(2π) −1 with λ as the wavelength of the incident light (ω i = cλ −1 ) that excites the SPP, n e is the effective index of the SPP and ∆n is a correction to this index induced by the perturbing 'geometric' potential V (σ) with: Here, 1 is the dielectric constant of the dielectric above the metallic surface and 2 = (n + iκ) 2 is the dielectric function of the metal, where n and κ are the refractive index and extinction coefficients of the metal substrate at the given frequency ω i .
For the standard time-independent Schrödinger equation (TISE) the eigenvalue is simply the energy of the state.Here, however, the energy of any SPP excitations is fixed at ω i by the incident light.Instead, ∆n acts to alter the effective wavelengths of the SPPs.This is seen in the z-direction Fourier component of the fields: e i(ne+∆n)z/λ , where e inez/λ is the Fourier component in the absence of the grating.In the same way that a Schrödinger particle has a Fourier time component of e iEt/ and as such an energy-eigenvalue TISE, the current system has a Fourier 'time' component of e i∆nz/λ and so an 'energy' of ∆n.
As such, the grating acts to create a non-zero ∆n and so modifies the effective wavelengths of the excited SPPs as λ SPP = λ(n e + ∆n) −1 .To be clear, a photon of wavelength λ that interacts with the metal-dielectric grating will do so to generate an SPP excitation of effective wavelength λ(n e + ∆n) −1 .Then, regardless of the sign of n e , ∆n is a correction to this refractive index and so is not restricted to take only positive values, as is the case when considering Schrödinger particles.∆n is free to both increase or decrease the effective index of the resultant SPP mode.As a result, the band-spectra as presented herein, which are in terms of this correction ∆n, are effectively dispersion spectra for the wavelength λ SPP of the excited modes.
In Della Valle and Longhi 1 , the potentials of the peaks and troughs have the same magnitude since the magnitude of the local radius of curvature, R(σ) = (∂ σ ϕ) −1 , never changes; only its sign.We instead consider a system in which the local radius of curvature varies between peaks and troughs in not only sign but also magnitude.
To mimic the SSH model, we take the radius of curvature of the peaks to be R δ = −a whilst those of the troughs are taken to be R v = a − t and R w = a + t alternately and periodically where t is a parameter that is free to vary within the range 1 − a ≤ t ≤ a − 1.In this way, by ultimately modifying a single parameter t, the length of the unit-cell is maintained to be constant.

III. THE PRIMITIVE UNIT-CELL SOLUTION
The solution within the bulk then proceeds through the standard scattering formalism that is used to solve the Kronig-Penney model 47 .The details of this calculation for the present system may be found in the supplementary material.The important result is that the simultaneous equations generated may be reduced to the following 2x2 unit-eigenvalue matrix equation: where the matrix is denoted S(k) and the matrix elements are given in the supplementary material.The transcendental equation that defines the energy bands may be found through det[S(k) − 2 ] = 0, which is solved numerically.Since the present system has been reduced to an emergent Kronig-Penney model that is governed by the Schrödinger equation, the topological information of the system resides in the Zak (one-dimensional Berry) phase of the wavefunction 15,[48][49][50][51] , which in the present system is in fact the amplitude F (σ): where u m k (σ) = e −ikσ F k (σ) is the unit-cell periodic amplitude, the summation is over all 'occupied' bands and the inner product signifies to integrate over the σ coordinate within the unit-cell.Since bosons are considered here, the meaning of 'occupied' is to say: all bands below the band gap.
In the finite system to be considered, two separate cases involving different physical parameters will be introduced.This is to best represent the topological edge/defect states.In both cases: λ = 0.8 µm, a = 8 µm, and the aperture angle is θ = 157°.
In one of the cases, the metal-dielectric is Ag-SiO 2 and so we take 1 = 3.9 whilst the metal substrate is assumed to be pristine, which is valid since it is covered by the silicon-based dielectric thereby eliminating oxidisation effects.Using data from Jiang et al. 52 we then have that: 2 = −27.6 + 0.919i.
In the other case of Au-air, we take 1 = 1 and the metal substrate is taken to be pristine also.Since gold does not readily oxidise, unlike silver, there is no issue with the dielectric above the substrate being air.Using data from Olmon et al. 53 As discussed in Ref. 1 , due to the fact that |Re( 2 )| |Im( 2 )| in both of these above cases, the imaginary part of 2 does not enter crucially into the asymptotic analysis and only appears, at leading order, as a term of the form e −ξ|z| with ξ ∈ R. As such, it does not affect the profile, F (σ), nor the behaviour of the SPP in the σ coordinate.
Since the asymptotic analysis used to derive Eq. ( 1) is based entirely upon the condition that |R| λ 1 , the present extension respects this.However, the aperture angles θ of both the convex peaks and concave troughs must be made to be identical to ensure that the angle ϕ(σ) does not vary discontinuously.
Plots of the bulk bands may be seen in Figs.3(a,b).Their similarity with those of the SSH model can be clearly seen.Indeed, the band gap closes at the transition point of t = 0. Taking the lower bulk band of either case in Figs.3(a,b), we may calculate the numerical Zak phase and find that θ Z = π for t < 0 and θ Z = 0 for t > 0. This step-wise change in the Zak phase invariant is also observed in the winding of the coefficient r(k).Those of r(k) for the lower band before and after the transition at t = 0 are shown in Figs.3(c,d).Thus the invariant that conforms with θ Z truly is then W = |1 − W r | where W r is the winding number of r(k) 41 .
We stress that the bands shown in Figs.3(a,b) have a different nature with respect to those normally encountered in non-bipartite metal gratings.In those systems, usually a pair of low-lying bands are observed which bare resemblance to bonding and anti-bonding orbitals.This is clearly seen by analyzing the eigenstates corresponding to energies around the gap.In non-bipartite gratings, these are centred at different points of the unit cell, namely the (only) peak and the corresponding trough 54 .This allows one to identify such bands as borne out of the linear superposition of both bound and excited states of the square potential wells generated by the grating.
On the contrary, in the bipartite case the eigenstates around the gap have equal weight on the two peaks in the unit cell.The bipartite bands of Figs.3(a,b) are therefore the result of the folding of the lowest band of a non-bipartite grating and are thus linear superpositions of only the bound states of square potential wells.A gap is opened by the unequal distances between peaks inside the cell and with peaks of neighbouring cells.When these become equal (i.e. at t = 0), the grating is in practice not bipartite anymore, and therefore the gap must vanish, as discussed above.

IV. THE FINITE AND EXTENDED UNIT-CELL SOLUTIONS
However, to find exponentially localised modes and/or observe any bulk-boundary correspondence, the finite system must be solved.The standard procedure is to terminate the edges of the lattice with hard walls such that the states are reflected back into the lattice at the edges.Such a hard termination may only be achieved perfectly with infinite-strength Dirac-delta potentials.
Physically, this may be approximated if the chain terminates at both ends with a trough that has both a large radius of curvature and a different dielectric deposited atop it, of constant 3 where 3 + 2 ∼ 0. This has the effect of making this final barrier (in the K-P model) have both a large width, as a result of a large R, and a tall height, as a result of 3 + 2 ∼ 0. Silicon could constitute such a material since its dielectric constant has a value of 3 ∼ 12.The physical effect at play here is that the SPP, which is excited at ω i < ω p (1 + 1 ) −1/2 , encounters the edge dielectric whereat ω i ω p (1 + 3 ) −1/2 and so has its wavevector 'shifted' along its dispersion.As a result, its group velocity tends to zero and thus it does not propagate meaningfully into the extended region.
Since, in the K-P model, the tunnelling probability of such a final barrier is many times smaller than those of the barriers of the bulk chain, it may be safely assumed that such a barrier may be approximated to have infinite strength.In a physical sense, the shear size of this final barrier acts to forbid the existence of SPPs within them and so it confines the SPPs within the region of interest.Theoretically, this may be realised effectively using hard walls.Then, due to the proximity of the final peaks to this effective infinitely positive potential, no states are permitted to reside therein.As a result, a system with M peaks will permit M − 2 states.
On the other hand, the opposite of the hard-wall would be an open boundary.In this case, the potential beyond the final well of the Kronig-Penney model ought to be zero, which is achieved by a flat surface in the physical system.Since a flat surface will always sustain SPPs 37 , such a boundary would cause the SPPs to leak out of the chain.Moreover, the abruptness of the change from the corrugation to the flat surface would no doubt introduce significant boundary effects at the resultant discontinuity.
As such, open boundary conditions, like those used in typical tight-binding models, are simply not applicable here.Finally, the finite barrier with ε 3 could be used explicitly rather than through an infinite one.However this introduces the Shockley state 55 , a separate entity, needlessly complicating the band spectra.
In Figs.4(a,b) are shown the band spectra for this very case involving hard walls.As may be seen, midgap and degenerate edge modes may be seen for t < 0 in conformity to the value of θ Z within the bulk when t < 0.Moreover, as may be seen in the behaviour of F (σ) for these degenerate edge states in Figs.4(c,d), they are confined to a single sublattice only.As such, the edge states are topological for both the reasons of degeneracy and single-sublattice-confinement 15 .
As such, since the states: (i) are degenerate, (ii) are confined to a single sublattice only, and (iii) conform with the bulk invariant θ Z , then a bulk-boundary correspondence is established thereby protecting the states against chiral-symmetry-preserving lattice perturbations.This is true without an a priori knowledge of the presence of chiral symmetry within the system.Upon the calculation, its presence may be inferred from the character of the edge states.It is possible to construct a tight-binding model for such a system in order to confirm this and is done so elsewhere 56 .
A further way to observe exponentially localised modes is to consider a defect within the chain that is sufficiently far away from any other defect and/or the edges.As such, we consider an extended unit cell that contains a single defect formed through the situation of two troughs of the same R next to one another.Then, by taking a large enough number of peaks along with Bloch's theorem for the entire unit-cell, a localised mode may be seen to exist with an exponentially decaying profile at the defect.
The band spectra of ∆n against t for such an extended unit cell may be seen in Figs.5(a,b).As may be observed, there exist mid-gap, and therefore exponentially localised, defect states.Within the simple SSH model, there is always a single mid-gap topological defect mode within the entire phase space guaranteed by the step-wise change in θ Z from π to 0 between the two chains 15 .
However, here an interesting and highly non-trivial feature is that, as t becomes more negative, the defect band tends towards the bulk bands.On the other hand, as t becomes more positive, another mode appears within the bulk gap in the Ag-SiO 2 system.This extra state is indeed defect-localised however it is non-topological since it does not appear at the topological transition at t = 0.Moreover, in this case, the defect state within the region t < 0 approaches and joins the bulk much more rapidly compared to the Au-air case.Finally, there is a defect state in the Ag-SiO 2 case above the upper band for t < 0. However, since it is not within the gap, it cannot be a topological state and this is shown (if plotted) by its weight upon both sublattices.
All these observations are to say that the defect is capable of being constructed in such a way that its localisation character is affected in a topologically trivial way.The local potential can become large enough, relative to the bulk, to forbid the state even though the system is topologically non-trivial.
In both systems, when t < 0, the potentials that neighbour the defect, V 3 ∝ (a − t) −1 , are reduced from their value at t = 0 whereas when t > 0 the potentials are greater.As such, the localisation of the defect is reduced when t < 0 and enhanced when t > 0.Then, the state within t < 0 'leaks out' of the defect to become a bulk mode whilst another state is allowed to exist within t > 0, which are respectively due to the decreased and increased apparent depths of the defect well in either case.This is an unfortunate state of affairs that may not be avoided with such a defect.
This indicates, rather nicely, a condition in which a symmetry-protected state is forbidden.The difference between the Ag-SiO 2 and Au-air systems is in the depths and heights of the wells and barriers of their respective Kronig-Penney models arising from their values for V 0,1,3 ∝ [−( 1 + 2 )] 1/2 .In the case of Au-air, the changing potential environment of the defect is less pronounced than in the Ag-SiO 2 case due to the smaller potential magnitudes in the former case compared to the latter case.Within a prototypical tight-binding model of such a system 56 these modified potential environments act to generate an on-site potential at the defect that is greater in magnitude to that of the bulk.Therefore the states are energetically forbidden from residing there and become bulk modes as a result.
If, instead, the defect is formed through a dielectric deposition technique akin to that which generates hard walls in the first considered case, then identical topological edge states that were observed there would exist here.They would be unlike standard defect states, however, as they would localise upon either side of the large defect trough in identical fashion to how the edge states in the first case localise to one side of the edge and not at all within the region defined by the extra deposited dielectric 3 .In analogy to the SSH model as applied to polyacetalene, this would be akin to forming a defect from a carbon-carbon triple bond.However, and crucially, the edge states, which arise when hard wall boundary conditions are imposed, are degenerate and have weights on a single sublattice only in both cases.As such, they are still symmetry-protected.It is only the defect mode in the Ag-SiO 2 case that loses its protection within a region of the phase space due to a topologically trivial boundary effect wherein the potential at the defect is greater than that within the bulk.Yet, within the other region of the phase space it still exhibits single-sublattice-confinement.

V. DISCUSSION
Within a tight-binding calculation, there are on-site potentials at each lattice site (peak) that transpire to be identical within the bulk 56 .Their strength not only depends on the depth of the potential well formed from the peak but also on the neighbouring potential barriers formed by the troughs.
In the case of the hard wall, the peak that precedes the infinite potential is neighboured to the right(left) by a trough of |R| = a + t and said infinite potential to the left(right).As such it has an on-site energy that is different, but crucially lesser, than in the bulk where each peak is neighboured by troughs with |R| = a − t and |R| = a + t.When t → −6 µm then this on-site energy starts to affect the edge state.Yet it remains degenerate, as seen in the bands, and confined to a single sublattice in the bulk of the chain away from the edges.It is only its edge behaviour that is affected.This is a testament to the robustness of the state.However, in the case of the defect state, it is neighboured by |R| = a + t on both sides and so has a different on-site potential than in the bulk.When t < 0 this potential is greater and transpires to force the state away from the defect and into the bulk.When t > 0 it is lesser so causes an increased localisation of the state to the defect and, moreover, the admittance of a further non-topological state when t → +6 µm.
The present theory is based upon an approximation: that of |R| λ.This holds for the parameters as used since we have that |R| is at most 14 µm and at least 2 µm, both of which are much larger than λ = 800 nm(2π) −1 ≈ 0.1 µm.In Ref. 1 , a full numerical calculation was conducted, which found that the approximation holds well for both 1 µm ≤ a ≤ 20 µm and 350 nm ≤ λ ≤ 1000 nm.

VI. CONCLUSION
Topological symmetry-protected localised SPP modes have been predicted to exist upon biharmonic (triharmonic?)metal-dielectric gratings by solving an emergent Kronig-Penney model found from Maxwell's equations.The necessary bulk-boundary correspondence was established between the Zak phase bulk invariant θ Z and the number of edge modes.The presence of chiral symmetry was then inferred from the solution thereby showing the states to be topologically protected.
However, interesting boundary/edge effects were observed that destroyed the protection of the defect states thereby indicating one such limit in which the bulk-boundary correspondence breaks down in one-dimension through a symmetry-breaking procedure, which has not been artificially introduced but instead arises naturally from the physical system under consideration.
Given the nature of the derivation of the Schrödingerlike equation through the asymptotic expansion and its applicability to such corrugated surfaces as presented, an extension of the work could be to consider a Moiré superlattice 57,58 .In this situation, the heights of the peaks and troughs vary in size through some sinusoidal envelope where the present approximation that |R| λ would remain applicable.In such systems, localised states may be observed and, due to their similarity to the present system and the applicability of the current model therein, such states may also be topological in nature.

VII. SUPPLEMENTARY MATERIAL A. The Generic Scattering Solution
The details of the bulk unit cell, as shown in Fig. 1, may be described.Firstly, it is constructed to be centred upon the trough that has R = a − t, thereby making the unit cell centro-symmetric, and so: Since the distances between the interfaces in σ are the arc-lengths of the curved surfaces between the edges of the peaks/troughs, the widths of the wells (peaks) in σ are both δ = a(π − θ) whilst the widths of the barriers (troughs) in σ are w = (a + t)(π − θ) and v = (a − t)(π − θ).Now we exactly solve the presented system by requiring that the amplitude F (σ) and its first derivative with respect to σ be continuous at the potential interfaces.The standard solution of the TISE is that the amplitude (wavefunction) is a superposition of right and left moving waves within each well and barrier.So, explicitly: , where N is the number of regions within the unit cell, which is five in our bulk unit cell, k is the Bloch wavevector that is defined through The amplitudes within each distinct region are given by: f j,k (σ) = C j e iqj (k)σ + D j e −iqj (k)σ , where: q j (k) = λ −1 [2n e (−∆n(k) − V j )] 1/2 , and we have that q 2 = q 4 ≡ q 0 and q 5 ≡ q 1 .The dependence of the wavevectors q j (k) upon the Bloch wavevector k is dropped from here on.
Due to the intricate relation of ∆n to k through the wavevectors q j , this may only be solved numerically.

B. Derivation of the Schrödinger-like Equation
The covariant formulation of Maxwell's equations in a general curved spacetime in the presence of dielectrics and the absence of free charges and currents is given by: where F µν is the contravariant electromagnetic field tensor and D µν is the covariant electromagnetic displacement tensor.In general, the Bianchi identity is in terms of covariant derivatives (; γ = ∇ γ ), however due to the symmetry and anti-symmetry of the electromagnetic field tensor and the Christofel symbols in their lower indices respectively, the covariant derivatives reduce to standard Minkowski derivatives, and so: On the other hand, the electromagnetic displacement tensor is defined in terms of the field tensor, the metric tensor, and the polarisation tensor as: where the polarisation tensor has been absorbed into F αβ to form Fαβ as a modified electromagnetic tensor with dielectric prefactors before the electric fields.As a result of this Lorentz-covariant formulation, Maxwell's equations may be recast in a completely covariant way in any coordinate system we choose so long as we know the metric tensor of said coordinate system.In the case of the corrugated surface considered herein, we introduce the set of curvilinear coordinates that are defined locally to the surface of (σ, η, z), where σ is the direction along the surface, i.e. the arc-length, and η is the coordinate perpendicular to the surface.At each point in σ along the surface, the curvilinear coordinates are related to the Cartesian coordinates by: where ϕ(σ) is shown in Fig. 1(a), is defined through e x • e σ = cos[ϕ(σ)], and depends explicitly upon σ.This dependence will be dropped from now on for brevity.The metric tensor, g µν , in a general spacetime satisfies: where ds is the infinitesimal spacetime arc-length and dx µ = (cdt, dx, dy, dz) is the infinitesimal 4-position, with a metric equivalent to the Minkowski metric (in the astrophysical convention) of η µν = diag(−1, 1, 1, 1), thus defining the arc-length in Minkowski space of: (ds Using the relationship between (x, y) and (σ, η) as in Eq. ( 16), we may see that: Thus, we see that: (ds This may be seen to simplify to: and in doing so observe that (σ, η) are indeed orthogonal to one another as they ought to be.Thus, the metric of the curved surface is g µν = diag(−1, h 2 1 , 1, 1) where: and R = (∂ σ ϕ) −1 is the local radius of curvature of the surface, which may be seen to follow from the definition of the mean curvature of the surface R −1 ≡ 2κ = −∇ • n where n is the unit normal of the surface.In the present case, this is none another than e η = − sin ϕe x + cos ϕe y and so: Through the chain rule: ∂ x,y ϕ = (∂ σ ϕ)(∂ x,y σ) and from the above matrix equation for x, y: ∂ x σ = cos ϕ and ∂ y σ = sin ϕ.Thus: Finally, the determinant is g = |g µν | = −h 2 1 and the contravariant metric is g µν = diag(−1, h −2 1 , 1, 1).Now, taking the equations as in Eq. ( 13) with this metric and: (since the covariant formulation allows for the simple renaming (x, y, z) → (σ, η, z)), where is the dielectric function that is free to vary in space, it may be shown through unilluminating algebra that: where: Gauss' laws in Eq. ( 25) do not help to determine the SPP fields from the incident fields.Rather, we consider Ampère's and Maxwell's equations Eq. (26).Evaluating the curls, rearranging the resultant expressions, assuming invariance in time, i.e. {E, B} ∝ e iωt , and measuring all the length scales within the problem in units of λ = λ/(2π) yields: where: At this point, the asymptotic expansion is implemented.Assuming that the surface is sufficiently smooth based upon the assumption that R λ, we make an expansion in a smallness parameter α such that R goes as α 2 .Since R = (∂ σ ϕ) −1 , ϕ is linear in σ.Thus, we expand σ in α simply as: σ = ασ 1 .Then, it is clear that R ∼ α −2 since we also now have that ∂ σ = α∂ σ1 .The fields are expanded in a like manner as u = u (0) + αu (1) however their amplitudes are also allowed to retain corrections in α to avoid secular growing terms within the analysis.As a result, we write: where Z(σ 1 , z) = Z 0 (σ 1 ) + αZ 1 (σ 1 , z) + α 2 Z 2 (σ 1 , z) + • • • with Z ∈ {F, G, H} are allowed to vary with z and σ 1 but we also have that q j = q (0) j + αq (1) j + α 2 q (2) j + O(α 3 ) since, at all orders, we must satisfy the interface boundary condition, n × (E 2 − E 1 ) = 0, at η = 0 whereby the dielectric jumps between 1 and 2 .
Taking the differential equations up to and including α 2 we see that for u: whilst for v we have that: where the fact that R = [α 2 ∂ σ1 ϕ(σ 1 )] −1 has been used.Thus at order α 0 , it is obvious to see that: in which case we arrive at: The solution to which is simply a = 1, b = p −1 j and c = −q (0) j −1 j .Thus, we see that: where the correct signs of q j must be chosen to ensure the correct behaviour of the function at infinity when in either of η > 0, wherein j = 1, or η < 0, wherein j = 2.We find p and q by making use of the interface condition upon the electromagnetic field at the η = 0 to see that q 1 .Thus, we find p and q (0) j through: as: Before continuing, the solvability condition for a differential equation must be introduced.This is the universal statement that, given an entirely general differential operator L x acting over a vector space x that describes the following differential equation: a solution may exist if and only if: where w L 0 (x) is the left-eigenvector that solves the homogeneous differential equation: and the inner product signifies to take an integral over the complete vector space x.At first order in α, we see that we have: which reduces to: Now, taking this equation to separate orders in α we see: and: and at order α which have the solutions , where the choice e = 0 is such that n • (B 2 − B 1 ) = 0 holds, and so the first correction to v is: Note that the determinant of the matrix on the left-hand side in the above is equal to zero.Thus, the solution is not unique; we could equally well have chosen any of the other components to be zero.However, we would have found a consistent result in either case.This may be seen by considering Gauss' law for magnetism ∇ • B = 0 to find that, given the above expansion, ∂ η B η = 0 showing that our choice was consistent.reduces to the following: (2) This equation must be satisfied when both η > 0 and η < 0 thus it follows that: 1 + iR −1 = q (0) 2 2q (2) thus we have two simultaneous equations for q (2) 1 and q (2) 2 , which we may solve to find that: or, in a simpler form: We now make use of the solvability condition again that must apply to this differential equation.Firstly, we find the left-eigenvector of: where w L 0 (x) = (r, s, t), as the following simultaneous equations: r + ps − q (0) j t = 0, −pr − j s = 0, q (0) j r − j t = 0, (69) whose solution is: r = 1, s = −p −1 j , t = q (0) j −1 j .Thus, we now make use of the solvability condition as: Thus, we find the following differential equation linking the amplitudes F 0 , F 2 : (2) which is the self-same differential equation as found from r − ps + q (0) j t = A j .Substituting our expression for q (2) j from the above into this equation yields: At this point we terminate the asymptotic expansion since we have arrived at the order that corresponds to the radius of curvature R −1 ∼ α 2 and so drop the expansion indices.Furthermore, we identify that p = n e is the effective index of the SPP and that q j = iγ j since 2 j /( 1 + 2 ) will always be a negative number at the allowed frequencies of SPP resonance.Finally, we reintroduce the dimensionality within the problem in terms of λ and so find: where: V (σ) = λn e 2R(σ) and R, being proportional ∂ σ ϕ(σ), varies step-wise along the surface even when the radii of curvature of the peaks and troughs varies in magnitude so long as the aperture angle is the same for all.
Then, if we take the SPP that 'forward' scatterers with a positive effective index of p = +n e , and assume translational invariance in z such that we may say F (σ, z) = F (σ)e i∆nz/λ , then we see: Conversely, if we choose the SPP that 'backward' scatterers with negative index p = −n e then we say that F (σ, z) = F (σ)e −i∆nz/λ and so: In other words, both values for p are described by the same Schrödinger-like equation for negative eigenvalues.The corrected index is then n = n e + ∆n in either case whose sign simply indicates the propagation direction along z.

FIG. 2 .
FIG. 2. (Colour online) Panel (a): The Fourier transform of a typical grating profile shows three dominant harmonic modes in its construction.Panel (b): (Not to scale) A plot of the function g(x) reveals that the bipartite corrugated structure is born from an underlying triharmonic pattern.

FIG. 3 .
FIG. 3. Panels (a,b): The band spectra for the bulk modes of the emergent Kronig-Penney model for the two cases with t = 2 µm.Panel (a) Ag-SiO 2 , panel (b), Au-air.Panels (b,c): The parametric windings of r(k) = rx(k)+iry(k) and the Zak phases for the lowest band of either case before and after the transition.Panel (c): t = −1 µm, panel (d): t = 1 µm.
FIG. 4. (Colour on-line) Panels (a,b): The band spectra of ∆n vs t for a chain of 44 combined peaks and troughs terminating with hard walls.Panel (a): Ag-SiO 2 , panel (b): Au-air.Although not clear, there are 20 bands in the plots.Panel (c): the 10th state for either system at t = −3 µm.Panel (d): the real-space magnetic field amplitude of this very mode.

FIG. 5 .
FIG. 5. (Colour on-line) Panels (a,b): The band spectra of ∆n vs t for a chain of 82 combined peaks and troughs within the extended defect-including unit-cell.Panel (a): Ag-SiO 2 , panel (b): Au-air.Panels (c,e): the 21st and 22nd states, respectively, of the Ag-SiO 2 system at t = 3 µm.Panels (d,f): the 21st state of the Au-air system for t = −3 µm and t = 3 µm, respectively.