## Abstract

A photonic crystal waveguide (PhC-WG) was reported to be usable as an optical sensor highly sensitive to various material parameters, which can be detected via changes in transmission through the PhC-WG caused by small changes of the refractive index of the medium filling its holes. To monitor these changes accurately, a precise optical model is required, for which the plane wave expansion (PWE) method is convenient. We here demonstrate the revision of the PWE method by employing the complex Fourier factorization approach, which enables the calculation of dispersion diagrams with fast convergence, i.e., with high precision in relatively short time. The PhC-WG is proposed as a line defect in a hexagonal array of cylindrical holes periodically arranged in bulk silicon, filled with a variable medium. The method of monitoring the refractive index changes is based on observing cutoff wavelengths in the PhC-WG dispersion diagrams. The PWE results are also compared with finite-difference time-domain calculations of transmittance carried out on a PhC-WG with finite dimensions.

© 2014 Optical Society of America

## 1. Introduction

#### 1.1. Photonic crystal sensing

Fast medical-diagnostic techniques have recently become of crucial importance, especially those consisting of miniaturized biological analyzing systems with the ability to detect single molecules without any assistance of special labels or tags (either radioactive or fluorescent). Potential real-time detection in small volumes with a minimum procedure of sample preparation is apt to considerably improve the medical care, biochemical studies, or drug development. A typical approach meeting these criteria is optical biosensing.

Various integrated optical biosensors have obtained considerable attention during the last two decades. A wide variety of optical devices and phenomena have been utilized to detect biological materials without necessity of labels. Among those, measurements of refractive index (RI) change (induced by biomolecular interactions) were found useful for the detection of DNA, proteins, antibody-antigen interactions, cells, and bacteria [1, 2].

Nowadays, commercially available biosensors are mainly based on the surface plasmon resonance (SPR) effect. Although the utilization of SPR in biosensing dates back to 1982 [3], new types of SPR sensors are still occasionally suggested [4, 5]. In principle, they sense RI changes near a metal surface by measuring the changes of reflectance due to the modified coupling of incident light to surface plasmons. Owing to the short decay length of evanescent electromagnetic fields, the probing depth of SPR sensors is limited to approximately 100 nm above the metal surface in aqueous media. Nevertheless, over the years the detection limit of SPR sensors was pushed down to 1 pg/mm^{2} [6]. They are, however, facing some critical issues such as intrinsic losses, low compactness, selectivity, and difficult fabrication and integration with other optical or electrical components. These issues have resulted in SPR commercial systems requiring very sophisticated technologies.

On the other hand, photonic crystal (PhC) biosensors treat these problems in an elegant manner. They provide extremely narrow resonant behavior, the resonant characteristics are readily tuned by adsorption of a biochemical material, they can be manufactured inexpensively, and they can be measured simply by illuminating with white light [7]. Therefore, a utilization of unique properties of mass-produced PhCs can lead to more efficient and accurate biological assays.

Various PhC biosensing devices such as plastic one-dimensional (1D) subwavelength gratings [8], two-dimensional (2D) silicon nitride arrays [9], 2D-patterned dielectric slabs [10], PhC fibers [11], and PhC waveguide (PhC-WG) systems [12] have been reported. They all exploit extremely sensitive resonant conditions for guided modes in the PhC with respect to RI changes in the ambient medium. On the other hand, some critical issues such as high polarization sensitivity or difficult light coupling into PhC structures must still be overcome. A proper PhC design for biosensing is therefore an essential task, which should be carefully handled to obtain the required sensing properties.

#### 1.2. Numerical troubles with plane wave expansion modeling

To propose and analyze PhC-based sensors and their functionality, accurate optical modeling techniques are required. We here use the plane wave expansion (PWE) method (whose effective implementation to a PhC-WG is the main subject of this paper) and the finite-difference time-domain (FDTD) method (for comparison).

The PWE method is particularly convenient for periodic structures such as PhCs since it is based on the expansion of periodic electromagnetic fields and permittivity functions into Fourier series, whose components are then treated by linear-algebraic numerical techniques. The same procedures have been carried out earlier in the frame of coupled-wave theory (or Fourier modal theory) of light diffraction by gratings (which can be referred to as PhC slabs). It has been reported that the numerical effectiveness of the PWE method, represented by the convergence rate in the truncated basis of the Fourier components, strongly depends on the specific implementation of the method, since there are various possibilities of the expansion of the material and field parameters [13, 14]. In particular, periodic discontinuities typical for PhC’s elements cause rather poor convergence and extremely time- and memory-consuming computer simulations, especially when the electric field vector has the nonzero perpendicular component to the permittivity discontinuities.

Fortunately, several authors have reported that a proper reformulation of the PWE equations can provide a high enhancement to the PWE method even for the troublesome polarization [15, 16, 17, 18]. Meanwhile, Li [19] showed that Laurent’s “direct” rule (adopted in conventional formulations to factorize the truncated Fourier series that corresponds to products of two periodic functions) yields bad convergence when the two functions of the product are simultaneously discontinuous. He suggested Fourier factorization (FF) rules (briefly summarized in Section B) and applied them to 1D gratings. This major breakthrough in the grating theory was soon applied by many authors to various grating structures with arbitrary periodic reliefs, anisotropic [20] and slanted [21] periodic systems, their various combinations [22, 23, 24], and other systems [25, 26, 27].

Later Li [28] applied the FF rules to 2D-periodic structures by using the “zigzag” Fourier expansion, which yielded an improvement for rectangular dots or holes. Then, Popov and Neviere [29] applied a coordinate transformation to treat individually the normal and tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the application of the correct rule for each field component and thus improved the convergence on relief gratings with round-shape elements. Later David et al. [30] utilized the normal–tangential field separation to 2D PhCs composed of circular or elliptical holes. Similarly, Schuster et al. [31] applied this method to 2D gratings, and also suggested more general distributions of polarization bases [32].

These approaches, always dealing with linear polarizations, enabled a significant improvement of the convergence properties. However, they ignored the fact that the transformation matrix between the Cartesian and the normal–tangential component bases of polarization became discontinuous at some points of the periodic cell, which slowed down the convergence. To overcome the discontinuities, a distribution of more complex (i.e., generally elliptic) polarization bases was suggested to improve optical simulations of 2D gratings and photonic crystals [33, 34], here referred to as complex Fourier factorization (CFF).

In this paper we employ the CFF approach into the PWE method to simulate a PhC-WG-based sensor and present fast convergence of our implementation compared to the standard approach. To increase the numerical performance, we determine the polarization basis distribution by calculating electrostatic field generated by hypothetical electric dipoles located at permittivity discontinuities (perpendicular to the discontinuities). After proper modification of the electrostatic field to generally complex-valued vectors, the polarization basis distribution becomes a continuous and smooth function. We define the PhC-WG as a line defect in a hexagonal array of cylindrical holes periodically arranged in a bulk silicon. The holes are then filled with a variable medium whose RI changes can be monitored via observing cutoff wavelengths in the PhC-WG dispersion diagrams. Finally we compare the PWE results with FDTD calculations of transmittance on a PhC-WG with finite dimensions.

We hope that this article will help sensor engineers in simple design and modeling optical sensors based on photonic crystals, or at least in understanding the principles of the PWE method employing CFF.

## 2. Plane wave expansion method

Eigenmodes and the optical response of PhCs or PhC-WGs can be calculated by means of the wave equation

where*ω*is the frequency,

*c*is the light velocity in vacuum,

**is the vector of electric field and**

*E**ε*(

*x*,

*y*) is the relative permittivity, being a 2D-periodic function of two variables throughout this paper. In the PWE method we expand each component

*j*∈ {

*x*,

*y*,

*z*} of the electric field vector into the pseudo-Fourier series, i.e.,

*p*= 2

*π*/Λ

*, and*

_{x}*q*= 2

*π*/Λ

*are the normalized reciprocal lattice vectors, with the periods Λ*

_{y}*and Λ*

_{x}*along the Cartesian axes. For brevity we have also defined*

_{y}**= [**

*k**k*;

_{x}*k*],

_{y}*p*=

_{m}*k*+

_{x}*mp*, and

*q*=

_{n}*k*+

_{y}*nq*. Analogously we could write the pseudo-Fourier expansion for any other field from Maxwell’s equations. Note that we assume the uniformity in the

*z*-axis direction for all the fields.

Similarly, we write the usual Fourier series for the relative permittivity or the impermittivity (*η* = 1/*ε*)

We here assume a hexagonal PhC made as 2D-periodic holes in a medium with a high RI. The geometry of the problem, the periodic cell, and the first Brillouin zone in the reciprocal space are depicted in Fig. 1(a). Using the same principle, we can simulate a PhC-WG made by omitting one row of holes as displayed in Fig. 1(b), where instead of a small periodic cell we have a large periodic supercell (with the size chosen so that two parallel waveguides do not influence each other). Assuming propagation within the *xy* plane and uniformity along the *z* axis (*∂ _{z}* = 0), Eq. (1) becomes

*E*-polarization and

**denotes a scaled electric displacement. In the case of E-polarization,**

*D̃**D̃*=

_{z}*εE*, we simply use the Laurent rule because the function

_{z}*E*is continuous everywhere, including the PhC discontinuities (the discontinuities of the function

_{z}*ε*).

The case of H-polarization (which is particularly important for a PhC-WG-based sensor),

**and**

*E***vectors by appropriate changes of the polarization bases at all points of the periodic cell. For this purpose we define a space dependent matrix**

*D̃***F**(

*x*,

*y*) of the polarization transformation, so that where

**,**

*x̂***are unit vectors pointing along the Cartesian coordinates, while**

*ŷ***,**

*û***denote space-dependent unit polarization vectors,**

*v̂***being always normal to all PhC discontinuities and**

*û***being tangential. In other words,**

*v̂**E*,

_{u}*E*are the components of the electric field in the new polarization bases at each point of the periodic cell, and analogously for the vector

_{v}**.**

*D̃*Taking advantage of generally elliptical polarizations, the elements *F _{ij}* of the matrix

**F**(now being complex numbers) can be chosen to be continuous functions of the whole space. Since we choose both basis vectors

**,**

*û***mutually orthogonal, they can be related**

*v̂**ξ*(

*x*,

*y*) and

*ζ*(

*x*,

*y*) are complex-valued continuous functions of the 2D space.

In the new polarization bases Eq. (13) becomes

*E*(tangential component of

_{v}**) and**

*E**D̃*(normal component of

_{u}**) are always continuous. Expressing this equation in the Cartesian basis yields**

*D̃*## 3. Distribution of the polarization basis

Before solving Eq. (27), the important point is determining the polarization transform **F** = [** û**,

**] or (equivalently) the functions**

*v̂**ξ*(

*x*,

*y*) and

*ζ*(

*x*,

*y*) so that

**is always linear at and perpendicular to all points of the permittivity discontinuities.**

*û*In this paper we assume a PhC-WG made as a line defect in a hexagonal PhC. The 2D-periodic supercell of the PhC-WG with the values of the permittivity is displayed in Fig. 2(a), where the red color stands for *ε* = 12 (approximately the permittivity of Si around the wavelength *λ* = 1700 nm), the blue color for *ε* = 1 (circular holes filled with vacuum), and the white circles denote the permittivity discontinuities. We set the period of the adjacent holes Λ* _{x}* = 370 nm and the diameter of the holes

*D*= 240 nm. Using this geometry (this 2D-periodic function

*ε*), we determine

**by the following procedure.**

*û*- We define a rectangular grid of
*N*-by-_{x}*N*vertices with coordinates [_{y}*x*,_{i}*y*] with_{j}*i*= 1, 2,...,*N*(do not confuse here with an imaginary unit),_{x}*j*= 1, 2,...,*N*,_{y}*x*=_{i}*i*Λ, and_{x}/N_{x}*y*=_{j}*j*Λ, so that we obtain a discrete permittivity_{y}/N_{y}*ε*. Below we calculate all space-dependent functions at these vertices. Since in our case ${\mathrm{\Lambda}}_{y}=4\sqrt{3}{\mathrm{\Lambda}}_{x}\approx 7{\mathrm{\Lambda}}_{x}$, we use_{i,j}*N*= 256 and_{x}*N*= 7 · 256 = 1792._{y} - We define
= [*g**∂*] (the gradient of the permittivity, calculated via the differences between adjacent vertices). We define the numerical gradient_{x}ε; ∂_{y}ε*g*= (_{x,ij}*ε*_{i}_{−1,}+_{j}*ε*_{i}_{−1,}_{j}_{−1}−*ε*−_{i,j}*ε*_{i,j}_{−1})/2 and*g*= (_{y,ij}*ε*_{i,j}_{−1}+*ε*_{i}_{−1,}_{j}_{−1}−*ε*−_{i,j}*ε*_{i}_{−1,})/2 assuming the unit size of the grid cell. The distribution of the gradient and other functions is displayed in Appendix C._{j} - To obtain a smoother gradient, we apply space averaging for its nonzero elements, i.e.,
*g′*= (_{ij}*g*+ ∑_{ij}*g*_{i}_{±1,}_{j}_{±1})/*N*, where_{ij}*N*is the number of the nonzero elements surrounding the [_{ij}*i*,*j*]-vertex (including the [*i*,*j*]-element). In our example we repeat the averaging 100 times. - We calculate the azimuth of the
vector,*g′**ϕ*_{0}= arg(*g′*+_{x}*ig′*), which we use to define a complex vector_{y}=*g″*exp(*g′**iϕ*_{0}). - Then we apply weighted space averaging for all elements,
= avg(*e*), precisely*g″**e*= [_{ij}*w*_{0}*g″*+_{i,j}*w*_{1}(*g″*_{i}_{+1,}+_{j}*g″*_{i}_{−1,}+_{j}*g″*_{i,j}_{+1}+*g″*_{i,j}_{−1}) +*w*_{11}(*g″*_{i}_{+1,}_{j}_{+1}+*g″*_{i}_{+1,}_{j}_{−1}+*g″*_{i}_{−1,}_{j}_{+1}+*g″*_{i}_{−1,}_{j}_{−1})]/(*w*_{0}+ 4*w*_{1}+ 4*w*_{11}). In our example we repeat this averaging,= avg(*e′*), 2000 times, but between each two iterations we keep the same value where*e**g*were nonzero. For effective averaging we set the weights_{i,j}*w*_{0}= 0.1,*w*_{1}= 1/6, and*w*_{11}= 1/12. Note that the ratio*w*_{1}/*w*_{11}= 2 is important for the numerical efficiency. - Finally we define the unit vector
= [*û**ξ; ζ*] =/|*e′*|, whose space distribution is displayed in Figs. 2(b) (its azimuth) and 2(c) (its angular ellipticity), where we can see that it is completely continuous and smooth.*e′*

Because the vector ** û** is complex (due to the point 4 above), we call this procedure the CFF method. If we omit the 4th point from the above list, we obtain only real-valued vectors including the vector

**. In this case we cal the procedure the normal vector method, because the vector**

*û***is normal to all the discontinuities of the permittivity. However, it has discontinuities at some points of the supercell, as visible from its azimuth distribution, displayed in Fig. 2(d). Note thas this method was introduced by David et al. [30], but with a more simple definition of the**

*û***distribution.**

*û*## 4. Results and discussion

To choose appropriate precision for calculation, we perform a brief numerical analysis of all the factorization methods, the standard method (named after Ho, who simply Fourier-expanded the permittivity function) represented by Eq. (10), the normal vector method represented by Eq. (27) (and using the real vector ** û**) and the CFF method represented also by Eq. (27) (and with complex

**) (where the polarization basis distributions**

*û***were described in the previous section). We calculate the 9th eigenfrequency of the PhC-WG described above at the X-point of the first Brillouin zone, which corresponds to the wave vector**

*û***=**

*k*

*x̂**k*=

_{x}

*x̂**π*/Λ

*(the reason for choosing this eigenvalue will be explained below) and we plot the results in Fig. 3.*

_{x}In PWE calculations we want the *M*-by-*N* ratio (the maximum numbers of the Fourier harmonics defined in Appendix A) same as the *N _{x}*-by-

*N*ratio of the supercell sizes in the direct space, in order to keep the same spacial resolution in both Cartesian directions. For this reason we choose

_{y}*N*= 7

*M*. According to the curves in Fig. 3, all the three methods are usable at the point

*M*= 4 for showing visual dispersion diagrams for which three effective digits are sufficient. However, the CFF method becomes considerably better for

*M*> 6, where it yields the precision with four effective digits. For this reason we here present dispersion relations with the precision corresponding to

*M*= 4 but the accurate frequencies or wavelengths with the precision corresponding to

*M*= 8.

The dispersion diagram of the PhC-WG as defined by Fig. 2(a) is displayed in Fig. 4(a). For its calculation, Eq. (27) was used, with the polarization basis distribution described in the previous section, with the numerical precision corresponding to *M* = 4 and *N* = 28. The dispersion diagram exhibits a band of forbidden frequencies between the 8th and the 9th bands. For better visualization, Fig. 4(b) shows a detail of the dispersion relation in the vicinity of the forbidden band (highlighted by two horizontal lines). The upper frequency of the forbidden band is denoted as the cutoff frequency, and its value was the one plotted in Fig. 3.

The purpose of the PhC-WG studied in this paper is sensing small changes of the RI of the medium inside holes via measuring the PhC-WG’s wavelength-dependent transmittance. Since the transmittance is approximately proportional to the density of propagation states of the structure, it is obviously zero in the range of the forbidden frequency band. The corresponding forbidden band of wavelengths starts at the cuttoff wavelength, which is *λ*_{cutoff} ≈ 1694 nm in the case of vacuum inside holes. The dependence of the cutoff wavelengths on the RI of the medium in holes [in other words, the function *λ*_{cutoff}(*n*_{in})] is displayed in Fig. 5, where the calculations were carried out with the precision *M* = 8 and *N* = 56.

Figure 5 shows that the function *λ*_{cutoff}(*n*_{in}) is nearly linear, so that the small change of the RI yields a proportional change of *λ*_{cutoff} with the proportionality factor given by the derivative d*λ*_{cutoff}/d*n*_{in} ≈ 78 nm.

In practice we can measure spectral transmittance of a realistic PhC-WG structure, which is a PhC-WG without periodic boundaries in the *y*-axis direction and made of only a finite number of periods in the propagation (*x*-axis) direction. For comparison, Fig. 6 displays the transmittance of the PhC-WG made as such a more realistic structure (assuming vacuum in holes), with 17 periods in the *x*-axis direction and 8 rows of holes on both sides of the cladding (as shown in the inset).

The transmittance spectrum was calculated by means of the FDTD software package MEEP [35]. The calculation was carried out in the 2D mode with perfectly matched layers at the edges of the calculation grid. The resolution of the calculation grid was 19 nm and the subpixel averaging was used. Gaussian pulse with the central frequency at 1423 nm and the pulse width of 925 nm located at the input of the PhC slab was employed as a radiation source for the calculation.

Figure 6 shows that the transmittance drops down at about 1702 nm, which is approximately in accord with the *λ*_{cutoff} result obtained by the PWE method. Similarly, Fig. 7 displays the spectral dependences of transmittance near the cutoff wavelengths for several values of RI (*n*_{in} = 1, 1.1,...,1.5) of the same finite PhC-WG structure, which again demonstrate good accordance with the cutoff wavelengths calculated by the PWE method.

The differences in the cutoff wavelengths yielded by both calculation methods are mainly due to the finite size of the PhC-WG supercell in the *y*-axis direction (or, equivalently, the nonzero size of the first Brillouin zone in the same direction). This finite size causes that part of light propagating along one waveguide (a missing row of holes) passes into another parallel waveguide, so that some propagating modes might also appear at the M-point of the first Brillouin zone, which corresponds to the wave vector ** k** = [

*k*,

_{x}*k*] = [

_{y}*π*/Λ

*,*

_{x}*π*/Λ

*].*

_{y}## 5. Conclusions

We have demonstrated an effective implementation of the PWE method by employing the CFF approach for modeling PhC-WG-based optical sensors. For the CFF approach we have developed a general procedure of determining a continuous and smooth polarization basis distribution for PWE calculations. Then we have shown that the design of the PhC-WG can be simply proposed by analyzing the dispersion relations of the PhC-WG, because the cutoff wavelengths of the transmittance spectra can be directly determined from the dispersion diagrams.

## A. Matrix notation

Suppose that *ε* is periodic and *f*, *d*, *g _{x}*, and

*g*are pseudoperiodic functions. Then the fundamental relations

_{y}*f*is a component of the electric field,

*d*of the electric displacement, and

*g*and

_{x}*g*of the magnetic field) get the form assuming the expansions

_{y}*d*(

*x*,

*y*) = ∑

_{m,n}*d*

_{mn}*e*

^{−i(pmx+qny)},

*g*(

_{x}*x*,

*y*) = ∑

_{m,n}*g*

_{x,mn}*e*

^{−i(pmx+qny)}, and

*g*(

_{y}*x*,

*y*) = ∑

_{m,n}*g*

_{y,mn}*e*

^{−i(pmx+qny)}.

Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the summation
${\sum}_{m=-M}^{+M}{\sum}_{n=-N}^{+N}$, we can renumber all the indices to replace the couple of two sets *m* ∈ {−*M*, −*M* + 1, ..., *M*} and *n* ∈ {−*N*, −*N* + 1, ..., *N*} by a single set of indices *α* ∈ {1, 2, ..., *α*_{max}}, with *α*_{max} = (2*M* + 1)(2*N* + 1), related

**div**” denotes the operation of integer division and “

**mod**” the remainder (the modulo operation). Then we can rewrite Eqs. (31)–(33) into the matrix relations where [

*f*], [

*d*], [

*g*], and [

_{x}*g*] are column vectors whose

_{y}*α*th elements are the Fourier [

*m*,

*n*] elements of the functions

*f*,

*d*,

*g*, and

_{x}*g*, indexed by

_{y}*α*(

*m*,

*n*) defined in Eq. (34), and where [[

*ε*]],

**p**, and

**q**are matrices whose elements are defined where the indices on the right hand parts are defined by Eqs. (35) and (36) and where

*δ*denotes the Kronecker delta. As a summary we can say that the multiplication by a function is in the reciprocal space represented by the matrix [[

_{αβ}*ε*]] (in the sense of the limit

*α*

_{max}→ ∞) and that the partial derivatives are represented by the diagonal matrices −

*i*

**p**and −

*i*

**q**.

## B. Rules of Fourier factorization

Although the rules were derived by Li [19] for 1D periodic functions, we here write them in the matrix formalism independent of the number of dimensions. Let *f*, *d*, and *ε* be piecewise-continuous functions with the same periodicity related

*f*], [

*d*], and [[

*ε*]] denote their matrices as defined in Appendix A.

*Laurent rule:* If *ε* and *f* have no concurrent discontinuities, then the Laurent rule applied to Eq. (43) converges uniformly on the whole period and hence

*Inverse rule:* If *ε* and *f* have one or more concurrent discontinuities but *d* is continuous, then Eq. (43) can be transformed into the case of the Laurent rule,

*ε*and

*f*have complementary discontinuities.

If none of the requirements for the two rules are satisfied, then none of them can be applied correctly because Eqs. (44) and (47) are no longer valid at the points of discontinuities, which considerably slows down the convergence. Therefore, we should carefully analyze the continuity of the functions and transform all the partial differential formulae to the two cases.

## C. Numerical details for the polarization basis distribution

Figure 8 displays the space distributions of several functions calculated during the procedure described in Section 3. Figures 8(a) and 8(b) display the distribution of the both Cartesian components of the vector ** g** determined by the point 2 of the procedure in Section 3. Figure 8(c) displays the smoothed function

*g′*according to the point 3 of the procedure. These three distributions correspond to both the CFF method and the normal vector method, whereas the remaining subfigures of Fig. 8 correspond only to the normal vector method, for which the point 4 of the procedure was omitted. Figures 8(d) and 8(e) display both Cartesian components of the real function

_{x}**according to the point 5, while Figs. 8(f) and 8(g) correspond to the real function**

*e′***according to the point 6 of the procedure. These two subfigures have their own color scale displayed on the right.**

*û*Figure 9 displays the space distribution of the function ** g″** determined according to the point 4 of the procedure, i.e., after making the gradient of the permittivity complex-valued. Figures 9(a) and 9(b) correspond to the real and imaginary part of

*g″*, while Figs. 9(c) and 9(d) correspond to the real and imaginary part of

_{x}*g″*, respectively.

_{y}Finally, Fig. 10 displays the distribution of the desired vector ** û** determined according to the last point 6 of the procedure. Analogously, Figures 10(a) and 10(b) correspond to the real and imaginary part of

*û*, while Figs. 10(c) and 10(d) correspond to the real and imaginary part of

_{x}*û*, respectively.

_{y}## Acknowledgments

The work was supported by the Grant Agency of the Czech Republic, projects P205/11/2137 and 13-30397S.

## References and links

**1. **A. Sassolas, B. D. Leca-Bouvier, and L. J. Blum, “Dna biosensors and microarrays,” Chem. Rev. **108**, 109–139 (2008). [CrossRef]

**2. **A. D. Taylor, J. Ladd, J. Homola, and S. Jiang, *Surface Plasmon Resonance (SPR) Sensors for the Detection of Bacterial Pathogens* (Springer, 2008), pp. 83–108.

**3. **C. Nylander, B. Liedberg, and T. Lind, “Gas-detection by means of surface-plasmon resonance,” Sens. Actuators **3**, 79–88 (1982). [CrossRef]

**4. **M. Piliarik, H. Sipova, P. Kvasnicka, N. Galler, J. R. Krenn, and J. Homola, “High-resolution biosensor based on localized surface plasmons,” Opt. Express **20**, 672–680 (2012). [CrossRef] [PubMed]

**5. **D. Regatos, D. Farina, A. Calle, A. Cebollada, B. Sepulveda, G. Armelles, and L. M. Lechuga, “Au/fe/au multilayer transducers for magneto-optic surface plasmon resonance sensing,” J. Appl. Phys. **108**, 054502 (2010). [CrossRef]

**6. **J. Homola, S. S. Yee, and G. Gauglitz, “Suraface plamson resonance sensors: review,” Sens. Actuators B **54**, 3–15 (1999). [CrossRef]

**7. **D. Threm, Y. Nazirizadeh, and M. Gerken, “Photonic crystal biosensors towards on-chip integration,” J. Biophotonics **5**, 601–616 (2012). [CrossRef]

**8. **B. Cunningham, B. Lin, J. Qiu, P. Li, J. Pepper, and B. Hugh, “A plastic colorimetric resonant optical biosensor for multiparallel detection of label-free biochemical interactions,” Sens. Actuators B Chem. **85**, 219–226 (2002). [CrossRef]

**9. **B. Cunningham, P. Li, B. Lin, and J. Pepper, “Colorimetric resonant reflection as a direct biochemical assay technique,” Sens. Actuators B Chem. **81**, 316–328 (2002). [CrossRef]

**10. **M. El Beheiry, V. Liu, S. Fan, and O. Levi, “Sensitivity enhancement in photonic crystal slab biosensors,” Opt. Express **18**, 22702–22714 (2010). [CrossRef] [PubMed]

**11. **J. Jensen, P. Hoiby, G. Emiliyanov, O. Bang, L. Pedersen, and A. Bjarklev, “Selective detection of antibodies in microstructured polymer optical fibers,” Opt. Express **13**, 5883–5889 (2005). [CrossRef] [PubMed]

**12. **N. Skivesen, A. Tetu, M. Kristensen, J. Kjems, L. H. Frandsen, and P. I. Borel, “Photonic-crystal waveguide biosensor,” Opt. Express **15**, 3169–3176 (2007). [CrossRef] [PubMed]

**13. **H. S. Sozuer, J. W. Haus, and R. Inguva, “Photonic bands: Convergence problems with the plane-wave method,” Phys. Rev. B **45**, 13962–13972 (1992). [CrossRef]

**14. **R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, “Accurate theoretical analysis of photonic band-gap materials,” Phys. Rev. B **48**, 8434–8437 (1993). [CrossRef]

**15. **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]

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

**17. **P. Lalanne, “Convergence performance of the coupled-wave and the differential methods for thin gratings,” J. Opt. Soc. Am. A **14**, 1583–1591 (1997). [CrossRef]

**18. **P Lalanne, “Improved formulation of the coupled-wave method for two-dimensional gratings,” J. Opt. Soc. Am. A **14**, 1592–1598 (1997). [CrossRef]

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

**20. **L. Li, “Reformulation of the Fourier modal method for surface-relief gratings made with anisotropic materials,” J. Mod. Opt. **45**, 1313–1334 (1998). [CrossRef]

**21. **B. Chernov, M. Neviere, and E. Popov, “Fast Fourier factorization method applied to modal analysis of slanted lamellar diffraction gratings in conical mountings,” Opt. Commun. **194**, 289–297 (2001). [CrossRef]

**22. **K. Watanabe, R. Petit, and M. Neviere, “Differential theory of gratings made of anisotropic materials,” J. Opt. Soc. Am. A **19**, 325–334 (2002). [CrossRef]

**23. **K. Watanabe, “Numerical integration schemes used on the differential theory for anisotropic gratings,” J. Opt. Soc. Am. A **19**, 2245–2252 (2002). [CrossRef]

**24. **L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A **5**, 345–355 (2003). [CrossRef]

**25. **P. Boyer, E. Popov, M. Neviere, and G. Tayeb, “Diffraction theory in TM polarization: application of the fast Fourier factorization method to cylindrical devices with arbitrary cross section,” J. Opt. Soc. Am. A **21**, 2146–2153 (2004). [CrossRef]

**26. **N. Bonod, E. Popov, and M. Neviere, “Light transmission through a subwavelength microstructured aperture: electromagnetic theory and applications,” Opt. Commun. **245**, 355–361 (2005). [CrossRef]

**27. **N. Bonod, E. Popov, and M. Neviere, “Fourier factorization of nonlinear Maxwell equations in periodic media: application to the optical Kerr effect,” Opt. Commun. **244**, 389–398 (2005). [CrossRef]

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

**29. **E. Popov and M. Neviere, “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]

**30. **A. David, H. Benisty, and C. Weisbuch, “Fast factorization rule and plane-wave expansion method for two-dimensional photonic crystals with arbitrary hole-shape,” Phys. Rev. B **73**, 075107 (2006). [CrossRef]

**31. **T. Schuster, J. Ruoff, N. Kerwien, S. Rafler, and W. Osten, “Normal vector method for convergence improvement using the RCWA for crossed gratings,” J. Opt. Soc. Am. A **24**, 2880–2890 (2007). [CrossRef]

**32. **P. Gotz, T. Schuster, K. Frenner, S. Rafler, and W. Osten, “Normal vector method for the RCWA with automated vector field generation,” Opt. Express **16**, 17295–17301 (2008). [CrossRef] [PubMed]

**33. **R. Antos, “Fourier factorization with complex polarization bases in modeling optics of discontinuous bi-periodic structures,” Opt. Express **17**, 7269–7274 (2009). [CrossRef] [PubMed]

**34. **R. Antos and M. Veis, “Fourier factorization with complex polarization bases in the plane-wave expansion method applied to two-dimensional photonic crystals,” Opt. Express **18**, 27511–27524 (2010). [CrossRef]

**35. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. **181**, 687–702 (2010). [CrossRef]