Open Access
Add to BibSonomyAbstract
We develop an efficient numerical method for computing defect modes in two dimensional photonic crystals based on the Dirichlet-to-Neumann (DtN) maps of the defect and normal unit cells. The DtN map of a unit cell is an operator that maps the wave field on the boundary of the cell to its normal derivative. The frequencies of the defect modes are solved from a condition that a small matrix is singular. The size of the matrix is related to the number of points used to discretize the boundary of the defect cell. The matrix is obtained by solving a linear system involving only discrete points on the edges of the unit cells in a truncated domain.
©2007 Optical Society of America
1. Introduction
In a photonic crystal (PhC) [1] with a point defect, defect modes [2, 3, 4] may exist for some frequencies in bandgaps. These localized eigenmodes of the wave field give rise to high quality-factor micro-cavities and have found applications in photonic circuits, light emitting devices, etc. In the ideal case, the defect modes are analyzed in an infinite PhC with a local perturbation of the otherwise perfectly periodic refractive index function. The frequencies of the defect modes are the eigenvalues of an eigenvalue problem formulated on the whole space. Since the eigenfunctions decay to zero away from the point defect, the computation domain is usually truncated. For highly extended defect modes that occur near band edges, the domain truncation approach is not efficient and the fictitious source superposition method [5] can be used. The defects modes can be calculated by the plane wave expansion method using the supercell approach that assumes periodic boundary conditions on the boundary of the truncated domain [6, 7]. In fact, the supercell approach allows us to use many existing methods for computing band structures [8, 9, 10, 11, 12, 13, 5, 14] to calculate the defect modes. Meanwhile, a simple zero boundary condition on the boundary of the truncated domain may be convenient, for example, in the finite element method [15]. The defect modes can also be calculated in time domain using the finite difference time domain (FDTD) method [16, 17, 18, 19, 20] or finite element time domain method [21, 15]. The eigenvalue problem of the defect modes usually leads to large matrices. In the plane wave expansion method [8, 9, 6, 7], this is caused by the large number of terms in Fourier series approximations to the wave field and the dielectric function. In the finite element method [10, 11, 15], a large number of elements are needed to discretize the truncated domain. The defect modes actually correspond to interior eigenvalues that are more difficult to solve than the extreme eigenvalues by existing numerical methods. Iterative methods such as the Krylov subspace methods are efficient only for a few largest and smallest eigenvalues. Furthermore, if the medium is dispersive, the eigenvalue problem becomes nonlinear and more difficult to solve.
Recently, the Dirichlet-to-Neumann (DtN) maps of unit cells have been used to develop efficient numerical methods for analyzing photonic crystal problems [23, 24, 25, 26, 27, 28]. The DtN map is an operator that maps the wave field on the boundary of the unit cell to its normal derivative. If we have K special solutions in the unit cell and choose K sampling points on the boundary of the unit cell, we can approximate the DtN map by a K × K matrix. For two-dimensional photonic crystals composed of parallel circular cylinders (air-holes or dielectric rods) in a lattice, the DtN map of a unit cell can be easily constructed using cylindrical waves as the special solutions. In the frequency domain, the DtN map allows us to reduce the computations to the edges of the unit cells for boundary value problems such as those related to finite PhCs [25, 26, 27]. It also allows us to formulate linear eigenvalue problems on edges of the unit cell for band structure problems even when the medium is dispersive [23, 24]. The formulations developed in [23, 24] assume that the frequency is a given parameter (not the eigenvalue) and calculate the dispersion relationships from linear eigenvalue problems, where the eigenvalue is related to one component of the Bloch wave vector. These formulations cannot be used to calculate defect modes, since these modes do not have a Bloch wave vector and the eigenvalue cannot be switched.
In this paper, we show that the concept of DtN map is still useful for computing defect modes. Using the DtN maps of the defect and normal unit cells, we formulate a nonlinear condition on the boundary of the defect cell for the eigenfrequency. The principle of our method is similar to boundary integral equation methods for eigenvalue problems [29, 30] and to mode matching method for computing waveguide modes [31]. These methods turn a linear eigenvalue problem to a nonlinear eigenvalue problem that involves a much smaller matrix. If the DtN maps are approximated by K × K matrices, the frequency of a defect mode is solved from the condition that a K × K matrix B is singular. The matrix B depends on the frequency and it is constructed from the DtN maps of the normal and defect cells. Although an iterative method is needed, our method is very efficient, since K is typically very small. Notice that the matrix in the standard formulation corresponds to a discretization of the supercell. Here, the matrix B corresponds to a discretization of the boundary of the defect cell. A large truncated domain is still used in our formulation, but only the wave field on edges of the unit cells are involved in the manipulation, because the DtN maps contain complete information about the wave field in the interior of the unit cell if the field on the edges is known. The matrix B is obtained when the wave field on all edges in the truncated domain, except the edges of the defect cell, are eliminated. We demonstrate our method by numerical examples.
2. Eigenvalue problems for defect modes
For two dimensional (2D) problems, the governing equation is
where k 0 = ω/c is the free space wavenumber, ω is the angular frequency, c is the speed of light in vacuum, n = n(x) is the refractive index function and x = (x,y). For the E polarization, u is the z-component of the electric field and ρ = 1. For the H polarization, u is the z-component of the magnetic field and ρ = n 2. For a 2D PhC [1], the refractive index function is periodic in two linearly independent directions. We have two primitive translation vectors a 1 and a 2, such that
where l 1 and l 2 are arbitrary integers. The xy plane can be divided into unit cells given periodically on the lattice defined by the two vectors a 1 and a 2. The unit cell can be chosen as the parallelogram defined by a 1 and a 2, but for the triangular lattice, a hexagon unit cell is often more convenient. The periodic variation of the refractive index gives rise to bandgaps which are intervals of the frequency, such that propagating Bloch waves of the form
where Φ satisfies the periodic condition (2), do not exist for any real Bloch wave vector k. For the defect mode problem, we assume that n(x) is given differently from (2) for x is a small bounded region specified by one or a few unit cells. An example is shown in Fig. 1, where the PhC is composed of a triangular lattice of parallel dielectric rods and the defect cell contains a rod with a different radius and/or different refractive index. The structure is supposed to be infinite and periodic except at the defect cell, although only four rings around the defect cell are shown in Fig. 1. Our problem is to find the frequencies in the bandgaps such that the homogeneous Helmholtz equation (1) has non-zero solutions that decay to zero at infinity. This is an eigenvalue problem formulated on the entire xy plane, where the eigenvalue is ω 2 or k 2 0.
Fig. 1. Point defect in a triangular lattice of dielectric rods. The defect cell contains a rod with a different radius and/or refractive index.
For practical numerical computations, the eigenvalue problem is solved on a finite domain S. In the supercell approach, S corresponds to a unit cell with translation vectors m a 1 and m a 2, where m is a positive integer. For the triangular lattice and m = 3, two different possible supercells are shown in Fig. 2. The parallelogram supercell contains m 2 = 9 elementary parallelogram unit cells (one defect cell and eight normal unit cells). The hexagon supercell contains one hexagon defect cell, six normal hexagon cells and six corner regions. Each of these corner regions corresponds to one third of a normal hexagon unit cell. Corresponding to duplicating the supercell on the lattice specified by the vectors m a 1 and m a 2, a periodic boundary condition is often used and it is convenient for the plane wave expansion method [3, 6, 7]. For other methods, such as the finite element method [15], a simple zero boundary condition can also be used. Since the eigenfunction decays to zero rapidly as the distance to the defect cell increases, the boundary condition makes little difference when the supercell is sufficiently large.
In our method, the xy plane is truncated to a domain that contains a few rings of normal cells around the defect cell. For the triangular lattice and using hexagon unit cells, the domain truncated to one ring contains exactly one defect cell and six normal cells. This domain is in fact different from the hexagon supercell for m = 3 as shown in Fig. 2, since the six corner regions are not included. The truncated domain with two rings around the defect cell is shown in Fig. 3. On the boundary of the truncated domain, we use a simple zero boundary condition.
If we know the frequency ω, the general solution of the Helmholtz equation (1) in a unit cell (normal or defect) including a circular cylinder can be written down analytically in terms of the cylindrical waves. The general solution allows us to find an operator (the DtN map) that maps the wave field on the boundary of the unit cell to its normal derivative there. In section 3, we use the DtN maps of the normal and defect unit cells to establish the condition
for the frequencies of the defect modes. Here, Γd is the boundary of the defect cell, B is an operator that depends on the frequency. In practice, with the general solution in a unit cell approximated by a finite sum of K special solutions, the DtN maps of the unit cells and the operator B are all approximated by K × K matrices. For a given ω, the matrix B can be calculated efficiently by solving a linear system involving only the wave field on the edges of the unit cells in the truncated domain. Since we are solving the frequencies of the defect modes, we need to search ω from the condition that the matrix B is singular. In section 4, we develop a modified version of the secant method to solve ω using the smallest singular value of B. Overall, we have turned the original linear eigenvalue problem (for a non-dispersive medium) on a two-dimensional domain to a nonlinear eigenvalue problem on a curve (the boundary of the defect cell). The method is iterative, but each iteration can be evaluated efficiently.
3. Formulation on the boundary of the defect cell
To reduce the eigenvalue problem to the boundary of the defect cell, we make use of the DtN maps of the normal and defect unit cells. Let Ω be a unit cell and Γ be the boundary of Ω, the DtN map of Ω is the operator Λ satisfying
where u is an arbitrary solution of the Helmholtz equation (1) and ν is a unit normal vector of Γ. Therefore, the DtN map Λ is the operator that maps the wave field on the boundary to its normal derivative. To find a matrix approximation to Λ, we follow the approach developed in [25, 23]. For an integer K, we choose K points on the boundary Γ, say x 1, x 2, ..., x K, and also assume that the general solution of the Helmholtz equation (1) in Ω can be approximated by a sum of K special solutions. That is
where ϕj satisfies (1). Evaluating the above at the K points on the boundary, we obtain
where Λ1 is a K × K matrix whose (i,j) entry is ϕj(x i). We can also evaluate the normal derivative at these K points. This leads to
where νi = ν(x i) is a unit normal vector of Γ at x i and Λ2 is a K × K matrix whose (i,j) entry is ∂νiϕj(x i). Therefore, we approximate the DtN map Λ by the matrix
If the unit cell contains a circular cylinder, the special solutions used in (5) can be chosen as the cylindrical waves which are available analytically [25, 23]. If the unit cell contains a cylinder with a more general cross section, we consider the scattering problem of this single cylinder on the entire plane (with the constant refractive index outside the cylinder extended to infinity) and let ϕj be the solution in connection with a plane incident wave. The different special solutions in (5) correspond to the same scatterer and different incident waves, and they can be solved together efficiently by a boundary integral equation method [28]. On the other hand, if the unit cell contains a few circular cylinders, we can use the multipole method [32] to find the special solutions [27]. Since the size of the unit cell is typical smaller than the free space wavelength, the integer K can be quite small. For a hexagon unit cell, we have K = 6N where N is the number of sampling points on each edge. Numerical examples in section 5 indicate that accurate solutions can be obtained for N ≤ 8.
The DtN maps of the unit cells allow us to write down one equation for each interior edge of the truncated domain S by matching the normal derivative of the wave field there. If N points are used on the edge, the wave field there will be approximated by a column vector of length N and the equation on that edge is in fact a system of N equations. To write down these equations, we first need to specify an ordering for the matrix approximation of the DtN map Λ. Consider the hexagon unit cell shown in Fig. 4, the wave field on the boundary Γ, i.e., u|Γ, is given as the column vector [u 1;u 2;...;u 6] (following the MATLAB notation) where u 1, u 2, ..., u 6 are column vectors of length N and they approximate the wave field u on the six edges following a counterclockwise ordering. To ease the matching on a common edge of two neighboring unit cells, we require that the sampling points on opposite edges follow the same order. More precisely, we have ordered the points x i, for 1 ≤ i ≤ 6N, such that on each edge the x coordinates increase as i is increased and if the x coordinates are constant (on the vertical edges) the y coordinates must increase. The case for N =3 is shown in Fig. 4. Such an ordering may seem to be less natural as the clockwise or counter-clockwise ordering, but it allows us to avoid a permutation when writing down the equation for each edge. We also assume that the normal vectors on opposite edges point to the same direction (one into the cell Ω and the other out). With these specifications, we write down the DtN map Λ as a 6 × 6 block matrix where each block is an N × N matrix. We have
To write down an equation for a given edge, we evaluate the normal derivative using the DtN maps of the two neighboring unit cells. To illustrate the procedure, we consider three hexagon unit cells with possibly different refractive index profiles. As shown in Fig. 5, the unit cells Ω1, Ω2 and Ω3 contain red, yellow and green cylinders, respectively. Let the DtN maps of the three unit cells Λ(1), Λ(2) and Λ(3) be partitioned as 6×6 block matrices accordingly. For edge 1 in Fig. 5, we have
Notice that edge 1 in the second unit cell Ω2 corresponds to the fourth edge, thus the fourth row of the DtN map Λ(2) appears in the second equation above. The above gives rise to one equation connecting the field on 11 edges of two unit cells Ω1 and Ω2. Similarly, for edge 2 in Fig. 5, we have
In the above, we have assumed that the three unit cells have the same background medium outside the cylinders. If this is not true, then the above equations are valid only for the E polarization. For the H polarization, we modify the above equations by matching n -2 ∂νu on the edges, where n is the refractive index of the background medium in each unit cell.
Fig. 4. Ordering of the edges and the sampling points on the edges of a hexagon unit cell for discretization of the DtN map.
Let us now order all interior edges of the truncated domain S so that the first six edges are the edges of the defect unit cell and the next 24 edges are the remaining edges of the six closest normal unit cell (i.e. the first ring). This implies that the wave field on all these edges is given as a column vector
where ui, is a column vector of length N for the field on the i-th edge, Ud = [u 1;u 2;...;u 6] is a column vector of length K = 6N for the boundary of the defect cell, Un = [u 7;u 8;...;u 30] is a column vector of length 24N for other edges in the first ring and Uf is the column vector for the remaining interior edges in the truncated domain S. Using the DtN maps of the normal and defect unit cells and following the procedure outlined above, we arrive at the following linear system:
where A 11 is K×K and A 22 is (24N)×(24N). The size of the matrix A is (MN)×(MN), where M is the total number of interior edges in the truncated domain S. In terms of N × N blocks, the matrix A contains at most 11 blocks in each row and column. Therefore, the matrix A is sparse with at most 11N non-zero entries in each row or column. Compared with a matrix that discretizes the truncated domain S or a related supercell, the matrix A is much smaller, since only the edges are involved. Nmerical examples in section 5 indicate that accurate solutions can be obtained with N as small as 5. This corresponds to 15 points for each unit cell, since each interior edge is shared by two unit cells. If the plane wave expansion method or a finite element method is used, a few hundreds points or basis functions are needed for each unit cell. Although the defect modes can be solved from the condition that A is singular, it is more efficient to use a reduced condition involving a much smaller matrix. We can solve Uf and Un in terms of Ud. In fact, we have two matrices D 1 and D 2 satisfying
then Un and Uf are given by
Therefore, the system (7) is reduced to
The above equation is the discrete version of (3). Since the DtN maps are constructed for a given ω, the matrix B depends on ω.
4. Searching defect mode frequencies
The frequencies of the defect modes can be determined from the condition that the matrix B in Eq. (10) is singular. Notice that the larger matrix A in Eq. (7) is also singular at these frequencies, but it is much easier to use the smaller matrix B. One way to find the defect mode is to solve ω from det B = 0. However, it is well known that the determinant is not a good indicator for near singularity of a matrix. Our approach is to solve the defect mode frequencies from
where s 1 is the smallest singular value of the matrix B, i.e., the square root of the smallest eigenvalue of B*B where B* the transpose and complex conjugate of B.
As a function of ω, s 1 is non-negative and behaves like the absolute value of another function. In particular, the derivative of s 1 with respect to ω is not continuous at its zeros. However, s 1 appears to smoothly connect to - s 1 as ω passes through a zero. This motivates us to develop a modified secant method for solving Eq. (11). In each step, we calculate two approximations. The first approximation is obtained from the standard secant method. Let ω 0 and ω 1 be two initial guesses satisfying s 1(ω 1) ≤ s 1(ω 0), then the first approximation is
To obtain the second approximation, we replace s 1(ω 1) by -s 1(ω 1). Thus
After that, we evaluate s 1 at both ω (1) 2 and ω (2) 2, and choose ω 2 to be the one that gives the smaller value of s 1. The process is then repeated using ω 1 and ω 2 to find ω 3, etc. The iterations can be terminated at some integer k, if s 1(ω k+1) > s 1(ωk) or |ω k+1 - ωk|/|ω k+1| < δ, where δ is a small number for error tolerance.
5. Numerical examples
To illustrate our method, we consider two numerical examples. The first example has been studied by Sakoda [17] and Rodríguez-Esquerre el al [15], and it is related to the experiment by Smith el al [4]. The PhC is a triangular lattice of dielectric rods in air. The refractive index and the radius of the rods are n 1 = 3 and R = 48a/127, respectively, where a is the lattice constant. The defect corresponds to one missing rod. Previous numerical results in [17] and [15] give the defect frequency in three significant digits as ωa/(2πc) = 0.468.
Using a truncated domain S with p = 6 rings around the defect cell and N = 8 sampling points on each edge, we calculate the defect mode frequency with initial guesses: ω 0 a/(2πc) = 0.46 and ω 1 a/(2πc) = 0.47. For the error tolerance δ = 10-12, the modified secant method gives us the solution ωa/(2πc) = 0.46798 in only four iterations. The fast convergence of the modified secant method can be explained by the graph of s 1 as a function of the frequency shown in Fig. 6. We observe that s 1 looks like the absolute value of a smooth function which is close to a straight line. For this problem, the defect mode belongs to the band gap given by 0.415 < ωa/(2πc) < 0.483. We have actually calculated the function s 1 on the entire band gap. It turns out that s 1 has only one zero near 0.468. Next, we check the convergence with respect to the number of sampling points on each edge, i.e, N. For that purpose, we fix p = 6 and δ = 10-12, then calculate the defect mode frequency using different values of N. Let ω (N) be the frequency obtained with N points on each edge, we compute a relative error using the solution obtained with N = 16 as the reference solution, that is
Fig. 6. The smallest singular value s 1 of the matrix B in Eq. (10) versus the normalized frequency ωa/(2πc) for a triangular lattice of dielectric rods with one missing rod in the defect cell.
As shown in Fig. 7, the relative errors appear to decay exponentially and the solutions are more accurate when N is odd. To check the convergence with respect to the number of rings p, we fix N = 8 and δ = 10-12, and calculate the defect mode for p = 3, 4,..., 12. Using the solution obtained with p = 12 as the reference solution, we calculate the relative error for different values of p. The results are shown in Fig. 8. As expected, the relative errors decay exponentially as p is increased. To have six significant digits in the solution, the relative error should be less than 10-6. From Fig. 7 and Fig. 8, it appears that we should choose N ≥ 7 and p ≥ 9. For N = 7 and p = 9, we have ωa/(2πc) = 0.467955.
Fig. 7. Relative errors of the defect mode frequency calculated using different values of N (the number of points on each edge of the unit cells).
After the defect mode frequency is obtained, we can calculate the eigenfunction in a few simple steps. First, we solve Ud from Eq. (10), where Ud represents the wave field of the defect mode on the six edges of the defect unit cell. At the defect mode frequency, the matrix B is singular, then Ud is the eigenvector for the zero eigenvalue of B. Next, we solve Un and Uf from Eq. (7). If the two matrices D 1 and D 2 are saved, we can use Eq. (9) to evaluate Un and Uf. This gives rise to the values of the eigenfunction on all edges of the unit cells. For each unit cell, we can find the coefficients {cj} in Eq. (5) and then evaluate the eigenfunction in its interior. For this numerical example, the eigenfunction is shown in Fig. 9.
Fig. 8. Relative errors of the defect mode frequency calculated using different values of p (the number of rings around the defect cell).
Fig. 9. Eigenfunction (z-component of the electric field) of the defect mode for a triangular lattice of dielectric rods with one missing rod (E polarization).
The second example was first analyzed by Feng and Arakawa [7] by the plane wave expansion method, but they only used a small 3×3 parallelogram supercell as shown in Fig. 2. Later, Sakoda [17] studied this problem with a finite difference time domain method and a larger supercell. The PhC is a triangular lattice of dielectric rods in air. The radius and the dielectric constant of the rods are R = 0.2a and ε 1 = n 2 1 = 13.0, respectively. The defect cell contains a circular rod with the same refractive index and a different radius Rd. The objective is to calculate the defect mode frequencies for various values of Rd/a. The PhC without the defect has a bandgap given by 0.2644 < ωa/(2πc) < 0.4345. In our calculations, we use a truncated domain S with p rings around the defect cell, where 6 ≤ p ≤ 9. When the defect mode frequency is close the band edges, the mode profile spreads out and a larger p is necessary. We also use a larger value of N to verify the calculations obtained with N = 5. For the modified secant method, a fixed error tolerance δ = 10-12 is chosen. Our results are shown in Fig. 10 and they agree well with those of Sakoda [17]. The analysis in [33] indicates that the curves in Fig. 10 should converge exponentially to the band edges. Since we have only calculated a limited number of points on each curve, the correct asymptotics near the band edges are not revealed in the figure.
Fig. 10. Normalized frequencies of the defect modes versus the ratio Rd/a between the radius of the defect rod and the lattice constant.
6. Conclusions
In this paper, we have developed an efficient numerical method for calculating defect modes for two-dimensional (2D) photonic crystals (PhCs) with defects. The method makes use of the Dirichlet-to-Neumann (DtN) maps of the unit cells and determines the defect mode frequencies from the condition that the matrix B in Eq. (10) is singular. The size of matrix B is K × K where K is the number of sampling points used on the boundary of the defect cell. A modification of the secant method is used to search the defect mode frequencies. In each iteration, it is necessary to calculate the matrix B by solving a linear system involving approximations only on the edges of unit cells in the truncated domain. Accurate results can be obtained with a small number of sampling points on each edge of the unit cells.
Our method was presented for a triangular lattice with a single defect cell, but it can be easily modified for square or rectangular lattices. The method can also be extended to PhCs with larger defect areas involving a few neighboring cells. For each edge in the truncated domain, we setup an equation as in section 3 and then reduce the system to the edges of the defect area by elimination. At the moment, the method is limited to pure 2D PhCs which are infinite and invariant in the third direction. The possibility of extending this method to three dimensional structures such as PhC slabs is being explored.
Acknowledgments
This research was partially supported by a grant from the Research Grants Council of Hong Kong Special Administrative Region, China (Project No. CityU 102207).
References and links
1. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light, Princeton University Press, Princeton, NJ. 1995.
2. S. L. McCall, P. M. Platzman, R. Dalichaouch, D. Smith, and S. Schultz, “Microwave Propagation in two-dimensional dielectric lattices,” Phys. Rev. Lett. 67, 2017–2020 (1991). [CrossRef] [PubMed]
3. E. Yablonovitch, T. J. Gmitter, R. D. Meade, A. M. Rappe, K. D. Brommer, and J. D. Joannopoulos, “Donor and acceptor modes in photonic band-structure,” Phys. Rev. Lett. 67, 3380–3383 (1991). [CrossRef] [PubMed]
4. D. R. Smith, R. Dalichaouch, N. Kroll, S. Schultz, S. L. McCall, and P. M. Platzman, “Photonic band structure and defects in one and two dimensions,” J. Opt. Soc. Am. B 10, 314–321 (1993). [CrossRef]
5. S. Wilcox, L. C. Botten, R. C. McPhedran, C. G. Poulton, and C. M. de Sterke, “Modeling of defect modes in photonic crystals using the fictitious source superposition method,” Phys. Rev. E 71, 056606 (2005). [CrossRef]
6. R. R. Villeneuve, S. H. Fan, and J. D. Joannopoulos, “Microcavities in photonic crystals: Mode symmetry, tunability, and coupling efficiency,” Phys. Rev. B 54, 7837–7842 (1996). [CrossRef]
7. X. P. Feng and Y. Arakawa, “Defect modes in two-dimensional triangular photonic crystals,” Japanese Journal of Applied Physics 36, L120–L123, (1997). [CrossRef]
8. K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. 65, 3152–3155 (1990). [CrossRef] [PubMed]
9. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001). [CrossRef] [PubMed]
10. D. C. Dobson, “An efficient method for band structure calculations in 2D photonic crystals,” J. Comput. Phys. 149, 363–376 (1999). [CrossRef]
11. W. Axmann and P. Kuchment, “An efficient finite element method for computing spectra of photonic and acoustic band-gap materials - I. Scalar case,” J. Comput. Phys. 150, 468–481 (1999). [CrossRef]
12. H. Y. D. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE Trans. Microwave Theory Tech. 44, 2688–2695 (1996). [CrossRef]
13. C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express 12, 1397–1408 (2004). [CrossRef] [PubMed]
14. P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75, 026703 (2007). [CrossRef]
15. V. F. Rodríguez-Esquerre, M. Koshiba, and H. E. Hernández-Figueroa, “Finite-element analysis of photonic crystal cavities: Time and frequency domains,” J. Lightw. Technol. 23, 1514–1521 (2005). [CrossRef]
16. K. Sakoda and H. Shiroma, “Numerical method for localized defect modes in photonic lattices,” Phys. Rev. B 56, 4830–4835 (1997). [CrossRef]
17. K. Sakoda, “Numerical study on localized defect modes in two-dimensional triangular photonic crystals,” Journal of Applied Physics , 84, 1210–1214 (1998). [CrossRef]
18. V. Kuzmiak and A. A. Maradudin, “Localized defect modes in a two-dimensional triangular photonic crystal,” Phys. Rev. B 57, 15242–15250 (1998). [CrossRef]
19. N. Stojíc, J. Glimm, Y. Deng, and J. W. Haus, “Transverse magnetic defect modes in two-dimensional triangular-lattice photonic crystals,” Phys. Rev. E 64, 056614 (2001). [CrossRef]
20. S. P. Guo and S. Albin, “Numerical techniques for excitation and analysis of defect modes in photonic crystals,” Opt. Express 11, 1080–1089 (2003). [CrossRef] [PubMed]
21. V. F. Rodríguez-Esquerre, M. Koshiba, and H. E. Hernández-Figueroa, “Finite-element time-domain analysis of 2-D photonic crystal resonant cavities,” IEEE Photon. Technol. Lett. 16, 816–818 (2004). [CrossRef]
22. R. Moussa, L. Salomon, F. de Fornel, and H. Aourag, “Numerical study on localized defect modes in two-dimensional lattices: a high Q-resonant cavity,” Physica B - Condensed Matter 338, 97–102 (2003). [CrossRef]
23. J. H. Yuan and Y. Y. Lu, “Photonic bandgap calculations using Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 23, 3217–3222 (2006). [CrossRef]
24. J. H. Yuan and Y. Y. Lu, “Computing photonic band structures by Dirichlet-to-Neumann maps: The triangular lattice,” Opt. Commun. 273, 114–120 (2007). [CrossRef]
25. Y. X. Huang and Y. Y. Lu, “Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps,” J. Lightw. Technol. 24, 3448–3453 (2006). [CrossRef]
26. Y. H. Huang and Y. Y. Lu, “Modeling photonic crystals with complex unit cells by Dirichlet-to-Neumann maps,” Journal of Computational Mathematics 25, 337–349 (2007).
27. S. J. Li and Y. Y. Lu, “Multipole Dirichlet-to-Neumann map method for photonic crystals with complex unit cells,” J. Opt. Soc. Am. A 24, 2438–2442 (2007). [CrossRef]
28. J. H. Yuan, Y. Y. Lu, and X. Antoine, “Modeling photonic crystals by boundary integral equations and Dirichlet-to-Neumann maps,” submitted for publication.
29. Y. Y. Lu and S.-T. Yau, “Eigenvalues of the Laplacian through boundary integral equations,” SIAM Journal on Matrix Analysis and Applications 12, 597–609 (1991). [CrossRef]
30. T. Lu and D. Yevick, “A vectorial boundary element method analysis of integrated optical waveguides,” J. Lightw. Technol. 21, 1793–1807 (2003). [CrossRef]
31. L. Prkna, M. Hubalek, and J. Ctyroky, “Vectorial eigenmode solver for bent waveguides based on mode matching,” IEEE Photon. Technol. Lett. 16, 2057–2059 (2004). [CrossRef]
32. D. Felbacq, G. Tayeb, and D. Maystre, “Scattering by a random set of parallel cylinders,” J. Opt. Soc. Am. A 11, 2526–2538 (1994). [CrossRef]
33. K. B. Dossou, R. C. McPhedran, L. C. Botten, A. A. Asatryan, and C. M. de Sterke, “Gap-edge asymptotics of defect modes in two-dimensional photonic crystals,” Opt. Express 15, 4753–4762 (2007). [CrossRef] [PubMed]
