Computing Photonic Crystal Defect Modes by Dirichlet-to-Neumann Maps

We develop an efficient numerical method for computing defect modes in two dimensional photonic crystals based on the Dirichletto-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 OCIS codes: (050.5298) Photonic crystals, (000.4430) Numerical approximation and analysis 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 twodimensional dielectric lattices,” Phys. Rev. Lett. 67, 2017-2020 (1991). 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). 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). 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). 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). 7. X. P. Feng and Y. Arakawa, “Defect modes in two-dimensional triangular photonic crystals,” Japanese Journal of Applied Physics 36, L120-L123, (1997). 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). 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). 10. D. C. Dobson, “An efficient method for band structure calculations in 2D photonic crystals,” J. Comput. Phys. 149, 363-376 (1999). 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). 12. H. Y. D. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE Trans. Microwave Theory Tech. 44, 2688-2695 (1996). 13. C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of twodimensional photonic crystals,” Opt. Express 12, 1397-1408 (2004). 14. P. J. Chiang and C. P. Yu and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75, 026703 (2007). 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). #87214 $15.00 USD Received 4 Sep 2007; revised 10 Oct 2007; accepted 17 Oct 2007; published 18 Oct 2007 (C) 2007 OSA 29 October 2007 / Vol. 15, No. 22 / OPTICS EXPRESS 14454 16. K. Sakoda and H. Shiroma, “Numerical method for localized defect modes in photonic lattices,” Phys. Rev. B 56, 4830-4835 (1997). 17. K. Sakoda, “Numerical study on localized defect modes in two-dimensional triangular photonic crystals,” Journal of Applied Physics, 84, 1210-1214 (1998). 18. V. Kuzmiak and A. A. Maradudin, “Localized defect modes in a two-dimensional triangular photonic crystal,” Phys. Rev. B 57, 15242-15250 (1998). 19. N. Stojı́c, J. Glimm, Y. Deng and J. W. Haus, “Transverse magnetic defect modes in two-dimensional triangularlattice photonic crystals,” Phys. Rev. E 64, 056614 (2001). 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). 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). 22. R. Moussa, L. Salomon, F. de Fornel and H. Aourag, “Numerical study on localized defect modes in twodimensional lattices: a high Q-resonant cavity,” Physica B Condensed Matter 338, 97-102 (2003). 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). 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). 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). 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). 28. J. H. Yuan, Y. Y. Lu, X. Antoine, “Modeling photonic crystals by boundary integral equations and Dirichlet-toNeumann 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). 30. T. Lu and D. Yevick, “A vectorial boundary element method analysis of integrated optical waveguides,” J. Lightw. Technol. 21, 1793-1807 (2003). 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). 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). 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).


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 qualityfactor 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 num-ber 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 twodimensional 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.

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 .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 ma 1 and ma 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 ma 1 and ma 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.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 twodimensional 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.

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 Λ = Λ 2 Λ −1 1 .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 equation connecting the field on 11 edges of two unit cells Ω 1 and Ω 2 .Similarly, for edge 2 in Fig. 5, we have 22 u 2 + Λ (1) 24 u 4 + Λ (1) 26 u 6 = Λ (3) 53 u 14 + Λ (3) 56 u 11 .
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.
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 u i is a column vector of length N for the field on the i-th edge, U d = [u 1 ; u 2 ; ...; u 6 ] is a column vector of length K = 6N for the boundary of the defect cell, U n = [u 7 ; u 8 ; ...; u 30 ] is a column vector of length 24N for other edges in the first ring and U f 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 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 U f and U n in terms of U d .
In fact, we have two matrices D 1 and D 2 satisfying then U n and U f 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 ω.

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 s 1 (B) = 0, (11) 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 ω (2) After that, we evaluate s 1 at both ω 2 and ω 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.

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 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  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.

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.

Fig. 1 .
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.

Fig. 3 .
Fig. 3. Truncated domain with two rings around the defect cell for a triangular lattice.

Fig. 4 .
Fig. 4. Ordering of the edges and the sampling points on the edges of a hexagon unit cell for discretization of the DtN map.

Fig. 5 .
Fig. 5. Three neighboring hexagon unit cells and a possible ordering of the edges.

Fig. 10 .
Fig. 10.Normalized frequencies of the defect modes versus the ratio R d /a between the radius of the defect rod and the lattice constant.