Investigation of band structures for 2 D non-diagonal anisotropic photonic crystals using a finite element method based eigenvalue algorithm

A finite element method based eigenvalue algorithm is developed for the analysis of band structures of two-dimensional non-diagonal anisotropic photonic crystals under the in-plane wave propagation. The characteristics of band structures for the square and triangular lattices consisting of anisotropic materials are examined in detail and the intrinsic effect of anisotropy on the construction of band structures is investigated. We discover some interesting relationships of band structures for certain directions of the wave vector in the first Brillouin zone and present a theoretical explanation for this phenomenon. The complete band structures can be conveniently constructed by means of this concept. © 2007 Optical Society of America OCIS codes: (230.3990) Microstructure devices; (160.3710) Liquid crystals; (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. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, Princeton, NJ, 1995). 4. J. C. Knight, T. A. Birks, P. St. J. Russell, and D. M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett. 21, 1547–1549 (1996). 5. S. G. Johnson, S. Fan, P. R. Villeneuve, J. D. Joannopoulos, and L. A. Kolodziejski, “Guided modes in photonic crystal slabs,” Phys. Rev. B 60, 5751–5758 (1999). 6. P. R. Villeneuve, S. Fan, and J. D. Joannopoulos, “Microcavities in photonic crystals: Mode symmetry, tunability, and coupling efficiency,” Phys. Rev. B 54, 7837–7842 (1996). 7. H. A. Haus, Waves and Fields in Optoelectronics (Prentice-Hall, Inc., Englewood Cliffs, NJ, 1984). 8. I. H. H. Zabel and D. Stroud, “Photonic band structures of optically anisotropic periodic arrays,” Phys. Rev. B 48, 5004–5012 (1993). 9. Z. Y. Li, B. Y. Gu, and G. Z. Yang, “Large absolute band gap in 2D anisotropic photonic crystals,” Phys. Rev. Lett. 81, 2574–2577 (1998). 10. C. Y. Liu and L. W. Chen, “Tunable band gap in a photonic crystal modulated by a nematic liquid crystal,” Phys. Rev. B 72, 045133 (2005). 11. 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). #79272 $15.00 USD Received 22 Jan 2007; revised 14 Apr 2007; accepted 17 Apr 2007; published 19 Apr 2007 (C) 2007 OSA 30 Apr 2007 / Vol. 15, No. 9 / OPTICS EXPRESS 5416 12. L. Zhang, N. G. Alexopoulos, D. Sievenpiper, and E. Yablonovitch, “An efficient finite-element method for the analysis of photonic band-gap materials,” in 1999 IEEE MTT-S Dig. 4, 1703–1706 (1999). 13. C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of twodimensional photonic crystals,” Opt. Express 12, 1397–1408 (2004). 14. P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75, 026703 (2007). 15. G. Alagappan, X. W. Sun, P. Shum, M. B. Yu, and D. den Engelsen, “Symmetry properties of two-dimensional anisotropic photonic crystals,” J. Opt. Soc. Am. A 23, 2002–2013 (2006). 16. J. Jin, The Finite Element Method in Electromagnetics (John Wiley and Sons, Inc., New York, 2002). 17. M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). 18. P. Yeh and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons, Inc., New York, 1999).


Introduction
Photonic crystals (PCs) are structures in which the materials are periodically organized to affect the propagation of electromagnetic (EM) waves.With the photonic band gaps (PBGs), the frequency ranges in which the EM waves are forbidden to propagate, resulting from the structure periodicity, PCs have inspired much interest and a lot of efforts have been devoted to their possible applications [1,2].Among varieties of PCs, two-dimensional (2D) PCs, characterized by easy fabrication, have attracted special attention and have been widely employed in various applications such as PC fibers, waveguides, and photonic cavities [3][4][5][6].Most of these applications make great use of the property of PBGs to achieve an extremely high degree of control over the propagation of EM waves.Consequently, accurate analyses of the band structures, and thus the PBGs, are quite significant for the design of these devices.
In addition to the spatial configuration of materials, material property itself plays a prominent role in the propagation of EM waves as well.The refractive index experienced by the EM wave propagating in isotropic materials is independent of the electric polarization, whereas that in anisotropic materials is strongly dependent on the electric polarization, which is further related to the wave propagation direction [7].Hence, by introducing anisotropic materials into PCs, the band structures change and some interesting phenomena have been observed and studied [8][9][10].However, the intrinsic effect on the concept of constructing the band structures caused by the introduction of anisotropic materials has not been discussed in any systematic way and is our main concern in this paper.
There are various kinds of numerical methods employed to calculate the band structures of 2D PCs.The plane-wave expansion (PWE) method [3] and the finite-difference time-domain (FDTD) method [11] are two usual schemes.In addition, the finite element method (FEM) [12] was utilized for this analysis likewise.Recently, the finite-difference frequency-domain (FDFD) method [13] and the pseudospectral method (PSM) [14] have also been demonstrated to provide alternative good approaches to obtain the band structures accurately.Most of these researches have focused on the structures made of either isotropic materials or anisotropic materials with diagonal permittivity and permeability tensors.But for many realistic applications, the principal coordinate system of the anisotropic materials may not always perfectly coincide with the coordinate system of the PCs [15].Consequently, the development of a more general numerical model which takes non-diagonal permittivity and permeability tensors into consideration becomes a highly valuable issue.Based on this kind of numerical model, it becomes possible to fully understand how the anisotropy influences the band structures and correctly take advantage of anisotropy on the utilization of PBGs.
In this paper, an FEM based eigenvalue algorithm for the analysis of band structures of 2D PCs formed by anisotropic materials is accomplished.The formulation is generalized to be able to deal with non-diagonal permittivity and permeability tensors in the cross-sectional plane of the PC for the in-plane wave propagation.For better accuracy and efficiency, the curved geom-etry in the structures is modeled by the curvilinear elements.The thorough derivation of this FEM based eigenvalue algorithm, including the adopted elements and the periodic boundary conditions (PBCs), is described in Section 2 in detail.Then the 2D anisotropic PCs with square and triangular lattices, composed of parallel silicon rods surrounded by liquid crystals (LCs) [10], are analyzed in Section 3. From these numerical results, we examine the relationships of band structures obtained from all possible directions of wave vector in the first Brillouin zone (BZ) and identify a guide for constructing the complete band structures of 2D anisotropic PCs.From these analyses we learn that the construction of band structures for anisotropic PCs have great differences from that for isotropic cases due to the intrinsic anisotropy which break the structure symmetry condition.The conclusion is summarized in Section 4.

Formulation
Under the source-free condition and with a time (t) dependence of the form exp( jωt) being implied, Maxwell's curl equations can be expressed as where ω is the angular frequency, μ 0 and ε 0 are the permeability and permittivity of free space, and [μ r ] and [ε r ] are, respectively, the relative permeability and permittivity tensors of the medium given by For a 2D PC which is uniform along the z direction and periodic in the x-y plane, the wave modes in the PC for the in-plane propagation can be decoupled into transverse-electric (TE) and transverse-magnetic (TM) to z modes if the tensor elements μ mn and ε mn ((m, n) = (x, z), (y, z), (z, x), and (z, y)) are set to be zero.

The TE and TM modes
For the TE mode which is composed of H z , E x , and E y components, the above curl equations can be reduced to From Eqs. ( 6) and ( 7) we can express E x and E y in terms of H z .Substituting these expressions into Eq.( 5), we can obtain the governing equation for the TE mode as where k 0 = ω √ μ 0 ε 0 is the wave vector in free space.
A similar procedure can be applied to the TM mode.For the TM mode comprising E z , H x , and H y components, Maxwell's curl equations can be simplified to From Eqs. ( 9) and (10) we can express H x and H y in terms of E z .Substituting these expressions into Eq.( 11), again we achieve the governing equation for the TM mode as For clarity we can express Eqs.( 8 [q] = ⎡ ⎣ q xx q xy q xz q yx q yy q yz q zx q zy q zz [q] = ⎡ ⎣ q xx q xy q xz q yx q yy q yz q zx q zy q zz for the TM mode (φ z = E z ).Note that Eq. (13) applies to both magnetic and dielectric materials with non-diagonal permeability and permittivity tensors provided that the TE and TM modes are decoupled.

Finite element discretization
For the scalar EM fields, H z or E z in Eq. ( 13), the quadratic triangular element [16,17] with six variables φ z1 -φ z6 , one at each of its three vertices and the other three at the middles of its three sides, as shown in Fig. 1(a), is extremely suitable for the discretization of the FEM.Moreover, in order to improve the accuracy and efficiency in the analysis, the curvilinear counterpart of Fig. 1(a), depicted in Fig. 1(b), is preferred and is adopted in our FEM based eigenvalue algorithm.
The six shape functions of the quadratic triangular element are given by provided that the nodes are numbered as shown in Fig. 1(a).In Eq. ( 18), L i (i = 1, 2, 3) are the simplex coordinates defined as with where x i and y i are the Cartesian coordinates of the ith vertex and Δ represents the area of the triangular element, that is, Here the index i in Eqs. ( 20)-( 23) assumes values 1, 2, 3, cyclically, so that if i = 3, then i + 1 = 1.
In the curvilinear element, the Cartesian coordinates, x and y, can be approximated as where x j and y j are the Cartesian coordinates of the jth node ( j = 1-6) within each element and N i is defined in Eq. ( 18).Noting that L 1 + L 2 + L 3 = 1 and selecting L 1 and L 2 as the independent variables, the relation of differentiation between the Cartesian coordinate system and the simplex coordinate system can be written as where [J] is the Jacobian matrix of the transformation given by Moreover, the integration of a function f (x, y) in the Cartesian coordinate system can be performed in the simplex coordinate system through e f (x, y)dxdy where |J| is the Jacobian, the determinant of [J].Hence, the required computation for element matrices associated with the governing equation can be obtained directly in the simplex coordinate system.Dividing the structure domain into a number of curvilinear triangular elements, the field φ z in each element can be expanded as φ z = {N} T {φ e z }, where {φ e z } is the variable vector for each element, and T denotes transpose.Applying Galerkin's method and assembling all element matrices, Eq. ( 13) can be transformed into the matrix form as where l denotes the contour or boundary enclosing the area, n is its outward normal vector, and ∑ e extends over all different elements.Note that Eq. ( 29) is not an eigenvalue matrix equation due to the existence of the vector {ψ} on the right-hand side.Fortunately, this term can be eliminated as a consequence of imposing the PBCs in our numerical model, as will be described in the next subsection.

Periodic boundary conditions for the 2D PCs
Due to the periodicity, PCs can be entirely described by the concept of unit cell along with the PBCs.Therefore, the correct setting of PBCs in the numerical model is extremely important.For a square unit cell, with the four sides labeled by I, II, III, and IV, respectively, as shown in Fig. 2(a), the PBCs can be expressed as where k x and k y are the wavenumbers in the x and y directions, respectively, and a is the lattice distance.To adequately incorporate these conditions into the numerical model, the whole boundary is grouped into three parts, as depicted in Fig. 2 29) are classified first according to the locations of the corresponding variables and Eq. ( 29) can be rewritten in the following form where {φ z } 0 stands for the vector of variables in the interior region, the subscripts, I, II, III, and IV, denote the sides at which the variables locate, and Through this procedure, the information at sides I and II is conveniently associated with that at sides III and IV, and the matrices and vectors in Eq. ( 29) are correctly modified to obtain the desired eigenvalue matrix equation of Eq. ( 44).Please note that, as every vertex belongs to a vertical and a horizontal side simultaneously, it is proper that the fields at the four vertices are linked together by the use of Eqs. ( 33)-(38), as has been automatically included in the above manipulation of Eq. ( 29).This idea is indicated by the arrows in Fig. 2(a).
As for a hexagonal unit cell, the PBCs can be written as PBCs.The fields at sides I and IV, II and V, and III and VI are, respectively, connected by Eqs. ( 45)-( 47), ( 48)-(50), and (51)-(53).Certainly, the six vertices should be considered separately because every vertex has two adjoining boundaries.Examining the relative positions of the six vertices carefully, the six vertices are partitioned into two independent parts, in each of which any vertex is at a distance a from the other two vertices.In each part, the fields at the three vertices can be linked together using Eqs.( 45)-( 53).The arrows in Fig. 2(b) illustrate this idea.Same as in the case of square unit cell, associating the information at sides I, II, and III with that at sides IV, V, and VI, respectively, will suitably modify the matrices and vectors in Eq. ( 29), and consequently the desired eigenvalue matrix equation can be derived.

Numerical results
For the analysis and discussion, we consider two 2D PCs formed by square-and trianglearranged parallel silicon rods surrounded by nematic LCs (5CB) [18].We choose the LCs as our anisotropic material for the convenient change of anisotropy by simply rotating the LC molecules.The components of the relative permittivity tensor [ε r ] of the nematic LCs we use are given as where n o = 1.5292 and n e = 1.7072 are, respectively, the ordinary and extraordinary refractive indices of the nematic LCs, θ c is the angle between the crystal c-axis and the z-axis, and φ c represents the angle between the projection of the crystal c-axis on the x-y plane and the x-axis, as defined in Fig. 3.For the decoupling of the TE and TM modes in the 2D condition.The angle, θ c , is set to be 90 • and thus the LC molecules will merely rotate in the x-y plane.
The band structures of these two kinds of lattices will be separately analyzed in the next two subsections.For the TM mode, the electric field only possesses z-component and will regard the anisotropic PC as an isotropic one.So we will focus our discussion on the band structures of the TE mode.From the numerical results we will explain the concept of constructing the band structures for anisotropic PCs clearly.

Square lattice
The first structure under consideration is the 2D PC with square lattice, as illustrated in Fig. 4(a), consisting of square-arranged parallel silicon rods with relative permittivity ε r = 11.56 and radius r = 0.2a in the background material of LCs.The unit cell of this structure is also indicated by the square region in Fig. 4(a), and the corresponding first BZ of the reciprocal lattice is depicted in Fig. 4(b).For isotropic PCs, a complete band structure can be conveniently constructed by simply considering all the possible directions of the wave vector k restricted in the irreducible BZ (IBZ) which is enclosed by points Γ, X, and M, as shown in Fig. 4(b).The reliability and completeness of the band structure obtained from this IBZ for isotropic PCs are fully based on the structure symmetry, including the rotation and reflection symmetry.For anisotropic PCs, the situation becomes more complicated because the circumstance experienced by the EM waves is highly direction-dependent.The anisotropy may break the structure symmetry, which in isotropic cases can be totally determined by the spatial configuration of materials.With this consideration, it is reasonable to inspect more band structures constructed from different directions of k in the first BZ.
To illustrate this phenomenon specifically, we take the square-arranged 2D PC of φ c = 30 • as an example.We compute the normalized frequencies for four distinct sub-zones in the first BZ marked by Γ-X-M, Γ-X -M, Γ-X -M , and Γ-X -M , as indicated in Fig. 4  comparisons for the four sub-zones are carefully examined in turn to figure out the relationships between them.In this figure, we use different colors to distinguish the four band structures obtained from the four distinct sub-zones defined before.In Fig. 5(a), the two band structures from Γ-X-M and Γ-X -M are dissimilar to each other except for the segment M-Γ.This can be easily understood because the segment M-Γ is shared by Γ-X-M and Γ-X -M.The consistency of results for segments Γ-X and M -Γ in Fig. 5(d) and (f), respectively, can be explained by this same reason.We then proceed to Fig. 5(c) where we can observe that the two band structures from Γ-X-M and Γ-X -M are exactly the same excluding the third part.The equivalence of the segments Γ-X and Γ-X can be visualized by the aid of Fig. 6(a).Although the direction of k 1 along segment Γ-X is opposite to that of k 2 along segment Γ-X , the rotation symmetry of the structure actually makes waves of these two wave vectors experience alike material property.Thus the band structure from these two segments are identical to each other.The reason for the equivalence of segments X-M and X -M is quite different because now k 1 along segment X-M does not point to the contrary direction of k 2 along segment X -M , and it cannot be attributed to the rotation symmetry.Besides, the appearance of LCs rotating at φ c = 30 • also prevents the reflection symmetry from a reasonable explanation.To explain what causes the identity of segments X-M and X -M , we appeal to the basic principle of periodic structures.For periodic structures, every lattice vector R and reciprocal lattice vector G must satisfy the requirement G • R = 2πN where N is an integer [3].Thus if k is incremented by G, the phase between cells is incremented by G • R, which is equal to 2πN and does not really affect the mode.For every pair of wave vectors, k 1 and k 2 , of segments X-M and X -M , the difference between them is exactly a reciprocal lattice vector G, as depicted in Fig. 6(b), so it is no surprise that the band structures of these two segments are the same.Next we examine segments X -M and X -M , another identical pair as shown in Fig. 5(d).The identity of this pair can be reasonably accepted by applying the two reasons just mentioned above simultaneously.This idea is illustrated in Fig. 6(c) which reveals that the difference between k 1 and k 2 is also exactly G.In Fig. 5(b) and (e), no conditions mentioned above occur, and, consequently, the whole band structures from these corresponding sub-zones are quite unlike.
When the LC molecules rotate to φ c = 45 • , another interesting behavior of the band structure can be observed.With LC molecules orienting to this special direction, the resultant structure, determined by the arrangement of silicon rods and the orientation of LC molecules, has reflection symmetry about the c-axis of LCs.According to this observation, it is absolutely reasonable to predict that the band structures from Γ-X-M and Γ-X -M are identical to those from Γ-X -M and Γ-X -M , respectively.The reflection symmetry of structures will vary with the rotation angle of LC molecules, but the primary idea discussed above works well and can be considered as a general principle applied to any anisotropic PCs.
We have demonstrated the possible relationships among the four distinct sub-zones in the first BZ clearly and comprehended that there really exists discrepancy among these sub-zones for anisotropic PCs.Consequently, to correctly construct the complete band structures of anisotropic PCs with square lattice, all the independent segments in these four distinct subzones should be considered collectively.We show the band structures of the TE mode for this anisotropic PC with square lattice of φ c = 0 • , 30 • , and 45 • in Fig. 7 from which we can see that there are no photonic band gaps for these different LC rotation angles.Notice that these band structures are constructed from the four sub-zones and the shadowed regions represent ignorable parts because they can be duplicated from unshadowed parts.

Triangular lattice
Fig. 8(a) shows the cross-section of the second structure to be analyzed, the 2D anisotropic PC with triangular lattice.Except for the lattice configuration, the structure parameters are the same as those of the square-arranged case.The unit cell for this 2D PC is a hexagon as indicated in Fig. 8(a) and the corresponding first BZ of the reciprocal lattice is also a hexagon, as plotted in Fig. 8(b).The analysis of band structures for this PC with triangular lattice is similar to that with square lattice discussed in the last subsection.For a hexagonal BZ, it is necessary to consider the six distinct sub-zones as marked in Fig. 8(b) for obtaining the complete band structures.The understanding about possible relationships among these six sub-zones can be conveniently accomplished in an analogous way as in the square BZ case.
We show the band structures of the TE mode for this PC of φ c = 0 • , 30 • , and 45 • in Fig. 9.To get an insight into the effect of the anisotropic material on the band gap, we concentrate on Fig. 9(a), the case of φ c = 0 • , first.The red region, if exists, in every one of the six distinct sub-zones stands for the band gap of the structure found from the corresponding sub-zone.For φ c = 0 • , Γ-M-K, Γ-M-K , and Γ-M -K are identical to Γ-M -K , Γ-M -K , and Γ-M -K , respectively.Accordingly, it is obvious that the band gaps corresponding to these respective pairs are exactly the same to each other.The main point worth emphasizing here is that the normalized frequencies between which the band gaps exist for Γ-M-K, Γ-M-K , and Γ-M -K are different.So to describe the band gap of this structure suitably, we should find the normalized frequency range appearing simultaneously in all distinct band gaps.Then this normalized frequent range will stand for the correct band gap of the structure.Specifically, the band gap found from Γ-M-K falls between the normalized frequencies of 0.636 and 0.656, whereas the

Fig. 2 .
Fig. 2. (a) A square unit cell and its corresponding PBCs.(b) A hexagonal unit cell and its corresponding PBCs.

Fig. 3 .
Fig. 3. Schematic definition of rotation angles for the LC molecule.

Fig. 4 .
Fig. 4. (a) The cross-section of a 2D PC with square lattice in the background of LCs.(b) The first BZ of (a).

Fig. 8 .
Fig. 8. (a) The cross-section of a 2D PC with triangular lattice in the background of LCs.(b) The first BZ of (a).
yx p yx p xy − p xx p yy Note that columns 2 and 4 as well as columns 3 and 5 of matrices [K] and [M] in Eq. (39) can be combined through Eqs.(33) and (36), respectively.Furthermore, if we divide row 2 in Eq. (39) by e − jk x a and add it to row 4, {ψ} I and {ψ} III will cancel each other by the use of Eqs.(34) and (35).Similarly, {ψ} II and {ψ} IV can cancel each other when row 3 in Eq. (39) is divided by e − jk y a and then added to row 5 with the help of Eqs.(37) and (38).By means of these operations, Eq. (39) becomes ⎡