## Abstract

In this paper we demonstrate that rigorous high-order perturbation of surfaces (HOPS) methods coupled with analytic continuation mechanisms are particularly well-suited for the assessment and design of nanoscale devices (e.g., biosensors) that operate based on surface plasmon resonances generated through the interaction of light with a periodic (metallic) grating. In this connection we explain that the characteristics of the latter are perfectly aligned with the optimal domain of applicability of HOPS schemes, as these procedures can be shown to be the methods of choice for low to moderate wavelengths of radiation and grating roughness that is representable by a few (e.g., tens of) Fourier coefficients. We argue that, in this context, the method can, for instance, produce full and precise reflectivity maps in computational times that are orders of magnitude faster than those of alternative numerical schemes (e.g., the popular “C-method,” finite differences, integral equations or finite elements). In this initial study we concentrate on the description of the basic principles that underlie the solution scheme, including those that relate to analytic continuation procedures. Within this framework, we explain how, in spite of conventional wisdom to the contrary, the resulting perturbative techniques can provide a most valuable tool for practical investigations in plasmonics. We demonstrate this with some examples that have been previously discussed in the literature (including treatments of the reflectivity and band gap structure of some simple geometries) and extend this to demonstrate the wider applicability of the proposed approach.

© 2013 Optical Society of America

## 1. INTRODUCTION

The interest in surface plasmons (SP), waves that propagate at the surface of a conductor [1,2], stems largely from their tight confinement of optical energy which allows for device miniaturization and sensitivities that are well beyond those achievable by current diffraction-limited technologies. The potential application of SP-based structures encompasses a variety of areas, including nanoscale photonic circuits [3], optical data storage [4], near-field imaging [5], photolithography [6], biosensing [7], etc. The research on and use of nanoplasmonic devices has seen a significant surge over the last twenty years as new phenomena have consistently been discovered and fabrication and characterization techniques have become ever more sophisticated and precise, to the point where surface roughness can be controlled to an accuracy below 1 nm (see e.g., [8] and the references therein). These advances, in turn, have provided an impetus for the development of accurate simulation tools that can resolve the effects of such small features and thus aid in the design process. A number of methodologies, both classical and specific to the relevant configurations, have been applied in this context, with varying degrees of success [9–12]. In this paper we show that high-order perturbation of surfaces (HOPS) methods enhanced with techniques of analytic continuation are well-suited for the analysis and design of nanoscale optical devices, and particularly so for those wherein the coupling of light to the metal is enabled through (periodic) corrugations on the surface of the latter.

The recent history of attempts at incorporating simulation work into the characterization of the devices that rely on the excitation of surface plasmon resonances (SPR) is very rich in number and variety of techniques. Indeed we find, for instance, that methods based on widely applicable numerical procedures for the modeling of scattering experiments, such as finite elements (see e.g., [13–17]), finite differences (e.g., [18–23]) and integral equations (e.g., [24–28]) have all been used in this context. As is the case in every other application and while, as we mentioned, these methods provide a wide domain of applicability, they also display significant challenges. Finite elements and finite differences, for instance, entail the (very fine, in most cases) discretization of entire domains (with the resulting computational expense) even if, as it is typically the case, the optical characteristics of the device are piecewise homogeneous. Moreover, and as a consequence of this, the computational domain must be artificially truncated, which gives rise to the need for the selection of appropriate “absorbing (or transparent) boundary conditions,” “perfectly matched layers” or other approximate treatment (see e.g., [29–31]). (Surface) integral equation based methods, on the other hand, avoid both of these impediments, as they rely on *surface* discretizations and their formulation intrinsically encodes the outgoing character of diffracted waves. These advantages, however, are attained at the expense of introducing *full*, i.e., not sparse, (impedance) matrices whose numerical inversion can render the methods uncompetitive, unless appropriate acceleration schemes are introduced (see, for instance, [24,28,32,33], and the references cited there). Furthermore, in the case of (infinite) periodic surfaces, a most significant challenge arises from the efficient evaluation of the appropriate Green function (see e.g., [34] and the references therein).

For the specific configurations provided by infinitely periodic (“rough”) surfaces, a number of alternative procedures have been devised (see e.g., the reviews in [35–37]). In this connection, a particularly popular method in nanoplasmonics is the “Chandezon Method” (or “C-method”) [38–40], which has been successfully applied in numerous studies (see e.g., [41–49] among many others). The method, and related schemes [40,50] (such as the “Rigorous Coupled Wave Analysis” (RCWA) and the “Fourier Modal Method” (FMM) [51–53]), relies on a simple coordinate transformation (or successive slicing, in the case of RCWA) that “flattens” the interfaces to simplify the imposition of boundary/transmission conditions. This, however, is done at the expense of foregoing the classical Rayleigh expansion [which is no longer representative of the field in the new coordinates; see Eq. (10) below], and it thus necessitates the derivation of new basis functions. This, in turn, can be posed as an eigenvalue problem, with the consequently high associated computational cost (see e.g., [50,54]).

Boundary perturbation schemes, such as the “small perturbation method” (SPM) [55], the “small slope approximation” (SSA) [56], or the “unified perturbation method” [57], on the other hand, typically rely on the explicitness of the Rayleigh expansion (or other representations, such as that provided by the $T$-matrix) to rapidly deliver scattering results. Most often, these techniques are used in low-order instances (see e.g., [58–60]) which restricts their domains of validity to very small deformations of a planar interface. Attempts at expanding these domains have been pursued through the derivation of higher-order approximations [61–65], but these can often result in only slight improvements. In spite of the long-held view that this is a consequence of the failure of the so-called “Rayleigh hypothesis” (RH) [36,66–68] for larger deformations, it has been rigorously shown that the latter plays no role in the intrinsic limitations of these procedures, as the fields can be proved to be *complex analytic* functions of a parameter measuring (analytic) boundary deformations for values of this parameter that include a complex neighborhood of *the entire real line* [69]. As a consequence, it follows that the limitations on these higher-order attempts are solely due to the finite radius of convergence of the resulting power series, that is, they arise as a result of the insistence on straightforward summation of the latter.

The observation that the field can, in fact, be analytically extended to a neighborhood of the real line, on the other hand, suggests that, while the totality of the information on the fields as the boundary perturbation parameter is varied is encoded in its Taylor series, a direct summation of the latter is not an optimal strategy. Based on these conclusions, a suite of novel numerical procedures was devised [70–73] that introduced mechanisms of analytic continuation (e.g., conformal maps [70], Padé approximation [71]) to significantly enlarge the domain of applicability of perturbative approaches. Indeed, it is now recognized that “… these results are of fundamental importance…” [74], and that the consequent algorithms constitute “… efficient numerical schemes…” that “… improve drastically the accuracy and the range of the SPM…” [36,58].

In what follows we show that the specific characteristics of these HOPS methods are well-suited to the evaluation of grating-mediated SPR generation. Indeed, as is well established, SPR arise in regimes wherein the height $h$ of the grating coupler is small compared to the period $d$, and the latter is comparable to the wavelength. Moreover, rich and varied reflectivity maps can be attained with the use of gratings that are representable by only a handful (tens) of Fourier modes, thus limiting their slope. As we explain below, in this regime the HOPS methods provide unparalleled performance and thus constitute excellent candidates to guide virtual design studies.

The rest of the paper is organized as follows. First, in Section 2 we set up the problem, describe our notation, and recall some basic facts on the origins and characteristics of surface plasmon polaritons (SPP). Section 3 is devoted to a brief review of the algorithms in the context of plasmonics, including a discussion of the materials and material properties that are relevant in the visible regime. The performance of the numerical procedure in this realm is analyzed in Section 4, where a variety of simulations are presented that demonstrate the effectiveness of the proposed procedures. We conclude in Section 5 with a summary and a survey of perspectives for the approach presented herein.

## 2. PRELIMINARIES: MAXWELL’S EQUATIONS AND SURFACE PLASMONS

To begin, we refer to the configuration in Fig. 1 where we display a generic (metallic) grating underlying a dielectric, and where we define our notation. The grating lines extend in the $z$ direction, and the profile is given by a function,

where $d$ is the period. The plane of incidence contains the grating wavevector, and the incidence angle is denoted by $\theta $; when applicable, the “line width” of a grating will be denoted by $w$.As is well known, in (“two-dimensional”) configurations such as that in Fig. 1, Maxwell’s equations decouple into two scalar problems corresponding to the transverse electric (TE) and transverse magnetic (TM) modes of polarization (see e.g., [71]). These are characterized by the alignment of the electric or magnetic fields with the grating lines, respectively. In the case of monochromatic waves oscillating at a single frequency $\omega $,

The incident wave, ${u}^{\text{inc}}$, on the other hand, will be assumed to be a monochromatic plane of unit amplitude wave,

where $\alpha ={k}^{+}\text{\hspace{0.17em}}\mathrm{sin}(\theta )$, $\beta ={k}^{+}\text{\hspace{0.17em}}\mathrm{cos}(\theta )$ and $\theta $ denotes the angle of incidence (see Fig. 1). The quasi-periodicity of the incident wave, implies a similar property for the scattered field, and this, in turn, allows for the (Rayleigh) expansion of the latter,*finite*number of reflected waves are propagating (namely, those corresponding to indices $r$ in the set $U=\{r:{\beta}_{r}^{+}>0\}$), while the rest decay in the far field.

The objective, then, is to design a scheme to determine the amplitudes ${B}_{r}^{\pm}$ of the reflected and transmitted waves that, in this case, will allow for the assessment and design of structures that support SPP. These are particular electromagnetic waves that propagate at the interface between a dielectric and a conductor, and are tightly confined in the direction of the normal to the surface. These surface waves arise through the coupling of the electromagnetic fields to oscillations of the conductor’s electron plasma [1], and they are closely related to (and can be deduced from) the Brewster mode [75,76].

The macroscopic study of SPP begins with the simple realization that, if the interface is flat, their dispersion relation can be explicitly derived to yield [1]

In terms of the overall expansions in Eq. (10), this behavior manifests itself through a transfer of incident radiation into modes that decay in the far field, in a manner so that the entire reflected field is significantly suppressed away from the dielectric/metal interface. More precisely, in this context, the condition for the occurrence of a SPR for a given angle of incidence can be identified with the presence of a drop in the reflectivity

or “normalized reflectivity,”In the next section we present an efficient high-order boundary perturbation method that allows for the rapid evaluation of the quantities in Eq. (15) as functions of the grating height and the wavelength of radiation and/or incidence angle.

## 3. A FAST, HIGH-ORDER BOUNDARY PERTURBATION

In this section we briefly review the HOPS method that was originally introduced in [70–72]; see also [73] for a more comprehensive review of subsequent improvements and applications. As indicated in Fig. 1, here we shall concentrate on “two-dimensional” configurations [71], where the lines of the gratings are assumed to have infinite extent; every development that follows, however, can be extended to fully three-dimensional (biperiodic) structures, at the expense of somewhat more complex derivations [72].

The basic observation underlying the procedure relates to the *explicit* nature of the solution in the case wherein the interface is entirely flat. This, in turn, suggests that this solution may be *analytically* continued to that of an undulated surface. More precisely, for a normalized profile $f(x)$ satisfying

As was shown in [71] (see also [69]), the coefficients ${B}_{r}(h)$ are *complex analytic* functions of the parameter $h$ (putting to rest a fifty-plus year controversy as to its validity as an actual convergent series; see e.g., [63,67,74,82] and the references therein), and as such these can be expressed in the form

*converges*for small values of the parameter $h$. More importantly, as shown in [69], these coefficients (and the field itself) are

*complex analytic in a whole neighborhood of the real line*; see the example in Fig. 2 (here and throughout, and whenever appropriate, the axes denote measures in

*nanometers*).

Arguing as in [71] it is straightforward to show that the coefficients ${d}_{n,r}^{\pm}$ satisfy the following *recursion*:

*convolution*nature of the inner sums (that can therefore be effected with fast Fourier transforms), and to the relatively low orders $n$ that suffice to attain accurate results throughout the domain of interest, if mechanisms of analytic continuation are incorporated. In relation with the latter, and following [71], we have chosen to rely on the

*Padé approximants*that can be readily evaluated from knowledge of the Taylor coefficients ${d}_{n,r}^{\pm}$ [Eq. (21)]. In this regard, throughout this paper we shall use the standard notation “[$L/M$]” to refer to a Padé approximant (i.e., a rational function approximation) with polynomials of degrees $L$ and $M$ in the numerator and denominator, respectively; we refer the interested reader to the treatise [83] for the mechanics of evaluation and fascinating properties of Padé approximation.

Clearly, the results are influenced by the dispersion model that is applicable in each instance. Figure 2(a), for example, was derived using a “classical” Lorentz model [84] that provides a characterization of the frequency-dependent permittivity in silver. Alternative models that depend (rather weakly) on the preparation of the sample include, for instance, the polynomial expressions of [85] and [86]. Models can also de derived from the well-known discretely sampled values of Johnson and Christy [87] or of Lynch and Hunter [88] (e.g., through suitable interpolation or by matching to an appropriately defined Lorentz model).

As is well-recognized [89] the appearance of classical SPR is confined to a rather well delineated subset in the parameter space consisting of the wavelength, angle of incidence, period, and the shape of the profile itself. Indeed, these grating-mediated SPR occur for small to moderate modulations that are much smaller than the wavelength of radiation which, in turn, is comparable to the period. The latter condition follows readily from the requirement that coupling occur at low ($0,\pm 1,\pm 2,\pm 3,\dots $) diffracted orders. Further increases in height, on the other hand, will typically conspire against efficient SPP coupling due to the concomitant changes in the phase of reradiated light [43,90]. It should be noted, however, that very deep gratings may still give rise to localized field enhancements, though this can be qualitatively distinguished as arising from waveguide and hybrid waveguide-SPP modes (see e.g., [85,91–93] and the references therein).

As we show next, and by its very nature, the HOPS method presented here is particularly advantageous for simulations of classical SPR, whose sharp features and high sensitivity to changes in the low amplitude harmonics of the interface make them ideally suited for applications that seek to quantify either the surface shape or the permittivities of the constituent media [43].

## 4. NUMERICAL RESULTS

In this section we shall present a variety of numerical results that demonstrate the accuracy and versatility of the proposed scheme for the analysis of grating-mediated SPR generation. As is well known, these resonances occur in the TM polarization [1], so we shall henceforth constrain our experiments to such illuminations [${C}_{0}={k}^{+}/{k}^{-}$ in Eq. (6)]. More precisely, we begin in Section 4.A with a demonstration of the benefits of incorporating mechanisms of analytic continuation into the HOPS methods within this context. Indeed, we show that these techniques can have a dramatic effect not only in expanding the domain of applicability of perturbative approaches [cf. Fig. 2(b)] but that they also effect an intrinsic acceleration of the convergence of the corresponding Taylor series even when the perturbations are sufficiently shallow to allow for its direct summation. Next, in Section 4.B we expand on the examples of the previous section to demonstrate the effectiveness of the approach in computing full reflectivity maps, as well as field representations. Finally, this demonstration is complemented in Section 4.C. where experiments on the (three-dimensional) dependence of the reflectivity on height, wavelength. and angle of incidence are shown to recover the appearance of “energy (photonic) band gaps.”

#### A. Analytic Continuation

In this section we present some results that demonstrate that analytic continuation mechanisms significantly expand on the applicability of boundary-perturbation methods and that they allow for a complete coverage of the range of SPP excitations, as described at the end of the previous section [allowing, in fact, for calculations that are well-beyond the regime of “shallow” modulations (see e.g., Fig. 5 below)]. In particular, as we show here, these procedures negate a variety of arguments that have been advanced over the last few decades (and that continue to be so) in connection with the suitability of perturbative approaches for SPP simulations (and, more generally, for simulations concerning diffractive nano-optical structures); see e.g., [41,94].

Our first example is motivated by that in [94] where it is argued that the results from perturbation theory display a behavior that ranges from differing “… from those of Maradudin’s and Rayleigh’s methods not only in magnitude but also in the angle…” for small perturbations (6 nm) to differences “… by 1 order of magnitude…” for slightly deeper profiles (10 nm). Unfortunately, it seems impossible to reproduce the specific results presented there with the description provided by the authors, as it *does not* appear feasible to excite a resonance at an incidence angle of 19.31° with a wave of energy 3.5 eV that illuminates a sinusoidal silver grating,

Under these conditions, in Fig. 3 we present, as in [94], approximations to the square of the absolute value ${|{B}_{-1}^{+}|}^{2}$ of the $-1$st order Rayleigh coefficient ${B}_{-1}^{+}$ [Eq. (10)] as garnered from the perturbative treatment described in Section 3 for heights $h=4$, 6 and 10 nm [Figs. 3(a)–3(c), respectively]. As in [94], the figure shows a clear deterioration of the results if a series with $N=5$ terms is used for increasing heights. In fact, the results corresponding to the deepest grating [Fig. 3(c)] show that such an approximation can differ *qualitatively* from more accurate representations corresponding to approximations of orders $N=7$, 9, and 11, which clearly converge in this case. More importantly, perhaps, the results also show that the limitations alluded to in [94] are solely due to the insistence on approximating the analytic continuation of ${B}_{-1}^{+}={B}_{-1}^{+}(h)$ through a Taylor series. Indeed, as the figure shows, a simple change to Padé approximation of order [2/2] (which uses exactly the same information as the Taylor series of order $N=5$) provides results with *several digits of accuracy*. Indeed, for instance, Fig. 4 shows the maximum *relative* error in the value of ${|{B}_{-1}^{+}|}^{2}$ for $h=10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ over the interval [19.2°,19.4°] for Taylor and Padé approximations of orders $N$ and $[(N-1)/2,(N-1)/2]$, respectively, when these are compared to the results from Padé approximants of order [12/12] ($N=25$), and it shows that the results of the [2/2] approximants have a maximum error of $3.3\times {10}^{-4}$, guaranteeing three-digit accuracy; similar calculations show that the accuracy of the [2/2] approximant for $h=4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and $h=6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ is five and four digits, respectively. These results constitute an instance of the property of analytic continuation mechanisms that we alluded to above, namely, the fact that these procedures *accelerate* convergence to the exact solution. This is clearly exemplified in Fig. 4 which shows, for instance, that a Taylor approximation of order $N=13$ still incurs an error of about $4.9\times {10}^{-3}$ while the corresponding [6/6] Padé approximant displays close to *full double precision accuracy*.

To complement the above, for our second example we retain the sinusoidal nature of the (silver) grating, but we explore a domain that lies well beyond that where low-order perturbative approaches are applicable. In fact, as we show, the fields corresponding to this configuration *do not admit* a high-order representation in terms of Taylor series, and Padé approximation (or other analytic continuation mechanisms) are essential to attain convergence. The configuration is taken from the work of Chen *et al.* [85] and is depicted in Fig. 5(a): a deep silver structure with period $d=258\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and height $h=124\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, illuminated at an angle of 30°; the permittivity of the metal is assumed to be given by the polynomial model on p. 1574 of [85]. Figure 5(b) shows the computed reflectivity at order zero, ${|{B}_{0}^{+}|}^{2}$, as a function of the wavelength throughout the visible range as computed from [32/32] Padé approximants; a dashed line marks the wavelength $\lambda =387\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ below which the $-1$st order is propagating. A numerical convergence analysis suggests that the displayed results have at least three digits of accuracy over the entire range, and at least four digits outside the interval from 405 to 415 nm.

In Fig. 6 we display results corresponding to the configuration in Fig. 5(a) as a function of the height $h$ for wavelengths $\lambda =410\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and $\lambda =495\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. In this connection it is worth noting that, in contrast with other simulation methodologies, the evaluation of the scattering amplitudes at different heights within the boundary perturbation framework is immediate from knowledge of the terms of the Taylor series or of its analytic continuation (as this simply entails a most straightforward and inexpensive evaluation of the approximation at specified points). Figures 6(a) and 6(c) show the reflectivity at order zero as a function of the height $h$: each figure displays various truncated Taylor series approximations as indicated, together with the corresponding Padé approximants (solid lines). In each case, the value ${h}_{p}$ corresponds to the absolute value of the closest (noncanceling) pole to the origin as approximated by the corresponding pole of the [32/32] Padé approximant ($N=65$). These values can be garnered from Figs. 6(b) and 6(d), where we show the complete set of poles and zeros in the complex $h$ plane of this latter approximation: the circles and crosses (square and stars) mark the location of the zeros and poles, respectively, that do (do not) cancel out. As it follows from these figures, the value of ${h}_{p}$ is strongly dependent on the frequency, and it provides a good approximation to the radius of convergence of the Taylor series. In particular, we see that in each case the series approximation is unable to capture the behavior for deep modulations such as that in Fig. 5(a), while the analytically continued representations (in the form of Padé approximants) do provide us with meaningful results.

Our final example in this section, entails a configuration that is far from sinusoidal as it corresponds to a (semi-) elliptical grating, which we assume to be optically described by the dispersion relation that follows from (cubic-spline) interpolation of the results in [88]. A Fourier representation with $F=25$ [Eq. (24)] was garnered from the exact description of the top portion of an ellipse with major axis of length $w=400\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and periodized with period $d=630\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$; the Fourier coefficients were summed in the Cesàro sense (i.e., using Fejér means) to minimize the oscillatory behavior at corners (see e.g., [95]).

Figure 7(a) shows the configuration to scale, with an amplitude $h=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$; the maximum slope, on the other hand is approximately 32 for this height, and 50 for $h=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. The results in (b)–(d) correspond to slight changes in the incidence angle, from $\theta =0\xb0$ (normal incidence), to $\theta =0.2\xb0$ to $\theta =1\xb0$ for heights of 20 and 30 nm. The results were obtained from [5/5] Padé approximants ($N=11$) and the error is no larger than $6\times {10}^{-7}$ and $5\times {10}^{-4}$ for $h=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and $h=30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, respectively. The corresponding Taylor sums display minimum errors of 25% and 100%, respectively, with $N=11$; increasing this to $N=65$ shows that, in fact, the series *diverges* for wavelengths in the vicinity of the resonances for heights that exceed approximately 15 nm.

It is worth noting, however, that in this case the amplitudes are sufficiently small so as to allow for a comparison with the resonances as predicted by Eq. (12) corresponding to a flat interface. Figure 8 provides a graphical display of these predictions. Indeed, for a grating of period $d$, Eq. (12) can be rewritten as

#### B. Reflection Data, Reflectivity Maps, and Electromagnetic Fields

In this section we expand on the results above and provide examples of entire “reflectivity maps” as functions of both wavelength and height. As we stated, the latter evaluation of the dependence of the reflectivity on the depth of the structure entails almost no computational effort within the HOPS method we review here, as this only demands the straightforward calculation of a Padé approximant (or other approximations to the analytic continuation of a Taylor series) at a number of points. In addition, we provide examples of actual field distributions as computed from the HOPS method that exemplify the “surface” character and subwavelength localization of the SPP.

To begin, in Fig. 9 we present the total reflectivity maps and an instance of the spatial distribution of the reflected (transverse) magnetic field for a configuration as in Fig. 5(a); similar results for the geometry in Fig. 7(a) are presented in Fig. 10.

In both cases the results correspond to simulations that rely on the HOPS method using Padé approximants of order [16/16] and spanning heights up to 125 nm. In the case of the elliptical grating, we further consider “negative” heights corresponding to elliptical *grooves* [rather than elliptical *bumps* as in Fig. 7(a)]. As in the examples above, we use the polynomial dispersion relation [85] for the sinusoidal profile, and the interpolated data from [88] in the case of the elliptical geometry. We note that, in the case of the sinusoidal structure, strong coupling to an SPP occurs for wavelengths that are slightly above that which marks the propagation of the $-1$st order mode and only for small to moderate modulations of less than 60 nm ($h/d\approx 0.23$), with truly negligible reflection for heights below 30 nm ($h/d\approx 0.12$). Indeed, Fig. 9(b) shows that for a height $h=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and a wavelength $\lambda =410\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, a six-fold enhancement of the field strength is attained in a subwavelength region near the interface, with nearly vanishing reflected fields a wavelength away. Similar results are obtained for the elliptical profile, where the period is significantly larger and so are the heights for which SPPs can be excited. Interestingly, for deeper profiles the longer wavelength branch that is present for small heights, and that can be identified with one that would be present on a flat interface [Fig. 8(d)], ceases to exist, as the resonant character of the structure begins to differ significantly from that which is derived from approximation by a planar surface. As shown in Fig. 10(b), at a height of 30 nm and a wavelength of 640 nm the local enhancement of the field is approximately ten-fold.

#### C. Photonic Energy Gaps (“band gaps”)

In this final section we present results that are similar in nature to those described in Section 4.B but where we interchange the variations of the height at a fixed incidence angle for variations of the incidence angle at a fixed height. In this manner, we are able to investigate the existence and dependencies of “band gaps” in the spectrum for SPP excitation (see [96] and the references therein). As with our previous examples, again here our results can be interpreted as countering the belief [41] that perturbation treatments lead to a “… limited range of applicability…” as “… shown by da Silva *et al.* [94]…” (see the first example in Section 4.A) and that perturbation methods that are “… applicable to larger amplitude gratings… must be in question…”.

For comparison purposes, the examples we present follow closely those in [41], where the effect of a second harmonic [$F=2$ in Eq. (24)] on the reflectivity of purely sinusoidal gratings is analyzed with the method of Chandezon *et al.* [39]. More precisely, we shall consider profiles of the form

*not*for a (slight) perturbation that includes a second harmonic (in fact, the Taylor series

*diverges*in this case; see Fig. 14 below). As shown in Fig. 13, on the other hand, this latter behavior can be overcome with the use of Padé approximants of order [2/2], even though, as we have alluded to before, this uses

*precisely*the same information as that encoded in the first five terms of the Taylor series.

As in [41] we find that for a purely sinusoidal structure the intersection of the SPR generated by the orders $\pm 1$ takes place at normal incidence, with little interaction among them as evidenced by the absence of a gap; see Fig. 13(a). In contrast, and again as in [41], Fig. 13(b) shows that a small additive perturbation in the second order harmonic allows for the coupling of these modes in a manner so as to generate a band gap. An approximate formula for the central wavelength $\overline{\lambda}$ of this gap is provided in Eq. (3b) in [96] (see also Eq. (5.52) in [41]), based on the neglect of higher order harmonics and the assumption that

in the formulation of the scattering problem as provided by the approach in [39]. This results in an expression,It is worth noting that, in contrast with the sentiment that “… the number of terms required to construct a solution to a given precision increases extremely rapidly as the depth of the grating rises…” [41], the results in Fig. 13 were obtained with only *five* terms and, thus, in the form of a [2/2] Padé approximant, this approximate solution should be equally amenable to an analytic study as the (truncated) $4\times 4$ “$T$-matrix” from [39] as used in [41]. Although we shall defer this analysis to a future publication, we remark here that a numerical investigation of the errors associated with the results in this figure suggests that the *maximum relative error* is less than 0.1% for Fig. 13(a) and about 3% for Fig. 13(b); a numerical convergence analysis is included in Fig. 14.

Our final example is related to the above, and it corresponds to that in Fig. 8 in [41]. In Fig. 15(a) we display the reflectivity corresponding to a sinusoidal profile, such as that in Fig. 12(a), but with a much deeper height of $h=60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$. To accurately deal with these depths we resort to a higher order version of the HOPS method using Padé approximants of order [8/8]. We see that, as a result of the change in height a gap opens in angle rather than in frequency. This transition is further detailed in Fig. 15(b), where the progression of reflectivities is shown from 10 to 60 nm in increments of 10 nm; the estimated maximum relative error for the entirety of the results in this latter figure is approximately $1\times {10}^{-5}$; that is, we estimate that each value of the reflectivity has at least 5 digits of accuracy (and increasing with decreasing height). We note that this figure is largely in qualitative agreement with Fig. 8 in [41], though some differences appear visible. On account of our estimated accuracy, we conjecture that these differences arise as a consequence of a lack of resolution in the results in [41] that arise from the use of the method in [39]. Finally, and for comparison purposes, we also include in Fig. 16 analogous results for profiles, such as that in Fig. 12(b) for varying heights $h=10\u201360\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$ and constant frequency content ratio $a/b=2.5$ [Eq. (29)]; the estimated maximum relative error in this case is about $3\times {10}^{-5}$.

## 5. CONCLUSIONS

In this paper we have introduced a *rigorous* high-order perturbation method (HOPS) to bear upon the area of plasmonics. Here, we have shown that the HOPS method [71], which couples ideas of regular perturbation theory and analytic continuation, is particularly well-suited to applications in this area, which typically relate to small and periodic deviations of flat surfaces (to couple light to SPP) in a regime where the period of the resulting structure is comparable to the wavelength of visible light. We have further explained how previous results and conclusions relating to boundary perturbation methods incorrectly attributed their failure to that of the RH and/or dismissed them as limited in applicability due to an unfounded insistence on summation of the resulting Taylor polynomials. In contrast, we have demonstrated that the limitations of HOPS methods are unrelated to the RH and, moreover, that these are significantly less restrictive for the investigation of SPR than previously assumed, provided that suitable mechanisms of analytic continuation (e.g., Padé approximation) are used.

As with other popular numerical simulation schemes (e.g., the C-method [38] or methods based on frequency domain integral equations [28]) the approach we elaborate upon here resolves the scattering problem at a single, fixed frequency and, thus, its repeated use is necessary to garner the complete spectral response of a given structure. In contrast with other techniques, however, the investigation of the dependence of the reflectivity upon the height of the grating coupler can be obtained through the HOPS discussed herein, with minimal additional computational cost. Moreover, we have explained how the very nature of the formulation allows for FFT-accelerated evaluations. The combination of these factors results in a scheme that can deliver full reflectivity maps at a fraction of the computational expense associated with alternative frequency-domain procedures. Comparisons with time-domain techniques [such as the popular finite-difference time-domain (FDTD) method], on the other hand, are left for future work; however, we anticipate that the extremely fine meshes that are known to be necessary to resolve plasmonic resonances (see e.g., [97]) and the extended time intervals that would allow for accurate reconstruction of the spectral response will render these uncompetitive (and, at times, simply impracticable) when compared to the schemes discussed above.

The developments presented here suggest a number of avenues for further investigation. In particular, as we mentioned, further quantitative comparison with and verification against state-of-the-art (e.g., efficient solvers based on surface integral equations) and/or popular (e.g., FDTD) alternative techniques will be forthcoming. In addition, the specifics related to modeling of SPR supports the idea of a search for adaptive procedures (e.g., ones that may take advantage of their “localized” nature in parametric space; see, e.g., Figs. 9 and 10) that may further lower the overall computational cost. This is particularly relevant in connection with an ultimate goal of incorporating mathematical optimization techniques that will allow for rapid virtual prototyping, particularly in the context of biosensing. In this connection, and due to the simplicity, accuracy and efficiency of the approach discussed here, we fully expect that thoughtful parallelization strategies (on CPUs and, particularly, on GPUs) will enable the development of a *real-time* simulator for use by practitioners. Finally, as we have explained, the extension of the present developments to allow for the treatment of fully three-dimensional geometries is rather straightforward [72]; slightly more complex, though still extremely efficient, is that to multi-layer structures (see e.g., [98]). The implications of both of these extensions in the study of SPR will be the subject of future work.

## ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation through grant number DMR-0941537. TWJ and SHO gratefully acknowledge support from the National Science Foundation through grant number CBET-1067681. The authors wish to thank the anonymous referees for their thoughtful and detailed comments.

## REFERENCES

**1. **S. A. Maier, *Plasmonics: Fundamentals and Applications* (Springer, 2007).

**2. **L. Novotny and B. Hecht, *Principles of Nano-Optics* (Cambridge University, 2012).

**3. **D. K. Gramotnev and S. I. Bozhevolnyi, “Plasmonics beyond the diffraction limit,” Nat. Photonics **4**, 83–91 (2010). [CrossRef]

**4. **M. Mansuripur, A. R. Zakharian, A. Lesuffleur, S.-H. Oh, R. J. Jones, N. C. Lindquist, H. Im, A. Kobyakov, and J. V. Moloney, “Plasmonic nano-structures for optical data storage,” Opt. Express **17**, 14001–14014 (2009). [CrossRef]

**5. **T. W. Johnson, Z. J. Lapin, R. Beams, N. C. Lindquist, S. G. Rodrigo, L. Novotny, and S.-H. Oh, “Highly reproducible near-field optical imaging with sub-20-nm resolution based on template-stripped gold pyramids,” ACS Nano **6**, 9168–9174 (2012). [CrossRef]

**6. **F. J. García-Vidal, L. Martin-Moreno, T. W. Ebbesen, and L. Kuipers, “Light passing through subwavelength apertures,” Rev. Mod. Phys. **82**, 729–787 (2010). [CrossRef]

**7. **J. Homola, “Surface plasmon resonance sensors for detection of chemical and biological species,” Chem. Rev. **108**, 462–493 (2008). [CrossRef]

**8. **N. C. Lindquist, P. Nagpal, K. M. McPeak, D. J. Norris, and S.-H. Oh, “Engineering metallic nanostructures for plasmonics and nanophotonics,” Rep. Prog. Phys. **75**, 036501 (2012). [CrossRef]

**9. **D. Barchiesi, B. Guizal, and T. Grosges, “Accuracy of local field enhancement models: toward predictive models?” Appl. Phys. B **84**, 55–60 (2006). [CrossRef]

**10. **G. Veronis and S. Fan, “Overview of simulation techniques for plasmonic devices,” in *Surface Plasmon Nanophotonics*, M. Brongersma and P. Kik, eds. Vol. 131 of Springer Series in Optical Sciences (Springer, 2007), pp. 169–182.

**11. **J. Hoffmann, C. Hafner, P. Leidenberger, J. Hesselbarth, and S. Burger, “Comparison of electromagnetic field solvers for the 3D analysis of plasmonic nano antennas,” Proc. SPIE **7390**, 73900J (2009).

**12. **J. Smajic, C. Hafner, L. Raguin, K. Tavzarashvili, and M. Mishrikey, “Comparison of numerical methods for the analysis of plasmonic structures,” J. Comput. Theor. Nanosci. **6**, 763–774 (2009). [CrossRef]

**13. **A. Schädle, L. Zschiedrich, S. Burger, R. Klose, and F. Schmidt, “Domain decomposition method for Maxwell’s equations: scattering off periodic structures,” J. Comput. Phys. **226**, 477–493 (2007). [CrossRef]

**14. **G. Demésy, F. Zolla, A. Nicolet, and M. Commandré, “Versatile full-vectorial finite element model for crossed gratings,” Opt. Lett. **34**, 2216–2218 (2009). [CrossRef]

**15. **M. Huber, J. Schöberl, A. Sinwel, and S. Zaglmayr, “Simulation of diffraction in periodic media with a coupled finite element and plane wave approach,” SIAM J. Sci. Comput. **31**, 1500–1517 (2009). [CrossRef]

**16. **K. Stannigel, M. König, J. Niegemann, and K. Busch, “Discontinuous Galerkin time-domain computations of metallic nanostructures,” Opt. Express **17**, 14934–14947 (2009). [CrossRef]

**17. **O. Tsilipakos, A. Pitilakis, A. Tasolamprou, T. Yioultsis, and E. Kriezis, “Computational techniques for the analysis and design of dielectric-loaded plasmonic circuitry,” Opt. Quantum Electron. **42**, 541–555 (2011). [CrossRef]

**18. **D. Christensen and D. Fowers, “Modeling SPR sensors with the finite-difference time-domain method,” Biosens. Bioelectron. **11**, 677–684 (1996). [CrossRef]

**19. **H. Sai, Y. Kanamori, K. Hane, and H. Yugami, “Numerical study on spectral properties of tungsten one-dimensional surface-relief gratings for spectrally selective devices,” J. Opt. Soc. Am. A **22**, 1805–1813 (2005). [CrossRef]

**20. **G. Veronis and S. Fan, “Modes of subwavelength plasmonic slot waveguides,” J. Lightwave Technol. **25**, 2511–2521 (2007). [CrossRef]

**21. **J. M. Montgomery, T.-W. Lee, and S. K. Gray, “Theory and modeling of light interactions with metallic nanostructures,” J. Phys. Condens. Matter **20**, 323201 (2008). [CrossRef]

**22. **N. C. Lindquist, T. W. Johnson, D. J. Norris, and S.-H. Oh, “Monolithic integration of continuously tunable plasmonic nanostructures,” Nano Lett. **11**, 3526–3530 (2011). [CrossRef]

**23. **A. Shahmansouri and B. Rashidian, “Comprehensive three-dimensional split-field finite-difference time-domain method for analysis of periodic plasmonic nanostructures: near- and far-field formulation,” J. Opt. Soc. Am. B **28**, 2690–2700 (2011). [CrossRef]

**24. **O. P. Bruno and M. C. Haslam, “Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences,” J. Opt. Soc. Am. A **26**, 658–668 (2009). [CrossRef]

**25. **A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A **26**, 732–740 (2009). [CrossRef]

**26. **B. Gallinet, A. M. Kern, and O. J. F. Martin, “Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral approach,” J. Opt. Soc. Am. A **27**, 2261–2271 (2010). [CrossRef]

**27. **J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araújo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A **28**, 1341–1348 (2011). [CrossRef]

**28. **H. Kurkcu, A. Ortan, and F. Reitich, “An efficient integral equation solver for two-dimensional simulations in nanoplasmonics,” in Proceedings of Waves 2013, Tunis, Tunisia (2013).

**29. **B. Alavikia and O. Ramahi, “Limitation of using absorbing boundary condition to solve the problem of scattering from a cavity in metallic screens,” in *Antennas and Propagation Society International Symposium (APSURSI)* (IEEE, 2010), pp. 1–4.

**30. **C.-C. Chao, S.-H. Tu, C.-M. Wang, H.-I. Huang, C.-C. Chen, and J.-Y. Chang, “Impedance-matching surface plasmon absorber for FDTD simulations,” Plasmonics **5**, 51–55 (2010). [CrossRef]

**31. **M. Wang, C. Engstrom, K. Schmidt, and C. Hafner, “On high-order FEM applied to canonical scattering problems in plasmonics,” J. Comput. Theor. Nanosci. **8**, 1564–1572 (2011). [CrossRef]

**32. **Y. Otani and N. Nishimura, “A periodic FMM for Maxwell’s equations in 3D and its applications to problems related to photonic crystals,” J. Comput. Phys. **227**, 4630–4652 (2008). [CrossRef]

**33. **A. D. Baczewski, N. C. Miller, and B. Shanker, “Rapid analysis of scattering from periodic dielectric structures using accelerated Cartesian expansions,” J. Opt. Soc. Am. A **29**, 531–540 (2012). [CrossRef]

**34. **H. Kurkcu and F. Reitich, “Stable and efficient evaluation of periodized Green’s functions for the Helmholtz equation at high frequencies,” J. Comput. Phys. **228**, 75–95 (2009). [CrossRef]

**35. **K. Warnick and W. Chew, “Numerical simulation methods for rough surface scattering,” Wave Random Media **11**, R1–R30 (2001). [CrossRef]

**36. **T. Elfouhaily and C. Guérin, “A critical survey of approximate scattering wave theories from random rough surfaces,” Wave Random Media **14**, R1–R40 (2004). [CrossRef]

**37. **J. A. deSanto, “Overview of rough surface scattering,” in *Light Scattering and Nanoscale Surface Roughness*, A. A. Maradudin and D. J. Lockwood, eds., Nanostructure Science and Technology (Springer, 2007), Chap. 8, pp. 211–235.

**38. **J. Chandezon, G. Raoult, and D. Maystre, “A new theoretical method for diffraction gratings and its numerical application,” J. Opt. **11**, 235–241 (1980). [CrossRef]

**39. **J. Chandezon, M. Dupuis, G. Cornet, and D. Maystre, “Multicoated gratings: a differential formalism applicable in the entire optical region,” J. Opt. Soc. Am. A **72**, 839–846 (1982). [CrossRef]

**40. **L. Li, J. Chandezon, G. Granet, and J.-P. Plumey, “Rigorous and efficient grating-analysis method made easy for optical engineers,” Appl. Opt. **38**, 304–313 (1999). [CrossRef]

**41. **W. L. Barnes, T. W. Preist, S. C. Kitson, and J. R. Sambles, “Physical origin of photonic energy gaps in the propagation of surface plasmons on gratings,” Phys. Rev. B **54**, 6227–6244 (1996). [CrossRef]

**42. **G. Granet, “Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: application to multilayer crossed gratings,” J. Opt. Soc. Am. A **15**, 1121–1131 (1998). [CrossRef]

**43. **R. A. Watts, A. P. Hibbins, and J. R. Sambles, “The influence of grating profile on surface plasmon polariton resonances recorded in different diffracted orders,” J. Mod. Opt. **46**, 2157–2186 (1999).

**44. **M. Kreiter, S. Mittler, W. Knoll, and J. R. Sambles, “Surface plasmon-related resonances on deep and asymmetric gold gratings,” Phys. Rev. B **65**, 125415 (2002). [CrossRef]

**45. **W.-C. Liu, “High sensitivity of surface plasmon of weakly-distorted metallic surfaces,” Opt. Express **13**, 9766–9773 (2005). [CrossRef]

**46. **S. Kahaly, S. K. Yadav, W. M. Wang, S. Sengupta, Z. M. Sheng, A. Das, P. K. Kaw, and G. R. Kumar, “Near-complete absorption of intense, ultrashort laser light by sub-*λ* gratings,” Phys. Rev. Lett. **101**, 145001 (2008). [CrossRef]

**47. **N. Bonod, E. Popov, L. Li, and B. Chernov, “Unidirectional excitation of surface plasmons by slanted gratings,” Opt. Express **15**, 11427–11432 (2007). [CrossRef]

**48. **I. S. Spevak, M. A. Timchenko, and A. V. Kats, “Design of specific gratings operating under surface plasmon–polariton resonance,” Opt. Lett. **36**, 1419–1421 (2011). [CrossRef]

**49. **L. Li and G. Granet, “Field singularities at lossless metal-dielectric right-angle edges and their ramifications to the numerical modeling of gratings,” J. Opt. Soc. Am. A **28**, 738–746 (2011). [CrossRef]

**50. **J. Bischoff and K. Hehl, “Perturbation approach applied to modal diffraction methods,” J. Opt. Soc. Am. A **28**, 859–867 (2011). [CrossRef]

**51. **M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. A **71**, 811–818 (1981). [CrossRef]

**52. **M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A **12**, 1077–1086 (1995). [CrossRef]

**53. **L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A **14**, 2758–2767 (1997). [CrossRef]

**54. **J. Bischoff, “Prospects and limits of the Rayleigh Fourier approach for diffraction modelling in scatterometry and lithography,” Proc. SPIE **7390**, 73901E (2009).

**55. **S. O. Rice, “Reflection of electromagnetic waves from slightly rough surfaces,” Commun. Pure Appl. Math. **4**, 351–378 (1951). [CrossRef]

**56. **A. G. Voronovich, “Small-slope approximation in wave scattering by rough surfaces,” Sov. Phys. J. Exp. Theor. Phys. **62**, 65–70 (1985).

**57. **E. Rodriguez and Y. Kim, “A unified perturbation expansion for surface scattering,” Radio Sci. **27**, 79–93 (1992). [CrossRef]

**58. **C.-A. Guérin and A. Sentenac, “Second-order perturbation theory for scattering from heterogeneous rough surfaces,” J. Opt. Soc. Am. A **21**, 1251–1260 (2004). [CrossRef]

**59. **D. Grieser, H. Uecker, S.-A. Biehs, O. Huth, F. Rüting, and M. Holthaus, “Perturbation theory for plasmonic eigenvalues,” Phys. Rev. B **80**, 245405 (2009). [CrossRef]

**60. **A. Trügler, J.-C. Tinguely, J. R. Krenn, A. Hohenau, and U. Hohenester, “Influence of surface roughness on the optical properties of plasmonic nanoparticles,” Phys. Rev. B **83**, 081412(R) (2011). [CrossRef]

**61. **A. A. Maradudin and E. R. Méndez, “Enhanced backscattering of light from weakly rough, random metal surfaces,” Appl. Opt. **32**, 3335–3343 (1993). [CrossRef]

**62. **K. A. O’Donnell, “High-order perturbation theory for light scattering from a rough metal surface,” J. Opt. Soc. Am. A **18**, 1507–1518 (2001). [CrossRef]

**63. **M. A. Demir and J. T. Johnson, “Fourth- and higher-order small-perturbation solution for scattering from dielectric rough surfaces,” J. Opt. Soc. Am. A **20**, 2330–2337 (2003). [CrossRef]

**64. **K. A. O’Donnell, “Small-amplitude perturbation theory for one-dimensionally rough surfaces,” in *Light Scattering and Nanoscale Surface Roughness*, A. A. Maradudin and D. J. Lockwood, eds., Nanostructure Science and Technology (Springer, 2007), Chap. 5, pp. 107–126.

**65. **H.-Y. Xie, M.-Y. Ng, and Y.-C. Chang, “Analytical solutions to light scattering by plasmonic nanoparticles with nearly spherical shape and nonlocal effect,” J. Opt. Soc. Am. A **27**, 2411–2422 (2010). [CrossRef]

**66. **J.-J. Greffet, “Scattering of electromagnetic waves by rough dielectric surfaces,” Phys. Rev. B **37**, 6436–6441 (1988). [CrossRef]

**67. **D. R. Jackson, D. P. Winebrenner, and A. Ishimaru, “Comparison of perturbation theories for rough-surface scattering,” J. Acoust. Soc. Am. **83**, 961–969 (1988). [CrossRef]

**68. **A. Soubret, G. Berginc, and C. Bourrely, “Backscattering enhancement of an electromagnetic wave scattered by two-dimensional rough layers,” J. Opt. Soc. Am. A **18**, 2778–2788 (2001). [CrossRef]

**69. **O. P. Bruno and F. Reitich, “Solution of a boundary value problem for the Helmholtz equation via variation of the boundary into the complex domain,” Proc. R. Soc. Edinburgh , Sect. A **122**, 317–340 (1992). [CrossRef]

**70. **O. P. Bruno and F. Reitich, “Numerical solution of diffraction problems: a method of variation of boundaries,” J. Opt. Soc. Am. A **10**, 1168–1175 (1993). [CrossRef]

**71. **O. P. Bruno and F. Reitich, “Numerical solution of diffraction problems: a method of variation of boundaries. II. Finitely conducting gratings, Padé approximants, and singularities,” J. Opt. Soc. Am. A **10**, 2307–2316 (1993). [CrossRef]

**72. **O. P. Bruno and F. Reitich, “Numerical solution of diffraction problems: a method of variation of boundaries. III. Doubly periodic gratings,” J. Opt. Soc. Am. A **10**, 2551–2562 (1993). [CrossRef]

**73. **O. Bruno and F. Reitich, “High-order boundary perturbation methods,” in *Mathematical Modeling in Optical Science*, G. Bao, L. Cowsar, and W. Masters, eds., Vol. 22 of Frontiers in Applied Mathematics (SIAM, 2001), Chap. 3, pp. 71–109.

**74. **L. Kazandjian, “A discussion of the properties of the Rayleigh perturbative solution in diffraction theory,” Wave Motion **42**, 169–176 (2005). [CrossRef]

**75. **J.-J. Greffet, “Introduction to surface plasmon theory,” in *Plasmonics: From Basics to Advanced Topics*, S. Enoch and N. Bonod, eds., Vol. 167 of Springer Series in Optical Sciences (Springer, 2012), Chap. 4, pp. 105–148.

**76. **D. Maystre, “Theory of Wood’s anomalies,” in *Plasmonics: From Basics to Advanced Topics*, S. Enoch and N. Bonod, eds., Vol. 167 of Springer Series in Optical Sciences (Springer, 2012), Chap. 2, pp. 39–83.

**77. **E. Kretschmann and H. Raether, “Radiative decay of non radiative surface plasmons excited by light,” Z. Naturforsch. A **23**, 2135–2136 (1968).

**78. **A. Otto, “Excitation of nonradiative surface plasma waves in silver by the method of frustrated total reflection,” Z. Phys. **216**, 398–410 (1968). [CrossRef]

**79. **B. Hecht, H. Bielefeldt, L. Novotny, Y. Inouye, and D. W. Pohl, “Local excitation, scattering, and interference of surface plasmons,” Phys. Rev. Lett. **77**, 1889–1892 (1996). [CrossRef]

**80. **H. Ditlbacher, J. R. Krenn, N. Felidj, B. Lamprecht, G. Schider, M. Salerno, A. Leitner, and F. R. Aussenegg, “Fluorescence imaging of surface plasmon fields,” Appl. Phys. Lett. **80**, 404–406 (2002). [CrossRef]

**81. **R. H. Ritchie, E. T. Arakawa, J. J. Cowan, and R. N. Hamm, “Surface-plasmon resonance effect in grating diffraction,” Phys. Rev. Lett. **21**, 1530–1533 (1968). [CrossRef]

**82. **J. Uretsky, “The scattering of plane waves from periodic surfaces,” Ann. Phys. **33**, 400–427 (1965). [CrossRef]

**83. **G. A. Baker and P. Graves-Morris, *Padé Approximants*, 2nd ed., Vol. 59 of Encyclopedia of Mathematics and its Applications (Cambridge University, 1996).

**84. **A. Rakic, A. Djurišic, J. Elazar, and M. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. **37**, 5271–5283 (1998). [CrossRef]

**85. **Z. Chen, I. R. Hooper, and J. R. Sambles, “Low dispersion surface plasmon–polaritons on deep silver gratings,” J. Mod. Opt. **53**, 1569–1576 (2006). [CrossRef]

**86. **D. J. Nash and J. R. Sambles, “Surface plasmon–polariton study of the optical dielectric function of silver,” J. Mod. Opt. **43**, 81–91 (1996).

**87. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**, 4370–4379 (1972). [CrossRef]

**88. **D. W. Lynch and W. R. Hunter, “Comments on the optical constants of metals and an introduction to the data for several metals,” in *Handbook of Optical Constants of Solids*, E. D. Palik, ed. (Academic, 1985), Vol. 1, pp. 275–367.

**89. **J. Homola, *Surface Plasmon Resonance Based Sensors*, Vol. 4 of Springer Series on Chemical Sensors and Biosensors (Springer, 2006).

**90. **Z. Chen, “Grating coupled surface plasmons in metallic structures,” Ph.D. dissertation (University of Exeter, 2007).

**91. **T. López-Rios, D. Mendoza, F. J. García-Vidal, J. Sánchez-Dehesa, and B. Pannetier, “Surface shape resonances in lamellar metallic gratings,” Phys. Rev. Lett. **81**, 665–668 (1998). [CrossRef]

**92. **F. García-Vidal, J. Sánchez-Dehesa, A. Dechelette, E. Bustarret, T. López-Rios, T. Fournier, and B. Pannetier, “Localized surface plasmons in lamellar metallic gratings,” J. Lightwave Technol. **17**, 2191–2195 (1999). [CrossRef]

**93. **I. R. Hooper and J. R. Sambles, “Dispersion of surface plasmon polaritons on short-pitch metal gratings,” Phys. Rev. B **65**, 165432 (2002). [CrossRef]

**94. **E. P. Da Silva, G. A. Farias, and A. A. Maradudin, “Analysis of three theories of scattering of electromagnetic radiation by gratings,” J. Opt. Soc. Am. A **4**, 2022–2024 (1987). [CrossRef]

**95. **A. J. Jerri, ed., *Advances in the Gibbs Phenomenon* (Σ Sampling, 2011).

**96. **W. L. Barnes, T. W. Preist, S. C. Kitson, J. R. Sambles, N. P. K. Cotter, and D. J. Nash, “Photonic gaps in the dispersion of surface plasmons on gratings,” Phys. Rev. B **51**, 11164–11167 (1995). [CrossRef]

**97. **A. Taflove and S. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain method* (Artech House, 2005).

**98. **A. Malcolm and D. P. Nicholls, “A field expansions method for scattering by periodic multilayered media,” J. Acoust. Soc. Am. **129**, 1783–1793 (2011). [CrossRef]