Photonic Band Gap Analysis Using Finite-Difference Frequency-Domain Method

This Article is brought to you for free and open access by the Electrical & Computer Engineering at ODU Digital Commons. It has been accepted for inclusion in Electrical & Computer Engineering Faculty Publications by an authorized administrator of ODU Digital Commons. For more information, please contact digitalcommons@odu.edu. Repository Citation Guo, Shangping; Wu, Feng; and Albin, Sacharia, "Photonic Band Gap Analysis Using Finite-Difference Frequency-Domain Method" (2004). Electrical & Computer Engineering Faculty Publications. 175. https://digitalcommons.odu.edu/ece_fac_pubs/175


Introduction
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][4][5] and the FDTD [6][7][8][9] method.The PWM is able to provide complete and accurate information.However, the algorithm complexity is O(N 3 ) 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 timedomain 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][11][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 [13].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.

Theory
We consider nonconductive materials under generalized coordinates denoted by three unit basis vectors u q (x,y,z) (q=1,2,3).The Maxwell's curl equations in complex form can be expressed as [9,15]: and the renormalized fields are: where k 0 is the wave vector in free space, Q i 's are the grid size along each direction.The εˆ and µ ˆ are respectively the effective relative permitivity and permeability constants which are 3x3 tensors under the generalized coordinate system: Q 0 is a constant introduced to be roughly equal to Qi's; 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 E ˆ (the magnetic field i H ˆ) 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:

I I
We use Yee's mesh [14] and finite difference to replace the derivates in Maxwell's curl equations [9,15] and formulate them in matrix form using the approach described in [10]: where U i and V i are coefficient matrices formed according to the boundary conditions, and they are proportional to An problem in frequency-domain is formed for either E ˆ or H ˆ by eliminating H ˆ or E ˆ 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 R l 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 u 1 and u 2 ; u 3 is coincident with the z direction.Since Q 3 is infinite, U 3 and V 3 in the equation (6-7) are zero, and simple eigenequations can be obtained.The lattice vector R l in Eq. ( 8) is chosen to be a q u q , and a q is the dimension of the unit cell or supercell along direction q.

Two-dimensional cases
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:  (10) where ( ) ( ) ) 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 U 1 V 1 and U 2 V 2 in the equations (10)(11)(12)(13) for TM modes.
According to Eq. ( 9), only 1 33 − ε is involved in TM mode.ε r is located at the same point as E 3 , so no averaging is needed for ε r3 .For TE mode, the ε r is half a grid away from E 1 and E 2 and the averaging is needed for ε r1 and ε r2 .The periodicity of ε r is applied for those grids outside the boundary.

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. [16].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 u 1 =[1 0 0], u 2 =[0 1 0], and u 3 =[0 0 1] and the metric [g] is a 3x3 identity matrix.For the triangular lattice, , and u 3 =[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. [16].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.The number of plane waves, 2: the number of grids along each direction.
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.

Fig. 1 .
Fig. 1.Yee's 2D mesh in general coordinates.The dotted components are at the boundaries.

Fig. 2 .Fig. 3 .
Fig. 2. The band structure for a 2D square lattice by FDFD (o) and PWM (-).441 plane waves are used for PWM and mesh resolution is a/80 for FDFD.Left: TM mode, Right: TE mode.

Table 1 .
Eigen-frequencies for the first five bands of TE wave (k=0) for a triangular lattice with air holes in dielectric materials.