OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 15, Iss. 22 — Oct. 29, 2007
  • pp: 14454–14466
« Show journal navigation

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

Shaojie Li and Ya Yan Lu  »View Author Affiliations


Optics Express, Vol. 15, Issue 22, pp. 14454-14466 (2007)
http://dx.doi.org/10.1364/OE.15.014454


View Full Text Article

Acrobat PDF (175 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

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

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

] with a point defect, defect modes [2

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

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

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]

] 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

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]

] 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

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

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]

]. In fact, the supercell approach allows us to use many existing methods for computing band structures [8

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

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

10. D. C. Dobson, “An efficient method for band structure calculations in 2D photonic crystals,” J. Comput. Phys. 149, 363–376 (1999). [CrossRef]

, 11

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

12. H. Y. D. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE Trans. Microwave Theory Tech. 44, 2688–2695 (1996). [CrossRef]

, 13

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]

, 5

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]

, 14

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]

] 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

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]

]. The defect modes can also be calculated in time domain using the finite difference time domain (FDTD) method [16

16. K. Sakoda and H. Shiroma, “Numerical method for localized defect modes in photonic lattices,” Phys. Rev. B 56, 4830–4835 (1997). [CrossRef]

, 17

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

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

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

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]

] or finite element time domain method [21

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]

, 15

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]

]. The eigenvalue problem of the defect modes usually leads to large matrices. In the plane wave expansion method [8

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

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]

, 6

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

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]

], 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

10. D. C. Dobson, “An efficient method for band structure calculations in 2D photonic crystals,” J. Comput. Phys. 149, 363–376 (1999). [CrossRef]

, 11

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]

, 15

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]

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

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

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

30. T. Lu and D. Yevick, “A vectorial boundary element method analysis of integrated optical waveguides,” J. Lightw. Technol. 21, 1793–1807 (2003). [CrossRef]

] and to mode matching method for computing waveguide modes [31

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]

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

ρx(1ρux)+ρy(1ρuy)+k02n2u=0,
(1)

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

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

], the refractive index function is periodic in two linearly independent directions. We have two primitive translation vectors a 1 and a 2, such that

n(x)=n(x+l1a1+l2a2)
(2)

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

u(x)=eikxΦ(x)

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

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]

, 6

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

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]

]. For other methods, such as the finite element method [15

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]

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

Fig. 2. Parallelogram and hexagon supercells for a triangular lattice and m = 3.

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.

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

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

B(ω)u|Γd=0
(3)

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

Λu|Γ=uν|Γ,
(4)

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

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]

, 23

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]

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

u(x)=j=1Kcjϕj(x),
(5)

where ϕj satisfies (1). Evaluating the above at the K points on the boundary, we obtain

[u(x1)u(x2)u(xK)]=Λ1[c1c2cK],

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

[ν1u(x1)ν2u(x2)νKu(xK)]=Λ2[c1c2cK],

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

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

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]

, 23

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]

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

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.

]. On the other hand, if the unit cell contains a few circular cylinders, we can use the multipole method [32

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]

] to find the special solutions [27

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]

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

Λ[u1u2u3u4u5u6]=[Λ11Λ12Λ13Λ14Λ15Λ16Λ21Λ22Λ23Λ24Λ25Λ26Λ31Λ32Λ33Λ34Λ35Λ36Λ41Λ42Λ43Λ44Λ45Λ46Λ51Λ52Λ53Λ54Λ55Λ56Λ61Λ62Λ63Λ64Λ65Λ66][u1u2u3u4u5u6]=[νu1νu2νu3νu4νu5νu6].
(6)

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

u1ν=Λ11(1)u1+Λ12(1)u2+Λ13(1)u3+Λ14(1)u4+Λ15(1)u5+Λ16(1)u6
=Λ41(2)u9+Λ42(2)u10+Λ43(2)u11+Λ44(2)u1+Λ45(2)u7+Λ46(2)u8.

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

u2ν=Λ21(1)u1+Λ22(1)u2+Λ23(1)u3+Λ24(1)u4+Λ25(1)u5+Λ26(1)u6
=Λ51(3)u12+Λ52(3)u13+Λ53(3)u14+Λ54(3)u15+Λ55(3)u2+Λ56(3)u11.

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.
Fig. 5. Three neighboring hexagon unit cells and a possible ordering of the edges.

U=[u1;u2;]=[Ud;Un;Uf],

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:

AU=[A11A12A21A22A23A32A33][UdUnUf]=0,
(7)

[A22A23A32A33][D1D2]=[A210],
(8)

then Un and Uf are given by

Un=D1Ud,Uf=D2Ud.
(9)

Therefore, the system (7) is reduced to

BUd=0forB=A11+A12D1.
(10)

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

s1(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

ω2(1)=ω1ω1ω0s1(ω1)s1(ω0)s1(ω1).

To obtain the second approximation, we replace s 1(ω 1) by -s 1(ω 1). Thus

ω2(2)=ω1ω1ω0s1(ω1)+s1(ω0)s1(ω1).

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

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

𝓔N=ω(N)ω(16)ω(16).
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

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]

] by the plane wave expansion method, but they only used a small 3×3 parallelogram supercell as shown in Fig. 2. Later, Sakoda [17

17. K. Sakoda, “Numerical study on localized defect modes in two-dimensional triangular photonic crystals,” Journal of Applied Physics , 84, 1210–1214 (1998). [CrossRef]

] 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

17. K. Sakoda, “Numerical study on localized defect modes in two-dimensional triangular photonic crystals,” Journal of Applied Physics , 84, 1210–1214 (1998). [CrossRef]

]. The analysis in [33

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]

] 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

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]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(050.5298) Diffraction and gratings : Photonic crystals

ToC Category:
Photonic Crystals

History
Original Manuscript: September 4, 2007
Revised Manuscript: October 10, 2007
Manuscript Accepted: October 17, 2007
Published: October 18, 2007

Citation
Shaojie Li and Ya Yan Lu, "Computing Photonic Crystal Defect Modes by Dirichlet-to-Neumann Maps," Opt. Express 15, 14454-14466 (2007)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-22-14454


Sort:  Year  |  Journal  |  Reset  

References

  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," Jpn. J. Appl. Phys. 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 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). [CrossRef]
  15. V. F. Rodr´ıguez-Esquerre, M. Koshiba and H. E. Hern´andez-Figueroa, "Finite-element analysis of photonic crystal cavities: Time and frequency domains," J. Lightwave 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," J. Appl. Phys. 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 twodimensional lattices: a high Q-resonant cavity," Physica B 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. Lightwave 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," J. Comput. Math. 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, 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 J. Matrix Anal. Appl. 12, 597-609 (1991). [CrossRef]
  30. T. Lu and D. Yevick, "A vectorial boundary element method analysis of integrated optical waveguides," J. Lightwave 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited