## Abstract

In this article we present a novel approach for determining the electromagnetic modes of photonic multilayer structures. We combine the plane wave expansion method with the method of lines resulting in a fast and accurate computational technique which we named the *plane wave admittance method*. In addition, we incorporate perfectly matched layers at the boundaries parallel to the multilayer surfaces which allow for easy determination of leaky modes. The convergence of the method is verified for the case of photonic crystal slab showing very good agreement with the results obtained with full three-dimensional plane wave expansion method while the numerical effort is largely reduced. The numerical implementation of the method will be soon available on the web.

©2005 Optical Society of America

## 1. Introduction

Recently, photonic crystals [1] (PC) have been attracting a lot of interest due to their novel optical properties such as the existence of photonic bandgaps, modified density of states, appearance of localized states of light, anomalous dispersion, suppressed or enhanced nonlinear effects [2] etc. Photonic crystals enabled the development of a range of novel devices, including waveguides and fibers with bandgap guidance, cavities with improved Q-factors, fiber lasers and semiconductor lasers [3] with low threshold and other novel properties, to mention only some. In this paper we introduce a novel modeling method for photonic crystal slabs (PCS) as well as layered PC-devices, which are particularly important for applications in integrated photonics. These are structures with two-dimensional periodicity in the horizontal plane (**x̂, ŷ**), and a uniform or layered composition in the orthogonal (**ẑ**) direction. In all cases, we assume the structures to have a finite thickness, and for simplicity we term them generally as PCS also when speaking of PC-multilayered stacks. The importance of finite thickness, except for technological aspects, relates to the the appearance of leaky modes and to the confinement mechanism which for single layered structures is the total internal reflection.

Numerical modeling of PC structures is difficult mainly due to the importance of vectorial effects that arise when high refractive index contrasts exist at sub-wavelength distances. This leads to the necessity of solving the fully vectorial wave equation and taking special care for the exact representation of the material properties near the boundaries. A large variety of methods have been proposed for PC modeling, both in time [4] and in frequency domains. Among the frequency domain modal methods, the plane-wave method [2, 5, 6] (PWM) in its numerous variants has been particularly widely used. The PWM naturally includes the periodic boundary conditions of the PC structures. It can be easily adapted to any desired form of the wave equation. The plane-wave representation of the permittivity and permeability has an analytical form for PC with most of regular geometries, however is not limited to them and can be easily determined numerically for arbitrary geometries as well [7]. Although the PWM is known to suffer from convergence problems [8], the nature of these problems is well understood and can be overcome by a proper choice of truncation rules [9] and by refractive index averaging [10], or by taking a sufficiently large number of basis functions. The PWM has been used for modal calculations as well as for bandgap determination in PC with full 3D periodicity [6], but also in a 2D variant for PC-fibers [11, 12] and PCS [13, 14]. As opposed to the case of PC-fibers, the proper analysis of PCS of finite thickness is a 3D problem. One way to reduce it into a 2D problem is to apply the effective index technique as in ref. [13, 15]. Another is to make use of the analytical form of Bloch solution of Maxwell equations, propagating within every layer along the direction **ẑ** of translational invariance. Then, the matching condition may follow from the transmission matrix or better from the scattering matrix approach as in refs. [14, 16]. The latter is particularly useful in the determination of the density of states of leaky modes in PCS [17].

Here, we propose to follow another matching condition based on the admittance transfer and developed by Pregla et al. under the name of Method of Lines [18, 19, 20] (MoL) or simply the admittance transfer method. This method is characterized by a very good numerical stability and has been successfully applied to layered structures with use of discretization by finite differences in the (**x̂**, **ŷ**) plane.

In this paper we propose to combine the admittance transfer method and PWM to achieve an efficient way of calculating the modes and bandgap structure of PCS. For the boundary conditions, we use the Bloch expansion in (**x̂, ŷ**) together with absorbing perfectly matched layer [21, 22] (PML) in **ẑ**, to account for the horizontal periodicity of the PCS and out-of-plane losses, respectively. The equations assume a diagonal and independent permittivity and permeability, and such media can be directly included in modeling. Clearly, adding a PML in the horizontal plane is also straightforward within the same material model. In the formulation of the PWM eigenequations for every layer, we take care not to put the frequency as the eigenvalue, therefore the treatment of dispersive and active materials is also direct.

The paper is organized as follows: In the second section we describe in details the proposed algorithm; we derive the necessary equations and suggest a way for their implementation. In the third section we apply our method to exemplary PCS structure and obtain the eigenmodes and their electromagnetic field distribution. Finally, in section 4 we demonstrate a full photonic band diagram computation and discuss the convergence and the accuracy of the results.

## 2. Computational method

#### 2.1. Generalized transmission line equations

In the method of lines (MoL) the analyzed structure is divided into several layers in which material parameters are constant along the *z* direction making it possible to obtain analytical solution [19]. In our work, the generalized transmission line (GTL) equations are used for this purpose. However, contrary to MoL, they are derived in plane-wave basis without using finite-differences.

The time independent Maxwell equations with diagonal anisotropic dielectric and magnetic tensors in Cartesian coordinates read

where *∂*_{ξ}
means partial derivative in the *ξ* direction. From the Eq. (3) and (6) it is possible to express the *z* components of the electric and magnetic vectors as

where *k*
_{0}=*ω/c*=*ω*(*µ*_{0}*ε*_{0}
)^{1/2} is the normalized frequency and *η*_{0}
=(*µ*
_{0}/*ε*
_{0})^{1/2} is the free space impedance. Substituting (7) and (8) into (1), (2), (4), and (5) we can obtain the GTL equations of the form

The partial derivatives in Eq. (9) and (10) have to be computed numerically. The traditional approach here is to use finite differences to discretize the right-hand sides of Eq. (9) and (10) leading to matrix equations

However, in order to ensure the satisfactory computational accuracy the number of the finite difference grid points must be large, approaching even the number of 10000 in the case of fully three-dimensional analysis [20]. This makes the computational procedure a very strenuous task, almost impossible to complete in a reasonable time, especially on a personal computer. For this reason we propose an alternative numerical representation of Eq. (9) and (10). As it will be shown later the size of the obtained matrices can be as small as 500 or even less still resulting in a satisfactory precision.

#### 2.2. Plane wave expansion of GTL equations

In general, GTL equations can be transformed into a finite problem by expanding the fields in some truncated basis {|*φG*〉}

where Ψ can be any of *E*_{x}*, -E*_{y}*, H*_{x}
or *H*_{y}
and {Ψ
^{G}
} is a set of expansion coefficients (To simplify the notation we expand -*E*_{y}
as -*E*_{y}
=∑
_{G}*E*^{G}*y*|*φG*〉.).

In our approach we propose a choice of such basis, namely the plane wave expansion, which is particularly useful in the case of periodic structures. For two-dimensional domain the basis functions take the form

where **r**‖ is the position in the *xy* plane and {**G**} forms a set of reciprocal lattice vectors. The plane wave basis fulfills the orthogonality condition, i.e.

where *δ*_{mn}
is the Kronecker delta and *G*_{m}*,G*_{n}
∈{**G**}.

The Bloch theorem states that in periodic structures each component of the electromagnetic field Ψ(**r**‖) can be expressed as

where **k**=*k*_{x}**x̂**+*k*_{y}
**ŷ** is the horizontal wavevector and Φ(**r**_{‖}) some periodic function. Hence the plane wave expansion of Eq. (13) and (14) becomes

Thanks to the orthogonality of the plane wave basis (Eq. (15)), Eq. (9) and (10) can be represented as

The matrix operators can be determined analytically by expanding the material parameters in plane wave basis as

and substituting this into Eq. (18) and (19) we obtain the final form of the plane-wave expansion of GTL equations

$$-i\sum _{\mathbf{G}\prime}\frac{{\eta}_{0}}{{k}_{0}}\left[\begin{array}{cc}-\left({G}_{y}+{k}_{y}\right)\left({G}_{y}^{\prime}+{k}_{y}\right){\kappa}^{\mathbf{G}-\mathbf{G}\prime}+{k}_{0}^{2}{\mu}_{x}^{\mathbf{G}-\mathbf{G}\text{'}}& \left({G}_{y}+{k}_{y}\right)\left({G}_{y}^{\prime}+{k}_{x}\right){\kappa}^{\mathbf{G}-\mathbf{G}\prime}\\ \left({G}_{x}+{k}_{x}\right)\left({G}_{y}^{\prime}+{k}_{y}\right){\kappa}^{\mathbf{G}-\mathbf{G}\prime}& -\left({G}_{x}+{k}_{x}\right)\left({G}_{x}^{\prime}+{k}_{x}\right){\kappa}^{\mathbf{G}-\mathbf{G}\prime}+{k}_{0}^{2}{\mu}_{y}^{\mathbf{G}-\mathbf{G}\prime}\end{array}\right]\phantom{\rule{.2em}{0ex}}\left[\begin{array}{c}{H}_{x}^{\mathbf{G}\prime}\\ {H}_{y}^{\mathbf{G}\prime}\end{array}\right],$$

$$-i\sum _{\mathbf{G}\prime}\frac{1}{{\eta}_{0}{k}_{0}}\left[\begin{array}{cc}-\left({G}_{x}+{k}_{x}\right)\left({G}_{x}^{\prime}+{k}_{x}\right){\gamma}^{\mathbf{G}-\mathbf{G}\prime}+{k}_{0}^{2}{\epsilon}_{y}^{\mathbf{G}-\mathbf{G}\text{'}}& -\left({G}_{x}+{k}_{x}\right)\left({G}_{y}^{\prime}+{k}_{y}\right){\gamma}^{\mathbf{G}-\mathbf{G}\prime}\\ -\left({G}_{y}+{k}_{y}\right)\left({G}_{x}^{\prime}+{k}_{x}\right){\gamma}^{\mathbf{G}-\mathbf{G}\prime}& -\left({G}_{y}+{k}_{y}\right)\left({G}_{y}^{\prime}+{k}_{y}\right){\gamma}^{\mathbf{G}-\mathbf{G}\prime}+{k}_{0}^{2}{\epsilon}_{x}^{\mathbf{G}-\mathbf{G}\prime}\end{array}\right]\phantom{\rule{.2em}{0ex}}\left[\begin{array}{c}{E}_{y}^{\mathbf{G}\prime}\\ {E}_{x}^{\mathbf{G}\prime}\end{array}\right],$$

where *G*_{x}*, G*_{y}*, k*_{x}
, and *k*_{y}
denote the *x* and *y* components of the **G** and **k** vectors respectively.

The equations (24) and (25) define the field vectors **Ē** and **H̄** and matrices **R̄****E** and **R̄**_{H}
from Eq. (11) and (12). They are necessary for analytical solution of Maxwell equations in z direction and determination of eigenmodes.

#### 2.3. Determination of eigenmodes

Consider Eq. (11) and (12). They can be combined to a second order differential equation for either the electric or the magnetic field

where **Q̄**_{E}
=**R̄**_{H}**R̄**_{E}
and **Q̄**_{H}**=R̄E**_{R}**̄**_{H}
. In general, these equations cannot be solved analytically, provided **Q̄ _{E}** and

**Q̄**are dense matrices. However, it is possible to overcome this difficulty by transforming the field to the principal axes [20]

_{H}so that

${\overline{\Gamma}}^{2}$
being a diagonal matrix. The transformation is performed by determining ${\overline{\Gamma}}^{2}$
as a set of the eigenvalues of the matrices **Q̄**_{E}
and **Q̄**_{H}
and **T**_{E}
and **T**_{H}
as their eigenvectors. It is worth noting that for uniform layers (e.g. claddings) **Q̄ _{E}**=

**Q̄**=${\overline{\Gamma}}^{2}$ which can further reduce the computation time.

_{H}In the transformed coordinates the Eq. (26) and (27) can be written as

which have the analytical solutions of the form

The parameters **A, B, C**, and **D** are determined by the admittance transfer procedure to satisfy the condition for fields continuity at layer boundaries in a multi-layered structure. This procedure is described in details in refs. [23, 19], so here we only mention its main points and the key equations. The main idea of admittance transfer is the evaluation, for each layer, of an admittance matrix **Ȳ**^{n}(*k*
_{0},**k**) so that

where *n* denotes the number of the layer (**Ē**
^{n}
and **H̄**^{n} are fields at the boundary between *n*-th and (*n*+1)-th layer). The admittance matrix can be computed using the iterative relation

where ${\overline{y}}_{1}^{n}$=-*i*
$\overline{\Gamma}$
^{n} tan^{-1}($\overline{\Gamma}$
^{n}
d
^{n}
) and ${\overline{y}}_{2}^{n}$=*i*
$\overline{\Gamma}$
^{n} sin^{-1}($\overline{\Gamma}$
^{n}*d*^{n}
), $\overline{\Gamma}$
^{n}
, **T**_{E}^{n} and **T**_{H}^{n} being the matrices defined in Eq. (30) for *n*-th layer with the thickness *d*^{n}
. To commence the iterative computation of the admittance matrices we need to assume the electric field at the boundary to be equal to zero **Ē**^{0}=0 and with this condition admittance matrix for the first layer can be written as

The above computation can be performed starting from both upper and lower boundary of the analyzed structure, divided into two parts by the matching interface called also the matching plane. At this interface the condition of fields continuity forms the following implicit eigenvalue problem

where **Ē** is the electric field at the interface. The solution of this eigenproblem gives *k*
_{0}(**k**) relation in the photonic crystal slab.

In most cases, especially in the case of leaky modes, the assumption that **Ē**=0 at the boundaries is unphysical. In order to alleviate it we introduce absorbing perfectly matched layers (PML) [4]. For our new form of the GTL equations (Eq. (24) and (25)) this is straightforward and does not increase the numerical effort.

## 3. Eigenmodes and field distribution in photonic crystal slab

In order to test the plane-wave admittance method we calculate the eigenmodes of an exemplary photonic crystal slab suspended in air, having the refractive index *n*_{R}
=3.5, with circular air holes arranged over a rectangular lattice. The slab thickness is *h*=0.6*a* and the hole radius *r*=0.3*a*, where *a* is an arbitrary lattice constant (see Fig. 1). Such structure has been analyzed in [24], which serves us as a benchmark. To remove the spurious air-confined modes and for accurate determination of the leaky modes radiating from the slab, we add a PML with thickness 0.5*a* at the distance 4.2*a* from the slab. The PML permittivity and permeability are diagonal tensors equal to *ε*
_{PML}=*µ*
_{PML}=diag([*s*_{z}*,s*_{z}, ${s}_{z}^{-1}$
]), where *s*_{z}
=1-2*i*, which is chosen to provide large absorption in the PMLs without too large number of spurious modes emerging.

As the structure possesses an inversion symmetry in the z direction we reduce our calculation to the half of it. The disadvantage of this procedure is the inability to determine odd (TM-like [25]) modes, as for such modes the vector **Ē** from Eq. (38) is equal to zero. However, the analysis of odd modes is straightforward and requires only different choice of matching interface, namely out of symmetry plane. The alternative approach would be the choice of the magnetic field instead of the electric one in Eq. (38).

The structure is divided into three layers uniform in the z direction—the slab (actually only a half of it is considered, as the matching plane is located in its center), the air layer and the PML. As the air and the PML are uniform layers, the computation is very efficient due to the abovementioned fact that for these layers the **Q**_{H}
matrices are already diagonal. As for the core layer the **R̄**_{E}
and **R̄**_{H}
matrices are computed using the Eq. (24) and (25) with the ${\epsilon}_{x}^{G}$, ${\epsilon}_{y}^{G}$*, κ*^{G}, ${\mathit{\mu}}_{x}^{\mathbf{G}}$, ${\mathit{\mu}}_{y}^{\mathbf{G}}$*, γ*^{G}
determined analytically as Fourier expansion coefficients [2]. In order to improve convergence the electric and magnetic constants are smoothed by convolution with Gaussian function

where σ is some small parameter which value has to be chosen empirically. For our calculations we have chosen its value as 0.001*a*
^{2}.

The solution of Eq. (38) is searched, for an arbitrary wavevector **k**, by the determination of the complex roots of a function *F*(*k*
_{0}), which is defined as

Figure 2 presents the plot of *F*(*k*
_{0}) as a function of the normalized frequency *ak*_{0}
/(2*π*) for **k**=[0.25*π*/*a*, 0]. As can be seen from this figure the function charval(*k*
_{0}) has significant minima, which values drop to zero for the eigenfrequencies (because of use of the absorbing PML the eigenfrequencies have some small imaginary part corresponding to absorption, which is not visible in the graph). The dotted lines denote the normalized frequencies of the following modes taken from the literature [24]. As can be seen from Fig. 2 our solutions match very well to them.

In the general case, the complex roots of the function *F*(*k*
_{0}) have to be determined by some root finding algorithm. We use for this purpose the browsing of the range of interest combined with the Newton-Raphson method. However, numerical problems can arise when the holes are very narrow, as one can see in the Fig. 2 (e.g. one around the frequency 0.08·2*π/a*) which may cause missing of some of the modes and therefore requires much denser sampling. This effect is especially strong for low frequencies.

After determining the eigenmodes the electric field at the matching interface can be obtained straightforward as the eigenvector from the Eq. (38). Consequently the magnetic field is computed using Eq. (35). The analytical solutions (33) and (34) allow to determine both fields in the whole three-dimensional domain.

Figure 3 represents the light intensity at the interface and its three dimensional distribution in the slab for one guided and two leaky modes. Additionally, the electric and magnetic field vectors are presented. As only even modes are analyzed, the magnetic field is perpendicular to the interface and the electric field is parallel to it.

By analyzing the spatial distribution of the electromagnetic field energy density we can easily distinguish between guided and leaky modes. For guided modes the energy density is exponentially dropping to zero with increasing the distance, while for the leaky modes it drops to a certain non-zero value which then remains constant through the whole air layer (Fig. 4). The computation of energy density can also help to select and reject spurious modes [24] that can emerge around the light-line for some wavevectors. These spurious modes are in large part confined into the PML, hence the ratio of *f*=*P*
_{PML}/*P*
_{core}, where *P*
_{PML} and *P*
_{core} is an average energy density for the plane in the middle of PML and the core respectively, is significantly larger for the spurious modes than for the real ones. For any plane located at the constant *z, P(z)* can be easily calculated from the formula

$$=\frac{1}{4}\sum _{\mathbf{G}}\sum _{\mathbf{G}\prime}\left[{\left({E}_{x}^{\mathbf{G}}\right)}^{*}{E}_{x}^{\mathbf{G}\prime}{\epsilon}_{x}^{\mathbf{G}-\mathbf{G}\prime}+{\left({E}_{y}^{\mathbf{G}}\right)}^{*}{E}_{y}^{\mathbf{G}\prime}{\epsilon}_{y}^{\mathbf{G}-\mathbf{G}\prime}+{\left({D}_{z}^{\mathbf{G}}\right)}^{*}{D}_{z}^{\mathbf{G}\prime}{\kappa}^{\mathbf{G}-\mathbf{G}\prime}\right]$$

$$+\frac{1}{4}\sum _{\mathbf{G}}\sum _{\mathbf{G}\prime}\left[{\left({H}_{x}^{\mathbf{G}}\right)}^{*}{H}_{x}^{\mathbf{G}\prime}{\mu}_{x}^{\mathbf{G}-\mathbf{G}\prime}+{\left({H}_{y}^{\mathbf{G}}\right)}^{*}{H}_{y}^{\mathbf{G}\prime}{\mu}_{y}^{\mathbf{G}-\mathbf{G}\prime}+{\left({B}_{z}^{\mathbf{G}}\right)}^{*}{B}_{z}^{\mathbf{G}\prime}{\gamma}^{\mathbf{G}-\mathbf{G}\prime}\right]$$

where ${D}_{z}^{\mathbf{G}}$ and ${B}_{z}^{\mathbf{G}}$ are determined from Eq. (7) and (8) and are equal to

## 4. Numerical results and discussion

We consider the same structure as in the previous section and compute the eigenfrequencies taking different numbers of planewaves for the two first modes at wavevector **k**=[0.25*π/a*, 0], The results, computation times on a personal computer (Athlon 2 GHz) and relative errors are listed in Tab. 1 and depicted in Fig. 5. They show very good convergence as the relative error

does not exceed 1% even for as small number of planewaves as 25 and drops below 0.1% for 225 plane waves in case of the guided first-order mode and for only 49 planewaves for the leaky second-order mode. For comparison the purely numerical plane wave expansion method requires more than 800 plane-waves for comparable convergence in a similar structure [24]. Thus the required computer memory, proportional to the square of the number of planewaves, is about 100 times larger than in case of our method.

The computed band diagram for even modes is presented in Fig. 6(a). In Fig. 6(b) we compare our results (blue dotted curves) with the band diagram obtained by full 3D plane-wave expansion method (red square curves) which we reproduce from ref. [24] (Fig. 8a). As one can see the results match very well for most of the bands. There are minor differences for the high order modes, however our results are in better agreement with the FDTD computations presented also in ref. [24].

In conclusion, we have presented a new method for determining the electromagnetic modes of photonic multilayer structures, which combines the plane wave expansion method with the method of lines. The additional application of perfectly matched layers as absorbing boundary conditions allowed to perform correct calculations for leaky modes and almost totally reduced spurious air-confined ones. The convergence of the method was improved by smoothing of the dielectric and magnetic tensors and proved to be very good even for small number of planewaves. Thus, the computations of band structure in photonic crystal slabs and general electromagnetic analysis of multi-layered structures can be performed with high efficiency.

## Acknowledgments

This work is the result of Short Term Scientific Mission supported by the COST Action P11 “Physics of linear, non-linear and active photonic crystals”. Additional contribution was provided by the Polish Ministry of Science in the SPB grant DWM 86 associated to COST P11 and the EC Network of Excellence on Micro-Optics “NEMO”.

We would like to thank Dr Tomasz Czyszanowski from Institute of Physics at the Technical University of Lodz for his help in completing this work.

## References and links

**1. **J. D. Joannopoulos, R. D. Meade, and a J. N. Winn, *Photonic Crystals: Molding the Flow of Light* (Princeton University Press, Princeton, 1995).

**2. **K. Sakoda, *Optical Properties of Photonic Crystals* (Springer-Verlag, Berlin2001).

**3. **N. Yokouchi, A. J. Danner, and K. D. Choquette, “Two-Dimensional Photonic Crystal Confined Vertical-Cavity Surface-Emitting Lasers,” IEEE J. Sel. Top. Quantum Electron. **9**, 1439–1445 (2003). [CrossRef]

**4. **A. Taflove, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Artec House Inc., Boston, 1995).

**5. **S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, “Photonic band gaps in periodic dielectric structures: The scalar-wave approximation,” Phys. Rev. B **46**, 10650–10656 (1992). [CrossRef]

**6. **S. Johnson and J. D. Joannopoulous, “Block Iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express **8**, 173–190 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [CrossRef] [PubMed]

**7. **S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express **11**, 167–175 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-2-167. [CrossRef] [PubMed]

**8. **H. S. Sözüer and J. W. Haus, “Photonic bands: convergence problems with the plane-wave method,” Phys. Rev. B **45**, 13962–13972 (1992). [CrossRef]

**9. **Ch. Sauvan, Ph. Lalanne, and J. P. Hugonin, “Truncation rules for modelling discontinuities with Galerkin method in electromagnetic theory,” Opt. Quantum Electron. **36**, 271–284, 2004. [CrossRef]

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

**11. **A. Ferrando, E. Silvestre, J. J. Miret, P. Andres, and M. V. Andres, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. **24**, 276–278 (1999). [CrossRef]

**12. **J. M. Pottage, D. Bird, T. D. Hedley, J. C. Knight, T. A. Birks, P. S. J Russell, and P. J. Roberts, “Robust photonic band gaps for hollow core guidance in PCF made from high index glass,” Opt. Express **11**, 2854–2861 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-22-2854. [CrossRef] [PubMed]

**13. **S. Shi, C. Chen, and D.W. Prather, “Revised plane wave method for dispersive material and its application to band structure calculations of photonic crystal slabs,” Appl. Phys. Lett. **86**, 043104 (2005), http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=VIRT01000011000004000041000001. [CrossRef]

**14. **P. Lalanne, “Electromagnetic Analysis of Photonic CrystalWaveguides Operating Above the Light Cone,” IEEE J. Quantum Electron. **38**, 800–804 (2002). [CrossRef]

**15. **M. Qiu, “Effective index method for heterostructure-slab-waveguide-based two-dimensional photonic crystals,” Appl. Phys. Lett. **81**1163–1165 (2002). [CrossRef]

**16. **P. Bienstman, “Two-stage mode finder for waveguides with a 2D cross-section,” Opt. Quantum Electron. **36**, 5–14, 2004. [CrossRef]

**17. **K. Ohtaka, J. Inoue, and S. Yamaguti, “Derivation of the density of states of leaky photonic bands,” Phys. Rev. B **70**035109 (2004). [CrossRef]

**18. **S. F. Helfert and R. Pregla, “Efficient Analysis of Periodic Structures,” J. Lightwave Technol. **16**, 1694–1702 (1998). [CrossRef]

**19. **O. Conradi, S. F. Helfert, and R. Pregla, “Comprehensive Modeling of Vertical-Cavity Laser-Diodes by the Method of Lines,” IEEE J. Quantum Electron. **37**, 928–935 (2001). [CrossRef]

**20. **S. F. Helfert, A. Barcz, and R. Pregla, “Three-dimensional vectorial analysis of waveguide structures with the method of lines,” Opt. Quantum Electron. **35**, 381–394 (2003). [CrossRef]

**21. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**22. **H. Derudder, F. Olyslager, D. De Zuter, and S. Van den Berghe, “Efficient Mode-Matcing Analysis of Discontinuities in Finit Planar Substrates Using Perfectly Matched Layers,” IEEE Trans. Antennas and Propagation **49**, 185–195 (2001). [CrossRef]

**23. **T. Czyszanowski, “Comparative Analysis of Validity Limits of Scalar and Vector Approaches to Optical Fields in Diode Lasers,” Ph.D. diss., Technical University of Lodz (2004).

**24. **S. Shi, C. Chen, and D. W. Prather, “Plane-wave expansion method for calculating band structure of photonic crystal slabs with perfectly matched layers,” J. Opt. Soc. Am. A **21**, 1769–1775 (2004) [CrossRef]

**25. **S. G. Johnson, S. Fan, P. R. Villeneuve, J. D. Joannopoulous, and L. A. Kolodziejski, “Guided modes in photonic crystal slabs,” Phys. Rev. B **60**, 5751–5758 (1999). [CrossRef]