Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals

A finite-difference frequency-domain method based on the Yee’s cell is utilized to analyze the band diagrams of two-dimensional photonic crystals with square or triangular lattice. The differential operator is replaced by the compact scheme and the index average scheme is introduced to deal with the curved dielectric interfaces in the unit cell. For the triangular lattice, the hexagonal unit cell is converted into a rectangular one for easier mesh generation. The band diagrams for both square and triangular lattices are obtained and the numerical convergence of computed eigen frequencies is examined and compared with other methods. © 2004 Optical Society of America OCIS codes: (230.3990) Microstructure devices; (260.2110) Electromagnetic theory; (999.9999) Photonic crystals. References and links 1. E. Yablonovitch, “Inhibited Spontaneous Emission in Solid-State Physics and Electronics,” Phys. Rev. Lett. 58, 2059–2062 (1987). 2. S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 2486– 2489 (1987). 3. K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existance of photonic gap in periodic dielectric structures,” Phys. Rev. Lett.65, 3152–3155 (1990). 4. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic crystal: modeling the flow of light(Princeton University Press., Princeton, NJ, 1995) 5. R. D. Meade, K. D. Brommer, A. M. Rappe, and J. D. Joannopoulous, “Existance of a photonic band gap in two dimensions,” Appl. Phys. Lett. 61, 495–497 (1992). 6. M. Plihal and A. A. Maradudin,“Photonic band structure of two-dimensional system: The triangular lattice,” Phys. Rev. B44, 8565–8571 (1991). 7. C. T. Chan, Q. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. 51, 16635–16642 (1995). 8. M. Qiu, and S. He, “A nonorthogonal finite-difference time-domain method for computing the band structure of a two-dimaensional photonic crystal with dielectric and metalic inclusions,” Appl. Phys. 87, 8268–8275 (1992). 9. M. J. Robertson, S. Ritchie, and P. Dayan, “Semiconductor waveguides: analysis of optical propagation in single rib structures and directional couplers,” Inst. Elec. Eng. Proc.-J. 132, 336–342 (1985). 10. M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” Inst. Elec. Eng. Proc.-J. 135, 56–63 (1988). 11. C. Vassallo, “Improvement of finite difference methods for step-index optical waveguides,” Inst. Elec. Eng. Proc.J.139, 137–142 (1992). 12. K. Bierwirth, N. Schulz, and F. Arndt, “Finite-difference analysis of rectangular dielectric waveguides by a new finite difference method,” J. Lightwave Technol. 34, 1104–1113 (1986). #3853 $15.00 US Received 20 February 2004; revised 25 March 2004; accepted 28 March 2004 (C) 2004 OSA 5 April 2004 / Vol. 12, No. 7 / OPTICS EXPRESS 1397 13. P. L̈usse, P. Stuwe, J. Sch üle, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–494 (1994). 14. G. R. Hadley and R. E. Smith, “Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions,” J. Lightwave Technol. 13, 465–469 (1994). 15. C. L. Xu, W. P. Huang, M. S. Stern, and S. K. Chaudhuri, “Full-vectorial mode calculations by finite difference method,” Inst. Elec. Eng. Proc.-J. 141, 281–286 (1994). 16. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite-difference vector beam propagation method. Analysis ans Assessment,” J. Lightwave Technol. 10, 295–305 (1992). 17. C. L. Xu W. P. Huang, and S. K. Chaudhuri, “Efficient and accurate vector mode calculations by beam propagation method,” J. Lightwave Technol. 11, 1209–1215 (1993). 18. K. S. Yee, “Numerical solution of initial boundary value problems involing Maxwell’s equations on isotropic media,” IEEE Trans. Antenna Propagat. 14, 302–307 (1966). 19. T. Ando, H. Nakayama, S. Numata, J. Yamauchi, and H. Nakano, “Eigenmode analysis of optical waveguides by a Yee-mesh-based imaginary-distance propagation method for an arbitrary dielectric interface,” J. Lightwave Technol.20, 1627–1634 (2002). 20. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10, 853–864 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853 21. H. Y. D. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE. Trans. Microwave Theory Tech. 34, 2688–2695 (1996). 22. L. Shen, S. He, and A. Xiao, ‘A finite-difference eigenvalue algorithm for calculating the band structure of a photonic crystal,” Comp. Phys. Commun. 143, 213–221 (2002). 23. A. Yefet and E. Turkel, “Fourth order compact implicit method for the Maxwell equations with discontinuous coefficients,” Appl. Numer. Math. 33, 125–134 (2000). 24. Y. P. Chiou, Y. C. Chiang, and H. C. Chang, “Improved three-point formulas considering the interface conditions in the finite-difference analysis of step-index optical devices,” J. Lightwave Technol. 18, 243–251 (2000). 25. Y. C. Chiang, Y. P. Chiou, and H. C. Chang, “Improved full-vectorial finite-difference mode solver for optical waveguides with step-index profiles,” J. Lightwave Technol. 20, 1609–1618 (2002). 26. 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), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173 27. R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulous, and O. L. Alerhand, “Accurate theoretical analysis of photonic band-gap materials,” Phys. Rev. B 48, 8434–8437 (1993).


Introduction
In recent years, photonic crystals (PCs) have inspired a lot of interest and many research efforts have been devoted to their possible applications.They are characterized by their photonic band gaps (PBGs) which result from structure periodicity [1][2][3][4][5][6][7][8].Two-dimensional (2-D) PCs, such as those composed of either dielectric rods or air columns, are widely employed in many applications such as new means of light waveguiding which offer a very high degree of control over the propagation of light.Even through sharp bends, a wide transmission band can be achieved by carefully designing the geometry of the waveguide.A variety of numerical methods have been utilized to calculate the band structures of 2-D PCs.Among these the most used are the planewave expansion (PWE) method [3][4][5][6] and the finite-difference time-domain (FDTD) method [7,8].
The finite difference method (FDM) has been widely adopted to deal with differential equations in numerical analysis.Due to the simplicity of its implementation and the sparsity of its resultant matrix, mode solvers based on the FDM have become attractive numerical methods for analyzing the propagation characteristics of the optical or dielectric waveguides, especially for those complex dielectric waveguides without analytical solutions.The FDM was first employed to solve the scalar waveguide modes under the weakly guiding approximation [9].For strongly guiding structures, semi-vectorial equations for optical waveguides with arbitrary index profiles were derived [10,11].To obtain more accurate mode fields, a full-vectorial finite difference scheme was then proposed [12][13][14].Each mesh point is surrounded by four subregions whose permittivities are allowed to be different.The sparse matrix is formed by dis-cretizing the Helmholtz equations and matching the boundary conditions for all the subregions.Another way is to discretize the wave equations derived from Maxwell's equations with the finite difference scheme [15][16][17], which is usually used in the beam propagation method (BPM) [16,17].Recently, full-vectorial finite-difference waveguide mode solvers based on the Yee's cell [18] were proposed.Ando et al. [19] derived a Yee-mesh-based imaginary-distance propagation method with modified FD formulae which takes the boundary conditions into account.Fast convergence was achieved for the analysis of step-index fibers and rib waveguides.Zhu and Brown [20] proposed a finite-difference frequency-domain (FDFD) method which directly discretizes Maxwell's equations using the central difference scheme based on a Yee's cell and the index average technique was adopted to improve numerical convergence and accuracy.
Besides the PWE method and the FDTD method, finite-difference eigenvalue algorithms for analyzing the band structures of 2-D PCs have been proposed.Yang [21] proposed an algorithm by discretizing the Helmholtz equation in terms of the magnetic field H z or the electric field E z for the transverse-electric (TE) and transverse-magnetic (TM) modes, using the central difference scheme in a similar way as that of [12] for waveguide analysis.Only square-lattice PCs were discussed and no curved dielectric interfaces were considered in the unit cell of the PC.Shen et al. [22] derived another algorithm from wave equations in terms of the parallel and transverse components of the magnetic field for the TE and TM modes, respectively, and employed an effective-medium technique for smoothing the field distribution across the dielectric interface.Again, only square PCs were considered.Fast convergent results for 2-D PCs involving dielectric rods were demonstrated.
In this paper, an FDFD eigenvalue algorithm based on the Yee's cell is derived in the similar way as in [20], and is utilized to analyze the band structures of 2-D PCs with either square or triangular lattice.Instead of discretizing the Helmholtz equations or the wave equations, we treat Maxwell's equations directly.To increase the accuracy of our FDFD algorithm, a formula based on a fourth-order accurate compact scheme [23] is adopted, and the index average scheme is employed to deal with the curved interface in the computing domain.The derivation of the formulation of this FDFD eigenvalue algorithm is described in Section 2 in detail, and some numerical results are presented and compared with other numerical methods in Section 3. We also discuss the convergence behaviors of our method and other FD methods in Section 4 which is followed by the conclusion.

The FDFD algorithm with compact scheme
The cross-section of a 2-D PC is shown in Figs.1(a) and (b) for the square lattice and the triangular lattice cases, respectively, with the respective unit cell illustrated by dashed lines, where a is the lattice distance and r is the radius of the circles which can be either dielectric rods or air columns.Since the 2-D PC is uniform along the z direction and periodic in the transverse plane, in the following discussion we only consider the in-plane propagation that has zero propagation constant in the z direction, β z = 0, so that the wave modes in the PC are either TE or TM to z mode.Assuming all the fields have the exp( jωt) dependence with ω being the angular frequency, Maxwell's curl equations can be expressed as where µ and ε are the relative permeability and permittivity of the medium, respectively, and µ 0 and ε 0 are the permeability and permittivity of free space.

The TE mode
For the TE mode, we have only H z , E x , and E y components and the above curl equations can be written as  The transverse plane is then discretized using Yee's mesh for the TE mode as shown in Fig. 2(a), and we can obtain where µ z , ε x , and ε y are the relative permeability or permittivity at the points of H z , E x , and E y .
Instead of using the central difference scheme, we introduce a fourth-order compact scheme [23] to replace the differential operator in the above equations to improve the accuracy.For example, ∂ H z ∂ x at points (i+1, j + 1 2 ), (i, j + 1 2 ), and (i−1, j + 1 2 ) can be expanded using the Taylor series expansion with respect to point (i where x is the grid size in the x direction and H.O.T. means higher order terms.Also, H z at point (i − 1 2 , j + 1 2 ) can be expressed as If we neglect the higher order terms (H.O.T.) in Eqs. ( 5)-( 8), we can obtain the following relationship between the ∂ H z /∂ x and H z values Thus, the term With Eqs. ( 9)-( 12), within the computing domain we can obtain the matrix equations where H z , E x , and E y are the vectors composed of the field components H z , E z , and E y , respectively, at the grid points, and where µ z , ε x , and ε y are diagonal matrices representing µ z , ε x , and ε y at the corresponding grid points.After some mathematical work, we can obtain an eigenvalue equation in terms of H z as where k 0 = ω(µ 0 ε 0 ) which, following the similar derivation using Yee's cell and the compact scheme, can be converted into the following matrix equation where the vectors E z , H x , and H y and the diagonal matrices ε z , µ x , and µ y have the similar definitions as in Eq. ( 13).The eigenvalue equation in terms of E z for the TM mode can similarly be obtained as Note that Eqs. ( 15) and ( 18) apply to both dielectric and magnetic materials.Therefore, the proposed algorithm is also suitable for calculating the band structures for 2-D PCs composed of magnetic materials.

Periodic boundary conditions for the 2-D PC
Due to its periodic geometry, in the analysis of the PCs, we only need to consider the unit cell along with the periodic boundary conditions (PBCs).For the PC with square lattice shown in Fig. 1(a), the PBCs can be expressed as where Ψ is either the electric or magnetic field in the unit cell, and k x and k y are the wavenumbers in the x and y directions, respectively.For the PC with triangular lattice shown in Fig. 1(b), we have the PBCs Please note that Eqs.(20a) and (20b) are sufficient for describing the BCs in an infinitely periodic lattice and Eq.(20c) can be derived from Eqs. (20a) and (20b).PBC1, PBC2, and PBC3 are indicated at the corresponding sides of the unit cell as shown in Fig. 3(a).Since the unit cell for the triangular lattice is hexagonal, it is not convenient to be fitted by Yee's rectangular mesh.To overcome this problem, an FDTD scheme in a nonorthogonal coordinate system was proposed in [8].By using the nonorthogonal mesh, the non-rectangular unit cell can be matched quite well but a coordinate transformation is needed and more complicated formulae need to be derived.In our method, we simply move the upperright triangle (the red one) to the lowerleft and the lowerright triangle (the blue one) to the upperleft, as shown in Fig. 3(b), generating a modified unit cell.The modified unit cell still obeys the same PBCs at its boundaries as indicated in Fig. 3(b), e.g., the upper portion of the left side and the lower portion of the right side, both labled as PBC1, are related by the PBC of Eq. (20a).Since it occupies a rectangular area, it is easier to be fitted by the rectangular mesh.

The average index scheme
As for the curved dielectric interface in the computing domain, the conventional treatment in the FDM is to use either the gradded-index [15][16][17] or the staircase index [12][13][14] approximation, and much finer mesh grids are needed to approach the original structure.Improved formulations for oblique or curved interfaces were derived in [24,25] using the Taylor series expansion and matching the interface conditions.Here, for simplicity, an averaging scheme which has been shown to be useful and efficient in increasing the analysis accuracy [20] is adopted to determine the values of µ and ε at each grid point.For example, in Fig. 2 where ε 1 and ε 2 are the relative permittivities of media 1 and 2, respectively, in the mesh and f is the filling fraction of medium 1 in the mesh.

Numerical results
We first consider the square-arranged 2-D PC, shown as the inset in the middle of Fig. 4(a), formed by parallel alumina rods with relative permittivity ε = 8.9 and radius r = 0.2a in the air, where a is the lattice distance.Each point along the boundary of the first Brillouin zone shown as the right inset in Fig. 4(a) determines k x and k y in the PBCs in Eq. ( 19), and thus the matrices in Eqs. ( 15) and (18).The eigen frequencies of this PC can then be obtained from solving the eigen values of Eqs. ( 15) and (18).The band diagrams of the TE and TM modes are plotted in Fig. 4(a) and (b), respectively.The circular dots are the results obtained using our compact FDFD algorithm with the index average scheme and the solid lines are obtained using the MIT Photonic-Bands (MPB) package [26] based on the PWE method with 128 × 128 resolution.It can been seen that the circular dots match quite well with the solid lines for both polarizations.For Fig. 4(a) and (b), 40 grid points are utilized in each lattice distance and it takes just a few minutes to obtain these two figures on Pentium IV 2.0 GHz personal computers.
To examine the efficiency of the index average (IA) scheme, we fix the wavenumber vector k at the M point in the first Brillouin zone, i.e., k = (π/a, π/a, 0), and search the eigen frequencies the TM modes as k is at the M point.The lines with circles and rectangles are the results obtained by our FDFD method without and with the index average scheme, respectively.
using various numbers of grid points with and without the IA scheme.Here and in the following figures, the number of grid points refers to the number of total grid points in the computing domain.In Fig. 5(a) and (b), the lines with circles represent the first eigen frequencies for the TE and TM modes, respectively, which are obtained using our FDFD algorithm without the IA scheme.One can see that it requires quite a lot of grid points to achieve numerical convergence.On the contrary, the lines with squares represent the results using the IA scheme and show quite excellent convergence for both TE and TM modes.
We then calculate other eigen frequencies for the same wavenumber vector k using our compact FDFD algorithm and three other eigenvalue algorithms for 2-D PCs, which are derived respectively following the FDFD mode solver [20], the FDM proposed by Lüsse et al. [13], and that by Xu et al. [15], with different numbers of total grid points.The relative error here refers to the difference between the result and a reference result obtained with the grid points being 2500.Figure 6(a) and (b) shows the relative errors of the first and second bands for the TE mode versus the number of grid points.All results from the above-mentioned algorithms are calculated with the IA scheme.It shows that as the number of total grid points is larger than 500, all the four algorithms can achieve a relative error less than 1%.Compared with the other three algorithms, our method can provide better convergence behavior thanks to the aid of the compact scheme.As for the TM mode, the relative errors of the first and second bands

Number of grid points
Relative error (%) TE first band this method method of [15] method of [13] method of [20]  TE second band this method method of [15] method of [13] method of [20] - TM first band this method method of [15] method of [13] method of [20]  TM second band this method method of [15] method of [13] method of [20] -3.0 are presented in Fig. 6(c) and (d), respectively.One can again see the convergence property of our method is good.As for the 2-D PC with triangular lattice, the band diagrams for both polarizations can be obtained similarly with our method using the modified unit cell.Figure 7(a) and (b) shows the band diagrams for the TE and TM modes of the structure consisting of dielectric cylinders with r = 0.2a and relative permittivity ε = 11.4 in the air.The triangles are our results using the compact FDFD method with the modified unit cell and the IA scheme, and the solid lines are obtained by the MPB package based on the PWE method [26].The circles represent the results from the FDTD method based on non-orthogonal coordinate [8] which was proposed to deal with such 2-D PC with non-rectangular unit cell and curved dielectric interfaces.One can see that there exists a visible difference between our results and those from the FDTD method, but fairly good agreement is observed between ours and those by the PWE method.To examine the convergence characteristics, we consider the k vector corresponding to the K point in the first Brillouin zone shown as the right inset in Fig. 7(a), i.e., k = (2π/ √ 3a, 2π/3a, 0).The convergent behaviors for the first and second bands of the TE and TM modes are illustrated in Fig. 8(a)-(d).Compared with the other three methods as in Fig. 6, our method can still provide better numerical convergence.
We discuss the order of numerical accuracy in the following.Although the compact FD scheme in space is a fourth-order one, as mentioned in [23], an arithmetic average of permittiv-  this method method of [15] method of [13] method of [20] 0 500 1000 1500 2000 2500

Number of grid points
Relative error (%) TE second band this method method of [15] method of [13] method of [20] - TM first band this method method of [15] method of [13] method of [20]  TM second band this method method of [15] method of [13] method of [20] ity at the dielectric interface could reduce the accuracy to second order.In our case, the accuracy of the proposed method is better than second order in the results of Fig. 6(a), (c), and (d) and Fig. 8(a), (c), and (d), whereas the accuracy is a little less than second order in the results of Fig. 6(b) and 8(b).Meade et al. [27] discussed numerical convergence in the PWE calculation of the first band at the zone edge for the square lattice of dielectric rods under different IA schemes.They found that the same IA scheme as in this paper showed very good convergence for the TM band, but not for the TE band.When using an inverse IA scheme involving a weighted average of the inverses of the permittivities, i.e., 1/ε z (i, j) = f /ε 1 + (1 − f )/ε 2 , the convergence performance became reversed between the TM and TE bands.They thus proposed a tensor dielectric scheme which performed well for both bands.In our FDFD calculation, the IA scheme performs well for both first bands as seen in Figs. 6 and 8.We have tried the inverse IA scheme in our FDFD analysis of the TE first and second bands.It reveals that the convergence becomes very good for the TE second band but much worse for the TE first band.Apparently, the effects of the IA schemes on the PWE method and on the FDFD method are different, which is worth further investigation.

Conclusion
We have demonstrated a new eigenvalue algorithm for solving the band diagrams of the 2-D PCs, which is based on the finite-difference frequency-domain method associated with the fourth-order accurate compact scheme [23].The index average scheme is also utilized to deal with the curved dielectric interfaces inside the unit cell.For the non-rectangular unit cell in triangular-lattice PCs, a modified unit cell of rectangular shape is introduced for easier numerical treatment.The band diagrams for both TE and TM modes in 2-D PCs are successfully obtained and good agreement with the PWE method calculation has been achieved.Our method is found to have better convergence property compared with other FD and FDFD algorithms.

Fig. 1 .
Fig. 1.The cross-sectional view of a 2-D PC and its unit cell with a being the lattice distance and r being the radius of the circles for (a) the square lattice and (b) the triangular lattice.

Fig. 3 .
Fig. 3. (a) The unit cell of the PC with triangular lattice and its corresponding PBCs.(b) The modified unit cell.

Fig. 4 .Fig. 5 .
Fig. 4. Band diagrams for the 2-D PC formed by square-arranged alumina rods with r/a = 0.2 and ε = 8.9 in the air.(a) TE mode and (b) TM mode.

Fig. 6 .
Fig. 6.The convergence properties of our method for the 2-D square-lattice PC compared with other methods.(a) TE first band; (b) TE second band; (c) TM first band; (d) TM second band.

Fig. 7 .
Fig. 7. Band diagrams of (a) the TE mode and (b) the TM mode for the 2-D PC formed by triangular-arranged dielectric cylinders with r/a = 0.2 and ε = 11.4 in the air.Our results (triangles) are compared with the results from the FDTD method (circles) and the PWE method (solid lines).

Fig. 8 .
Fig. 8.The convergence properties of our method for the 2-D triangular-lattice PC compared with other methods.(a) TE first band; (b) TE second band; (c) TM first band; (d) TM second band.