## Abstract

We present a detailed description of a computationally efficient, semi-analytical method (SAM) to calculate the electomagnetic field distribution in a 1D-periodic, subwavelength-structured metal film placed between dielectric substrates. The method is roughly three orders of magnitude faster than the finite-element method (FEM). SAM is used to study the resonant transmission of light through nanoplasmonic structures, and to analyze the role of fundamental and higher-order Bloch surface plasmons in transmission enhancement. The method is also suitable for solving the eigenvalue problem and finding modes of the structure. Results obtained with SAM, FEM, and the finite-difference time-domain method show very good agreement for various parameters of the structure.

© 2008 Optical Society of America

## 1. Introduction

A surface plasmon polariton (SPP) is a fundamental electromagnetic excitation which is localized at the interface between media with positive (dielectric) and negative (conductor) permittivity [1, 2]. Excitation of the SPP is a resonant phenomenon. It is very sensitive to variations of geometrical and material parameters of a structure prompting application of SPPs in bio- and chemical sensing [3]. Recent advances in nanofabrication have sparked enormous interest in surface waves on the nanoscale. The term “plasmonics” [4] (or “nanoplasmonics” [5] for structures with a characteristic size of ~100 nm) has been recently introduced and collectively refers to various optical phenomena in subwavelength-size metal-dielectric media, such as metal films with subwavelength slits or apertures, nanoresonators, arrays of metal nanoparticles, or nanoantennas.

The subwavelength periodicity in a metallic medium may be sufficient to induce localized excitations on the surface. This is due to the excitation of localized Bloch modes, which are counterparts of SPPs in periodic media. However, unlike SPPs, Bloch surface plasmons (BSPs) [6] may have zero momentum, and therefore will not require a prism for their excitation. Similar to periodic dielectric structures, such as Bragg gratings and photonic crystals, periodicity in plasmonic structures is responsible for the emergence of novel optical effects. For example, BSPs were shown to contribute to resonantly enhanced transmission through a periodically structured metal film [7]–[14]. Enhanced transmission plays an important role in photonic applications, because being a resonant phenomenon, it is very sensitive to variation of optical properties of the dielectric substrates. Periodic nanoplasmonic structures promise exciting applications in nonlinear optics, silicon photonics, subwavelength lithography, sensing, and near-field microscopy [6, 15].

Modeling plays an ever increasing role in studying periodic nanoplasmonic structures. Widely used space- and time-domain numerical methods require large computation time or/and significant memory resources because of nanometer feature sizes, but typically micrometer lengths along other dimensions. The situation is particularly aggravated if higher-order (HO) plasmonic modes, which need a finer computational grid for accurate resolution, are excited [16]. Therefore, methods such as FEM [17] or finite-difference time-domain (FDTD) [18], typically require at least a multi-processor workstation. These methods also may be prone to numerical artifacts, and extra care must be taken in proper application of boundary conditions. A fast, reliable, and easy to implement method is therefore highly desirable for studying and optimizing nanophotonic structures.

Due to the complexity of the problem, purely analytical methods are usually restricted to a perturbation-like periodicity. For example, an analytical approach proposed in [10] is limited to a weak, sinusoidal modulation of the film’s permittivity. A solution to Maxwell’s equations was only found for a symmetric structure with identical substrate and superstrate. This approach could not be extended to the case of, e.g., step-like modulation of the film. Within the framework of Fourier modal methods (FMMs) [19]–[23], one can use a periodic ansatz for the permittivity of the film, and consequently for the electromagnetic field in the whole structure. Boundary conditions provide a set of linear algebraic equations for harmonic amplitudes of the field components. For a small-amplitude, single-harmonic (i.e. sinusoidal) modulation, one can calculate resonant frequencies of the structure analytically [9]. The same approach was then applied to calculate the transmittance of a thin metal film with the same type of periodicity [24]. Small-amplitude, single-harmonic modulation was a clear limitation of the approach, since high-contrast, step-profile metal films and slit arrays, or even a single-harmonic phase grating with a large modulation depth, require a much larger number of harmonic terms *N* in the field expansion. The study of HO BSPs also calls for a more complete field expansion. An extension of the method [24] to *N*>1 harmonics has been conceptually outlined in [25], but no implementation of the method has been proposed. As a result, all reported calculations have been performed for the case of *N*=1, raising concerns over applicability of the method to true lamellar gratings [26].

In this paper, we discuss in detail a generalization of the method that was introduced in [24] to an arbitrary number of field harmonics. We also propose a robust, computationally efficient implementation of the algorithm. Our approach can be viewed as a variation of FMMs (sometimes called the rigorous coupled-wave analysis - RCWA) with several numerical steps replaced by analytical calculations. Only the polynomial eigenvalue problem and two sets of linear algebraic equations are solved numerically. All other calculations are performed analytically. We thus term this method as the semi-analytical method (SAM). Our computations are orders of magnitude faster than FEM and FDTD. Due to reduced number of matrix operations we expect the SAMbe faster than the RCWA. The SAMallows for obtaining the resonant frequency of the fundamental and the HO BSP, the electromagnetic field distribution in the structure, and transmission and reflection coefficients. It also gives a clear physical interpretation of the scattering process, revealing the role of the BSPs in the resonant scattering. We give a step-by-step recipe for numerical implementation of the algorithm. With minimum programming efforts, it should allow an efficient study of BSPs in a certain class of nanoplasmonic structures. The SAM can also help verify convergence of other numerical techniques.

The paper is organized as follows. In Section 2 we describe the structure that we studied. The details of the SAM and instructions on its numerical implementation are given in Section 3. Fundamental and HO BSPs and associated resonant transmission through the structure are discussed in Section 4. Here, we also compare results obtained with SAM, FEM, and FDTD for various parameters of the structure, demonstrating the high accuracy of our method. Section 5 contains concluding remarks.

## 2. 1D-periodic nanoplasmonic structure

As a practical example of the grating, we assume a rectangular profile of the complex-valued permittivity (Fig. 1)

with period Λ, duty cycle *ρ*, thickness *h*, and the contrast Δ*ε*=ℜ(*ε*
_{T} - *ε*
_{B})>0. In this work we mostly focus on planar bimetallic gratings with ℜ(*ε*
_{T})<0. However, we also present some encouraging results when applying the SAM to a true slit geometry with *ℜ*
_{T}=1. Unlike metal-dielectric (true slit) structures, where both evanescent and waveguiding modes contribute to the transmitted light [7, 13, 14], bimetallic gratings only support evanescent (plasmonic) modes. As a result, much narrower resonant peaks (especially for small Δ*ε*) appear in the transmitted spectrum as compared to the spectra of slit arrays. This feature is essential for nonlinear optical and sensing applications. The described modulation of the metal’s permittivity can be achieved, e.g., by alternating stripes of two metals that are optically different in some wavelength range, such as gold and silver in the visible. For other possible combinations, see Fig. 2. Other options to create bimetallic gratings may include geometric structuring to modify the effective permittivity of a stripe and using alloys or metal-dielectric composites having a prescribed effective permittivity [27].

Several closed-form analytic expressions for dispersion in metals are available (see, e.g., [28]–[30]). Since our results were found to be independent of details of the dispersion model, without loss of generality, both “steps” of the grating are assumed to have the same functional dependence

where λ is the wavelength of the incoming light. For example, the plasma wavelength λ_{p}=145 nm and *ε*
^{∞}
_{B}=1.53 correspond to gold [29] and both real and imaginary parts of permittivity depend on the wavelength. To observe transmission features associated with the BSPs more clearly, we assume artificially low loss in both segments of the grating. In our simulations in Section 4, the loss parameter γ_{p} varies between 10^{7} nm and 5×10^{5} nm, where the latter value corresponds to ℑ(*ε*
_{T,B})≈0.1 at 900 nm. We also show several results for *γ*
_{p}=4×10^{4} nm, which gives realistic loss in metal.

To keep the length of our description reasonable, we restrict ourselves to a normal incidence of light (Fig. 1). The results for the oblique incidence will be discussed elsewhere. In what follows, roman subscripts “x”, “y”, and “z” denote the respective component of the electromagnetic field, while the running indices are italicized. Superscripts “+” and “-” refer to the superstrate and the substrate, respectively. Variables denoting matrices and vectors are typed in bold face.

## 3. The semi-analytical method: details

Below we introduce the equations for the electromagnetic field components. This is followed by analysis of the general solution for the field in the structured metal film, the superstrate, and the substrate. We show how boundary conditions lead to a system of linear algebraic equations for the unknown coefficients of the periodic ansatz, and give explicit formulas for the distribution of the electromagnetic field in the structure. Also derived are the expressions for transmission and reflection coefficients. The key computational steps of the algorithm are listed in the end of this section.

#### 3.1. Maxwell’s equations for the TM-mode and the periodic ansatz

The periodic permittivity of the metal film can be represented by a truncated Fourier series

where *ε*
_{0} is the average permittivity of the film (not to be confused with the free-space permittivity) and *g*=2*π*/Λ is the magnitude of the wave vector of the grating. Without loss of generality, we only consider the even harmonics, since the origin of *x*-axis can be chosen arbitrarily (Fig. 1). For the stepwise function (1), the Fourier coefficients are

Since our main focus is on low-loss metals, for bimetallic gratings it is sufficient to take ℑ(*ε*
_{0}) as the average loss of the film and set ℑ(*ε _{m}*)=0. A TE mode does not participate in the resonant interaction with the structure [24]. Therefore, we only consider a TM plane wave with the field components

**E**=(

*E*

_{x},0,

*E*)

_{z}*e*+c.c. and

^{-iwt}**H**=(0,

*H*

_{y},0)

*e*+c.c. After substitution in Maxwell’s equations

_{-iwt}they give a system of two coupled PDEs for the electric field components *E*
_{x} and *E*
_{z}

where *K*=2*π*/λ. The *H*
_{y} component of the field is related to *E*
_{x} and *E*
_{z} by

where *Z*=376.6 Ohms is the free space impedance.

We are looking for a solution of (5) as a superposition of harmonics of the structure’s periodicity. Therefore, a general ansatz for the solution inside the film takes the form

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}+\left({a}_{\sigma}+\sum _{n=1}^{N}\left[{p}_{\sigma n}\mathrm{cos}\left(\mathrm{ngx}\right)+{q}_{\sigma n}\mathrm{sin}\left(\mathrm{ngx}\right)\right]\right){e}^{-k\left(z+h\right)},$$

where σ is either x or z. The first and the second term in (7) show the contribution of the top (*z*=0) and the bottom (*z*=*-h*) interfaces, respectively. The unknown coefficients *A*
_{σ}, *P*
_{σn}, *Q*
_{σn}, *a*
_{σ}, *p*
_{σn}, *q*
_{σn} and values of the wave parameter *k* are to be found from (5) and from the electromagnetic boundary conditions. For consistency, the number of harmonics retained in the ansatz (7) should be no less than the number of harmonics in the representation of the modulated medium (3), i.e. *N*≥*M*. Unless otherwise stated, we assume *M*=*N*. For the superstrate (substrate), only the first (second) term in (7) represents the general solution.

#### 3.2. Solution in the structured metal film, -h<z<0

Substituting (7) into (5) and equating coefficients of cos(*ngx*), sin(*ngx*), and the free terms, we arrive at two decoupled systems of linear algebraic equations. In a matrix form, these systems can be written as

where **U**
_{N} and **V**
_{N} are square matrices of size 2*N*+1, ** ξ**=(

*A*

_{x},

*P*

_{x1},

*Q*

_{z1},

*P*

_{x2},

*Q*

_{z2}, …,

*P*,

_{xN}*Q*)

_{zN}^{T},

**ζ**=(

*A*

_{z},

*P*

_{z1},

*Q*

_{x1},

*P*

_{z2},

*Q*

_{x2}, …,

*P*

_{zN},

*Q*

_{xN})

^{T}. The analysis is simplified for the normal incidence of light because the incoming plane wave lacks any

*z*-component of the field. Therefore, none of the components

*A*

_{z},

*P*

_{zn}, and

*Q*

_{xn}are excited within the structure, resulting in

which reduces the the number of unknown coefficients in (7) by half.

Matrix **U**
_{N} is a sum of two matrices **Y**
_{N} and **u**
_{N}:

which originate, respectively, from the product of *ε*(*x*) with the field, and from the terms having the second-order mixed derivatives in Maxwell’s Eq. (5).

The nonzero elements of **Y**
_{N} are

where *n*,*m*=1,2…,*N* and

The upper and lower sign in (13) applies to *f*
^{A}
_{m} and *f*
^{B}
_{m}, respectively; Θ(*x*)=1 for *x*≥0 and Θ(*x*)=0 otherwise. The other matrix in (11), **u**
_{N}, has the following nonzero elements

For example, for *N*=2 we obtain
${\mathbf{U}}_{2}=\left(\begin{array}{ccccc}{k}^{2}+{\epsilon}_{0}{K}^{2}& {\epsilon}_{1}{K}^{2}& 0& {\epsilon}_{2}{K}^{2}& 0\\ 2{\epsilon}_{1}{K}^{2}& {k}^{2}+\left({\epsilon}_{0}+{\epsilon}_{2}\right){K}^{2}& -\mathrm{kg}& {\epsilon}_{1}{K}^{2}& 0\\ 0& \mathrm{kg}& -{g}^{2}+\left({\epsilon}_{0}-{\epsilon}_{2}\right){K}^{2}& 0& {\epsilon}_{1}{K}^{2}\\ 2{\epsilon}_{2}{K}^{2}& {\epsilon}_{1}{K}^{2}& 0& {k}^{2}+{\epsilon}_{0}{K}^{2}& -2\mathrm{kg}\\ 0& 0& {\epsilon}_{1}{K}^{2}& 2\mathrm{kg}& -4{g}^{2}+{\epsilon}_{0}{K}^{2}\end{array}\right),$
For Eq. (8) to have a solution, condition

should be satisfied. This equation determines the allowable values of the wave number *k* along the *z*-axis, resulting in an (*N*+1)th order polynomial for *k*
^{2}, whose zeros are ±*k _{j}*,

*j*=1, …,

*N*+1. For a more efficient numerical treatment, one can recast (15) in the form

which is known as the quadratic eigenvalue problem. As follows from (11), (12), and (14), the only nonzero elements of matrices **A** and **B** are *a*
_{11}=*a*
_{2n,2n}=1 and *b*
_{2n+1,2n}=-*b*
_{2n,2n+1}=*g*, respectively. Matrix **C** equals *K*
^{2}
**Y**
_{N} with some extra terms added on the main diagonal, namely *c*
_{2n+1,2n+1}=*K*
^{2}
*y*
_{2n+1,2n+1}-*n*
^{2}
*g*
^{2}. An elaborate discussion about the methods for solving (16) can be found in [31]. As a ready implementation, one can use, for example, function polyeig in MATLAB^{®} [32].

Each eigenvalue *k*=±*κ _{j}* of

**U**

_{N}corresponds to an independent solution of (7). A linear superposition of all solutions gives the general form of the electric field distribution in the film. For a bounded medium (-

*h*<

*z*<0), any sign of the eigenvalue results in an acceptable solution. We keep the “+” sign for ℜ(

*κ*)>0 and, accounting for (10), obtain from (7)

_{j}The superscript in parentheses refers to a particular eigensolution, while the subscript *n* shows the order of the harmonic. The amplitudes *A*
^{(j)}
_{x}, *a*
^{(j)}
_{x}, *P*
^{(j)}
_{xn}, *p*
^{(j)}
_{xn}, *Q*
^{(j)}
_{zn}, *q*
^{(j)}
_{zn} are to be determined by the boundary conditions.

#### 3.3. Solution in the superstrate, z>0

By taking the first term in (7) with a reduction imposed by (10), one can look for the solution in the superstrate *z*>0. Matrix **U**
_{N} is considerably simplified in the absence of periodic variation of permittivity in the dielectric medium. With *ε*
_{0} replaced by *ε*
^{+} and any other *ε _{m}* ≡ 0, it becomes block-structured, i.e. it is reduced to

*N*+1 independent equations for each harmonic. The corresponding eigenvalues are given by

for *n*=0,1, …,*N*. The zero-order eigenvalue -*ε*
^{+}
*K*
^{2} corresponds to plane wave solutions for the *E*
^{+}
_{x} component of the field, namely
${R}_{x}{e}^{-i\sqrt{{\epsilon}^{+}}\mathrm{Kz}}$
as a reflected wave and
${A}_{x}^{\mathrm{in}}{e}^{i\sqrt{{\epsilon}^{+}}\mathrm{Kz}}$
as an incident wave. The electric field in the superstrate can then be written as

with *R*
_{x}, *P*
^{+}
_{xn}, and *Q*
^{+}
_{zn} being the unknown amplitudes to be determined by the boundary conditions.

In (19), *s*±=±1 is the sign of the eigenvalue. It must be chosen properly to obtain physically meaningful solutions. For a lossless dielectric, *n*
^{2}
*g*
^{2}-*ε*
^{+}
*K*
^{2} can be either positive (evanescent harmonic) or negative (propagating harmonic). In the superstrate, the amplitude of the evanescent wave should decay for *z*→+∞, which imposes *s*+=-1. To ensure the correct sign of the phase for reflected harmonics (i.e. make the phase front propagate in the direction of the wave vector), one has to take *s*+=1 if *n*
^{2}
*g*
^{2}-*ε*^{+}*K*
^{2}<0. Therefore, for the purely real *ε*
^{+}

This alternation of the sign *s*
^{±} occurs due to discontinuity of the square root function across its branch-cut along the negative real axis. If the loss in the dielectric is taken into account, i.e. ℑ(*ε*
^{+})>0, then (*η*
^{+}
_{n})^{2} is complex-valued. For this case, *s*
^{+}= -1 for both propagating and evanescent waves will result in the correct behavior of both evanescent and propagating harmonics (Fig. 3, left plot).

Substitution of *k*=*η _{n}* in matrix

**U**

_{N}for the superstrate, with a subsequent solution of (15), gives a relation between the harmonic amplitudes

#### 3.4. Solution in the substrate, z<-h

By following the same steps as for derivation of (18), we obtain the eigenvalues

and the solution for the field components in the substrate as

where *T _{x}* is the amplitude of the zero-order diffracted wave. For a lossless dielectric, i.e. with ℑ(

*ε*

^{-})=0 we have

where the same reasoning as in deriving (20) was used. If material loss is accounted for in the substrate, then *s*
^{-} ≡1 (Fig. 3, right plot). Just as in (21), the following relation between the harmonic amplitudes holds for the substrate

#### 3.5. Boundary conditions

We have reduced a generic form (7) of the solution to Maxwell’s equations to more specific expressions, namely (17), (19), and (23). We have also determined a discrete set of wave numbers *κ _{j}* (for the film) and

*η*

^{±}

_{n}(for dielectric media). Next, we apply boundary conditions to uniquely determine the unknown amplitudes. The electromagnetic boundary conditions at

*z*=0 are

Applying the first boundary condition (26a) to (17a) and (19a), we obtain the following relations for the unknown amplitudes

The second boundary condition (26b) combined with (17b) and (19b), results in

As discussed in [9], only one equation (for the zero-order term) is needed from the third boundary condition (26c)

A similar set of equations can be obtained from the boundary conditions at *z*=-*h*:

#### 3.6. The main matrix equation

The number of the unknown harmonic amplitudes can be reduced by using the compatibility of system (8). For example, one can use the eigenvectors of matrix **U**
_{N} as follows. Because of (15), the equations of (8) are linearly dependent. Therefore, one equation in (8), e.g., first, can be dropped. Then, the remaining system can be solved to find the harmonic amplitudes *P*
^{(j)}
_{xn} and *Q*
^{(j)}
_{zn} expressed through the zero-order amplitude *A*
^{(j)}
_{x}. In other words, we define *α _{nj}* and

*β*as

_{nj}and rewrite (8) as

where the right-hand side vector **r**=-2*K*
^{2}(*ε*
_{1}, 0, *ε*
_{2}, 0, …, *ε _{N}*, 0)

^{T}is obtained from the first column of matrix

**U**

_{N},

**=(**

*ρ′**α*

_{1j},

*β*

_{1j},

*α*

_{2j},

*β*

_{2j}, …,

*α*,

_{Nj}*β*)

_{Nj}^{T}and

**W**

^{(j)}

_{N}(

*κ*) is obtained from

_{j}**U**

_{N}by removing the first row and the first column and substituting

*k*=

*k*. For example, for

_{j}*N*=2 ${\mathbf{W}}_{2}^{\left(j\right)}\left({\kappa}_{j}\right)=\left(\begin{array}{cccc}{\kappa}_{j}^{2}+\left({\epsilon}_{0}+{\epsilon}_{2}\right){K}^{2}& -{\kappa}_{j}g& {\epsilon}_{1}{K}^{2}& 0\\ {\kappa}_{j}g& -{g}^{2}+\left({\epsilon}_{0}-{\epsilon}_{2}\right){K}^{2}& 0& {\epsilon}_{1}{K}^{2}\\ {\epsilon}_{1}{K}^{2}& 0& {\kappa}_{j}^{2}+{\epsilon}_{0}{K}^{2}& -2{\kappa}_{j}g\\ 0& {\epsilon}_{1}{K}^{2}& 2{\kappa}_{j}g& -{4g}^{2}+{\epsilon}_{0}{K}^{2}\end{array}\right).$

Similarly, we can replace *κ _{j}* in (34) with -

*κ*to solve the system for

_{j}**ρ″**=(

*$\tilde{\alpha}$*

_{1j},

*$\tilde{\beta}$*

_{1j},

*$\tilde{\alpha}$*

_{2j},

*$\tilde{\beta}$*_{2j}, …,

*$\tilde{\alpha}$*,

_{Nj}*$\tilde{\beta}$*)

_{Nj}^{T}where

Now, when ** ρ′** and

**are known, we can express all unknown amplitudes in (17) through 2**

*ρ″**N*+2 zero-order harmonics that form a vector

Combining (27b), (28), and (21) we obtain equality

where we have also used definitions (33) and (35). Formula (37) represents *N* equations (for each *n*), where unknown are the elements of vector **a**. A similar relation can be obtained from (30b), (31), and (25), which gives *N* additional equations for the total of 2*N*+2 elements of vector **a**. To complete the set, the remaining two equations are obtained by substituting *T*
_{x} from (30a) in (32) to obtain

and by substituting *R*
_{x} from (27a) in (29) to obtain

The resulting system of 2*N*+2 equations for the elements of vector **a** can be written in matrix form as

where the only nonzero element of the right-hand side vector **v** is
${v}_{2N+2}=2\sqrt{{\epsilon}^{+}}K{A}_{x}^{\mathrm{in}}$
and matrix **M** is

The dimension of each of four submatrices **m** is *N*×*N*+1, with *p* and *q* being the row and the column indices, respectively, *p*=1, …, *N*, *q*=1, …, *N*+1. The submatrices have the following elements

Equations (42a) and (42b) originate from (37). They are obtained by collecting coefficients of *A*
^{(j)}
_{x} and of *a*
^{(j)}
_{x}, respectively. The last two rows of (41) originate from (38) and (39), respectively. Numerical solution of (40) essentially completes the computational algorithm of the SAM.

#### 3.7. Electromagnetic field distribution

Once elements of vector **a** are known, one can easily calculate the electromagnetic field for the whole structure. As follows from (17), (33), and (35), the field components in the metal film (-*h*≤*z*≤0) expressed through the elements of vector **a** are

The field components in the superstrate (*z*>0) are then

where we have used (19), (27b), (21) together with (33), (35) and *R*
_{x} is given by (27a). Similarly, (23), (30b), (25) result in the field components in the substrate (*z*<-*h*):

where *T*
_{x} is given by (30a). The magnetic field components *H*
_{y} are obtained by differentiation of (43), (44), and (45) according to (6).

#### 3.8. Transmission and reflection coefficients

If the wavelength of the incoming light satisfies the condition

and the superstrate and substrate are lossless, only the zeroth-order, plane wave harmonics *T*
_{x} and *R*
_{x} contribute to the total transmitted and reflected power. Because *η*
^{±}
_{n} in (19) and (23) are purely imaginary for all *n*, all higher-order harmonics are evanescent and do not contribute to the total field at *z*=±∞. In this case, the zeroth-order transmitted amplitude is given by (30a), while the zeroth-order reflected amplitude is obtained from (27a) as

The corresponding transmission and reflection coefficients are then

If for particular system parameters, higher-order harmonics become propagating, they will start contributing to the far field. Then, transmittance should be calculated as the ratio of the *z*-component of the Poynting vector integrated over the period ∫_{Λ}
**S**
_{z}d*x*|_{z=-∞} to the same quantity that is calculated for the incident wave at *z*=+∞.

Upon substitution of (45) and (6) in the expression for the *z*-component of the Poynting vector
${S}_{z}=\frac{1}{2}\Re \left({E}_{x}{H}_{y}^{*}\right)$
we arrive at the expression for the transmission coefficient of the structure

where condition *ε*-≥(*n*λ/Λ)^{2} is satisfied for *n*=1,2,…,*L*, i.e. *L*≥1 harmonics in the substrate are propagating. In obtaining (49), we have used that ∫_{Λ}Σ_{n} cos(*ngx*)Σ_{m} cos(*mgx*)d*x* equals Λ/2 if *m*=*n* and zero if *m*≠*n*. By using (25) and (22), formula (49) can be simplified to

where *T*
_{x} and *P*
^{-}
_{xn} are given by (30a) and (30b), respectively. Analogously, we obtain the reflection coefficient as

where *R*
_{x} and *P*
^{+}
_{xn} are given by (47) and (27b), respectively.

#### 3.9. Computation steps of the SAM: a summary

As a short summary, we list the key steps one follows to calculate the electromagnetic field in the structure along with its transmission and reflection coefficients:

## 4. Results and discussion

First, we demonstrate several test cases to verify that calculations with the SAM agree with results of the FEM and the FDTD. A commercial code from COMSOL^{®} [33] was used as FEM implementation. It has been run on a 2-processor 4 GB workstation. FDTD calculations were performed on an 8-processor, 32 GB workstation while the SAM code runs on a standard laptop with 1 GB memory and 1.8 GHz processor. This section also features fundamental and higher-order Bloch surface plasmons and the SAM modeling of a true slit grating.

#### 4.1. Validation of the SAM with FEM and FDTD

Figure 4 shows the results of all three methods for bimetallic gratings of various thicknesses. The parameters of the grating and the range of wavelengths are chosen to ensure excitation of the first-order or fundamental BSP (FBSP). FBSPs are discussed in more detail in Section 4.2. Excitation of the FBSP corresponds to transmission maxima in Fig. 4. For thicker films, all three methods agree almost perfectly (Fig. 4a). For thinner films, the interaction between even and odd FBSP becomes stronger and regions of highly suppressed transmission appear in the spectrum. These regions correspond to excitation of a linear combination of symmetric and antisymmetric FBSP. Thus, depending on the excitation wavelength and film thickness, FBSPs can both enhance and suppress transmission through the film. The fact that in both scenarios FBSPs are involved was a source of controversy over understanding the role of surface plasmons in enhanced transmission. For example, suppressive role was assigned to surface plasmons in resonant transmission through metallic gratings with narrow slits [34]. In our example, for *h*=72 nm the transmission drops to nearly -60 dB (i.e. 10^{-6}), which is four orders of magnitude below the tunneling rate. Such regimes require a high-resolution computation grid and are quite difficult to resolve, especially with the FDTD method. As a result, the FDTD solutions of Figs. 4(b),(c) are not fully convergent to the solution. Nevertheless, good agreement between the SAM and the FEM is seen in Figs. 4(b)(c), even for very low transmittance levels.

The number of terms *N* one needs to keep in the field expansion depends on the system parameters. We have found that for low contrast bimetallic gratings, *N* can be quite low. In fact, for Δ*ε*=5 just one harmonic is sufficient (Fig. 5). This leads us to a non-intuitive conclusion that the response of a true rectangular but *low-contrast* bimetallic grating with *ε*
_{T,B}<0 can be accurately modeled by response of a sinusoidal (phase) grating. On the contrary, even the phase grating itself (*M*≡1) might require *N*>1 terms in the filed expansion if the amplitude of modulation *ε*
_{1} is high (Fig. 6).

#### 4.2. Fundamental Bloch surface plasmons

As was shown in [10, 9, 24], at some resonant wavelength an FBSP is excited. An approximate value of that wavelength λ^{res}
_{1} can be estimated from the dispersion relation for the SPP on a smooth metal surface. It is, therefore, given by the solution of the transcendental equation λ=Λ(*ε*
_{0}(λ)*ε*
^{±}(λ)/[*ε*
_{0}(λ)+*ε*
^{±}(λ)])^{1/2}. The exact value of λ^{res}
_{1} is determined by the modal equation det **M**=0, which for a relatively thin film has two distinct solutions (Fig. 7(a). The solutions are complex-valued even in the absence of any loss in metal [9]. The two solutions correspond to an odd and an even BSP which are excited if the wavelength of the incident light coincides with the resonant frequency, i.e. λ=ℜ(λ^{res}
_{1}). This excitation is accompanied by a transmission enhancement of the structure (Fig. 7(b). Note that both values of λ^{res}
_{1} exactly coincide with the maxima of transmittance. The nonzero imaginary part of λ^{res}
_{1} reflects the finite lifetime *τ*=ℜ(λ^{res})^{2}/[2*πc*Im(λ^{res})] of each BSP mode. Since the lower-frequency mode (λ=918.2 nm) has a shorter lifetime, the right resonance peak on Fig. 7(b) is broader than the left peak. For this simulation, we have taken a close to realistic loss coefficient in the metal with γ_{p}=4×10^{4} nm, which corresponds to ℑ(ε_{T,B})=0.81. Excellent agreement with the FEM can be seen for the calculated transmittance. On the computational side, when loss in a system is high, with a relatively large size of matrix **M**, the absolute value of det **M** becomes very large making it difficult to accurately find zeros of the eigenvalue equation.

Excitation of FBSPs accompanied by resonantly enhanced light transmission can be clearly seen from the animation FBSP.avi. A snapshop of a movie for a wavelength near an FBSP resonance is shown in Fig. 8. Strong enhancement of the field (especially its *E*
_{z} component) can be observed near both resonant wavelengths. Also notable is the behavior of the Poynting vector **S** in the wavelength region near resonances. There, the the energy flow changes direction along the *x*-axis and vector **S** forms circular areas. If the FBSP is not excited, the nonresonant tunneling becomes the only mechanism of light transmission. It is accompanied by a strong reflection from the film, visible as a pronounced interference pattern for the *E*
_{x} component of the field.

#### 4.3. Higher-order Bloch surface plasmons

According to the Bloch theorem, the same plasmonic structure that reveals the FBSP can support HO modes at the higher resonant frequencies. To observe the effect in the wavelength range with low loss in the metal, one can decrease the grating period or increase the permittivity of the substrate or/and the superstrate. For a grating period of ~500 nm, third-order BSPs can be excited in a silicon-based structure [16]. Excitation of the fourth-order BSPs and corresponding resonantly enhanced transmission of light can be observed by watching animation 4order-BSP.avi. Figure 9 shows a frame of the movie for a wavelength corresponding to resonance with the even fourth-order BSP. The number of lobes of the mode over half a period equals the order number of the BSP. The HO modes are much more tightly localized compared to the FLSP. For example, a single peak of the mode in Fig. 9 measures about 50 nm along the *x*-axis. The light confinement in the transverse direction is therefore ~λ/16, which is quite promising for the subwavelength lithography. The penetration depth of the localized field is even less (~20 nm). To minimize the spot size for a fixed period of the structure, one has to increase the index of the substrates.

Computations related to HO BSPs pose more difficulties because of amuch denser mesh used for proper resolution of the mode. For example, in computing transmittance of the structure near the third order resonance (Fig. 10) we were not able to achieve perfect agreement with FEM because of memory restrictions of our workstation. The FDTD grid size of 0.25 nm was nearly sufficient but required several days of computation time. The corresponding runtime for the FEM code was several hours compared to less than a minute for the SAM.

#### 4.4. The true slit geometry

In Fig. 11, we show results on light transmission through a film with true slits, i.e. for *ε*
_{T}=1 in (1). Except for the regime of nearly complete transparency, the quantitative agreement of FEM and FDTD with the SAM is very good. The wavelength of maximum transmittance can be predicted quite accurately with a more accurate ansatz in the SAM. However, the SAM underestimates the magnitude of maximum transmittance 𝒯_{max} by about 5 dB. The reason for this discrepancy is an incomplete account for loss of metal in the current implementation of the SAM. Taking only the zero-order term in the expansion of the imaginary part of *ε*(*x*) is not sufficient to achieve a perfect result. Further development of the method is needed to improve its stability when dealing with the unequal loss of the grating’s components. Highly conductive gratings under TM incidence are known to suffer from numerical instabilities [21, 20, 22, 23]. A variety of techniques developed to improve the stability of FMMs (see, e.g., [35] for an overview) can be combined with the SAM. For example, the filtering method [36] proposed recently can be used to accurately treat a true slit array within the SAM. We note that the artificially decreased loss coefficient in the metal part of the grating decreases the discrepancy with 𝒯_{max} but gives a slightly offset resonant wavelength. The computation time with the SAM, however, remains within several minutes even for metal-dielectric structures with a small duty cycle.

## 5. Conclusions

We have presented a computationally efficient semi-analytical method for calculating electromagnetic field in a nanoplasmonic structure consisting of a thin, subwavelength-structured metal film sandwiched between two dielectric layers. The method is easy to implement, and it does not require significant computational resources. It is particularly suitable for studying fundamental and higher-order Bloch surface plasmons and their role in resonant transmission of light. It can also be used for optimizing nanoplasmonic structures for a targeted nonlinear-optical or optoelectronic application. Predictions of the SAM agree very well with calculations done with finite-element and FDTD methods, with the SAM being orders of magnitude faster. Bimetallic gratings, i.e. planar metal films with both segments having negative permittivity, are handled by this method very accurately. We have found that a low-contrast, step-like modulation can be accurately represented by just a few harmonics in the Fourier series of the grating profile. Some work needs to be done on the proper description of unequal material loss in different segments of the grating. However, qualitative agreement with FEM and FDTD was achieved even for a true slit, lamellar metal-dielectric structure. This method can be extended to multilayer structures. Extension of the SAM to 2D geometries also seems feasible.

## Acknowledgments

The authors are grateful to Kevin Sparks for useful remarks on the manuscript and to Nikolay Panin for help with FEM simulations of Fig. 11.

## References and links

**1. **V. M. Agranovich and D. L. Mills, eds., *Surface Polaritons* (North Holland, Amsterdam, 1982).

**2. **A. D. Boardman, ed., *Electromagnetic Surface Modes* (Wiley, 1982).

**3. **J. Homola, ed., *Surface Plasmon Resonance Based Sensors* (Springer, 2006). [CrossRef]

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

**5. **S. Kawata and H. Masuhara, eds., *Nanoplasmonics. From Fundamentals to Applications* (Elsevier, 2006).

**6. **A. V. Zayats, I. I. Smolyaninov, and A. A. Maradudin, “Nano-optics of surface plasmon polaritons,” Phys. Rep. **408**, 131–314 (2005). [CrossRef]

**7. **J. A. Porto, F. J. García-Vidal, and J. B. Pendry, “Transmission resonances on metallic gratings with very narrow slits,” Phys. Rev. Lett. **83**, 2845–2848 (1999). [CrossRef]

**8. **A. Krishnan, T. Thio, T. J. Kim, H. Lezec, T. W. Ebbesen, P. A. Wolff, J. Pendry, L. Martin-Moreno, and F. J. Garcia-Vidal, “Evanescently coupled resonance in surface plasmon enhanced transmission,” Opt. Commun. **200**, 1–7 (2001). [CrossRef]

**9. **S. A. Darmanyan and A. V. Zayats, “Light tunneling via resonant surface plasmon polariton states and the enhanced transmission of periodically nanostructured metal films: An analytical study,” Phys. Rev. B **67**, 035424 (2003). [CrossRef]

**10. **A. M. Dykhne, A. K. Sarychev, and V. M. Shalaev, “Resonant transmission through metal films with fabricated and light-induced modulation,” Phys. Rev. B **67**, 195402 (2003). [CrossRef]

**11. **N. Bonod, S. Enoch, L. Li, E. Popov, and M. Nevière, “Resonant optical transmission through thin metallic films with and without holes,” Opt. Express **11**, 482–490 (2003). [CrossRef]

**12. **D. Gérard, L. Salomon, F.
de Fornel, and A. V. Zayats, “Ridge-enhanced optical transmission through a continuous metal film,” Phys. Rev. B **69**, 113405 (2004). [CrossRef]

**13. **Y. Xie, A. R. Zakharian, J. V. Moloney, and M. Mansuripur, “Transmission of light through a periodic array of slits in a thick metallic film,” Opt. Express **13**, 4485–4491 (2005). [CrossRef]

**14. **N. Garcia and M. Nieto-Vesperinas, “Theory of electromagnetic wave transmission through metallic gratings of subwavelength slits,” J. Opt. A: Pure Appl. Opt. **9**, 490–495 (2007). [CrossRef]

**15. **K. Busch, G. von Freymann, S. Linden, S. F. Mingaleev, L. Tkeshelashvili, and M. Wegener, “Periodic nanostructures for photonics,” Phys. Rep. **444**, 101–202 (2007). [CrossRef]

**16. **A. Kobyakov, A. Mafi, A. R. Zakharian, and S. A. Darmanyan, “Fundamental and higher-order Bloch surface plasmons in planar bimetallic gratings on silicon and glass substrates,” J. Opt. Soc. Am. B (submitted).

**17. **J. L. Volakis, A. Chatterjee, and J. L. Kempel, *Finite Element Method for Electromagnetics* (IEEE Press, New York, 1998). [CrossRef]

**18. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas and Prop. **14**, 302–307 (1966). [CrossRef]

**19. **M. G. Moharam, E. B. Grann, D. A. Pommet, and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A **12**, 1068–1076 (1995). [CrossRef]

**20. **P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A **13**, 779–784 (1996). [CrossRef]

**21. **L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A **13**, 1870–1876 (1996). [CrossRef]

**22. **G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization,” J. Opt. Soc. Am. A **13**, 1019–1023 (1996). [CrossRef]

**23. **E. Popov and M. Nevière, “Grating theory: new equations in Fourier space leading to fast converging results for TM polarization,” J. Opt. Soc. Am. A **17**, 1773–1784 (2000). [CrossRef]

**24. **S. A. Darmanyan, M. Nevière, and A. V. Zayats, “Analytical theory of optical transmission through periodically structured metal films via tunnel-coupled surface polariton modes,” Phys. Rev. B **70**, 075103 (2004). [CrossRef]

**25. **A. Benabbas, V. Halté, and J.-Y. Bigot, “Analytical model of the optical response of periodically structured metallic films,” Opt. Express **13**, 8730–8745 (2005). [CrossRef]

**26. **E. Popov and M. Nevière, “Analytical model of the optical response of periodically structured metallic films: Comment,” Opt. Express **14**, 6583–6587 (2006). [CrossRef]

**27. **W. Cai, D. A. Genov, and V.M. Shalaev, “Superlens based metal-dielectric composites,” Phys. Rev. B **72**, 193101 (2005). [CrossRef]

**28. **M. I. Markovic and A. D. Rakic, “Determination of the reflection coefficients of laser light of wavelengths λ (0.22 *µ*m, 200 *µ*m) from the surface of aluminium using the Lorentz-Drude model,” Appl. Opt. **29**, 3479–3483 (1990). [CrossRef]

**29. **P. G. Etchegoin, E. C. L. Ru, and M. Meyer, “An analytic model for the optical properties of gold,” J. Chem. Phys. **125**, 164705 (2006). [CrossRef]

**30. **E. Palik and G. Ghosh, eds., *The Electronic Handbook of Optical Constants of Solids* (Academic, New York, 1999).

**31. **F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. **43**, 235–286 (2001). [CrossRef]

**34. **Q. Cao and P. Lalanne, “Negative role of surface plasmons in the transmission of metallic gratings with very narrow slits,” Phys. Rev. Lett. **88**, 057403 (2002). [CrossRef]

**35. **M. Nevière and E. Popov, *Light Propagation in Periodic Media: Differential Theory and Design* (Marcel Dekker, 2003).

**36. **N. M. Lyndin, O. Parriaux, and A. V. Tishchenko, “Modal analysis and suppression of the Fourier modal method instabilities in highly conductive gratings,” J. Opt. Soc. Am. A **24**, 3781–3788 (2007). [CrossRef]