Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy

Sen-ming Hsu; Hung-chun Chang

doi:10.1364/OE.15.015797

1. Introduction

Photonic crystals (PCs) are structures formed by periodically-organized materials that may create some frequency ranges, named the photonic band gaps (PBGs), in which the electromagnetic (EM) waves are forbidden to propagate. With this peculiar property, PCs are especially characterized by the capability of achieving an extremely high degree of control over the propagation of EM waves. This significant character makes PCs quite different from other structures and has drawn much attention for them. Therefore, since the pioneering works in [1,2], many efforts have been devoted to the potential applications of PCs, especially the two-dimensional (2D) PCs [3–6] which can be fabricated easily. This paper concerns in-plane wave propagations in 2D PCs.

The fundamental characteristics of PCs are highly related to the configurations and properties of materials of the unit cell. Various kinds of numerical methods have been employed to accomplish the analysis of band structures for 2D PCs. The plane-wave expansion (PWE) method [3] and the finite-difference time-domain (FDTD) method [7] are two familiar schemes. In addition, some famous numerical schemes, such as the finite element method (FEM) [8], the finite-difference frequency-domain (FDFD) method [9], and the pseudospectral method (PSM) [10], have been utilized for this analysis likewise. To obtain the complete band structures for 2D PCs from these numerical models, the wave modes are often decoupled into transverse-electric (TE) and transverse-magnetic (TM) modes, and then the eigen frequencies of these two sets of modes are separately solved. In other words, these numerical models are based on a scalar algorithm in which the scalar equations for the TE and TM modes are manipulated apart to construct the band structures for 2D PCs.

If the 2D PCs are composed of isotropic materials, the above-mentioned numerical models are perfectly sufficient to perform the analysis of band structures. However, sometimes the anisotropic materials are intentionally introduced into the 2D PC structures for the purpose of altering the patterns of the band structures, and thus controlling the behaviors of the PBGs [11–13]. For the 2D PCs made of diagonal anisotropic materials, these numerical models may still be capable of obtaining the complete band structures provided that the completeness of the irreducible Brillouin zone (IBZ) for constructing the band structures is carefully modified [14]. Nevertheless, in situations where the permittivity and permeability tensors of the anisotropic PCs are not in the diagonal form, the correct band structures can not be obtained using these numerical models through separating TE and TM modes. Recently some more general models have been developed to investigate 2D non-diagonal anisotropic PCs [14,15]. But in those researches, the anisotropy was still restricted for the guarantee that the wave modes can simply be decoupled into TE and TM modes and the scalar algorithm can still be applied.

For the analysis of band structures of 2D PCs with arbitrary 3D anisotropy, the traditional scalar algorithm is obviously insufficient, and a different algorithm beyond the scope of the scalar algorithm is absolutely necessary. For example, a finite element-integral equation (FEIE) method [16] has been formulated for treating planar PBG structures with finite thickness [17,18], which were essentially 3D problems. However, the matrix eigenvalue equation formulated was not a standard one with the coefficient matrix elements depending on the eigenvalue to be searched, and some iteration procedure was needed for finding the eigenvalues [17]. In this paper, a full-vectorial FEM based eigenvalue algorithm, which is capable of analyzing the band structures of 2D PCs with arbitrary permittivity and permeability tensors, is proposed. To deal with the most general situation of 2D problems with 3D material anisotropy, we simultaneously take all the components of the electric or magnetic fields into consideration, and successfully achieve a generalized full-vectorial matrix eigenvalue equation for the analysis of band structures. The matrix equation obtained is in the form of a standard matrix eigenvalue equation. The thorough derivation of this full-vectorial FEM based eigenvalue algorithm is described in Section 2 in detail. Then several 2D PCs with square and triangular lattices, including the isotropic and anisotropic cases, are analyzed in Section 3. From these numerical results and corresponding theoretical explanations, the performance of this algorithm is examined, and the correctness and usefulness of this algorithm are demonstrated. The conclusion is summarized in Section 4.

2. Formulation

2.1. The governing equation

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

\nabla \times E = - j ω μ_{0} [μ_{r}] H

\nabla \times H = j ω ε_{0} [ε_{r}] E

where ω is the angular frequency, µ ₀ and ε ₀ 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

[μ_{r}] = [\begin{matrix} μ_{xx} & μ_{xy} & μ_{xz} \\ μ_{yx} & μ_{yy} & μ_{yz} \\ μ_{zx} & μ_{zy} & μ_{zz} \end{matrix}]

[ε_{r}] = [\begin{matrix} ε_{xx} & ε_{xy} & ε_{xz} \\ ε_{yx} & ε_{yy} & ε_{yz} \\ ε_{zx} & ε_{zy} & ε_{zz} \end{matrix}] .

From Eqs. (1) and (2), we can derive the vectorial wave equation as

\nabla \times ([p] \nabla \times Φ) - k_{0}^{2} [q] Φ = 0

where $k_{0} = ω \sqrt{μ_{0} ε_{0}}$ is the wavenumber in free space, Φ is either the electric field E or the magnetic field H, and the tensors [p] and [q] are, respectively, given by

[p] = [\begin{matrix} p_{xx} & p_{xy} & p_{xz} \\ p_{yx} & p_{yy} & p_{yz} \\ p_{zx} & p_{zy} & p_{zz} \end{matrix}] = {[\begin{matrix} μ_{xx} & μ_{xy} & μ_{xz} \\ μ_{yx} & μ_{yy} & μ_{yz} \\ μ_{zx} & μ_{zy} & μ_{zz} \end{matrix}]}^{- 1}

[q] = [\begin{matrix} q_{xx} & q_{xy} & q_{xz} \\ q_{yx} & q_{yy} & q_{yz} \\ q_{zx} & q_{zy} & q_{zz} \end{matrix}] = [\begin{matrix} ε_{xx} & ε_{xy} & ε_{xz} \\ ε_{yx} & ε_{yy} & ε_{yz} \\ ε_{zx} & ε_{zy} & ε_{zz} \end{matrix}]

for Φ=E and

[p] = [\begin{matrix} p_{xx} & p_{xy} & p_{xz} \\ p_{yx} & p_{yy} & p_{yz} \\ p_{zx} & p_{zy} & p_{zz} \end{matrix}] = {[\begin{matrix} ε_{xx} & ε_{xy} & ε_{xz} \\ ε_{yx} & ε_{yy} & ε_{yz} \\ ε_{zx} & ε_{zy} & ε_{zz} \end{matrix}]}^{- 1}

[q] = [\begin{matrix} q_{xx} & q_{xy} & q_{xz} \\ q_{yx} & q_{yy} & q_{yz} \\ q_{zx} & q_{zy} & q_{zz} \end{matrix}] = [\begin{matrix} μ_{xx} & μ_{xy} & μ_{xz} \\ μ_{yx} & μ_{yy} & μ_{yz} \\ μ_{zx} & μ_{zy} & μ_{zz} \end{matrix}]

for Φ=H.

Fig. 1. The cross-section of a 2D PC and the corresponding unit cell for (a) the square lattice and (b) the triangular lattice.

Download Full Size | PDF

In this paper, we will consider 2D PCs with two kinds of spatial configuration, the square- and triangle-arranged lattices, as shown in Fig. 1(a) and (b), respectively, where a is the lattice distance and r is the radius of the parallel cylinders. The square region, enclosed by sides I, II, III, and VI, in Fig. 1(a) indicates the unit cell of the 2D PC with square lattice, while the hexagonal region, surrounded by sides I, II, …, VI, in Fig. 1(b) denotes the unit cell of the 2D PC with triangular lattice. For these 2D PCs which are uniform along the z direction and periodic in the x-y plane, the field distribution Φ of the wave modes in the PCs for the in-plane propagation can be expressed as

Φ (x, y, z) = Φ_{t} (x, y) + \hat{z} Φ_{z} (x, y)

where Φ_t and Φ_z, both assumed to be function of x and y, are the transverse and longitudinal field components of Φ, respectively. Substituting Eq. (10) into Eq. (5) and using

\nabla = \nabla_{t} + \nabla_{z}

where ∇_t and ∇_z are, respectively, the transverse and longitudinal parts of the ∇ operator, Eq. (5) can be separated into its transverse component

\nabla_{t} \times ([\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & p_{zz} \end{matrix}] \nabla_{t} \times Φ_{t} + [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ p_{zx} & p_{zy} & 0 \end{matrix}] \nabla_{t} \times \hat{z} Φ_{z})

and its longitudinal component

\nabla_{t} \times ([\begin{matrix} 0 & 0 & p_{xz} \\ 0 & 0 & p_{yz} \\ 0 & 0 & 0 \end{matrix}] \nabla_{t} \times Φ_{t} + [\begin{matrix} p_{xx} & p_{xy} & 0 \\ p_{yx} & p_{yy} & 0 \\ 0 & 0 & 0 \end{matrix}] \nabla_{t} \times \hat{z} Φ_{z})

Equations (12) and (13) can be used to construct a full-vectorial eigenvalue algorithm for the analysis of band structures of 2D anisotropic PCs with arbitrary permeability and permittivity tensors under the in-plane wave propagation.

2.2. The FEM based matrix eigenvalue equation

Notice that all field components are involved simultaneously in Eqs. (12) and (13). Therefore, to develop an FEM based matrix eigenvalue equation according to Eqs. (12) and (13), the simple edge or nodal elements are apparently insufficient. In addition, to achieve better accuracy and efficiency, a higher order element which can model the curved geometry in the structures better is preferred. Consequently, the curvilinear hybrid edge/nodal element [19,20], as shown in Fig. 2, is adopted for the discretization of the FEM. An edge element with eight variables, ϕ _t1-ϕ _t8, based on linear tangential and quadratic normal (LT/QN) vector basis functions and a quadratic nodal element with six variables, ϕ _z1-ϕ _z6, are employed for the transverse and longitudinal field components, respectively. The vector-based shape function for the curvilinear edge element, x̂{U}+ŷ{V}, and the scalar-based shape function for the curvilinear nodal element, {N}, are listed in Table 1, where {ϕ ^e _t} and {ϕ ^e _z} are, respectively, the edge- and nodal-variable vectors for each element, L _i (i=1,2,3) are the simplex coordinates, |J| is the Jacobian, the determinant of the Jacobian matrix, [J], defined by

[J] = [\begin{matrix} \frac{\partial x}{\partial L_{1}} & \frac{\partial y}{\partial L_{1}} \\ \frac{\partial x}{\partial L_{2}} & \frac{\partial y}{\partial L_{2}} \end{matrix}]

and |∇_t L _i|_j and |J|_j are, respectively, the values of |∇_t L _i| (i=1,2,3) and |J| at the nodal point j (j=1,2, …, 6).

The relation of differentiation between the Cartesian coordinate system and the simplex coordinate system can be written as

[\begin{matrix} \frac{\partial}{\partial L_{1}} \\ \frac{\partial}{\partial L_{2}} \end{matrix}] = [J] [\begin{matrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{matrix}]

and the integration of a function f (x,y) in the Cartesian coordinate system can be performed in the simplex coordinate system through

\int \int_{e} f (x, y) dxdy = \int_{0}^{1} [\int_{0}^{1 - L_{1}} f (L_{1}, L_{2}, L_{3}) ∣ J (L_{1}, L_{2}, L_{3}) ∣ d L_{2}] d L_{1} .

Hence, the required computation for element matrices associated with the governing equation can be obtained directly in the simplex coordinate system.

Fig. 2. Curvilinear hybrid edge/nodal element.

Download Full Size | PDF

Dividing the structure domain into a number of curvilinear hybrid edge/nodal elements, the field Φ in each element can be expanded as

Φ = Φ_{t} + \hat{z} Φ_{z} = [\begin{matrix} {U}^{T} {ϕ_{t}^{e}} \\ {V}^{T} {ϕ_{t}^{e}} \\ {N}^{T} {ϕ_{z}^{e}} \end{matrix}]

where T denotes transpose. Applying Galerkin’s method, assembling all element matrices, and incorporating the periodic boundary conditions (PBCs) [14] given by

Table 1. Vector- and scalar-based shape functions

View Table

{Φ ∣}_{I} = e^{- j k_{x} a} {Φ ∣}_{III}

{\frac{\partial Φ}{\partial x} ∣}_{I} = e^{- j k_{x} a} {\frac{\partial Φ}{\partial x} ∣}_{III}

{\frac{\partial Φ}{\partial y} ∣}_{I} = e^{- j k_{x} a} {\frac{\partial Φ}{\partial y} ∣}_{III}

{Φ ∣}_{II} = e^{- j k_{y} a} {Φ ∣}_{IV}

{\frac{\partial Φ}{\partial x} ∣}_{II} = e^{- j k_{y} a} {\frac{\partial Φ}{\partial x} ∣}_{IV}

{\frac{\partial Φ}{\partial y} ∣}_{II} = e^{- j k_{y} a} {\frac{\partial Φ}{\partial y} ∣}_{IV}

for square-arranged lattice and

{Φ ∣}_{I} = e^{- j (k_{x} \frac{\sqrt{3} a}{2} + k_{y} \frac{a}{2})} {Φ ∣}_{IV}

{\frac{\partial Φ}{\partial x} ∣}_{I} = e^{- j (k_{x} \frac{\sqrt{3} a}{2} + k_{y} \frac{a}{2})} {\frac{\partial Φ}{\partial x} ∣}_{IV}

{\frac{\partial Φ}{\partial y} ∣}_{I} = e^{- j (k_{x} \frac{\sqrt{3} a}{2} + k_{y} \frac{a}{2})} {\frac{\partial Φ}{\partial y} ∣}_{IV}

{Φ ∣}_{II} = e^{- j k_{y} a} {Φ ∣}_{V}

{\frac{\partial Φ}{\partial x} ∣}_{II} = e^{- j k_{y} a} {\frac{\partial Φ}{\partial x} ∣}_{V}

{\frac{\partial Φ}{\partial y} ∣}_{II} = e^{- j k_{y} a} {\frac{\partial Φ}{\partial y} ∣}_{V}

{Φ ∣}_{III} = e^{j (k_{x} \frac{\sqrt{3} a}{2} - k_{y} \frac{a}{2})} {Φ ∣}_{VI}

{\frac{\partial Φ}{\partial x} ∣}_{III} = e^{j (k_{x} \frac{\sqrt{3} a}{2} - k_{y} \frac{a}{2})} {\frac{\partial Φ}{\partial x} ∣}_{VI}

{\frac{\partial Φ}{\partial y} ∣}_{III} = e^{j (k_{x} \frac{\sqrt{3} a}{2} - k_{y} \frac{a}{2})} {\frac{\partial Φ}{\partial y} ∣}_{VI}

for triangle-arranged lattice, where k _x and k _y are, respectively, the wavenumbers in the x and y directions, and the subscript denotes the side label of the unit cell, as defined in Fig. 1(a) and (b), Eqs. (12) and (13) can be transformed into the matrix form as

[K] {ϕ} - k_{0}^{2} [M] {ϕ} = {0}

where {0} is a null vector, and the variable vector {ϕ} and the matrices, [K] and [M], are given by

{ϕ} = [\begin{matrix} {ϕ_{t}} \\ {ϕ_{z}} \end{matrix}]

[K] = [\begin{matrix} [K_{tt}] & [K_{tz}] \\ [K_{zt}] & [K_{zz}] \end{matrix}]

[M] = [\begin{matrix} [M_{tt}] & [M_{tz}] \\ [M_{zt}] & [M_{zz}] \end{matrix}]

with

[K_{tt}] = \sum_{e} \iint [p_{zz} \frac{\partial {V}}{\partial x} \frac{{\partial {V}}^{T}}{\partial x} - p_{zz} \frac{\partial {V}}{\partial x} \frac{{\partial {U}}^{T}}{\partial y}

[K_{tz}] = \sum_{e} \iint [p_{zx} \frac{\partial {V}}{\partial x} \frac{{\partial {N}}^{T}}{\partial x} - p_{zy} \frac{\partial {V}}{\partial x} \frac{{\partial {N}}^{T}}{\partial y}

[K_{zt}] = \sum_{e} \iint [p_{xz} \frac{\partial {N}}{\partial x} \frac{{\partial {V}}^{T}}{\partial x} - p_{xz} \frac{\partial {N}}{\partial x} \frac{{\partial {U}}^{T}}{\partial y}

[K_{zz}] = \sum_{e} \iint [p_{xx} \frac{\partial {N}}{\partial x} \frac{{\partial {N}}^{T}}{\partial x} - p_{xy} \frac{\partial {N}}{\partial x} \frac{{\partial {N}}^{T}}{\partial y}

[M_{tt}] = \sum_{e} \iint [q_{xx} {U} {U}^{T} + q_{xy} {U} {V}^{T}

[M_{tz}] = \sum_{e} \iint [q_{xz} {U} {N}^{T} + q_{yz} {V} {N}^{T}] dxdy

[M_{zt}] = \sum_{e} \iint [q_{zx} {N} {U}^{T} + q_{zy} {N} {V}^{T}] dxdy

[M_{zz}] = \sum_{e} \iint [q_{zz} {N} {N}^{T}] dxdy .

where Σ_e extends over all different elements. Equation (33) is the desired matrix eigenvalue equation in the proposed full-vectorial algorithm for the analysis of 2D anisotropic PCs with arbitrary permittivity and permeability tensors.

3. Numerical results

In the following, the band structures of several 2D PCs will be analyzed using the proposed full-vectorial algorithm. For the cases in which the wave modes can be decoupled into TE and TM modes, we also show the results obtained from the scalar algorithm in [14] for comparison and discussions.

3.1. Isotropic PCs

To understand the characteristics of this full-vectorial algorithm, we start from the band-structure analysis of two familiar 2D isotropic PCs. The first is formed by square-arranged parallel alumina rods with relative permittivity ε=8.9 and radius r=0.2a in the air, and the second consists of triangle-arranged dielectric cylinders with relative permittivity ε=11.4 and radius r=0.2a in the air. The first BZ for the 2D PC with square lattice is shown in Fig. 3(a), in which four sub-zones are marked. For isotropic PCs, considering all the possible directions of the wave vector k in any single sub-zone, identical to the IBZ, is sufficient for the construction of complete band structures. Similarly, the six sub-zones in the first BZ for the 2D PC with triangular lattice, as shown in Fig. 3(b), are exactly the same as isotropic PCs. Consequently, the band structures of the two isotropic PCs considered here can be completely constructed from sub-zones Γ-X-M and Γ-M-K, respectively. The band structures of these two PCs obtained from the full-vectorial algorithm are, respectively, shown in Fig. 4(a) and (b), in which the band structures constructed from the scalar algorithm [14] are also plotted for comparison. The blue solid lines denote the results obtained using the full-vectorial algorithm, and the red and green circular dots represent, respectively, the results of the TE and TM modes calculated from the scalar algorithm. With the validation using the scalar algorithm, it is believed that the full-vectorial algorithm can correctly construct the band structures of the TE and TM modes for 2D isotropic PCs.

Fig. 3. The first BZ of a 2D PC with (a) square lattice and (b) triangular lattice.

Download Full Size | PDF

We further examine the convergence behaviors of the full-vectorial algorithm, compared with the scalar algorithm. For the 2D PC with square lattice, k is fixed at the X point in the first BZ and various numbers of elements are used to search the eigen frequencies of several bands for both TE and TM modes. The convergence behaviors of the first and third bands for the TE mode are, respectively, shown in Fig. 5(a) and (b), and those for the TM mode are, respectively, shown in Fig. 5(c) and (d). In Fig. 5(a) and (b), the lines with blue triangles, blue squares, and red circles represent the results calculated from the E-formulation full-vectorial algorithm, the H-formulation full-vectorial algorithm, and the TE scalar algorithm, respectively. It can be seen that the eigen frequencies obtained from the E- and H-formulation full-vectorial algorithms will converge to the same value with enough number of elements used, although their convergence behaviors are obviously different from each other. On the other hand, the agreement of the convergence behaviors between the H-formulation full-vectorial algorithm and the TE scalar algorithm provides a powerful support for the correctness of this full-vectorial algorithm. According to the definitions in [14], the TE mode is composed of H _z, E _x, and E _y components, and the unknown field component to be solved in the scalar algorithm for the TE mode is H _z. Therefore, when the H-formulation is employed in the full-vectorial algorithm, H _z will appear in the unknown field vector explicitly, and it is quite reasonable that the results from these two algorithms would be consistent with each other.

Fig. 4. The band structures for the 2D isotropic PC formed by (a) square-arranged alumina rods with ε=8.9 and r=0.2a in the air and (b) triangle-arranged dielectric cylinders with ε=11.4 and r=0.2a in the air.

Download Full Size | PDF

Fig. 5. The convergence behaviors of the full-vectorial algorithm for the 2D isotropic PC with square lattice compared with the scalar algorithm as k is fixed at the X point. (a) TE first band; (b) TE third band; (c) TM first band; (d) TM third band.

Download Full Size | PDF

In Fig. 5(c) and (d), the lines with green circles represent the results calculated from the TM scalar algorithm, and the lines with blue triangles and blue squares stand for the same meanings as in Fig. 5(a) and (b). Because the unknown field component to be solved in the scalar algorithm for the TM mode is E _z, it is no surprise that its convergence behaviors agree well with those of the E-formulation full-vectorial algorithm.

As for the 2D PC with triangular lattice, k is fixed at theMpoint in the first BZ, and the counterparts of Fig. 5(a), (b), (c), and (d) are displayed in Fig. 6(a), (b), (c), and (d), respectively. We can observe that the convergence behaviors of the E- and H-formulation full-vectorial algorithms, and the scalar algorithms for the TE and TM modes for triangular lattice are similar to those for square lattice. Hence, the discussions for the square lattice can be applied to the triangular lattice as well, and the correctness of the proposed full-vectorial algorithm for solving the band structures of isotropic 2D PCs is verified again.

Fig. 6. The convergence behaviors of the full-vectorial algorithm for the 2D isotropic PC with triangular lattice compared with the scalar algorithm as k is fixed at the M point. (a) TE first band; (b) TE third band; (c) TM first band; (d) TM third band.

Download Full Size | PDF

3.2. Anisotropic PCs

We then study the case of 2D anisotropic PCs. To provide a structure for real applications, the materials used to comprise the PCs are chosen carefully. The anisotropic material Te, a kind of positive uniaxial crystal with principal indices of n _o=4.8 and n _e=6.2, can be used to provide a chance of obtaining the absolute band gap [21]. In addition, liquid crystals (LCs) are widely used in the research of tunable PCs for the convenient change of the anisotropy by simply rotating the LC molecules, and the tunable PBGs can be achieved [13]. To utilize the advantages of these two materials collectively, we consider a 2D anisotropic PC which is composed of triangle-arranged LC columns in the background of Te, as depicted in Fig. 7(a). The radius r of the LC columns is 0.45a and the components of the relative permittivity tensor [ε _r] of the nematic LCs (5CB) [22] we use are given as

ε_{xx} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \cos^{2} ϕ_{c}

ε_{xx} = ε_{yx} = (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \sin ϕ_{c} \cos ϕ_{c}

ε_{xz} = ε_{zx} = (n_{e}^{2} - n_{o}^{2}) \sin θ_{c} \cos ϕ_{c} \cos ϕ_{c}

ε_{yy} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \sin^{2} ϕ_{c}

ε_{yz} = ε_{zy} = (n_{e}^{2} - n_{o}^{2}) \sin θ_{c} \cos θ_{c} \sin ϕ_{c}

ε_{zz} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \cos^{2} θ_{c}

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. 7(b).

Fig. 7. (a) Cross-section of the 2D anisotropic PC with triangular lattice. (b) Schematic definition of rotation angles for the LC molecule.

Download Full Size | PDF

With the introduction of anisotropic materials into the PC structures, the four sub-zones in the first BZ for the 2D PC with square lattice, as shown in Fig. 3(a), and the six sub-zones in the first BZ for the 2D PC with triangular lattice, as shown in Fig. 3(b), must be taken into consideration for the construction of complete band structures [14]. Fig. 8 shows the complete band structures for this 2D PC of θc=0°, 30°, 45°, 60°, and 90°, and ϕ _c=30°. Notice that for the special cases, i.e. θ _c=0° or 90°, the scalar algorithm in [14] can be employed to obtain the band structures as well. Consequently, in Fig. 8(a) and (e), for the cases of θ _c=0° and 90°, respectively, not only the results from the full-vectorial algorithm, represented by the blue solid lines, but also those from the scalar algorithm, represented by the red and green circular dots, are displayed. It can be seen that, even for anisotropic PCs, the results of these two algorithms still agree with each other very well provided that the wave modes can be decoupled into the TE and TM modes. The significance of the full-vectorial algorithm is conspicuously revealed from Fig. 8(b), (c), and (d). These band structures can not be obtained from the scalar algorithm because the wave modes in the PC are hybrid ones when θ _c is not exactly 0° or 90°. However, by taking advantage of the full-vectorial algorithm, we can construct these band structures smartly.

In Fig. 8, the red region stands for the absolute band gap for all directions in the x-y plane of this PC, and the shadowed regions represent ignorable parts as a result of the symmetry conditions and the basic principle of periodic structures [14]. We can see that there always exist absolute band gaps for this PC of different θ _c’s when ϕ _c=30°. To systematically understand the influences of θc and ϕc on the absolute band gaps and clearly demonstrate how the anisotropy affects the wave modes in the PC structures, we analyze the band structures of more different θ _c’s and ϕ _c’s. Fig. 9 shows the normalized frequency range of the absolute band gap versus ϕc for five different θ _c’s from 0° to 90°. We can find out that the absolute band gap can be noticeably tuned by changing θ _c regardless of the value of ϕ _c, and the upper limit of the absolute band gap is much dependent on θ _c than the lower limit for this PC. On the other hand, it is seen that, for every fixed θ _c, the normalized frequency range of the absolute band gap almost keeps constant when ϕ _c increases from 0° to 90°.

For θ _c=0° or 90°, the above phenomenon is absolutely reasonable and can be directly understood by the aid of Fig. 8(a) and (e), from which we can observe that the normalized frequency range of the absolute band gap is determined by the TM modes which are independent of the

value of ϕ _c since the electric field is perpendicular to the x-y plane. However, for θ _c=30°, 45°, and 60°, this ϕ _c-independent phenomenon seems unreasonable. Theoretically, in these conditions all the wave modes, including those to determine the absolute band gap, in the PC are hybrid ones and are expected to be dependent on the anisotropy of the PC structure, which is related to the orientation of the LC molecules. Consequently, rotating the LC molecules around the z-axis should influence the absolute band gap in some way. To confirm this theoretical argument concretely, the numerical data are presented in another way which can obviously show the variation of the absolute band gap. Specifically, for every fixed θ _c, the calculated absolute band gap limits for any ϕ _c are normalized through the transformation defined as

Fig. 8. The band structures of the 2D anisotropic PC with triangular lattice of (a) θ _c=0° and ϕ _c=30°, (b) θ _c=30° and ϕ _c=30°, (c) θ _c=45° and ϕ _c=30°, (d) θ _c=60° and ϕ _c=30°, and (e) θ _c=90° and ϕ _c=30°.

Download Full Size | PDF

Fig. 9. The normalized frequency range of the absolute band gap for the 2D anisotropic PC with triangular lattice versus ϕ _c for five different θ _c’s from 0° to 90°.

Download Full Size | PDF

Fig. 10. The variation of the absolute band gap limits for the 2D anisotropic PC with triangular lattice versus ϕ _c for five different θ _c’s from 0° to 90°.

Download Full Size | PDF

{\frac{ω a}{2 π c} ∣}_{normalized} = \frac{ω a}{2 π c} - \frac{{\frac{ω a}{2 π c} ∣}_{max} + {\frac{ω a}{2 π c} ∣}_{min}}{2}

where ωa/2πc|_max and ωa/2πc|_min are, respectively, the maximum and minimum values of the limits for ϕ _c from 0° to 90°. The results after the transformation of Eq. (51) are shown in Fig. 10. It is seen that the absolute band gap limits indeed vary with ϕ _c for θ _c=30°, 45°, and 60°. We also observe that the upper limit is still much dependent on ϕ _c than the lower limit for this PC, just as their dependence on θ _c. Another interesting phenomenon is that the variations of the upper and lower limits are opposite, i.e. when the upper limit reaches the maximum value, the lower limit will reach the minimum value, and vice versa. Finally, it is noticed that the variation of the limits has a period of 60° over ϕ _c. This phenomenon can be reasonably accepted because the unit cell of a PC with triangular lattice is a hexagon and rotating the LC molecules around the z-axis by 60° dose not really change the PC structure due to the rotation symmetry. Based on these numerical results and theoretical discussions, the correctness and usefulness of the proposed full-vectorial algorithm have been clearly demonstrated.

4. Conclusion

We have successfully developed a full-vectorial FEM based eigenvalue algorithm for the analysis of band structures of 2D PCs with arbitrary 3D anisotropy for in-planewave propagations. In this full-vectorial algorithm, we attempt an idea of considering all the field components simultaneously and formulate a matrix eigenvalue equation to deal with the most general situation of 2D anisotropic PC problems. The relationship of the traditional scalar algorithm and this full-vectorial algorithm is described through the analysis of isotropic PCs first. Then, we examine how the rotation of the LC molecules in the considered anisotropic PC structure influences the wave modes and the absolute band gap. The effects of θ _c and ϕ _c on the behaviors of the absolute band gap limits are discussed in detail. By considering the numerical results carefully and explaining the observed phenomena theoretically, we demonstrate the performance and significance of this full-vectorial algorithm, and see that a deeper insight about 2D anisotropic PCs can be obtained through this full-wave analysis.

Acknowledgments

This workwas supported in part by the National Science Council of the Republic of China under grant NSC95-2221-E-002-328, in part by the Ministry of Economic Affairs of the Republic of China under grant 95-EC-17-A-08-S1-0006, in part by the Excellent Research Projects of National Taiwan University under grant 95R0062-AE00-06, and in part by the Ministry of Education of the Republic of China under “The Aim of Top University Plan” grant.

References and links

1. E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. 58, 2059–2062 (1987). [CrossRef] [PubMed]

2. S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 2486–2489 (1987). [CrossRef] [PubMed]

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). [CrossRef] [PubMed]

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). [CrossRef]

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). [CrossRef]

7. 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]

8. 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).

9. C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express12, 1397–1408 (2004), http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-7-1397 [CrossRef] [PubMed]

10. 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). [CrossRef]

11. I. H. H. Zabel and D. Stroud, “Photonic band structures of optically anisotropic periodic arrays,” Phys. Rev. B 48, 5004–5012 (1993). [CrossRef]

12. 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). [CrossRef]

13. 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). [CrossRef]

14. S. M. Hsu, M. M. Chen, and H. C. Chang, “Investigation of band structures for 2D non-diagonal anisotropic photonic crystals using a finite element method based eigenvalue algorithm,” Opt. Express15, 5416–5430 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-9-5416 [CrossRef] [PubMed]

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). [CrossRef]

16. G. E. Antilla and N. G. Alexopoulos, “Scattering from complex three-dimensional geometries by a curvilinear hybrid finite-element-integral equation approach,” J. Opt. Soc. Am. A 11, 1445–1457 (1994). [CrossRef]

17. L. Zhang and N. G. Alexopoulos, “Finite-element based techniques for the modeling of PBG materials,” Electromagnetics 19, 225–239 (1999). [CrossRef]

18. D. Sievenpiper, L. Zhang, R. F. J. Broas, N. G. Alexopoulos, and E. Yablonovitch, “High-impedance electromagnetic surfaces with a forbidden frequency band,” IEEE Trans. Microwave Theory Tech. 47, 2059–2074 (1999). [CrossRef]

19. J. Jin, The Finite Element Method in Electromagnetics (John Wiley and Sons, Inc., New York, 2002).

20. M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems,” J. Lightwave Technol. 18, 737–743 (2000). [CrossRef]

21. 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). [CrossRef]

22. P. Yeh and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons, Inc., New York, 1999).

edge	{ϕ ^e _t }	x̂{U}+yŷ{V}	nodal	{ϕ ^e _z }	{N}
	ϕ _t1	\|J\|₁\|∇ _t L ₃\|₁(L ₁∇ _t L ₂)		ϕ _z1	L ₁(2L ₁-1)
	ϕ _t2	\|J\|₂\|∇ _t L ₁\|₂(L ₂∇_t L ₃)		ϕ _z2	L ₂(2L ₂-1)
	ϕ _t3	\|J\|₃\|∇ _t L ₂\|₃(L ₃∇ _t L ₁)		ϕ _z3	L ₃(2L ₃-1)
LT/QN	ϕ _t4	\|J\|₂\|∇ _t L ₃\|₂(L ₂∇ _t L ₁)	quadratic	ϕ _z4	4L ₁ L ₂
	ϕ _t5	\|J\|₃\|∇ _t L ₁\|₃(L ₃∇ _t L ₂)		ϕ _z5	4L ₂ L ₃
	ϕ _t6	\|J\|₁\|∇ _t L ₂\|₁(L ₁∇ _t L ₃)		ϕ _z6	4L ₃ L ₁
	ϕ _t7	4L ₁ L ₂∇ _t L ₃/\|∇ _t L ₃\|₄
	ϕ _t8	4L ₂ L ₃∇ _t L ₁/\|∇ _t L ₁\|₅

Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy

Abstract

1. Introduction

2. Formulation

2.1. The governing equation

2.2. The FEM based matrix eigenvalue equation

3. Numerical results

3.1. Isotropic PCs

3.2. Anisotropic PCs

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (10)

Tables (1)

Equations (58)

Optics Express