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,

 

Fig. 1. Generic grating configuration for SPR coupling.

Download Full Size | PDF

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,

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

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

 

Fig. 2. (a) (Normalized) reflectivity [“NR” in Eq. (16)] map for a sinusoidal silver profile with period $d=650\mathrm{nm}$ at normal incidence in TM polarization; (b) poles (open markers) and zeros (filled markers) of (the [32/32] Padé approximant of) the coefficient $B_0(h)$ for complex $h$ and different values of the wavelength of radiation $\lambda$ ($660\le \lambda \le 700$). Note that no poles are present in a complex neighborhood of the real line (shaded area). The zeros indicated with an arrow correspond to close-to-zero reflectivity, indicating the presence of a plasmon resonance.

Download Full Size | PDF

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

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\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\mathrm{nm}$ and $h=6\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.

 

Fig. 3. (a) Square of the absolute value ${|B_{-1}^{+}|}^2$ of the $-1$st order Rayleigh coefficient $B_{-1}^{+}$ in Eq. (10) (“reflectivity at order $-1$”) for a sinusoidal silver grating of pitch $d=400\mathrm{nm}$ and height $h=4\mathrm{nm}$, illuminated with a plane wave of energy 2.26 eV, corresponding to a wavelength $\lambda =548.60\mathrm{nm}$, as a function of incidence angle; (b) same as in (a) but with height $h=6\mathrm{nm}$; (c) same as in (a) but with height $h=10\mathrm{nm}$.

Download Full Size | PDF

 

Fig. 4. Logarithm (base 10) of the maximum relative error in ${|B_{-1}^{+}|}^2$ (circles) and best linear fit (dashed line) corresponding to incidences in the interval [19.2°,19.4°] for the case in Fig. 3(c) ($h=10\mathrm{nm}$). The optical constants used in the simulations are from [87].

Download Full Size | PDF

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\mathrm{nm}$ and height $h=124\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\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.

 

Fig. 5. (a) Configuration corresponding to a deep silver grating of height 124 nm and period 258 nm illuminated with a plane wave with an incidence angle of 30°; (b) reflectivity at order 0, ${|B_0|}^2$ [Eq. (10)], as a function of wavelength for a polynomial permittivity [85] as computed from a [32/32] Padé approximant. The dashed vertical line marks the wavelength below which the $-1$st Rayleigh mode becomes propagating.

Download Full Size | PDF

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\mathrm{nm}$ and $\lambda =495\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.

 

Fig. 6. (a) Zeroth order reflectivity as a function of height $h$ for the example in Fig. 5 at a wavelength $\lambda =410\mathrm{nm}$; the solid lines display the values of the Padé approximations corresponding to each value of $N$ for which the Taylor partial sums are shown. The vertical line at $h=h_p$ depicts 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$). (b) Poles and zeros in the complex $h$ plane of the [32/32] Padé approximant for the configuration in (a); the circles and crosses (square and stars) mark the location of the zeros and poles, respectively, that do (do not) cancel out. (c) Same as (a) but for a wavelength $\lambda =495\mathrm{nm}$. (d) Same as (b) for the configuration in (c).

Download Full Size | PDF

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\mathrm{nm}$ and periodized with period $d=630\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\mathrm{nm}$; the maximum slope, on the other hand is approximately 32 for this height, and 50 for $h=30\mathrm{nm}$. The results in (b)–(d) correspond to slight changes in the incidence angle, from $\theta =0°$ (normal incidence), to $\theta =0.2°$ to $\theta =1°$ 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\mathrm{nm}$ and $h=30\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.

 

Fig. 7. (Total) Reflectivity for an elliptical grating configuration, as shown in (a) (period $d=630\mathrm{nm}$, linewidth $w=400\mathrm{nm}$), for different incidence angles. The vertical lines mark the wavelength below which the $-1$st order becomes propagating.

Download Full Size | PDF

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

 

Fig. 8. Left- (horizontal lines) and right- (curves) hand sides of Eq. (28). Their intersections, marked by vertical lines, define the resonant wavelengths that would be excited for a silver interface of zero height and period $d=630$, as depicted in (a), under the dispersion derived from (cubic-spline) interpolation of the results in [88].

Download Full Size | PDF

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.

 

Fig. 9. Results for a sinusoidal profile of period $d=258\mathrm{nm}$ illuminated with an incidence angle $\theta =30°$. (a) Total reflectivity as a function of height and wavelength; the thick dashed vertical line marks the wavelength below which the $-1$st order becomes propagating. (b) Spatial distribution of the (transverse) reflected magnetic field, $|u^{+,\text{scat}}|=|H_z^{+,\text{scat}}|$ [Eq. (3)], corresponding to the intersection of the dotted–dashed lines in (a): height $h=20\mathrm{nm}$ and wavelength $\lambda =410\mathrm{nm}$.

Download Full Size | PDF

 

Fig. 10. Same as Fig. 9 for the configuration of Fig. 7(a) with incidence angle $\theta =1°$ [Fig. 7(d)]. (a) Reflectivity map; (b) spatial (transverse) magnetic field distribution corresponding to the intersection of the dotted–dashed lines in (a): height $h=30\mathrm{nm}$ and wavelength $\lambda =640\mathrm{nm}$.

Download Full Size | PDF

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\mathrm{nm}$ and a wavelength $\lambda =410\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

 

Fig. 11. (a) A purely sinusoidal profile with amplitude $a=5\mathrm{nm}$; and (b) a modulated profile with a second harmonic component as in Eq. (29) ($a=5\mathrm{nm}$, $b=2\mathrm{nm}$).

Download Full Size | PDF

 

Fig. 12. Reflectivities as a function of incidence angle and wavelength for: (a) the purely sinusoidal profile in Fig. 11(a); and (b) for the modulated profile Fig. 11(b). In both cases the simulations correspond to the HOPS approach using a truncated Taylor series of order $N=5$. The white area in the middle of the latter figure (corresponding to values that exceed the maximum in the plotted scale) indicates the region where the five-term Taylor series is already inaccurate.

Download Full Size | PDF

 

Fig. 13. Reflectivities as a function of incidence angle and wavelength for: (a) the purely sinusoidal profile in Fig. 11(a); and (b) a modulated profile with a second harmonic component [Eq. (29)] as depicted in Fig. 11(b). For the latter, the line $\lambda =\overline{\lambda }$ marks the approximate location of the center of the band gap as predicted in Eq. (3b) in [96]. In both cases the simulations correspond to the HOPS approach using [2/2] Padé approximants ($N=5$).

Download Full Size | PDF

 

Fig. 14. (a) Numerical convergence analysis (and best linear fit) for the examples in Figs. 12(a) and 13(a); and (b) same for those in Figs. 12(b) and 13(b). Logarithm of the relative maximum error in the range of wavelengths between 635 and 670 nm, and angles between $-1.5°$ and 1.5°. Note the slow convergence of the Taylor series in (a), and its divergence in (b).

Download Full Size | PDF

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

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\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–60\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}$.

 

Fig. 15. (a) Reflectivity for a sinusoidal profile [$b=0$ in Eq. (29)] with period $d=634\mathrm{nm}$ and height $h=2a=60\mathrm{nm}$. (b) Reflectivities corresponding to the same profile for varying heights $h=10–60\mathrm{nm}$.

Download Full Size | PDF

 

Fig. 16. (a) Reflectivity for a Fourier grating [as in Eq. (29)] with period $d=634\mathrm{nm}$, height $h=60\mathrm{nm}$, and $a/b=2.5$. (b) Reflectivities corresponding to the same profile for varying heights $h=10–60\mathrm{nm}$.

Download Full Size | PDF

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.