A finite-difference frequency-domain (FDFD) method is applied for photonic band gap calculations. The Maxwell’s equations under generalized coordinates are solved for both orthogonal and non-orthogonal lattice geometries. Complete and accurate band gap information is obtained by using this FDFD approach. Numerical results for 2D TE/TM modes in square and triangular lattices are in excellent agreements with results from plane wave method (PWM). The accuracy, convergence and computation time of this method are also discussed.
©2004 Optical Society of America
Photonic band gap materials and devices have been under intense research for over a decade following the seminal papers [1–2]. There are several methods for band structure analysis, such as the plane wave method (PWM) [3–5] and the FDTD [6–9] method. The PWM is able to provide complete and accurate information. However, the algorithm complexity is O(N3) and the computation is heavy for large problems. The order-N method based on FDTD can effectively reduce computation. It solves the Maxwell’s equations within the unit cell in time-domain by applying an initial field that covers all the possible symmetries; the eigen-modes are identified as the spectral peaks from the Fourier transform of the time-variant fields. The drawback of this method is that the accuracy depends on the number of iterations in time. There is also a possibility of losing true eigen-mode if the corresponding peak is too small, or resolution is too low. Moreover, spurious modes may arise from spectral noise. The FDFD method has been proposed for optical waveguide analysis [10–12], which is accurate and stable. In this paper, we show that this technique can be applied in photonic band gap analysis and we note that an FDFD approach using Helmholtz equation has been shown in . First we describe the derivation of the FDFD algorithm under generalized coordinate system and then apply the algorithm on 2D photonic crystals with two different geometries. The accuracy, convergence, and computation time in the FDFD method are compared with those of PWM.
and the renormalized fields are:
where k0 is the wave vector in free space, Qi’s are the grid size along each direction. The and are respectively the effective relative permitivity and permeability constants which are 3×3 tensors under the generalized coordinate system:
Q0 is a constant introduced to be roughly equal to Qi′s; |u 1·u 2×u 3| is the volume of the unit cell, εri (µri) is the relative dielectric constant (the relative permeability constant) at the position where the electric field Êi (the magnetic field Ĥi ) is located. g is the metric tensor that can be obtained using the three unit vectors,
and the length in generalized coordinate can be calculated using:
where Ui and Vi are coefficient matrices formed according to the boundary conditions, and they are proportional to 1/Qi .
An eigen-value problem in frequency-domain is formed for either Ê or Ĥ by eliminating Ĥ or Ê in Eqs. (6–7). For a given wave vector k, all the referred components outside the unit cell boundary can be obtained using Bloch’s periodic boundary condition:
where Rl can be an arbitrary lattice vector, and here it is limited in the unit cell or supercell.
A 2D Yee’s mesh under generalized coordinate system is shown in Fig. 1 for both TE and TM modes. E and H are arranged along two basis vectors u1 and u2; u3 is coincident with the z direction. Since Q3 is infinite, U3 and V3 in the equation (6–7) are zero, and simple eigen-equations can be obtained. The lattice vector Rl in Eq. (8) is chosen to be aquq, and aq is the dimension of the unit cell or supercell along direction q.
For TM modes, the eigen-equation is shown as follows:
The fields E and H in the two dimensional grids are arranged row by row into column vectors. Subsequently the Bloch boundary conditions are applied to get the matrices U and V:
For TE modes, the eigen-equation is:
The U and V matrices for TE mode are similar to those for TM modes and can be obtained by doing the exchange U1𒇔V1 and U2𒇔V2 in the equations (10–13) for TM modes.
According to Eq. (9), only is involved in TM mode. εr is located at the same point as E3, so no averaging is needed for εr3. For TE mode, the εr is half a grid away from E1 and E2 and the averaging is needed for εr1 and εr2. The periodicity of εr is applied for those grids outside the boundary.
3. Numerical results
All matrices involved are sparse; hence we can apply sparse matrix techniques to save computation time and memory. We implemented the algorithm using MATLAB since it provides convenient tools for sparse matrix operation with minimal programming efforts. Here we show a few numerical examples using our FDFD method and compare against the PWM using the program similar to that in Ref. . The first example is a 2D square lattice with dielectric alumina rods in air: εa=8.9, εb=1.0 and radius of the rod R=0.20a (a is the lattice constant). The second example is a 2D triangular lattice with air holes in dielectric GaAs materials: εa=1.0, εb=13.0 and hole radius R=0.20a. For the square lattice u1=[1 0 0], u2=[0 1 0], and u3=[0 0 1] and the metric [g] is a 3×3 identity matrix. For the triangular lattice, u1=[1 0 0], u2=[0 1/2 √3/2], and u3=[0 0 1] and [g] is calculated by Eq. (4). We used sub-cell averaging technique to smooth the transition at the interface and Eq. (5) is used for the rod profile transformation. The details of discretization and transformation of the cylindrical rod are shown in Ref. .
The TE/TM band gap for 2D square lattice and triangular lattice of the above examples are shown in Fig. 2 and Fig. 3 respectively. The FDFD results are indicated by ‘o’ and PWM results are plotted as solid lines. The results from the two methods show excellent agreements.
In Table 1 we list the first five bands at k=0 for the 2D triangular lattice shown above using the two methods in order to compare their accuracy and computation time. The computation time is measured on a 2.4GHz mobile Celeron® notebook with 256MB memory. From the Table we see that FDFD can reach the same accuracy as PWM in a shorter time.
A convergence curve for the eigen-frequency of band 5 at k=0 is shown in Fig. 4 versus the number of grids used along each direction. The computation time is also presented in the figure. The eigen-values converge to the accurate value at a moderate mesh size, for example, a/80. The computation time is highly dependent on the memory available on the computer.
When the unit cell or supercell has symmetry properties, computation time could be saved by using part of the unit cell under proper boundary conditions .
Next, we show a defect mode analysis using FDFD for the 2D square lattice of alumina rods in air as in the first example. A 5×5 supercell is selected and 200 grids are used along each direction. In this case, only the defect frequency is of interest since the band gap information is already known. Therefore we only have to find a certain number of eigen-frequencies of interest and the computation time is effectively reduced. The eigen-frequency obtained by FDFD is 0.3930. The mode field is shown in Fig. 5. Both results agree well with those by PWM and FDTD [16–17].
In conclusion, we have presented a FDFD method for photonic band gap calculations. This method is able to provide complete and accurate information about the band structure of a photonic crystal. The results of 2D TE/TM modes for two different geometries are compared with those obtained using plane wave method, and excellent agreement is achieved. By using a generalized coordinate system, various lattice geometries can be analyzed in the same manner.
The research at Old Dominion University is supported by NASA Langley Research Center through NASA-University Photonics Education and Research Consortium (NUPERC). We thank the reviewer for bringing references [11–13] to our attention.
References and links
4. R. D. Meade and A. M. Rappe et al., “Accurate theoretical analysis of photonic band gap materials,” Phys. Rev. B 48, 8434–8437 (1993). [CrossRef]
6. M. Qiu and S. He, “A nonorthogonal finite-difference time-domain method for computing the band structure of a two-dimensional photonic crystal with dielectric and metallic inclusions,” J. Appl. Phys. 87, 8268–8275 (2000). [CrossRef]
7. C. T Chan, Y. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. B 51, 16635–16642 (1995). [CrossRef]
8. J. Arriaga, A. J. Ward, and J. B. Pendry, “Order-N photonic band structures for metals and other dispersive materials,” Phys. Rev. B 59, 1874–1877 (1999). [CrossRef]
9. A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. 43, 773–793 (1996). [CrossRef]
10. 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 [CrossRef] [PubMed]
11. K. Bierwith, N. Schulz, and F. Arndt, “Finite-difference analysis of rectangular dielectric waveguide structures,” IEEE Trans. Microwave Theory Tech. 34, 1104–1113 (1986). [CrossRef]
12. P. Lusse, P. Stuwe, J Schule, 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). [CrossRef]
13. H Y D Yang, “Finite-difference analysis of 2D photonic crystals,” IEEE Trans. Microwave Theory Tech. 44, 2688–2695 (1996). [CrossRef]
14. K.S Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302–307 (1966). [CrossRef]
15. A. J. Ward, “Order-N program documentation,” http://www.sst.ph.ic.ac.uk/photonics/ONYX/orderN.html
16. S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express 11, 167–175 (2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-2-167 [CrossRef] [PubMed]
17. S. Guo and S. Albin, “Numerical techniques for excitation and analysisof defect modes in photonic crystals,” Opt. Express 11, 1080–1089 (2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-9-1080 [CrossRef] [PubMed]
18. P R McIssac, “Symmetry induced modal characteristics of uniform waveguides-I: Summary of results,” IEEE Trans. MTT 23, 421–429 (1975). [CrossRef]