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² [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

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

 

Fig. 1 Geometry of the hexagonal PhC (a) and PhC-WG made by omitting one row (b).

Download Full Size | PDF

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

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 û, v̂ mutually orthogonal, they can be related

In the new polarization bases Eq. (13) becomes

3. Distribution of the polarization basis

Before solving Eq. (27), the important point is determining the polarization transform F = [û, v̂] or (equivalently) the functions ξ(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.

 

Fig. 2 Space distribution of the permittivity function (a) and of the vector û of the polarization basis, its azimuth (b) and ellipticity (c) in degrees, plotted within the PhC-WG supercell. The white circles denote the permittivity discontinuities (the hole edges). Note that û is always linear at and perpendicular to the hole edges. The azimuth distribution corresponding the normal vector method (d) is similar to the azimuth in (b).

Download Full Size | PDF

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̂π/Λ_x (the reason for choosing this eigenvalue will be explained below) and we plot the results in Fig. 3.

 

Fig. 3 Convergence properties of the standard factorization (Ho) method (using Eq. (10), the dashed curve with squares), the normal vector method (using Eq. (27) and real û), the dash-dotted curve with triangles), and the complex Fourier factorization method (also using Eq. (27), but with complex û, the solid curve with circles). The 9th eigenfrequency at the X-point (the cutoff frequency) is displayed according to the number of the retained Fourier harmonics M (in our PhC-WG M corresponds to N = 7M).

Download Full Size | PDF

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_y 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 N = 7M. 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.

 

Fig. 4 (a) Dispersion diagram of the PhC-WG. Normalized eigenfrequencies are plotted according to the normalized k_x (in the units of 2π/Λ_x). (b) Detail of the dispersion diagram showing the PhC-WG’s forbidden band and the position of the cutoff frequency.

Download Full Size | PDF

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.

 

Fig. 5 Calculated cutoff wavelengths according to the RI of the medium inside holes in the PhC-WG.

Download Full Size | PDF

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/dn_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).

 

Fig. 6 Transmission spectrum of a finite PhC-WG with vacuum inside holes. The inset displays the PhC-WG geometry.

Download Full Size | PDF

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.

 

Fig. 7 Details of the transmission spectra for various RIs of media inside holes.

Download Full Size | PDF

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.

Appendix

A. Matrix notation

Suppose that ε is periodic and f, d, g_x, and g_y are pseudoperiodic functions. Then the fundamental relations

(28)$$d(x,y)=\epsilon (x,y)f(x,y),$$

(29)$$g_x(x,y)={\partial }_{x}f(x,y),$$

(30)$$g_y(x,y)={\partial }_{y}f(x,y)$$

(typically, f is a component of the electric field, d of the electric displacement, and g_x and g_y of the magnetic field) get the form

(31)$$d_{mn}=\sum _{k,l=-\infty }^{+\infty }{\epsilon }_{m-k,n-l}{f}_{kl},$$

(32)$$g_{x,mn}=-i{p}_m{f}_{mn},$$

(33)$$g_{y,mn}=-i{q}_n{f}_{mn}.$$

assuming the expansions d(x, y) = ∑_m,n d_mn e^{−i(p_mx+q_ny)}, g_x(x, y) = ∑_m,n g_x,mn e^{−i(p_mx+q_ny)}, and g_y(x, y) = ∑_m,n g_y,mn e^{−i(p_mx+q_ny)}.

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 = (2M + 1)(2N + 1), related

(34)$$\alpha (m,n)=m+M+1+(n+N)(2M+1),$$

(35)$$n(\alpha )=(\alpha -1)\mathbf{div}(2M+1)-N,$$

(36)$$m(\alpha )=(\alpha -1)\mathbf{mod}(2M+1)-M,$$

where “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

(37)$$[d]=[[\epsilon ]][f],$$

(38)$$[g_x]=-i\mathbf{p}[f],$$

(39)$$[g_y]=-i\mathbf{q}[f],$$

where [f], [d], [g_x], and [g_y] are column vectors whose αth elements are the Fourier [m, n] elements of the functions f, d, g_x, and g_y, indexed by α(m, n) defined in Eq. (34), and where [[ε]], p, and q are matrices whose elements are defined

(40)$${[[\epsilon ]]}_{\alpha \beta }={\epsilon }_{m(\alpha )-m(\beta ),n(\alpha )-n(\beta )},$$

(41)$$p_{\alpha \beta }=p_{m(\alpha )}{\delta }_{\alpha \beta },$$

(42)$$q_{\alpha \beta }=q_{n(\alpha )}{\delta }_{\alpha \beta },$$

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 −ip and −iq.

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

(43)$$d=\epsilon f,$$

and let [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

(44)$$[d]=[[\epsilon ]][f]$$

can be applied with fast convergence.

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,

(45)$$f=\frac{1}{\epsilon }d,$$

and hence

(46)$$[f]=\left[\left[\frac{1}{\epsilon }\right]\right][d].$$

Accordingly, we can state

(47)$$[d]={\left[\left[\frac{1}{\epsilon }\right]\right]}^{-1}[f],$$

referred to as the inverse rule. We say that the functions ε 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′_x 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 e′ according to the point 5, while Figs. 8(f) and 8(g) correspond to the real function û according to the point 6 of the procedure. These two subfigures have their own color scale displayed on the right.

 

Fig. 8 Space distribution of the functions g_x (a), g_y (b), g′_x (c), e′_x (d), e′_y (e), û_x (f), and û_y (g). The distributions (a)–(c) are same for the CFF method and normal vector method. The distributions (d)–(g) correspond only to the normal vector method.

Download Full Size | PDF

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″_x, while Figs. 9(c) and 9(d) correspond to the real and imaginary part of g″_y, respectively.

 

Fig. 9 Space distribution of the functions Re(g″_x) (a), Im(g″_x) (b), Re(g″_y) (c), and Im(g″_y) (d), all corresponding the CFF method.

Download Full Size | PDF

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 û_x, while Figs. 10(c) and 10(d) correspond to the real and imaginary part of û_y, respectively.

 

Fig. 10 Space distribution of the functions Re(û_x) (a), Im(û_x) (b), Re(û_y) (c), and Im(û_y) (d), all corresponding the CFF method.

Download Full Size | PDF

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]