## Abstract

In-plane light propagation in two-dimensional (2D) photonic crystals (PCs) has been investigated by using the finite element method (FEM) in frequency domain. Conventionally, the band structures of 2D PCs were calculated by either the plane-wave expansion method (PWEM) or the finite difference time domain method. Here, we solve the eigenvalue equations for the band structures of the 2D PCs using the adaptive FEM in real space. We have carefully examined the convergence of this approach for the desired accuracy and efficiency. The calculated results show some discrepancies when compared to the results calculated by the PWEM. This may be due to the accuracy of the PWEM limited by the discontinuous nature of the dielectric functions. After acquiring the whole information of the dispersion relations within the irreducible Brillouin zone of the 2D PCs, the in-plane photon density of states for both the transverse electric (TE) and transverse magnetic (TM) modes can be calculated, accurately. For the case, the width of the complete band gap predicted by the FEM is much smaller, only about 65 % of that calculated by the PWEM. Therefore, the discrepancy in the prediction of complete band gaps between these two methods can be quite large, although the difference in band structure calculations is only a few percent. These results are relevant to the spontaneous emission by an atom, or to dipole radiation in two-dimensional periodic structures.

©2007 Optical Society of America

## 1. Introduction

In the past two decades, photonic crystals (PCs) have attracted much attention [1, 2]. To control the optical properties of materials has become a key issue in material engineering. Photonic crystals are periodic dielectric materials characterized by photonic band gaps (PBGs). A PBG can prohibit the propagation of electromagnetic (EM) waves whose frequencies fall within the band gap region. These materials are expected to find many applications in optoelectronics and optical communications. It was proposed that the emission of electromagnetic radiation can be modified by the environment [3, 4]. Several environments such as metallic cavities [5], dielectric cavities [6], and superlattices [7–12] have been studied. The environmental effects have been described by the photon density of states (PDOS) which is related to the transition rate of Fermi golden rule. In principle, a complete band gap along all dimensions in space can be best realized in a three-dimensional (3D) system. However, the difficulty in fabricating such 3D crystals with band gaps in the optical regime prohibits the progression of many applications. On the other hand, two-dimensional (2D) PCs have been extensively studied because they provide the possibility to control the propagation of light yet remain comparatively easy to fabricate [13–15].

Many numerical methods have been developed and applied to the analysis and investigation of photonic crystals including the plane-wave expansion method (PWEM) [16–22], the finite-difference frequency-domain/finite-difference time-domain (FDTD) method [23–25], the multiple-scattering method [26, 27], and the finite-element frequency-domain/finite-element time-domain method [28–35], and others [36–42]. In spite of successful computations, there are several problems with Fourier-based methods. First of all, the dielectric function is discontinuous, so Fourier-type expansions converge slowly. Many precautions must be taken in order to ensure that the calculated spectra are correct. For instance, it was found that the discontinuous nature of the dielectric function severely limits the accuracy of the PWEM [22]. The other is the field localization due to the complexity of geometry. The difference relations in the FDTD method generate errors that are due to numerical dispersion, which limits the overall size of the structure. In addition, when the shape of the grid does not conform to the shape of a real material such as a circular rod, additional errors can be produced [43]. As experimental techniques to fabricate PCs are improved, a more rigorous method, which is free from these drawbacks, will be required. On the other hand, the finite-element method (FEM) has proved to be a flexible and efficient numerical tool with which to design various types of microwave components with inhomogeneous and complex structures. Because the geometry is defined by elements in the FEM, not by points as in the finite-difference method, there is no aliasing for the description of the geometry. This characteristic provides a flexible way to model highly inhomogeneous structures [44].

In this work, in-plane light propagation in a 2D photonic crystal is investigated by using an adaptive finite element method in frequency domain. The polarization characteristics including both the transverse electric (TE) and transverse magnetic (TM) modes are considered in our simulation model. The finite-element method is employed to discretize the dielectric function profile of the PCs and the in-plane band structures are calculated by solving eigenvalue equations with proper periodic boundary conditions following the Bloch theorem [45, 46]. We have carefully examined the convergence of this approach for the desired accuracy and efficiency. The calculated results show some discrepancies when compared to the results calculated by the PWEM [21]. Based on the finite-element analysis of the in-plane band structures of the 2D PCs in the irreducible Brillouin zone, the in-plane photon density of states of the 2D PCs for the TE and TM modes can be accurately calculated, respectively. The model is formulated in Sec. 2. The calculated results and discussion are presented in Sec. 3. The conclusions are given in Sec. 4.

## 2. Formulation

#### 2.1. In-plane band structure

The propagation of light in a photonic crystal can be studied by solving the Maxwell’s equations. For time-harmonic fields it is convenient to use a phasor notation. The Maxwell’s equations lead to the wave equations, or the master equations:

and

where *ε*(*r*) and *μ*(*r*) are the permittivity and permeability functions of the PCs, respectively, and *ω* is the angular eigen-frequency. In a 2D periodic system, the dielectric function is a periodic function of *x* and *y*. We assume that the materials are linear, homogeneous, isotropic, lossless, and nonmagnetic. We have

where *ε _{r}*(

*x*,

*y*) is the dielectric function profile, and

*ε*and

_{a}*ε*are the dielectric constants of the air and dielectric regions, respectively. The two master equations are reduced to two homogeneous Helmholtz’s equations for the air (dielectric) region

_{d}and

A two-dimensional photonic crystal is periodic in two directions (*x*,*y*) and homogeneous in the third (*z*). For light propagating in the *xy*-plane, we can separate the modes into two independent polarizations, TM and TE modes, and consider the band structure and photon density of states of each. The propagation properties of TM and TE modes can be characterized by the field components parallel to the rods, *E _{z}*(

*x*,

*y*) and

*H*(

_{z}*x*,

*y*), respectively. For the in-plane propagation (

*k*= 0), the Helmholtz’s equations, Eqs. (4) and (5), for the air (dielectric) region can be rearranged as:

_{z}and

As the system has discrete translational symmetry in the *xy*-plane, we have *ε*(**r**
_{∥}) = *ε*(**r**
_{∥} + **R**), where **r**
_{∥} is the in-plane position vector and **R** is the linear combination of the primitive lattice vectors. By applying Bloch theorem, we can focus our attention on the values of the in-plane wave vector, **k**
_{∥} = *k _{x}x*̂ +

*k*, that are in the first Brillouin zone. We can calculate the in-plane band structures of 2D PCs by solving Eqs. (6) and (7) with the following boundary conditions,

_{y}ŷand

By symmetry, we further restrict the value of *k _{x}* and

*k*in the irreducible Brillouin zone.

_{y}For illustrations, we consider two symmetric lattices, square and triangular ones, as shown in Figure 1. The cross-sectional views of a square lattice of dielectric columns, with a radius *r* = 0.38 *a* and dielectric constant *ε _{d}* = 9, and the corresponding reciprocal lattice are shown in Fig. 1(a) and Fig. 1(b), respectively, where

*a*is the lattice constant and Γ: (0,0),

*X*: (

*π*/

*a*(1,0)), and

*M*: [

*π*/

*a*(1, 1)) are the high-symmetry points at the corners of the irreducible Brillouin zone (1/8 of the first Brillouin zone, in shaded region). The cross-sectional views of a triangular array of air columns with a radius

*r*= 0.4297

*a*drilled in a dielectric substrate of dielectric constant

*ε*= 11.9, and the corresponding reciprocal lattice are shown in Fig. 1(c) and Fig. 1(d), respectively. The corresponding high-symmetry points at the corners of the irreducible Brillouin zone (1/12 of the first Brillouin zone, in shaded region) are Γ : (0,0),

_{d}*M*: (

*π*/

*a*(1, -1/ √3)), and

*K*:(

*π*/

*a*(4/3,0)).

For the band diagram of the square lattice, the dimensionless frequency (*ωa*/2*πc*) is calculated and plotted as a function of Bloch’s vector (**k**
_{∥}) as the boundary Γ*XM*Γ of the irreducible Brillouin zone. For that of the triangular lattice, the corresponding path is Γ*MK*Γ. Table 1 shows the reciprocal lattice paths and the corresponding ranges of *k*
_{∥} and phase changes, **k**
_{∥} ∙ **R**
_{1} and **k**
_{∥} ∙ **R**
_{2}, for both the square and triangular lattices. Note that the direct lattice basis vectors **R**
_{1} and **R**
_{2} are (*a*, 0) and (0, *a*), respectively, for the square lattice and those of the triangular lattice are (*a*,0) and (*a*/2, √3*a*/2).

#### 2.2. Finite element method

In any periodic structure only a single cell needs to be considered. The unit cells of the square and triangular lattices we chose are also shown in Fig. 1(a) and Fig. 1(c), respectively. With the information of the phase changes mentioned above, the 2D Helmholtz’s equations, Eqs. (6) and (7), and the corresponding equations of boundary conditions, Eqs. (8) and (9), form two sets of eigenvalue problems for the TM and TE polarizations, respectively. To solve the eigenvalue equations, we employ an eigensolver for partial differential equations (PDEs) based on an adaptive finite element method. The finite element method can be adapted to problems of great complexity and unusual geometry, it is an extremely powerful tool in the solution of important problems in heat transfer, fluid mechanics, and mechanical systems. In addition, the method has been demonstrated for the successful design and analysis of a wide variety of electromagnetic components. For the adaptive FEM, a 2D, 6-node triangular finite element with a 5th-order basis function providing continuous derivatives between elements is used and the subdomains are partitioned into triangles, or mesh elements. The adaptive mesh generation will identify the regions where high resolution is needed and produces an appropriate mesh, which may be nonuniform, as shown in Figure 2. Therefore, in the adaptive finite element analysis, the spatial domain will be discretized into smaller and smaller triangles for satisfying the Helmholtz’s equations and the imposed boundary conditions until the desired accuracy is achieved. One should note that the shapes of the meshes conform the geometry of the PCs very well, as one can see in Figure 2.

#### 2.3. In-plane photon density of states

The dispersion relations of not only the boundaries but the interior of the irreducible Brillouin zone are needed for the calculation of photon density of states. For the dispersion relation of the square lattice, the dimensionless frequency is calculated as a function of Bloch’s vector as in the interior of the triangle Γ*XM* of the irreducible Brillouin zone. For that of the triangular lattice, the corresponding *k*-points are distributed inside the triangle Γ*MK*. Then, the dispersion relations can be plotted as 3D (*k _{x}* -

*k*-

_{y}*ω*) diagrams and it will be a surface for each band. To perform the 2D PDOS calculation, we construct two equifrequecy lines

*ω*(

*k*,

_{x}*k*) =

_{y}*ω*and

*ω*(

*k*,

_{x}*k*) =

_{y}*ω*+ Δ

*ω*), where

*ω*) is an arbitrary value of the frequency and Δ

*ω*is an infinitesimal increment. The differential surface element in

*k*-space within the lines is

The area of the total phase space surface contributing to the frequency range (*ω*, *ω* + *dω*), by definition, is

Then, the number of states within the range (*ω*, *ω* + *dω*) is obtained from dividing the surface calculated in Eq. (11) by the surface corresponding to one mode, (2*π*)^{2}/*S*, where *S* is the surface of the system. By definition, we have *dN*(*ω*) ≡ *D*(*ω*)*dω*. Finally, we arrive at the expression for the PDOS as

## 3. Results and discussion

We have carefully examined the convergence of the FEM approach for the desired accuracy and efficiency. Figure 3 shows our convergence tests of the first 15 bands at the *M* point for the TE and TM modes of a triangular array. As one can see in the figure, the fractional errors (*δω*/*ω*) of all eigenvalues are smaller than 10^{-4} when the mesh number increases to 16,340. The accuracy and convergence of FEM have been demonstrated. The dispersion relations of 2D photonic crystals are basically controlled by the dielectric contrast, the lattice type, and the filling ratio. For illustrations, two cases have been calculated by this finite element based approach. One is an example in the textbook [13], a square array of dielectric columns with a radius *r* = 0.38 *a* and dielectric constant *ε _{d}* = 9. The band structures for the TE and TM modes of the square lattice are plotted in Fig. 4(a). The in-plane wave vectors go along the edges of the irreducible Brillouin zone, from Γ to

*X*,

*X*to

*M*, and

*M*to Γ, as shown in the inset of Fig. 4 (a). Clearly, the lattice of dielectric cylinders has a gap for the TM modes but not for the TE modes. Our calculated results are similar to those in Ref [13]. The other case is a triangular array of air columns with a radius

*r*= 0.4297

*a*drilled in a silicon substrate (

*ε*≈ 11.9) in Ref. [21]. Figure 4 (b) shows the calculated band structures for the TE and TM modes of the triangular lattice following the paths, from Γ to

_{d}*M*,

*M*to

*K*, and

*K*to Γ, as shown in the inset of Fig. 4 (b). A complete photonic band gap opens at

*ωa*/2

*πc*~ 0.4. The calculated results are quite similar to those in Ref [21].

However, some discrepancies have been found after we carefully compared the data from our approach with those calculated by the PWEM in Ref. [21], where 720 plane waves were used to produce the results that deviated from convergent results by less than 0.5 %. Figure 5 shows the comparisons of band structure of the TE and TM modes calculated by the FEM and the PWEM. The discrepancy is getting larger as the frequency is going to higher regime. To make it clear, the frequency shift of PWEM from FEM corresponding to some symmetric *k*-points (Γ,*M*,*K*) for both the TE and TM modes is plotted in Figure 6. As one can see in the figure, the frequency shift approximately increases as the frequency increases. For the case, the discrepancy in band structure calculations is up to about 4 %. This may be due to the accuracy of the PWEM limited by the discontinuous nature of the dielectric functions [22]. As the contrast of the dielectric constant is high, the step-like dielectric function is usually approximated by limited number of Fourier basis in the PWEM. Therefore, the convergency of the PWEM to realistic cases will be slow and sometimes not so good.

After acquiring the whole information of the dispersion relations, the photon density of states can be calculated. The eigenfrequencies for 861 *k*-points uniformly distributed in the irreducible Brillouin zone have been calculated for the two cases via the FEM approach. It means that, due to the symmetry consideration, we effectively discretize the first Brillouin zone into 6,241 and 9,852 grid points for the square array and the triangular array, respectively. As Wang *et al*. found that a grid number of 256,000 is enough for the 3D case [47], the grid numbers we employed are enough to ensure the convergency accuracy to be better than 1 % for the 2D PDOS calculations. The dispersion relations corresponding to the TE and TM modes of the square lattice are plotted in Figs. 7(a) and 7(b), respectively. Those corresponding to the TE and TM modes of the triangular lattice are plotted in Figs. 7(c) and 7(d), respectively. For clarity, we only show the lowest four dispersion surfaces for each case. The in-plane PDOS for the TE and TM modes of both the square lattice and the triangular one have been calculated and plotted, as shown in Figure 8, where a frequency increment of 0.001 is used. It took nearly 18 h for each of the PDOS calculations to be performed on a PC equipped with a AMD Athlon 64 3000+ CPU and 1 GB RAM. The results of in-plane PDOS calculated by the PWEM in Ref. [21] are also given in Figs. 8(c) and 8(d) for comparison. For the case, the complete band gap, i.e., the overlap of the gaps of TE and TM modes, predicted by our method is much smaller and shifts to lower frequency, as one can see in the inset of Fig. 8(d). The frequency range of the complete band gap calculated by the PWEM method is from 0.383 to 0.411 and that predicted by our method is from 0.382 to 0.400. The corresponding bandwidth of the former case is 0.028 and that of the latter case is 0.018. The width of the band gap predicted by the FEM is only about 65 % of that calculated by the PWEM. Therefore, the discrepancy in the prediction of complete band gaps between these two methods can be quite large, although the difference in band structure calculations is only a few percent.

## 4. Conclusion

We have investigated the in-plane light propagation in 2D photonic crystals by using the finite element method in frequency domain. Conventionally, the band structures of the 2D PCs were calculated by either the plane-wave expansion method or the finite difference time domain method. The polarization characteristics including both the TE and TM modes are considered in our simulation model. The adaptive finite-element method is employed to discretize the dielectric function profile of the PCs and the in-plane band structures are calculated by solving eigenvalue equations with proper periodic boundary conditions following the Bloch theorem. We have carefully examined the convergence of this approach for the desired accuracy and efficiency. The calculated results show some discrepancies when compared to the results calculated by the PWEM. This may be due to the accuracy of the PWEM limited by the discontinuous nature of the dielectric functions. The FEM can be easily adapted to solve problems of great complexity and unusual geometry. The eigenvalues can be accurately and efficiently calculated no matter how complex the geometric structures are. Based on the finite-element analysis of the in-plane dispersion relations of the 2D PCs in the irreducible Brillouin zone, the in-plane photon density of states for both the TE and TM modes can be calculated accurately. For the case, the width of the complete band gap predicted by the FEM is much smaller, only about 65 % of that calculated by the PWEM. Therefore, the discrepancy in the prediction of complete band gaps between these two methods can be quite large, although the difference in band structure calculations is only a few percent. These results are relevant to the spontaneous emission by an atom, or to dipole radiation in two-dimensional periodic structures.

## Acknowledgement

The authors would like to thank Prof. B. Y. Gu at Institute of Physics, CAS and Prof. C. T. Chan at Department of Physics, HKUST for the helpful comments, discussions, and encouragement. This work was partially supported by the National Science Council, Taiwan, R. O. C. under Grant No. NSC 95-2112-M- 030-002, National Center for Theoretical Sciences, and National Center for High-performance Computing, Taiwan, R. O. C. which provides the computing resources.

## References and links

**1. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic Crystals* (Princeton University Press, Princeton, New Jersey1995).

**2. **J. D. Joannopoulos, P. R. Villeneuve, and S. Fan, “Photonic crystals: putting a new twist on light,” Nature (London) **386**,143–149 (1997). [CrossRef]

**3. **E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev. **69**,681 (1946).

**4. **D. Kleppner, “Inhibited Spontaneous Emission,” Phys. Rev. Lett. **47**,233–236 (1981). [CrossRef]

**5. **A. O. Barut and J. P. Dowling, “Quantum electrodynamics based on self-energy: Spontaneous emission in cavities,” Phys. Rev. A **36**,649–654 (1987). [CrossRef] [PubMed]

**6. **H. Rigneault and S. Monneret, “Modal analysis of spontaneous emission in a planar microcavity,” Phys. Rev. A **54**,2356–2368 (1996). [CrossRef] [PubMed]

**7. **J. P. Dowling and C. M. Bowden, “Atomic emission rates in inhomogeneous media with applications to photonic band structures,” Phys. Rev. A **46**,612–622 (1992). [CrossRef] [PubMed]

**8. **T. Suzuki and P. K. L. Yu, “Emission power of an electric dipole in the photonic band structure of the fcc lattice,” J. Opt. Soc. Am. B **12**,570–582 (1995). [CrossRef]

**9. **A. Kamli, M. Babiker, A. Al-Hajry, and N. Enfati, “Dipole relaxation in dispersive photonic band-gap structures,” Phys. Rev. A **55**,1454–1461 (1997). [CrossRef]

**10. **A. S. Sánchez and P. Halevi, “Spontaneous emission in one-dimensional photonic crystals,” Phys. Rev. E **72**,056609 (2005). [CrossRef]

**11. **P. Halevi and A. S. Sánchez, “Spontaneous emission in a high-contrast one-dimensional photonic crystal,” Opt. Commun. **251**,109–114 (2005). [CrossRef]

**12. **M. C. Lin and R. F. Jao, “Quantitative analysis of photon density of states for a realistic superlattice with omnidirectional light propagation,” Phys. Rev. E **74**,046613 (2006). [CrossRef]

**13. **K. Sakoda, *Optical Properties of Photonic Crystals* (Springer, Berlin2001).

**14. **M. M. Sigalas, R. Biswas, K. M. Ho, and C. M. Soukoulis, “Theoretical investigation of off-plane propagation of electromagnetic waves in two-dimensional photonic crystals,” Phys. Rev. B **58**,6791–6794 (1998). [CrossRef]

**15. **Z. Y. Li and Y. Xia, “Omnidirectional absolute band gaps in two-dimensional photonic crystals,” Phys. Rev. B **64**,153108 (2001). [CrossRef]

**16. **K. Sakoda, “Optical transmittance of a two-dimensional triangular photonic lattice,” Phys. Rev. B **51**,4672–4675 (1995). [CrossRef]

**17. **K. Sakoda, “Transmittance and Bragg reflectivity of two-dimensional photonic lattices,” Phys. Rev. B **52**,8992–9002 (1995). [CrossRef]

**18. **R. Hillebrand, W. Hergert, and W. Harms, “Theoretical band gap studies of two-dimensional photonic crystals with varying column roundness,” Phys. stat. sol. (b) **217**,981–989 (2000). [CrossRef]

**19. **M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: The triangular lattice,” Phys. Rev. B **44**,8565–8571 (1991). [CrossRef]

**20. **O. J. F. Martin, C. Girard, D. R. Smith, and S. Schultz, “Generalized field propagator for arbitrary finite-size photonic band gap structures,” Phys. Rev. Lett. **82**,315–318 (1999). [CrossRef]

**21. **K. Busch and S. John, “Photonic band gap formation in certain self-organizing systems,” Phys. Rev. E **58**,3896–3908 (1998). [CrossRef]

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

**23. **A. J. Ward and J. B. Pendry, “Calculating photonic Green’s functions using a nonorthogonal finite-difference time-domain method,” Phys. Rev. B **58**,7252–7259 (1998). [CrossRef]

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

**25. **C. T. Chan, Q. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. B **51**,16635–16642 (1995). [CrossRef]

**26. **G. Tayeb and D. Maystre, “Rigorous theoretical study of finite-size two-dimensional photonic crystals doped by microcavities,” J. Opt. Soc. Am. A **14**,3323–3332 (1997). [CrossRef]

**27. **W. Zhang, C. T. Chan, and P. Sheng, “Multiple scattering theory and its application to photonic band gap systems consisting of coated spheres,” Opt. Express **8**,203–208 (2001). [CrossRef] [PubMed]

**28. **J. K. Hwang, S. B. Hyun, H. Y. Ryu, and Y. H. Lee, “Resonant modes of two-dimensional photonic bandgap cavities determined by the finite-element method and by use of the anisotropic perfectly matched layer boundary condition,” Opt. Soc. Am. B **15**,2316–2324 (1998). [CrossRef]

**29. **W. Axmann and P. Kuchment, “An efficient finite element method for computing spectra of photonic and acoustic band-gap materials,” J. Comput. Phys. **150**,468–481 (1999). [CrossRef]

**30. **M. Koshiba, Y. Tsuji, and M. Hikari, “Time-domain beam propagation method and its application to photonic crystal circuits,” J. Lightwave Tech. **18**,102–110 (2000). [CrossRef]

**31. **G. Pelosi, A. Cocchi, and A. Monorchio, “A hybrid FEM-based procedure for the scattering from photonic crystals illuminated by a Gaussian beam,” IEEE Trans. Antennas Propag. **48**,973–980 (2000). [CrossRef]

**32. **D. C. Dobson, J. Gopalakrishnan, and J. E. Pasciak, “An efficient method for band structure calculations in 3D photonic crystals,” J. Comput. Phys. **161**,668–679 (2000). [CrossRef]

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

**34. **B. P. Hiett, J. M. Generowicz, S. J. Cox, M. Molinari, D. H. Beckett, and K. S. Thomas, “Application of finite element methods to photonic crystal modeling,” IEE Proc. -Sci. Meas. Technal. **149**,293–296 (2002). [CrossRef]

**35. **W. J. Kim and J. D. O’Brien, “Optimization of a two-dimensional photonic crystal waveguide branch by simulated annealing and the finite-element method,” J. Opt. Soc. Am. B **21**,289–295 (2004). [CrossRef]

**36. **A. Figotin and Y. A. Godin, “The Computation of Spectra of Some 2D Photonic Crystals,” J. Comput. Phys. **136**,585–598, 1997. [CrossRef]

**37. **R. C. McPhedran, L. C. Botten, J. McOrist, A. A. Asatryan, C. M. de Sterke, and N. A. Nicorovici, “Density of states functions for photonic crystals,” Phys. Rev. E **69**,016609, 2004. [CrossRef]

**38. **E. Moreno, D. Erni, and C. Hafner, “Band structure computations of metallic photonic crystals with the multiple multipole method,” Phys. Rev. B **65**,155120, 2002. [CrossRef]

**39. **J. M. Elson and P. Tran, “Coupled-mode calculation with the R-matrix propagator for the dispersion of surface waves on a truncated photonic crystal,” Phys. Rev. B **54**,1711–1715, 1996. [CrossRef]

**40. **L. C. Botten, N. A. Nicorovici, R. C. McPhedran, C. Martijn de Sterke, and A. A. Asatryan, “Photonic band structure calculations using scattering matrices,” Phys. Rev. E **64**,046603 (2001). [CrossRef]

**41. **D. Hermann, M. Frank, K. Busch, and P. wölfle, “Photonic band structure computations,” Opt. Express **8**,167–172 (2001). [CrossRef] [PubMed]

**42. **J. B. Pendry and A. MacKinnon, “Calculation of photon dispersion relations,” Phys. Rev. Lett. **69**,2772–2775 (1992). [CrossRef] [PubMed]

**43. **A. Taflove, *Computational Electrodynamics: The Finite- Difference Time-Domain Method*, Artech, Boston, Mass. (1995).

**44. **J. Jin, *The Finite Element Method in Electromagnetics, 2nd ed*. (Wiley-Interscience, New York2002).

**45. **A. Yariv and P. Yeh, *Optical Waves in Crystals* (Wiley, New York1984).

**46. **C. Kittel, *Introduction to Solid State Physics* (Wiley, New York1976).

**47. **X. H. Wang, R. Wang, B. Y. Gu, and G. Z. Yang, “Decay distribution of spontaneous emission from an assembly of atoms in photonic crystals with psudogaps,” Phys. Rev. Lett. **88**,093902 (2002) [CrossRef] [PubMed]