A Novel Numerical Method for the Analysis of 2 D Photonic Crystals : The Cell Method

Abstract: The Cell Method, a new efficient numerical method suitable for working with periodic structures having anisotropic inhomogeneous media with curved shapes, is proposed in order to calculate the band gap of 2D photonic crystals for in-plane propagation of TM and TE waves. Moreover some numerical comparisons with other numerical methods will be provided. c © 2002 Optical Society of America OCIS codes: (000) General (000.3860) Mathematical methods in physics


Introduction
In the study of the properties of photonic band gap materials is important to find if some bands of frequencies exist (photonic band gap) where the propagation of the electromagnetic waves is prohibited in all the directions.In order to find these frequencies a lot of methods have been proposed up to now.The Plane Wave Expansion Method (PWE) [1] works with inhomogeneous media having curved shapes and deals with large not sparse hermitian matrices.The Finite Elements Method (FEM) works with inhomogeneous media having curved shapes as well and a generalized eigenvalues problem with large sparse (quasi band diagonal) matrices has to be solved.The Finite Difference Method (FD) works with inhomogeneous media without curved shapes and deals with large sparse (quasi band diagonal) hermitian matrices [2] in non generalized eigenvalues problems.In this paper we will use a new method, The Cell Method (CM) [3][4] [7], that works with inhomogeneous media having curved shapes, like FEM, and deals with large sparse (quasi band diagonal) hermitian matrices in non generalized eigenvalues problems, like FD.

The Cell Method
The Cell Method is based on the followings tools: • Two kinds of oriented space elements (space cells): the space elements endowed with an inner orientation that we will call primal space cells (primal points P, lines L, surfaces S and volumes V) and the space elements endowed with an outer orientation that we will call dual space cells (dual points P, lines L, surfaces S and volumes Ṽ) (Fig. 1a).
• Two oriented cell complexes and a relationship of duality between them.A space complex, synonymous of three dimensional grid, is a structured collection of points, lines, surfaces and volumes (cells).The primal complex is made by inner oriented space cells, the dual complex is made by outer oriented space cells and the relationship of duality means that we need a biunivocal matching respectively between the primal points P and the dual volumes Ṽ, the primal lines L and the dual surfaces S, the primal surfaces S and the dual lines L and between the primal volumes V and the dual points P.  • The use of the global (integral) variables to represent physical variables and their physically coherent association with the oriented space cells (Fig. 1b) [3][4]: The electromagnetic laws, that are links between global variables, can be divided into two categories, according to the type of variables that are linked: • Topological equations (Field equations).These equations link global variables associated with space cells belonging to the same kind of complex (either the primal or the dual one).Since they do not concern metric concepts, they can be enforced on the cell complexes in a discrete not approximated form by using appropriate incidence matrices.If we denote with G, C, D and G = D T , C = C T , D = −G T the incidence matrices [8] related to the primal and dual cell complexes respectively, which are the discrete counterparts of the differential operators gradient, curl and divergence, the topological equations in the frequency domain can be expressed as follows: - where V, Φ, F, Ψ, I, Q c are scalar arrays.
• Constitutive relations.These equations link global variables associated with space cells belonging to different kind of cell complexes (i.e. the electric constitutive relation links the electric voltage V , associated with the primal lines L, and the electric flux Ψ, associated with the dual surfaces S).They concern metric concepts and material parameters and they can be enforced on the cell complexes only approximatively in a discrete form by using suitable constitutive matrices: If an orthogonality occurs (L⊥ S and L⊥S, as in the Fig. 1b), then the matrices M ε ,M µ can be built diagonal [8] and assure, at the same time, the convergence of the method.In general, both for orthogonal or non-orthogonal grids, it is possible to build the matrices M ε ,M µ by the Microcell Interpolation Scheme (MIS) [4].

Analysis of 2d periodic structures with the Cell Method
The 2D photonic crystals are structures which have a transverse periodicity (xy-plane) and a longitudinal uniformity (z-axes).In order to study their properties it is enough to deal with a unit cell, the atomic part of the periodic structure (Fig. 2), with appropriate periodic boundary conditions [1][2].In case of in-plane propagation, which is the propagation only in the transverse plane, it is possible to analyze separately the TMz modes from the TEz modes.

Two dimensional TMz case
In CM for the TMz case we have to work on very simple 3D cell complexes that are, in particular, 2D extrude cell complexes along the z-axes (Fig. 3).This fact allows us to work like in a pure 2D case (Fig. 4a) by projecting the 3D cell complexes on the xy-plane, where: • The 2D unit cell is discretized by a triangular primal cell complex and a barycentric dual cell complex (Fig. 4a).
• • The constitutive relations are: From the previous equations and since C = −G T in the 2D cases, we obtain the matricial equation: that is the Helmholtz equation in CM for the TMz case, where V is the array of the N (=number of nodes=number of primal points P) electric voltages.If we build the matrices M ν and M ε with the MIS, for triangular grids the G T M ν GV is a symmetric positive semidefinite band diagonal matrix and the M ε is a diagonal positive definite matrix.For a rectangular unit cell of sizes [a, b] (Fig. 2a) the appropriate boundary condition is [2]: where β x and β y are the phase constants in the x and y directions.Since some voltages on the boundary are dependent of other voltages because the periodic conditions (Fig. 2c), it is possible to reduce the dimension of the problem by means of the following relation: where V is a vector of the independent voltages, T(β x , β y ) is a sparse matrix whose entries T ij are 1 if V i = V j , either e −jβxa or e −jβyb or e −jβxa−jβ y b according to the boundary condition between V i and V j , 0 otherwise.Thus the Helmholtz equation with the periodic boundary conditions reduces to: where A = T H G T M ν GT is a hermitian semidefinite positive matrix and B = T H M ε T is a diagonal real positive definite matrix.This generalized eigenvalues problem can be transformed into a simple eigenvalues problem: where The matrix M T M is still sparse hermitian positive semidefinite and, by an appropriate renumbering of the nodes, the structure of this matrix becomes quasi band diagonal.Thus, in order to find its eigenvalues, can be used the same efficient method developed for the FD method in [2].

Two dimensional TEz case
In CM for the TEz case it is possible to do the same considerations as the TMz cases by swaping the role of the primal and dual cell complexes (Fig. 4b): • The 2D unit cell is discretized by a barycentric primal cell complex and a triangular dual cell complex.
• The constitutive relations are: From the previous equations and since GT = C in the 2D cases, we obtain the matricial equation: that is the Helmholtz equation in CM for the TEz modes, where F is the array of the N (=number of nodes=number of dual points P ) magnetic voltages.For a rectangular unit cell of sizes [a, b], by the boundary condition: a similar eigenvalue problem like (5) has to be solved: In order to value the reliability of CM, some tests have been performed.The first test was the calculus of the photonic band structure for T M and T E waves in a square lattice (lattice constant a) of dielectric rods ( r = 12, radius=0.2a).The results, compared with PWE [5], show the existence of a TM photonic band gap (Fig. 5a).In (Fig. 5b) the normalized frequencies of the first T M mode in M and the second T M mode in X, for the first test, calculated by CM and FEM, are compared for different number of grid points.The second test was the calculus of the photonic band structure for T M and T E waves in a triangular lattice (lattice constant a) of air columns (radius=0.48a) in a dielectric substrate ( r = 13).The results, in excellent agreement with [1] where a PWE has been used, exhibit a complete band gap (Fig. 6a).The third test was the calculus of the photonic band structure for T M and T E waves in a square lattice of anisotropic dielectric rods ( xx = yy = 23.04,zz = 38.44,Filling Fraction=0.4) (Fig. 6b).The results, in good agreement with [6] where a PWE has been used, exhibit a complete band gap between (0.219-0.254) 2πc/a.In conclusion, the good agreement of the CM results with the results of other methods shows the effectiveness of the present approach.Moreover, in the test 1, CM exhibits more accuracy than FEM for the same number of grid points.

Conclusion
In this paper an efficient and accurate method for the analysys of the band gap of 2D photonic crystals for in-plane propagation of TM and TE modes with anisotropic media has been proposed and some comparisons with the results of other methodologies have been provided.Further work is however necessary to extend the method to the out of plane propagation cases.

Fig. 1 .
Fig. 1. (a)Inner and outer orientations of the space cells in the 3D space (b)Primal and dual complexes and global variables associated with their primal and dual space cells respectively.

Fig. 2 .
Fig. 2. Examples of 2D photonic crystals and periodicity of the unit cells in the xyplane (a) for rectangular unit cells (b) for quadrangular unit cells.(c) Dependence of some voltages on the boundary of a rectangular unit cell by the periodic boundary condition.

•
The association global variables → space cells is: Electric voltage V → Primal point P, Magnetic flux Φ → Primal line L Magnetic voltage F → Dual line L, Electric flux Ψ → Dual surface S The topological equations are: Faraday-Neumann law: GV = −jωΦ, Maxwell-Ampère law: CF = jωΨ.

FEM
Second TM mode frequency in X FEM First TM mode frequency in M CM First TM mode frequency in M CM Second TM mode frequency in X

Fig. 5 .Fig. 6 .
Fig. 5. (a) Results of the test 1 (b) A comparison of some CM and FEM results in the test 1.