**1. Introduction**

28 Photonic Crystals Book 1

318 Photonic Crystals – Introduction, Applications and Theory

Yablonovitch, E. (1987). Inhibited spontaneous emission in solid-state physics and electronics,

Yablonovitch, E. (1993). Photonic band-gap crystals, *Journal of Physics: Condensed Matter*

Zhang, Z., Yoshie, T., Zhu, X., Xu, J. & Scherer, A. (2006). Visible two-dimensional photonic

crystal slab laser, *Applied Physics Letters* 89(7): 071102 –071102–3.

*Physical Review Letters* 58(20): 2059–2062.

5(16): 2443–2460.

Photonic crystals are modern artificially designed periodical systems capable to affect the motion of photons in a similar way that the periodicity of the atomic potential in a semiconductor crystal affects the motion of electrons. The physical properties of light in a photonic crystal resemble those of electrons in atomic crystals, leading to forbidden propagation of electromagnetic modes at certain frequencies, as demonstrated by Yablonovitch (1987), Ho et al. (1990), Joannopoulos et al. (1995; 1997), etc. The existence of the optical band gap (which is the part of the spectrum for which the wave propagation is not possible) makes photonic crystals broadly interesting from many viewpoints of fundamental research and applications. Recent studies are motivated by promising applications such as purely optical integrated circuits [e.g., by Lin et al. (1998), Noda (2006), or Hugonin et al. (2007)], artificial metamaterials with high tunability [Datta et al. (1993); Genereux et al. (2001); Krokhin et al. (2002); Reyes et al. (2005)], high-sensitivity photonic biosensors [Block et al. (2008); Skivesen et al. (2007)], or devices based on phenomena not accessible in conventional media [Benisty (2009); Kosaka et al. (1998); Krokhin & Reyes (2004)]. It is also necessary for accounting for structural colors of wings of butterflies or beetles, feathers of birds, or iridescent plants [Kinoshita & Yoshioka (2005); Vukusic & Sambles (2003)].

Since the optical properties of photonic crystals strongly depend on their geometrical structure, used materials, etc., their proper design is crucial for the correct device functionality. Thorough theoretical analysis therefore takes place in the development, using various numerical methods for calculation including methods of finite difference in the time or frequency domain or finite element methods [Joannopoulos et al. (1995)]. One of the most common calculation techniques applied to photonic crystals is the plane wave expansion method, which is a frequency-domain approach based on the expansion of the fields and material parameters into the Fourier (or reciprocal) space. The components of this expansion represent the definite-frequency states. After some necessary truncation of the complete basis (plane waves with a finite cutoff), the partial differential equations are then solved as a linear-algebraic problem. However, the convergence rate of this method strongly depends on the implementation of Maxwell's equations in the truncated plane-wave basis [Meade et al. (1993); Sozuer et al. (1992)]. In the case of periodic discontinuities (typical for photonic crystals) the convergence is rather poor so that the computer calculations might become extremely time- and memory-consuming.

and 7 will propose how to factorize anisotropic and three-dimensional (3D) photonic crystals,

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 321

The modes of photonic crystals are in principle the eigensolutions of the wave equation with an inhomogeneous, periodic relative permittivity *ε*(r). One possible version of the wave

∇ × (∇ × <sup>E</sup>) = *<sup>ω</sup>*<sup>2</sup>

<sup>−</sup>*ik*0(*mpx*+*nqy*+*lsz*) <sup>=</sup> ∑

is its eigenfunction, which has the form of a pseudoperiodic Floquet–Bloch function. Here *p* = 2*π*/Λ*x*, *q* = 2*π*/Λ*y*, and *s* = 2*π*/Λ*<sup>z</sup>* are the normalized reciprocal lattice vectors. For simplicity we assume the periods Λ*<sup>j</sup>* along the Cartesian axes throughout this chapter. For brevity we have also defined k = *k*0[*p*0, *q*0,*s*0], *pm* = *p*<sup>0</sup> + *mp*, *qn* = *q*<sup>0</sup> + *nq*, and *sl* = *s*<sup>0</sup> + *ls*. (Analogously we could write the wave equation for an unknown magnetic field H or any

Owing to the periodicity of the problem, the plane wave expansion method is the reference method for the mode calculation. It is based on the Fourier expansion of the field such as in Equation 2 and on the Fourier expansion of a material function, either the permittivity or the

*εmnle*

*ηmnle*

The rules for choosing the most appropriate material parameter and the most appropriate field for the Fourier expansion are governed by various methods of Fourier factorization. In the past, the E method (*η* and the electric displacement D were expanded), H method (*η* and H were expanded), and Ho method (*ε* and E were expanded) were the typical choices.

Now we carry out the transformation of the partial differential equations into matrix equations in order to solve the eigenproblem by linear-algebraic methods. For simplicity we limit ourselves to 1D and 2D photonic crystals, and always choose the direction of propagation

*εmn e*

*fmn e*

−*i*(*mpx*+*nqy*)

−*i*(*pmx*+*qny*)

+∞ <sup>∑</sup>*m*,*n*=−<sup>∞</sup>

+∞ <sup>∑</sup>*m*,*n*=−<sup>∞</sup> *m*,*n*,*l*

*emnle*

where *ω*2/*c*<sup>2</sup> is its eigenvalue (*ω* is the frequency and *c* the light velocity in vacuum) and

*<sup>c</sup>*<sup>2</sup> <sup>E</sup>, (1)

<sup>−</sup>*ik*0(*mpx*+*nqy*+*lsz*) (3)

<sup>−</sup>*ik*0(*mpx*+*nqy*+*lsz*) (4)

, (5)

. (6)

<sup>−</sup>*ik*0(*pmx*+*qny*+*sl <sup>z</sup>*) (2)

equation is the equation for the unknown electric field with an unknown frequency,

1 *ε*

*emnle*

*ε*(r) = ∑ *m*,*n*,*l*

*η*(r) = ∑

in the *xy* plane, so that *∂<sup>z</sup>* = 0. With these restrictions we write

*ε*(*x*, *y*) =

*f*(*x*, *y*) =

*m*,*n*,*l*

respectively.

**2. Plane wave expansion**

E(r) = *e*

other field from Maxwell's equations.)

impermittivity *η*(r) = 1/*ε*(r),

**2.2 Matrix notation**

<sup>−</sup>*i*k·<sup>r</sup> ∑ *m*,*n*,*l*

**2.1 General remarks**

Because the underlying physical phenomenon in the optical behavior of photonic crystals is based on diffraction (therefore the lattice constant of a periodic structure has to be in the same length-scale as half the wavelength of the electromagnetic wave), several conclusions can be advantageously adopted from the classical coupled-wave theory, one of the most effective methods for modeling diffraction of electromagnetic waves by periodic gratings, which was developed during several past decades. This can provide a high enhancement to the plane wave expansion method, resulting in the reduction of the computation resources.

In the mid 1990s Li (1996) showed that Laurent's "direct" rule, which had always been adopted in conventional formulations to factorize the truncated Fourier series that corresponds to products of two periodic functions, presents bad convergence when the two functions of the product are simultaneously discontinuous. He suggested three Fourier factorization rules (briefly summarized in Section 2.3) and applied them to one-dimensional (1D) diffraction gratings. This major breakthrough in the grating theory (called "fast Fourier factorization") was soon applied by many authors to various grating structures with arbitrary periodic reliefs, anisotropic [Li (1998)] and slanted [Chernov et al. (2001)] periodic systems, their various combinations [Li (2003); Watanabe (2002); Watanabe et al. (2002)], and other systems [Bonod et al. (2005a;b); Boyer et al. (2004)].

Later Li (1997) applied the factorization rules to two-dimensional (2D) periodic structures treated by "zigzag" Fourier expansion, which yielded an improvement for rectangular dots or holes. However, Popov & Neviere (2000) have pointed out that the staircase approximation (of the coupled wave theory using the slicing of relief profiles) in combination with the traditional formulation of differential equation within one slice violates Li's factorization rules. This was a major complication for the analysis of the periodic systems made of rounded elements. Therefore, they applied a coordinate transform to treat individually the normal and tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the application of the correct rule for each field component and thus improved the convergence.

Later David et al. (2006) utilized the normal–tangential field separation to 2D photonic crystals composed of circular or elliptical holes. Similarly, Schuster et al. (2007) applied this method to 2D gratings, and also suggested more general distributions of polarization bases [Gotz et al. (2008)]. These approaches, always dealing with linear polarizations, enabled a significant improvement of the convergence properties, but ignored the fact that the transformation matrix between the Cartesian and the normal–tangential component bases of polarization became discontinuous at the center and along the boundaries of the periodic cell, which slowed down the resulting convergence. To overcome these discontinuities, a distribution of more complex (i.e., generally elliptic) polarization bases was recently suggested to improve optical simulations of 2D gratings and photonic crystals [Antos (2009); Antos & Veis (2010)].

Our chapter will describe in detail the application of the Fourier factorization rules to the plane wave expansion method for numerical analysis of general photonic crystals. Section 2 will introduce the principle of the plane wave expansion together with the notation of matrices and factorization theorems. Section 3 will refer to 1D photonic structures made as periodic stratified media. The consistency of the correct factorization rules with classical theory of Yeh et al. (1977) and Yariv & Yeh (1977) will be shown, pointing to the correct boundary conditions of the tangential components of the electric and magnetic field on multilayer interfaces. Section 4 will repeat our previously described methodology for 2D photonic crystals made of circular elements, and Section 5 will generalize it to elements of other shapes. Sections 6 and 7 will propose how to factorize anisotropic and three-dimensional (3D) photonic crystals, respectively.

#### **2. Plane wave expansion**

#### **2.1 General remarks**

2 Will-be-set-by-IN-TECH

Because the underlying physical phenomenon in the optical behavior of photonic crystals is based on diffraction (therefore the lattice constant of a periodic structure has to be in the same length-scale as half the wavelength of the electromagnetic wave), several conclusions can be advantageously adopted from the classical coupled-wave theory, one of the most effective methods for modeling diffraction of electromagnetic waves by periodic gratings, which was developed during several past decades. This can provide a high enhancement to the plane

In the mid 1990s Li (1996) showed that Laurent's "direct" rule, which had always been adopted in conventional formulations to factorize the truncated Fourier series that corresponds to products of two periodic functions, presents bad convergence when the two functions of the product are simultaneously discontinuous. He suggested three Fourier factorization rules (briefly summarized in Section 2.3) and applied them to one-dimensional (1D) diffraction gratings. This major breakthrough in the grating theory (called "fast Fourier factorization") was soon applied by many authors to various grating structures with arbitrary periodic reliefs, anisotropic [Li (1998)] and slanted [Chernov et al. (2001)] periodic systems, their various combinations [Li (2003); Watanabe (2002); Watanabe et al. (2002)], and other

Later Li (1997) applied the factorization rules to two-dimensional (2D) periodic structures treated by "zigzag" Fourier expansion, which yielded an improvement for rectangular dots or holes. However, Popov & Neviere (2000) have pointed out that the staircase approximation (of the coupled wave theory using the slicing of relief profiles) in combination with the traditional formulation of differential equation within one slice violates Li's factorization rules. This was a major complication for the analysis of the periodic systems made of rounded elements. Therefore, they applied a coordinate transform to treat individually the normal and tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the application of the correct rule for each field component and thus improved the convergence. Later David et al. (2006) utilized the normal–tangential field separation to 2D photonic crystals composed of circular or elliptical holes. Similarly, Schuster et al. (2007) applied this method to 2D gratings, and also suggested more general distributions of polarization bases [Gotz et al. (2008)]. These approaches, always dealing with linear polarizations, enabled a significant improvement of the convergence properties, but ignored the fact that the transformation matrix between the Cartesian and the normal–tangential component bases of polarization became discontinuous at the center and along the boundaries of the periodic cell, which slowed down the resulting convergence. To overcome these discontinuities, a distribution of more complex (i.e., generally elliptic) polarization bases was recently suggested to improve optical simulations of 2D gratings and photonic crystals [Antos (2009); Antos & Veis (2010)]. Our chapter will describe in detail the application of the Fourier factorization rules to the plane wave expansion method for numerical analysis of general photonic crystals. Section 2 will introduce the principle of the plane wave expansion together with the notation of matrices and factorization theorems. Section 3 will refer to 1D photonic structures made as periodic stratified media. The consistency of the correct factorization rules with classical theory of Yeh et al. (1977) and Yariv & Yeh (1977) will be shown, pointing to the correct boundary conditions of the tangential components of the electric and magnetic field on multilayer interfaces. Section 4 will repeat our previously described methodology for 2D photonic crystals made of circular elements, and Section 5 will generalize it to elements of other shapes. Sections 6

wave expansion method, resulting in the reduction of the computation resources.

systems [Bonod et al. (2005a;b); Boyer et al. (2004)].

The modes of photonic crystals are in principle the eigensolutions of the wave equation with an inhomogeneous, periodic relative permittivity *ε*(r). One possible version of the wave equation is the equation for the unknown electric field with an unknown frequency,

$$\frac{1}{c}\nabla \times (\nabla \times \mathbf{E}) = \frac{\omega^2}{c^2} \mathbf{E},\tag{1}$$

where *ω*2/*c*<sup>2</sup> is its eigenvalue (*ω* is the frequency and *c* the light velocity in vacuum) and

$$E(\mathbf{r}) = e^{-i\mathbf{k}\cdot\mathbf{r}} \sum\_{m,n,l} e\_{mnl} e^{-i\mathbf{k}\_0(m\mathbf{p}\times+n\mathbf{q}\mathbf{y}+l\mathbf{s}z)} = \sum\_{m,n,l} e\_{mnl} e^{-i\mathbf{k}\_0(p\_mx+q\_ny+s\_lz)}\tag{2}$$

is its eigenfunction, which has the form of a pseudoperiodic Floquet–Bloch function. Here *p* = 2*π*/Λ*x*, *q* = 2*π*/Λ*y*, and *s* = 2*π*/Λ*<sup>z</sup>* are the normalized reciprocal lattice vectors. For simplicity we assume the periods Λ*<sup>j</sup>* along the Cartesian axes throughout this chapter. For brevity we have also defined k = *k*0[*p*0, *q*0,*s*0], *pm* = *p*<sup>0</sup> + *mp*, *qn* = *q*<sup>0</sup> + *nq*, and *sl* = *s*<sup>0</sup> + *ls*. (Analogously we could write the wave equation for an unknown magnetic field H or any other field from Maxwell's equations.)

Owing to the periodicity of the problem, the plane wave expansion method is the reference method for the mode calculation. It is based on the Fourier expansion of the field such as in Equation 2 and on the Fourier expansion of a material function, either the permittivity or the impermittivity *η*(r) = 1/*ε*(r),

$$\varepsilon(\mathbf{r}) = \sum\_{m,n,l} \varepsilon\_{mnl} e^{-i\mathbf{k}\_0(mpx+nqy+lsz)} \tag{3}$$

$$\eta(\mathbf{r}) = \sum\_{m,n,l} \eta\_{mnl} e^{-i\mathbf{k}\_0(mpx+nqy+l\mathbf{s}z)} \tag{4}$$

The rules for choosing the most appropriate material parameter and the most appropriate field for the Fourier expansion are governed by various methods of Fourier factorization. In the past, the E method (*η* and the electric displacement D were expanded), H method (*η* and H were expanded), and Ho method (*ε* and E were expanded) were the typical choices.

#### **2.2 Matrix notation**

Now we carry out the transformation of the partial differential equations into matrix equations in order to solve the eigenproblem by linear-algebraic methods. For simplicity we limit ourselves to 1D and 2D photonic crystals, and always choose the direction of propagation in the *xy* plane, so that *∂<sup>z</sup>* = 0. With these restrictions we write

$$\varepsilon(x,y) = \sum\_{m,n=-\infty}^{+\infty} \varepsilon\_{mn} \, e^{-i(mpx+nqy)},\tag{5}$$

$$f(x,y) = \sum\_{m,n=-\infty}^{+\infty} f\_{mn} e^{-i(p\_mx+q\_ny)}.\tag{6}$$

For 1D periodicity we choose the inhomogeneity along the *y* axis. In this case the *α* index is

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 323

+∞ ∑ *n*=−∞

+∞ ∑ *n*=−∞

Although the theorems were derived by Li (1996) for 1D periodic functions, we here summarize them in the matrix formalism independent of the number of dimensions. Let *f* , *h*,

**Theorem 1.** If *ε* and *f* have no concurrent discontinuities, then the Laurent rule applied to

**Theorem 2.** If *ε* and *f* have one or more concurrent discontinuities but *h* is continuous, then

*<sup>f</sup>* <sup>=</sup> <sup>1</sup> *ε*

−<sup>1</sup>

referred to as the inverse rule. We say that the functions *ε* and *f* have complementary

**Theorem 3.** If none of the requirements of the first two theorems are satisfied, then none of the rules can be applied correctly because Equations 27 and 30 are no longer valid at the points of discontinuities, which considerably slows down the convergence. Therefore, we should carefully analyze the continuity of the functions and transform all the partial differential

[ *<sup>f</sup>* ] = 1 *ε* 

[*h*] = 1 *ε* *εn e*

*fn e*

<sup>−</sup>*inqy*, (22)

<sup>−</sup>*iqny*, (23)

[[*ε*]]*nk* = *<sup>ε</sup>n*−*k*, (24) *qnk* = *qnδnk*. (25)

*h* = *εf* , (26)

[*h*] = [[*ε*]][ *f* ] (27)

*h*, (28)

[*h*]. (29)

[ *f* ] (30)

*ε*(*y*) =

*f*(*y*) =

and *ε* be piecewise-continuous functions with the same periodicity related

and let [ *f* ], [*h*], and [[*ε*]] denote their matrices as defined in Section 2.2.

Equation 26 converges uniformly on the whole period and hence

Equation 26 can be transformed into the case of Theorem 1,

Here [[*ε*]] is a Toeplitz matrix (a matrix with constant diagonals).

not necessary because the *n* index is sufficient,

**2.3 Simplified theorems of Fourier factorization**

can be applied with fast convergence.

and hence

discontinuities.

Accordingly, we can state

formulae to the first two cases.

where *f* is a component of the electric field. We will now derive the matrix expressions of the fundamental relations

$$h(\mathbf{x}, y) = \varepsilon(\mathbf{x}, y) f(\mathbf{x}, y), \tag{7}$$

$$\log\_x(\mathbf{x}, y) = \partial\_{\mathbf{x}} f(\mathbf{x}, y), \tag{8}$$

$$
\partial\_y(\mathbf{x}, y) = \partial\_y f(\mathbf{x}, y),
\tag{9}
$$

i.e., the relations of multiplication by a function and applying partial derivatives. Assuming the expansions of the new functions *h* = ∑*m*,*<sup>n</sup> hmn e*−*i*(*pmx*+*qny*), *gx* = ∑*m*,*<sup>n</sup> gx*,*mn e*−*i*(*pmx*+*qny*), and *gy* = ∑*m*,*<sup>n</sup> gy*,*mn e*−*i*(*pmx*+*qny*), we rewrite Equations 7–9 using the convolution rule, and applying the partial derivatives as follows:

$$h\_{mn} = \sum\_{k,l=-\infty}^{+\infty} \varepsilon\_{m-k,n-l} f\_{kl} \tag{10}$$

$$\mathbf{g}\_{\mathbf{x},mn} = -\mathbf{i}p\_{m}\mathbf{f}\_{mn} \tag{11}$$

$$\mathbf{g}\_{y,mn} = -\mathbf{i}q\_nf\_{mn}.\tag{12}$$

Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the summation ∑+*<sup>M</sup> <sup>m</sup>*=−*<sup>M</sup>* <sup>∑</sup>+*<sup>N</sup> <sup>n</sup>*=−*N*, we can renumber all the indices to replace the couple of two sets *m* ∈ {−*M*, −*M* + 1, . . . , *M*} and *n* ∈ {−*N*, −*N* + 1, . . . , *N*} by a single set of indices *α* ∈ {1, 2, . . . , *α*max}, with *α*max = (2*M* + 1)(2*N* + 1), related

$$\alpha(m, n) = m + M + 1 + (n + N)(2M + 1),\tag{13}$$

$$n(\mathfrak{a}) = (\mathfrak{a} - 1)\mathbf{div}(2M + 1) - N\_\prime \tag{14}$$

$$m(a) = (a-1)\mathbf{mod}(2M+1) - M,\tag{15}$$

where "**div**" denotes the operation of integer division and "**mod**" the remainder (the modulo operation). Then we can rewrite Equations 10–12 into the matrix relations

$$[h] = [\![\varepsilon]][f]. \tag{16}$$

$$\mathbf{p}\left[\mathbf{g}\_{\mathbf{x}}\right] = -i\mathbf{p}\left[f\right]\_{\prime} \tag{17}$$

$$\begin{bmatrix} \mathbf{g}\_{\mathbf{y}} \end{bmatrix} = -i\mathbf{q} \begin{bmatrix} f \end{bmatrix} \tag{18}$$

where [ *f* ], [*h*], [*gx*], and [*gy*] are column vectors whose *α*th elements are the Fourier [*m*, *n*] elements of the functions *f* , *h*, *gx*, and *gy*, indexed by *α*(*m*, *n*) defined in Equation 13, and where [[*ε*]], **p**, and **q** are matrices whose elements are defined

$$\|\boldsymbol{\varepsilon}\|\_{a\mathcal{B}} = \varepsilon\_{m(a) - m(\boldsymbol{\beta}), n(a) - n(\boldsymbol{\beta})} \tag{19}$$

$$p\_{a\mathcal{B}} = p\_{m(a)} \delta\_{a\mathcal{B}'} \tag{20}$$

$$
\mathfrak{q}\_{a\mathfrak{B}} = \mathfrak{q}\_{n(a)} \delta\_{a\mathfrak{B}'} \tag{21}
$$

where the indices on the right hand parts are defined by Equations 14 and 15 and where *δαβ* denotes the Kronecker delta. As a summary we can say that the multiplication by a function is in the reciprocal space represented by the matrix [[*ε*]] (in the sense of the limit *α*max → ∞) and that the partial derivatives are represented by the diagonal matrices −*i***p** and −*i***q**.

For 1D periodicity we choose the inhomogeneity along the *y* axis. In this case the *α* index is not necessary because the *n* index is sufficient,

$$\varepsilon(y) = \sum\_{n=-\infty}^{+\infty} \varepsilon\_n e^{-inqy} \,\_{\prime} \tag{22}$$

$$f(y) = \sum\_{n=-\infty}^{+\infty} f\_n e^{-iq\_ny} \,\tag{23}$$

$$\|\varepsilon\|\_{nk} = \varepsilon\_{n-k\prime} \tag{24}$$

$$
\eta\_{nk} = \eta\_n \delta\_{nk}.\tag{25}
$$

Here [[*ε*]] is a Toeplitz matrix (a matrix with constant diagonals).

#### **2.3 Simplified theorems of Fourier factorization**

Although the theorems were derived by Li (1996) for 1D periodic functions, we here summarize them in the matrix formalism independent of the number of dimensions. Let *f* , *h*, and *ε* be piecewise-continuous functions with the same periodicity related

$$h = \mathfrak{e}f,\tag{26}$$

and let [ *f* ], [*h*], and [[*ε*]] denote their matrices as defined in Section 2.2.

**Theorem 1.** If *ε* and *f* have no concurrent discontinuities, then the Laurent rule applied to Equation 26 converges uniformly on the whole period and hence

$$[\mathfrak{h}] = [\mathfrak{e}][f] \tag{27}$$

can be applied with fast convergence.

**Theorem 2.** If *ε* and *f* have one or more concurrent discontinuities but *h* is continuous, then Equation 26 can be transformed into the case of Theorem 1,

$$f = \frac{1}{\varepsilon}h,\tag{28}$$

and hence

4 Will-be-set-by-IN-TECH

where *f* is a component of the electric field. We will now derive the matrix expressions of the

i.e., the relations of multiplication by a function and applying partial derivatives. Assuming the expansions of the new functions *h* = ∑*m*,*<sup>n</sup> hmn e*−*i*(*pmx*+*qny*), *gx* = ∑*m*,*<sup>n</sup> gx*,*mn e*−*i*(*pmx*+*qny*), and *gy* = ∑*m*,*<sup>n</sup> gy*,*mn e*−*i*(*pmx*+*qny*), we rewrite Equations 7–9 using the convolution rule, and

> +∞ ∑ *k*,*l*=−∞

Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the

sets *m* ∈ {−*M*, −*M* + 1, . . . , *M*} and *n* ∈ {−*N*, −*N* + 1, . . . , *N*} by a single set of indices

where "**div**" denotes the operation of integer division and "**mod**" the remainder (the modulo

where [ *f* ], [*h*], [*gx*], and [*gy*] are column vectors whose *α*th elements are the Fourier [*m*, *n*] elements of the functions *f* , *h*, *gx*, and *gy*, indexed by *α*(*m*, *n*) defined in Equation 13, and

where the indices on the right hand parts are defined by Equations 14 and 15 and where *δαβ* denotes the Kronecker delta. As a summary we can say that the multiplication by a function is in the reciprocal space represented by the matrix [[*ε*]] (in the sense of the limit *α*max → ∞) and that the partial derivatives are represented by the diagonal matrices −*i***p** and −*i***q**.

*hmn* =

operation). Then we can rewrite Equations 10–12 into the matrix relations

 *gy* 

*α* ∈ {1, 2, . . . , *α*max}, with *α*max = (2*M* + 1)(2*N* + 1), related

where [[*ε*]], **p**, and **q** are matrices whose elements are defined

*h*(*x*, *y*) = *ε*(*x*, *y*)*f*(*x*, *y*), (7) *gx*(*x*, *y*) = *∂<sup>x</sup> f*(*x*, *y*), (8) *gy*(*x*, *y*) = *∂<sup>y</sup> f*(*x*, *y*), (9)

*gx*,*mn* = −*ipm fmn*, (11) *gy*,*mn* = −*iqn fmn*. (12)

*<sup>n</sup>*=−*N*, we can renumber all the indices to replace the couple of two

*α*(*m*, *n*) = *m* + *M* + 1 + (*n* + *N*)(2*M* + 1), (13) *n*(*α*)=(*α* − 1)**div**(2*M* + 1) − *N*, (14) *m*(*α*)=(*α* − 1)**mod**(2*M* + 1) − *M*, (15)

> [*h*] = [[*ε*]][ *f* ], (16) [*gx*] = −*i***p** [ *f* ], (17)

[[*ε*]]*αβ* <sup>=</sup> *<sup>ε</sup>m*(*α*)−*m*(*β*),*n*(*α*)−*n*(*β*), (19) *pαβ* = *pm*(*α*)*δαβ*, (20) *qαβ* = *qn*(*α*)*δαβ*, (21)

= −*i***q** [ *f* ], (18)

*<sup>ε</sup>m*−*k*,*n*−*<sup>l</sup> fkl*, (10)

fundamental relations

summation ∑+*<sup>M</sup>*

applying the partial derivatives as follows:

*<sup>m</sup>*=−*<sup>M</sup>* <sup>∑</sup>+*<sup>N</sup>*

[ *<sup>f</sup>* ] = 1 *ε* [*h*]. (29)

Accordingly, we can state

$$[h] = \left[\frac{1}{\varepsilon}\right]^{-1}[f] \tag{30}$$

referred to as the inverse rule. We say that the functions *ε* and *f* have complementary discontinuities.

**Theorem 3.** If none of the requirements of the first two theorems are satisfied, then none of the rules can be applied correctly because Equations 27 and 30 are no longer valid at the points of discontinuities, which considerably slows down the convergence. Therefore, we should carefully analyze the continuity of the functions and transform all the partial differential formulae to the first two cases.

where the first and the second row corresponds to the TE and the TM polarization, respectively. The formal labels of the relative permittivity, (L) and (I), means that the multiplication with the corresponding component of the electric field will be (according to the factorization rules) treated by either the Laurent rule (L) or the inverse rule (I). This is due to the fact that *Ez* and *Ex* are continuous functions (because they are tangential to the discontinuities of *ε*), whereas *Ey* is the component normal to the discontinuities and hence discontinuous (while the product *εEz* is continuous). By separating one field in each case, we

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 325

*Ez* = *<sup>ω</sup>*<sup>2</sup>

or, after expanding the permittivity or impermittivity and the fields into the Fourier series, the

For the small-period limit the validity of these rules can be analytically verified by treating the periodic structure as alternation of two homogeneous layers, where we use the boundary conditions for the continuity of the tangential electric and magnetic fields on all interfaces. Now the function *ε* is assumed constant (*ε<sup>I</sup>* or *εII*) within each layer of the thickness *d*<sup>1</sup> or *d*2.

> *aj bj*

Applying the small-period approximation (*Pj*)±<sup>1</sup> <sup>≈</sup> <sup>1</sup> <sup>∓</sup> *ik*0*qjdj* to the problem of propagation

(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

, **Ω** = **C**<sup>2</sup>

Ω<sup>±</sup> = 1 ± *ik*0Λ

Ω<sup>±</sup> = 1 ± *ik*0Λ

[*Ez*] = *<sup>ω</sup>*<sup>2</sup>

[*Hz*] = *<sup>ω</sup>*<sup>2</sup>

<sup>−</sup>*ik*0*qj*(*y*−*yj*) + *bje*

<sup>0</sup>)1/2, with *<sup>ε</sup><sup>j</sup>* <sup>=</sup> *<sup>ε</sup><sup>I</sup>* for odd *<sup>j</sup>* and *<sup>ε</sup><sup>j</sup>* <sup>=</sup> *<sup>ε</sup>II* for even *<sup>j</sup>*. It is coupled with the

, **C***<sup>j</sup>* =

<sup>2</sup> (<sup>1</sup> <sup>−</sup> *<sup>ε</sup>jqj*+1/*εj*+1*qj*) for the TM polarization, and *Pj* <sup>=</sup> *<sup>e</sup>*−*ik*0*qjdj* for both

 **C**1

 *P*<sup>2</sup> 0 0 <sup>1</sup> *P*2

*<sup>ε</sup>*<sup>0</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

 *α<sup>j</sup> β<sup>j</sup> β<sup>j</sup> α<sup>j</sup>*

<sup>2</sup> (<sup>1</sup> <sup>−</sup> *qj*/*qj*<sup>+</sup>1) for the TE polarization, or *<sup>α</sup><sup>j</sup>* <sup>=</sup> <sup>1</sup>

 *P*<sup>1</sup> 0 0 <sup>1</sup> *P*1 <sup>0</sup>, (TE) (40)

<sup>0</sup>*η*0)*ε*0. (TM) (41)

<sup>0</sup>)1/2, assuming the *<sup>e</sup>*−*ik*<sup>0</sup> *<sup>p</sup>*<sup>0</sup> *<sup>x</sup>* factor of all fields.

*ik*0*qj*(*y*−*yj*)

*Hz* = *<sup>ω</sup>*<sup>2</sup>

1 *ε*(L) *∂y*  *<sup>c</sup>*<sup>2</sup> *ε*(L)*Ez*, (TE) (33)

*<sup>c</sup>*<sup>2</sup> *Hz* (TM) (34)

*<sup>c</sup>*<sup>2</sup> [*Ez*], (TE) (35)

*<sup>c</sup>*<sup>2</sup> [*Hz*]. (TM) (36)

], (37)

, (38)

, (39)

<sup>2</sup> (1 +

obtain the wave equations

matrix formulae

where *qj* = (*ε<sup>j</sup>* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

where *α<sup>j</sup>* = <sup>1</sup>

*εjqj*+1/*εj*+1*qj*), *β<sup>j</sup>* = <sup>1</sup>

through the whole period,

− *∂*2 *<sup>x</sup>* + *∂*<sup>2</sup> *y* 

[[*ε*]]−<sup>1</sup> *p*2 <sup>0</sup> <sup>+</sup> **<sup>q</sup>**<sup>2</sup>

**3.3 Consistency with Yeh's theory for the small-period limit**

*<sup>z</sup>* (*x*, *y*) = *e*

 = **C***<sup>j</sup> Pj* 0 0 <sup>1</sup> *Pj*

<sup>0</sup> <sup>+</sup> **<sup>q</sup>**[[*ε*]]−1**<sup>q</sup>**

According to the geometry in Fig. 1, the field in the *j*th layer has the dependence

<sup>−</sup>*ik*<sup>0</sup> *<sup>p</sup>*<sup>0</sup> *<sup>x</sup>*[*aje*

− <sup>1</sup> *ε*(I) *∂*2 *<sup>x</sup>* + *∂<sup>y</sup>*

 1 *ε p*2

*E*(*j*)

field in the next layer by the matrix equation

polarizations. Obviously, *qj* = (*ε<sup>j</sup>* <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

 *a*3 *b*3 = **Ω** *a*1 *b*1 

yields the eigenvalues of the **Ω** operator

 *aj*<sup>+</sup><sup>1</sup> *bj*<sup>+</sup><sup>1</sup>

<sup>2</sup> (<sup>1</sup> <sup>+</sup> *qj*/*qj*<sup>+</sup>1), *<sup>β</sup><sup>j</sup>* <sup>=</sup> <sup>1</sup>

#### **3. One-dimensional photonic crystals**

#### **3.1 Geometry of the problem**

Fig. 1 shows the geometrical configuration of a 1D photonic crystal made as periodic alternation of two different layers with relative permittivities *ε<sup>I</sup>* and *εII*, thicknesses *d*<sup>1</sup> and *d*2, with periodicity Λ along the *y* axis. The relative thickness of the first layer (with respect to the period) is denoted *w*; the relative thickness of the second layer is then 1 − *w*. The coordinate system is chosen to get the uniform problem along the *z* axis. This means that the plane of incidence (plane defined by the vector of periodicity and the wave vector of incidence k = *k*0[*p*0, *q*0, 0], here only hypothetical since the photonic crystal is infinite) coincides with the *xy* plane.

Fig. 1. Geometry of a 1D photonic crystal made as periodic alternation of two layers

Then we distinguish between two polarizations of electromagnetic fields which can be treated independently. The transverse electric (TE) polarization has E perpendicular to the plane of incidence (*Ez*, *Hx*, and *Hy* are nonzero). The transverse magnetic (TM) polarization has H perpendicualar to the plane of incidence (*Hz*, *Ex*, and *Ey* are nonzero).

#### **3.2 Application of Fourier factorization**

Propagation in a 1D periodic medium, whose inhomogeneity along the *y*-axis is described by the relative permittivity function *ε*(*y*), is governed by Maxwell's equations (choosing the coordinate system uniform along the *z*-axis, i.e., *∂<sup>z</sup>* = 0)

$$
\partial\_Y E\_{\overline{z}} = -i\omega\mu\_0 H\_{\overline{x}}, \quad \partial\_{\overline{x}} E\_{\overline{z}} = i\omega\mu\_0 H\_{\overline{y}}, \quad \partial\_{\overline{x}} H\_{\overline{y}} - \partial\_{\overline{y}} H\_{\overline{x}} = i\omega\varepsilon\_0 \varepsilon\_{(\mathcal{L})} E\_{\overline{z}}.\tag{31}
$$

$$
\partial\_y H\_z = i\omega \varepsilon\_0 \varepsilon\_{(\mathcal{L})} E\_{\rm x} \quad \partial\_x H\_z = -i\omega \varepsilon\_0 \varepsilon\_{(\mathcal{I})} E\_{\rm y} \quad \partial\_x E\_{\rm y} - \partial\_y E\_{\rm x} = -i\omega \mu\_0 H\_{\rm z} \tag{32}
$$

6 Will-be-set-by-IN-TECH

Fig. 1 shows the geometrical configuration of a 1D photonic crystal made as periodic alternation of two different layers with relative permittivities *ε<sup>I</sup>* and *εII*, thicknesses *d*<sup>1</sup> and *d*2, with periodicity Λ along the *y* axis. The relative thickness of the first layer (with respect to the period) is denoted *w*; the relative thickness of the second layer is then 1 − *w*. The coordinate system is chosen to get the uniform problem along the *z* axis. This means that the plane of incidence (plane defined by the vector of periodicity and the wave vector of incidence k = *k*0[*p*0, *q*0, 0], here only hypothetical since the photonic crystal is infinite) coincides with

plane of incidence

*y*

**3. One-dimensional photonic crystals**

**3.1 Geometry of the problem**

*I* 

*x*

*d*<sup>1</sup> = *w*-

**3.2 Application of Fourier factorization**

( ) <sup>1</sup> <sup>0</sup> <sup>0</sup> *ik <sup>p</sup> <sup>x</sup> <sup>q</sup> yI a e* <sup>+</sup>

( ) <sup>1</sup> <sup>0</sup> <sup>0</sup> *ik <sup>p</sup> <sup>x</sup> <sup>q</sup> yI b e*

*II* 

0 -

*z*

*I* 

*d*<sup>1</sup> *d*<sup>2</sup>

Then we distinguish between two polarizations of electromagnetic fields which can be treated independently. The transverse electric (TE) polarization has E perpendicular to the plane of incidence (*Ez*, *Hx*, and *Hy* are nonzero). The transverse magnetic (TM) polarization has H

Propagation in a 1D periodic medium, whose inhomogeneity along the *y*-axis is described by the relative permittivity function *ε*(*y*), is governed by Maxwell's equations (choosing the

> *∂yEz* = −*iωμ*0*Hx*, *∂xEz* = *iωμ*0*Hy*, *∂xHy* − *∂yHx* = *iωε*0*ε*(L)*Ez*, (31) *∂yHz* = *iωε*0*ε*(L)*Ex*, *∂xHz* = −*iωε*0*ε*(I)*Ey*, *∂xEy* − *∂yEx* = −*iωμ*0*Hz*, (32)

Fig. 1. Geometry of a 1D photonic crystal made as periodic alternation of two layers

perpendicualar to the plane of incidence (*Hz*, *Ex*, and *Ey* are nonzero).

*II* 

(4)

*d*2 = (1–*w*)-

coordinate system uniform along the *z*-axis, i.e., *∂<sup>z</sup>* = 0)

(1) (2) (3)

( ) <sup>2</sup> <sup>0</sup> <sup>0</sup> *ik <sup>p</sup> <sup>x</sup> <sup>q</sup> <sup>y</sup> II a e* <sup>+</sup>

( ) <sup>2</sup> <sup>0</sup> <sup>0</sup> *ik <sup>p</sup> <sup>x</sup> <sup>q</sup> <sup>y</sup> II b e*

the *xy* plane.

where the first and the second row corresponds to the TE and the TM polarization, respectively. The formal labels of the relative permittivity, (L) and (I), means that the multiplication with the corresponding component of the electric field will be (according to the factorization rules) treated by either the Laurent rule (L) or the inverse rule (I). This is due to the fact that *Ez* and *Ex* are continuous functions (because they are tangential to the discontinuities of *ε*), whereas *Ey* is the component normal to the discontinuities and hence discontinuous (while the product *εEz* is continuous). By separating one field in each case, we obtain the wave equations

$$-\left(\partial\_x^2 + \partial\_y^2\right) E\_z = \frac{\omega^2}{c^2} \varepsilon\_{(\mathcal{L})} E\_{z\prime} \quad \text{(TE)}\tag{33}$$

$$-\left(\frac{1}{\varepsilon\_{(\mathcal{I})}}\partial\_{\mathbf{x}}^{2} + \partial\_{\mathbf{y}}\frac{1}{\varepsilon\_{(\mathcal{L})}}\partial\_{\mathbf{y}}\right)H\_{\mathbf{z}} = \frac{\omega^{2}}{c^{2}}H\_{\mathbf{z}} \quad \text{(TM)}\tag{34}$$

or, after expanding the permittivity or impermittivity and the fields into the Fourier series, the matrix formulae

$$\left[\left[\varepsilon\right]^{-1}\left(p\_0^2 + \mathbf{q}^2\right)\left[E\_z\right] = \frac{\omega^2}{c^2}\left[E\_z\right], \quad \text{(TE)}\tag{35}$$

$$\mathbf{q}\left(\left[\frac{1}{\varepsilon}\right]p\_0^2 + \mathbf{q}\left[\varepsilon\right]^{-1}\mathbf{q}\right)[H\_z] = \frac{\omega^2}{c^2}[H\_z].\tag{36}$$

#### **3.3 Consistency with Yeh's theory for the small-period limit**

For the small-period limit the validity of these rules can be analytically verified by treating the periodic structure as alternation of two homogeneous layers, where we use the boundary conditions for the continuity of the tangential electric and magnetic fields on all interfaces. Now the function *ε* is assumed constant (*ε<sup>I</sup>* or *εII*) within each layer of the thickness *d*<sup>1</sup> or *d*2. According to the geometry in Fig. 1, the field in the *j*th layer has the dependence

$$E\_z^{(j)}(x,y) = e^{-ik\_0 p\_0 x} [a\_j e^{-ik\_0 q\_j (y-y\_j)} + b\_j e^{ik\_0 q\_j (y-y\_j)}],\tag{37}$$

where *qj* = (*ε<sup>j</sup>* <sup>−</sup> *<sup>p</sup>*<sup>2</sup> <sup>0</sup>)1/2, with *<sup>ε</sup><sup>j</sup>* <sup>=</sup> *<sup>ε</sup><sup>I</sup>* for odd *<sup>j</sup>* and *<sup>ε</sup><sup>j</sup>* <sup>=</sup> *<sup>ε</sup>II* for even *<sup>j</sup>*. It is coupled with the field in the next layer by the matrix equation

$$
\begin{bmatrix} a\_{j+1} \\ b\_{j+1} \end{bmatrix} = \mathbf{C}\_{j} \begin{bmatrix} P\_{j} & 0 \\ 0 & \frac{1}{\mathcal{V}\_{j}} \end{bmatrix} \begin{bmatrix} a\_{j} \\ b\_{j} \end{bmatrix}, \quad \mathbf{C}\_{j} = \begin{bmatrix} \alpha\_{j} & \beta\_{j} \\ \beta\_{j} & \alpha\_{j} \end{bmatrix}, \tag{38}
$$

where *α<sup>j</sup>* = <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>+</sup> *qj*/*qj*<sup>+</sup>1), *<sup>β</sup><sup>j</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>−</sup> *qj*/*qj*<sup>+</sup>1) for the TE polarization, or *<sup>α</sup><sup>j</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> (1 + *εjqj*+1/*εj*+1*qj*), *β<sup>j</sup>* = <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>−</sup> *<sup>ε</sup>jqj*+1/*εj*+1*qj*) for the TM polarization, and *Pj* <sup>=</sup> *<sup>e</sup>*−*ik*0*qjdj* for both polarizations. Obviously, *qj* = (*ε<sup>j</sup>* <sup>−</sup> *<sup>p</sup>*<sup>2</sup> <sup>0</sup>)1/2, assuming the *<sup>e</sup>*−*ik*<sup>0</sup> *<sup>p</sup>*<sup>0</sup> *<sup>x</sup>* factor of all fields.

Applying the small-period approximation (*Pj*)±<sup>1</sup> <sup>≈</sup> <sup>1</sup> <sup>∓</sup> *ik*0*qjdj* to the problem of propagation through the whole period,

$$
\begin{bmatrix} a\_3 \\ b\_3 \end{bmatrix} = \boldsymbol{\Omega} \begin{bmatrix} a\_1 \\ b\_1 \end{bmatrix}, \quad \boldsymbol{\Omega} = \mathbf{C}\_2 \begin{bmatrix} P\_2 & 0 \\ 0 & \frac{1}{T\_2} \end{bmatrix} \mathbf{C}\_1 \begin{bmatrix} P\_1 & 0 \\ 0 & \frac{1}{T\_1} \end{bmatrix}, \tag{39}
$$

yields the eigenvalues of the **Ω** operator

$$
\Omega\_{\pm} = 1 \pm ik\_0 \Lambda \sqrt{\varepsilon\_0 - p\_{0'}^2} \quad \text{(TE)} \tag{40}
$$

$$
\Omega\_{\pm} = 1 \pm ik\_0 \Lambda \sqrt{(1 - p\_0^2 \eta\_0) \epsilon\_0}. \quad \text{(TM)} \tag{41}
$$

*y*

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 327

2*R* -

*-*

 

 - 

(a) Square periodicity

denoted Γ, X, M, or K according to the corresponding periodicity.

 

*<sup>y</sup>*

*<sup>x</sup>*

**-**

> *-*

Fig. 4. Two examples of 2D periodic arrangements; below are the first Brillouin zones in the

As an example, in Fig. 3 we have chosen the *xz* plane as the plane of incidence, but we could choose any plane parallel with the *z* axis provided that we want to treat the TE and TM polarizations independently. In general, the plane of incidence is determined by the k vector in the reciprocal space as displayed in Fig. 4, whose particular symmetry points are

In this section we study a 2D photonic crystal composed of infinite cylinders with a circular cross-section with either square [Fig. 4(a)] or hexagonal [Fig. 4(b)] periodicity. For the square symmetry the unit cell has the dimensions Λ*<sup>x</sup>* = Λ*<sup>y</sup>* = Λ, and for the hexagonal symmetry

*x* 

*Hz*, *Ex*, *Ey*), while the TM polarization has E along the *z* axis (nonzero *Ez*, *Hx*, *Hy*), which is




 

*-*

*-* **-**

 

(b) Hexagonal periodicity

Fig. 3. Geometry of a 2D photonic crystal made as bi-periodic alternation of rods or holes

*z*

plane of periodicity

now the simple case.

reciprocal space

*x*

plane of incidence

*y* 

with

$$
\varepsilon\_0 = w \varepsilon\_I + (1 - w) \varepsilon\_{\Pi'} \tag{42}
$$

$$
\eta\_0 = w \eta\_I + (1 - w) \eta\_{II'} \tag{43}
$$

Assuming a Floquet mode with the eigenvalues of **<sup>Ω</sup>** being <sup>Ω</sup><sup>±</sup> <sup>=</sup> *<sup>e</sup>*±*ik*0*q*0<sup>Λ</sup> <sup>≈</sup> <sup>1</sup> <sup>±</sup> *ik*0*q*0Λ, we see that in the small-period limit the periodic structure behaves as a homogeneous anisotropic medium with the *<sup>y</sup>*-component of the normalized wave vector *<sup>q</sup>*<sup>0</sup> = (*ε*<sup>0</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup> <sup>0</sup>)1/2 for the TE polarization and *<sup>q</sup>*<sup>0</sup> = [(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup> <sup>0</sup>*η*0)*ε*0] 1/2 for the TM polarization. These formulae are identical with Equations 35 and 36 if we retain only the 0th element of all matrices. This is a very interesting disclosure that the results obtained by Yeh and coauthors already in 1970s are consistent with the extensive, more general research carried out in 1990.

#### **3.4 Numerical example**

Example of the comparison of applying the correct Fourier factorization rules with applying the opposite ones is shown in Fig. 2 for both polarizations. The normalized eigenfrequency *ω*Λ/2*πc* of the first band is displayed according to *N*; the structure is made as two alternating layers of the equal thicknesses 500 nm (Λ = 1000 nm) and permittivities *ε<sup>I</sup>* = 3 and *εII* = 1. The wave vector is chosen k = (0.5*π*/Λ)[1, 1, 0].

#### **4. Two-dimensional photonic crystals with circular elements**

#### **4.1 Geometry of the problem**

Fig. 3 displays the geometrical arrangement of a 2D photonic crystal made as bi-periodic alternation of rods or holes with a cylindrical cross-section. Instead of a single vector of periodicity we now have the plane of periodicity (determined by two vectors of periodicities defining a unit cell), which here coincides with the *xy* plane. For simplicity we choose the incidence direction in the plane of periodicity, so that we can again distinguish between two independent polarizations. The plane of incidence is now determined along the propagation wave vector k = *k*0[*p*0, *q*0, 0] and perpendicular to the plane of periodicity. The TE polarization has now H along the *z* axis, which is now the more difficult case (with nonzero 8 Will-be-set-by-IN-TECH

Assuming a Floquet mode with the eigenvalues of **<sup>Ω</sup>** being <sup>Ω</sup><sup>±</sup> <sup>=</sup> *<sup>e</sup>*±*ik*0*q*0<sup>Λ</sup> <sup>≈</sup> <sup>1</sup> <sup>±</sup> *ik*0*q*0Λ, we see that in the small-period limit the periodic structure behaves as a homogeneous anisotropic

with Equations 35 and 36 if we retain only the 0th element of all matrices. This is a very interesting disclosure that the results obtained by Yeh and coauthors already in 1970s are

Example of the comparison of applying the correct Fourier factorization rules with applying the opposite ones is shown in Fig. 2 for both polarizations. The normalized eigenfrequency *ω*Λ/2*πc* of the first band is displayed according to *N*; the structure is made as two alternating layers of the equal thicknesses 500 nm (Λ = 1000 nm) and permittivities *ε<sup>I</sup>* = 3 and *εII* = 1.

> 0.2694 0.2695 0.2695 0.2696 0.2696 0.2697 0.2697 0.2698 0.2698 0.2699 0.2699

normalized frequency

Fig. 2. Convergence properties of the correct factorization compared with the opposite one

Fig. 3 displays the geometrical arrangement of a 2D photonic crystal made as bi-periodic alternation of rods or holes with a cylindrical cross-section. Instead of a single vector of periodicity we now have the plane of periodicity (determined by two vectors of periodicities defining a unit cell), which here coincides with the *xy* plane. For simplicity we choose the incidence direction in the plane of periodicity, so that we can again distinguish between two independent polarizations. The plane of incidence is now determined along the propagation wave vector k = *k*0[*p*0, *q*0, 0] and perpendicular to the plane of periodicity. The TE polarization has now H along the *z* axis, which is now the more difficult case (with nonzero

medium with the *<sup>y</sup>*-component of the normalized wave vector *<sup>q</sup>*<sup>0</sup> = (*ε*<sup>0</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

<sup>0</sup>*η*0)*ε*0]

correct factorization opposite factorization

**4. Two-dimensional photonic crystals with circular elements**

The wave vector is chosen k = (0.5*π*/Λ)[1, 1, 0].

0 10 20 30 40 50

*N*

(a) TE polarization

consistent with the extensive, more general research carried out in 1990.

*ε*<sup>0</sup> = *wε<sup>I</sup>* + (1 − *w*)*εII*, (42) *η*<sup>0</sup> = *wη<sup>I</sup>* + (1 − *w*)*ηII*, (43)

1/2 for the TM polarization. These formulae are identical

0 10 20 30 40 50

correct factorization opposite factorization

*N*

(b) TM polarization

<sup>0</sup>)1/2 for the TE

with

polarization and *<sup>q</sup>*<sup>0</sup> = [(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup>

**3.4 Numerical example**

0.246 0.2462 0.2464 0.2466 0.2468 0.247 0.2472

**4.1 Geometry of the problem**

normalized frequency

Fig. 3. Geometry of a 2D photonic crystal made as bi-periodic alternation of rods or holes

*Hz*, *Ex*, *Ey*), while the TM polarization has E along the *z* axis (nonzero *Ez*, *Hx*, *Hy*), which is now the simple case.

Fig. 4. Two examples of 2D periodic arrangements; below are the first Brillouin zones in the reciprocal space

As an example, in Fig. 3 we have chosen the *xz* plane as the plane of incidence, but we could choose any plane parallel with the *z* axis provided that we want to treat the TE and TM polarizations independently. In general, the plane of incidence is determined by the k vector in the reciprocal space as displayed in Fig. 4, whose particular symmetry points are denoted Γ, X, M, or K according to the corresponding periodicity.

In this section we study a 2D photonic crystal composed of infinite cylinders with a circular cross-section with either square [Fig. 4(a)] or hexagonal [Fig. 4(b)] periodicity. For the square symmetry the unit cell has the dimensions Λ*<sup>x</sup>* = Λ*<sup>y</sup>* = Λ, and for the hexagonal symmetry

**4.2.2 Normal vector method (Model B)**

(a) Model A


(c) Model C

Jones matrix transform **F**),

correct rules,

*r* 

Fig. 5. Distribution of the basis polarization vector u for the factorization models

 *Ex Ey* = **F** *Eu Ev* 

According to the factorization theorems, neither the Laurent rule nor the inverse rule is correct for both products in Equations 50 and 51, because both pairs of functions have concurrent discontinuities and both products *D*˜ *<sup>x</sup>* and *D*˜ *<sup>y</sup>* are discontinuous as well. On the other hand, by an appropriate change of the polarization bases at all points (using a space-dependent

we can treat independently the normal (*u*) and tangential (*v*) components of the fields by the


(b) Model B

*R* 

/2

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 329


(d) Model C'

[*D*˜ *<sup>u</sup>*] = [[1/*ε*]]−1[*Eu*], (54) [*D*˜ *<sup>v</sup>*] = [[*ε*]][*Ev*]. (55)

*r* 

*r* 

/2

, (53)

*R* 

it can be chosen as an Λ-by-Λ <sup>√</sup>3 rectangle. The corresponding first Brillouin zones of the reciprocal space are depicted in the bottom part of Fig. 4.

Maxwel's equation are again

$$
\partial\_y E\_z = -i\omega\mu\_0 H\_{z\prime} \quad \partial\_x E\_z = i\omega\mu\_0 H\_{y\prime} \quad \partial\_x H\_y - \partial\_y H\_x = i\omega\varepsilon\_0 \varepsilon E\_{z\prime} \tag{44}
$$

$$
\partial\_y H\_z = \mathrm{i}\omega \varepsilon\_0 \varepsilon E\_x, \quad \partial\_x H\_z = -\mathrm{i}\omega \varepsilon\_0 \varepsilon E\_y, \quad \partial\_x E\_y - \partial\_y E\_x = -\mathrm{i}\omega \mu\_0 H\_z. \tag{45}
$$

However, now we cannot put there the labels (L) or (I) for the factorization type as easily as above, because the discontinuities of the permittivity function are now mixed among the components of E.

Assuming a hypothetical anisotropy of the relative permittivity function, we define a scaled electrical displacement D˜ ,

$$
\begin{bmatrix} \tilde{D}\_{\mathbf{X}} \\ \tilde{D}\_{\mathbf{y}} \end{bmatrix} = \boldsymbol{\varepsilon} \begin{bmatrix} E\_{\mathbf{x}} \\ E\_{\mathbf{y}} \end{bmatrix} = \begin{bmatrix} \varepsilon\_{\mathbf{xx}} \ \varepsilon\_{\mathbf{xy}} \\ \varepsilon\_{\mathbf{yx}} \ \varepsilon\_{\mathbf{yy}} \end{bmatrix} \begin{bmatrix} E\_{\mathbf{x}} \\ E\_{\mathbf{y}} \end{bmatrix} \tag{15}
$$

$$
\tilde{D}\_z = \varepsilon\_{zz} E\_z = \varepsilon\_{(\mathcal{L})} E\_{\mathcal{Z}\prime} \quad \text{(TM)} \tag{47}
$$

where *εzz* is obviously the only component of the permittivity tensor for which we can use the Laurent rule.

Defining also a 2-by-2 matrix of electrical impermittivity η = ε−<sup>1</sup> helps in the formulation of the wave equations

$$(-\partial\_{\mathbf{y}}\eta\_{\mathbf{xx}}\partial\_{\mathbf{y}} + \partial\_{\mathbf{x}}\eta\_{\mathbf{yx}}\partial\_{\mathbf{y}} + \partial\_{\mathbf{y}}\eta\_{\mathbf{xy}}\partial\_{\mathbf{x}} - \partial\_{\mathbf{x}}\eta\_{\mathbf{yy}}\partial\_{\mathbf{x}})H\_{\mathbf{z}} = \frac{\omega^2}{c^2}H\_{\mathbf{z}} \tag{48}$$

$$-\left(\partial\_x^2 + \partial\_y^2\right) E\_z = \frac{\omega^2}{c^2} \varepsilon\_{(\mathcal{L})} E\_{z\prime} \quad \text{(TM)}\tag{49}$$

where *ηjk* are the components of the electrical impermittivity. For the simplicity of the TM polarization case we below focus our attention only to the TE polarization.

#### **4.2 Methods of Fourier factorization**

In this section we compare several models corresponding to different factorization approaches.

#### **4.2.1 Elementary (Cartesian) method (Model A)**

First, Model A assumes the solution in the basis of the ˆx and ˆy polarizations uniform within the periodic cell, where in accordance with Ho et al. (1990) we choose the Laurent rules

$$\begin{bmatrix} \tilde{D}\_{\mathbf{x}} \end{bmatrix} = \begin{bmatrix} \varepsilon \mathbf{E}\_{\mathbf{x}} \end{bmatrix} = \begin{bmatrix} \mathbb{E} \end{bmatrix} \begin{bmatrix} \mathbf{E}\_{\mathbf{x}} \end{bmatrix} \tag{50}$$

$$[\tilde{D}\_{\mathcal{Y}}] = [\varepsilon E\_{\mathcal{Y}}] = [\mathbb{E}][E\_{\mathcal{Y}}].\tag{51}$$

The components of the electric impermittivity in Equation 48 then becomes [[*ε*]]−<sup>1</sup> for the cases of *ηxx*, *ηyy*, and zero for the cases of *ηxy*, *ηyx*, or

$$\begin{aligned} \, \|\boldsymbol{\eta}\|\!\_{\mathcal{A}} &= \begin{bmatrix} \|\boldsymbol{\varepsilon}\|\!\_{\boldsymbol{\varepsilon}}^{-1} & \begin{bmatrix} \mathbf{0} \end{bmatrix} \\ \|\!\_{\boldsymbol{\varepsilon}}\|\!\_{\boldsymbol{\varepsilon}}^{-1} & \begin{bmatrix} \boldsymbol{\varepsilon} \end{bmatrix} \end{bmatrix} . \end{aligned} \tag{52}$$

For illustration we show the distribution of the first basis polarization vector (identical with the constant vector ˆx) in Fig. 5(a), where the black circle denotes the element boundary (the permittivity discontinuity).

#### **4.2.2 Normal vector method (Model B)**

10 Will-be-set-by-IN-TECH

However, now we cannot put there the labels (L) or (I) for the factorization type as easily as above, because the discontinuities of the permittivity function are now mixed among the

Assuming a hypothetical anisotropy of the relative permittivity function, we define a scaled

where *εzz* is obviously the only component of the permittivity tensor for which we can use the

Defining also a 2-by-2 matrix of electrical impermittivity η = ε−<sup>1</sup> helps in the formulation of

*Ez* = *<sup>ω</sup>*<sup>2</sup>

where *ηjk* are the components of the electrical impermittivity. For the simplicity of the TM

In this section we compare several models corresponding to different factorization

First, Model A assumes the solution in the basis of the ˆx and ˆy polarizations uniform within the periodic cell, where in accordance with Ho et al. (1990) we choose the Laurent rules

The components of the electric impermittivity in Equation 48 then becomes [[*ε*]]−<sup>1</sup> for the cases

For illustration we show the distribution of the first basis polarization vector (identical with the constant vector ˆx) in Fig. 5(a), where the black circle denotes the element boundary (the

[[*ε*]]−<sup>1</sup> [[0]] [[0]] [[*ε*]]−<sup>1</sup> [[η]]<sup>A</sup> =

(−*∂yηxx∂<sup>y</sup>* <sup>+</sup> *<sup>∂</sup>xηyx∂<sup>y</sup>* <sup>+</sup> *<sup>∂</sup>yηxy∂<sup>x</sup>* <sup>−</sup> *<sup>∂</sup>xηyy∂x*)*Hz* <sup>=</sup> *<sup>ω</sup>*<sup>2</sup>

polarization case we below focus our attention only to the TE polarization.

*∂yEz* = −*iωμ*0*Hz*, *∂xEz* = *iωμ*0*Hy*, *∂xHy* − *∂yHx* = *iωε*0*εEz*, (44) *∂yHz* = *iωε*0*εEx*, *∂xHz* = −*iωε*0*εEy*, *∂xEy* − *∂yEx* = −*iωμ*0*Hz*. (45)

> *Ex Ey*

*<sup>D</sup>*˜ *<sup>z</sup>* = *<sup>ε</sup>zzEz* = *<sup>ε</sup>*(L)*Ez*, (TM) (47)

[*D*˜ *<sup>x</sup>*]=[*εEx*] = [[*ε*]][*Ex*], (50) [*D*˜ *<sup>y</sup>*]=[*εEy*] = [[*ε*]][*Ey*]. (51)

. (52)

, (TE) (46)

*<sup>c</sup>*<sup>2</sup> *Hz*, (TE) (48)

*<sup>c</sup>*<sup>2</sup> *ε*(L)*Ez*, (TM) (49)

<sup>√</sup>3 rectangle. The corresponding first Brillouin zones of the

it can be chosen as an Λ-by-Λ

Maxwel's equation are again

components of E.

Laurent rule.

approaches.

the wave equations

**4.2 Methods of Fourier factorization**

**4.2.1 Elementary (Cartesian) method (Model A)**

of *ηxx*, *ηyy*, and zero for the cases of *ηxy*, *ηyx*, or

permittivity discontinuity).

electrical displacement D˜ ,

reciprocal space are depicted in the bottom part of Fig. 4.

 *D*˜ *x D*˜ *y* = ε *Ex Ey* = *εxx εxy εyx εyy*

> − *∂*2 *<sup>x</sup>* + *∂*<sup>2</sup> *y*

Fig. 5. Distribution of the basis polarization vector u for the factorization models

According to the factorization theorems, neither the Laurent rule nor the inverse rule is correct for both products in Equations 50 and 51, because both pairs of functions have concurrent discontinuities and both products *D*˜ *<sup>x</sup>* and *D*˜ *<sup>y</sup>* are discontinuous as well. On the other hand, by an appropriate change of the polarization bases at all points (using a space-dependent Jones matrix transform **F**),

$$
\begin{bmatrix} E\_{\chi} \\ E\_{\mathcal{Y}} \end{bmatrix} = \mathbf{F} \begin{bmatrix} E\_{\mathcal{U}} \\ E\_{\mathcal{U}} \end{bmatrix} \tag{53}
$$

we can treat independently the normal (*u*) and tangential (*v*) components of the fields by the correct rules,

$$[\tilde{D}\_{\mathfrak{u}}] = [\mathbbm{1}/\varepsilon]^{-1}[E\_{\mathfrak{u}}]\_{\prime} \tag{54}$$

$$\mathbb{E}\left[\mathcal{D}\_{\upsilon}\right] = \left\llbracket \varepsilon \right\rrbracket \left[E\_{\upsilon}\right]. \tag{55}$$

(still orthogonal). By means of rotation *θ* and ellipticity *�* we define the first basis vector

1 + cos *<sup>π</sup><sup>r</sup> R*

*<sup>D</sup>*(*φ*) = <sup>Λ</sup>/2

u = lim *r*→*D*(*φ*)

is the distance from the cell's center to its edge. In Equation 62 the Jones vector on the right represents a polarization ellipse (with ellipticity *�*) oriented along the *x* coordinate, the matrix in the middle rotates this polarization by the azimuth *θ*, and the factor *ei<sup>θ</sup>* preserves the continuity of the phase at the center and along the boundaries of the cell. This continuity

The distribution of the basis polarization vector u within the periodic cell is shown in Fig. 5(c). Here the azimuth of the polarization ellipse is constant along the lines coming from the cell's center, which is similar to Model B. However, the ellipticity is now zero (corresponding to linear polarization) only on the boundaries of the circular element, has the maximum value (*π*/4 for circular polarization) at the cell's center and along its boundaries, and continuously varies (with a smooth sine dependence) in the intermediate ranges. Thus we obtain a smooth and completely continuous matrix function **F**(*x*, *y*), which is analogously used to calculate the

[[*ζζ*∗]], [[ <sup>1</sup>

[[*ξ*∗*ζ*]], [[ <sup>1</sup>

In the case of the hexagonal periodicity we define u and the other periodic quantities inside one hexagon (half the area of the rectangular unit cell) where we can use formally the same

which is now the distance from the hexagon's center to its edge. Here Λ is the hexagon's

*<sup>D</sup>*(*φ*) = <sup>Λ</sup>/2 max *n*=0,...,5 cos

*<sup>ε</sup>* ]][[*ξζ*∗]] <sup>−</sup> [[*ε*]]−<sup>1</sup>

*<sup>ε</sup>* ]][[*ζζ*∗]] + [[*ε*]]−<sup>1</sup>

*<sup>φ</sup>* <sup>−</sup> *<sup>n</sup><sup>π</sup>* 3

[[*ξζ*∗]]

[[*ξξ*∗]]

, (68)

. (67)

<sup>u</sup> <sup>=</sup> <sup>1</sup> √2 1 *i* 

1 + cos *<sup>π</sup>*[*r*+*D*(*φ*)−2*R*] *D*(*φ*)−*R*

cos *θ* − sin *θ* sin *θ* cos *θ*

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 331

 cos *� i* sin *�*

*θ*(*r*, *φ*) = *φ*, (63)

(*<sup>r</sup>* <sup>≤</sup> *<sup>R</sup>*)

max(<sup>|</sup> cos *<sup>φ</sup>*|, <sup>|</sup> sin *<sup>φ</sup>*|) (65)

, (62)

(*<sup>r</sup>* <sup>&</sup>gt; *<sup>R</sup>*). (64)

, (66)

u = *e iθ* 

> 8

*π* 8 

lim *r*→0

which is the vector of left circular polarization (independent of *φ*).

*<sup>ε</sup>* ]][[*ξξ*∗]] + [[*ε*]]−<sup>1</sup>

*<sup>ε</sup>* ]][[*ξ*∗*ζ*]] <sup>−</sup> [[*ε*]]−<sup>1</sup>

*�*(*r*, *<sup>φ</sup>*) = *<sup>π</sup>*

Here *R* denotes the radius of the circular element and

can be easily checked by evaluating the limits

impermittivity in the reciprocal space

equations as above, except for

[[η]]<sup>C</sup> =

 [[ 1

[[ 1

shortest width (equal to the width Λ*<sup>x</sup>* of the rectangular cell).

where

The field components *Eu*, *D*˜ *<sup>u</sup>* are normal to the discontinuities of the relative permittivity function, while *Ev*, *D*˜ *<sup>v</sup>* are tangential. The factorization rules used in Equations 54 and 55 are justified simply because *Ev* and *D*˜ *<sup>u</sup>* are continuous.

A suitable distribution of the matrix **F** within the periodic cell can obviously be the rotation

$$\mathbf{F} = \begin{bmatrix} \cos \phi - \sin \phi \\ \sin \phi & \cos \phi \end{bmatrix} \text{ \textit{\textsuperscript{0}}} \tag{56}$$

where the polar angle *φ*(*x*, *y*) is in the first cell distributed according to the polar coordinates *rei<sup>φ</sup>* = *x* + *iy*, and then periodically repeated over the entire 2D space. This enables defining the matrices [[*c*]], [[*s*]] from the corresponding 2D-periodic functions *c* = cos *φ*, *s* = sin *φ*.

Let u and v be the two columns of the matrix **F**, both being mutually orthogonal basis vectors of linear polarization. From the above definitions we see that u is a polarization vector normal to the structure discontinuities, whereas v is tangential. In Fig. 5(b) we show the distribution of u within the periodic cell. The basis polarization vectors are constant along the lines of the constant azimuth (*φ* = const) and rotate as *φ* increases. It is obvious that the matrix function **F**(*x*, *y*) has no discontinuities concurrent with the electric field, so that we can use both Laurent and inverse rules for the transformation of polarization, e.g.,

$$
\begin{bmatrix} \begin{bmatrix} E\_x \end{bmatrix} \\ \begin{bmatrix} E\_y \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} E\_{\mathcal{U}} \end{bmatrix} \\ \begin{bmatrix} E\_{\mathcal{U}} \end{bmatrix} \end{bmatrix} \tag{57}
$$

$$\begin{bmatrix} \mathbf{F} \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} \mathbf{c} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} -s \end{bmatrix} \\ \begin{bmatrix} \mathbf{c} \end{bmatrix} \end{bmatrix} \end{bmatrix} \tag{58}$$

Combining Equations 54, 55, and 57 yields

$$
\begin{bmatrix} \begin{bmatrix} E\_X \\ E\_Y \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \\ \begin{bmatrix} \mathbf{0} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{F}^{-1} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \tilde{D}\_x \\ \tilde{D}\_y \end{bmatrix} \end{bmatrix} \tag{59}
$$

from where we derive the electric impermittivity in the reciprocal space (corresponding to Model B)

$$\begin{aligned} \begin{bmatrix} \mathfrak{h} \end{bmatrix} \mathbf{B} &= \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \\ \begin{bmatrix} \left[\mathbf{0}\right] \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix}^{-1} \end{bmatrix} \begin{bmatrix} \mathbf{F}^{-1} \end{bmatrix} \\ &= \begin{bmatrix} \left[\frac{1}{\varepsilon}\right] \begin{bmatrix} \varepsilon^{2} \end{bmatrix} \begin{bmatrix} \varepsilon^{2} \end{bmatrix} + \begin{bmatrix} \varepsilon \end{bmatrix}^{-1} \begin{bmatrix} \varepsilon^{2} \end{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \begin{bmatrix} \left[\varepsilon s \right] - \left[\varepsilon \right]^{-1} \left[\varepsilon s\right] \\ \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \begin{bmatrix} \left[\varepsilon s \right] - \left[\varepsilon \right]^{-1} \left[\varepsilon s\right] \end{bmatrix} \end{aligned} \tag{60} \end{aligned} \tag{61}$$

whose components are immediately applicable to Equation 48.

#### **4.2.3 Method with elliptical polarization bases (Model C)**

The above approach (Model B) only deals with linear polarizations and thus suffers from the fact that the matrix function **F**(*x*, *y*) has a singularity at the central point of the periodic cell and other discontinuities along the cell boundaries. This slows down the convergence of the numerical implementation, as will be evidenced below. On the other hand, we can make **F** continuous by using complex functions *ξ* and *ζ* or, in other words, by defining u, v as complex vectors corresponding to generally elliptic polarizations,

$$\begin{aligned} \mathbf{u} = \begin{bmatrix} \mathfrak{f} \\ \zeta \end{bmatrix}, \quad \mathbf{v} = \begin{bmatrix} -\zeta^\* \\ \tilde{\xi}^\* \end{bmatrix} \end{aligned} \tag{61}$$

(still orthogonal). By means of rotation *θ* and ellipticity *�* we define the first basis vector

$$\mathbf{u} = e^{i\theta} \begin{bmatrix} \cos\theta - \sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} \cos\epsilon \\ i\sin\epsilon \end{bmatrix} \tag{62}$$

where

12 Will-be-set-by-IN-TECH

The field components *Eu*, *D*˜ *<sup>u</sup>* are normal to the discontinuities of the relative permittivity function, while *Ev*, *D*˜ *<sup>v</sup>* are tangential. The factorization rules used in Equations 54 and 55 are

A suitable distribution of the matrix **F** within the periodic cell can obviously be the rotation

cos *φ* − sin *φ* sin *φ* cos *φ*

where the polar angle *φ*(*x*, *y*) is in the first cell distributed according to the polar coordinates *rei<sup>φ</sup>* = *x* + *iy*, and then periodically repeated over the entire 2D space. This enables defining the matrices [[*c*]], [[*s*]] from the corresponding 2D-periodic functions *c* = cos *φ*, *s* = sin *φ*.

Let u and v be the two columns of the matrix **F**, both being mutually orthogonal basis vectors of linear polarization. From the above definitions we see that u is a polarization vector normal to the structure discontinuities, whereas v is tangential. In Fig. 5(b) we show the distribution of u within the periodic cell. The basis polarization vectors are constant along the lines of the constant azimuth (*φ* = const) and rotate as *φ* increases. It is obvious that the matrix function **F**(*x*, *y*) has no discontinuities concurrent with the electric field, so that we can use both Laurent

> [[*c*]] [[−*s*]] [[*s*]] [[*c*]]

from where we derive the electric impermittivity in the reciprocal space (corresponding to

The above approach (Model B) only deals with linear polarizations and thus suffers from the fact that the matrix function **F**(*x*, *y*) has a singularity at the central point of the periodic cell and other discontinuities along the cell boundaries. This slows down the convergence of the numerical implementation, as will be evidenced below. On the other hand, we can make **F** continuous by using complex functions *ξ* and *ζ* or, in other words, by defining u, v as complex

, v =

 −*ζ*<sup>∗</sup> *ξ*∗ 

 [[**F**−1]]

*<sup>ε</sup>* ]][[*c*2]] + [[*ε*]]−1[[*s*2]], [[ <sup>1</sup>

*<sup>ε</sup>* ]][[*cs*]] <sup>−</sup> [[*ε*]]−1[[*cs*]], [[ <sup>1</sup>

 [[**F**−1]] [*D*˜ *<sup>x</sup>*] [*D*˜ *<sup>y</sup>*] 

*<sup>ε</sup>* ]][[*cs*]] <sup>−</sup> [[*ε*]]−1[[*cs*]]

*<sup>ε</sup>* ]][[*s*2]] + [[*ε*]]−1[[*c*2]]

, (56)

, (57)

. (58)

, (59)

, (60)

(61)

**F** = 

justified simply because *Ev* and *D*˜ *<sup>u</sup>* are continuous.

and inverse rules for the transformation of polarization, e.g.,

Combining Equations 54, 55, and 57 yields

Model B)

 [*Ex*] [*Ey*] = [[**F**]] [[ 1 *<sup>ε</sup>* ]] [[0]] [[0]] [[*ε*]]−<sup>1</sup>

[[η]]<sup>B</sup> = [[**F**]]

= [[ 1

 [[ 1 *<sup>ε</sup>* ]] [[0]] [[0]] [[*ε*]]−<sup>1</sup>

[[ 1

whose components are immediately applicable to Equation 48.

**4.2.3 Method with elliptical polarization bases (Model C)**

vectors corresponding to generally elliptic polarizations,

u = *ξ ζ* 

 [*Ex*] [*Ey*] = [[**F**]] [*Eu*] [*Ev*] 

[[**F**]] = 

$$
\theta(r,\phi) = \phi,\tag{63}
$$

$$\epsilon(r,\phi) = \begin{cases} \frac{\pi}{8} \left( 1 + \cos\frac{\pi r}{\mathbb{R}} \right) & (r \le R) \\\frac{\pi}{8} \left\{ 1 + \cos\frac{\pi \left[ r + D(\phi) - 2R \right]}{D(\phi) - R} \right\} \text{ (\$r > R)}. \end{cases} \tag{64}$$

Here *R* denotes the radius of the circular element and

$$D(\phi) = \frac{\Lambda/2}{\max(|\cos \phi|, |\sin \phi|)}\tag{65}$$

is the distance from the cell's center to its edge. In Equation 62 the Jones vector on the right represents a polarization ellipse (with ellipticity *�*) oriented along the *x* coordinate, the matrix in the middle rotates this polarization by the azimuth *θ*, and the factor *ei<sup>θ</sup>* preserves the continuity of the phase at the center and along the boundaries of the cell. This continuity can be easily checked by evaluating the limits

$$\lim\_{r \to 0} u = \lim\_{r \to D(\phi)} u = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ i \end{bmatrix} . \tag{66}$$

which is the vector of left circular polarization (independent of *φ*).

The distribution of the basis polarization vector u within the periodic cell is shown in Fig. 5(c). Here the azimuth of the polarization ellipse is constant along the lines coming from the cell's center, which is similar to Model B. However, the ellipticity is now zero (corresponding to linear polarization) only on the boundaries of the circular element, has the maximum value (*π*/4 for circular polarization) at the cell's center and along its boundaries, and continuously varies (with a smooth sine dependence) in the intermediate ranges. Thus we obtain a smooth and completely continuous matrix function **F**(*x*, *y*), which is analogously used to calculate the impermittivity in the reciprocal space

$$\begin{split} \left[\boldsymbol{\mathfrak{n}}\right]\_{\mathsf{C}} = \begin{bmatrix} \left[\frac{1}{\varepsilon}\right] \left[\left[\boldsymbol{\xi}\boldsymbol{\xi}^{\ast}\right] + \left[\boldsymbol{\varepsilon}\right]^{-1} \left[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\ast}\right], \left[\frac{1}{\varepsilon}\right] \left[\left[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\ast}\right] - \left[\boldsymbol{\varepsilon}\right]^{-1} \left[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\ast}\right] \right] \\ \left[\frac{1}{\varepsilon}\right] \left[\left[\boldsymbol{\zeta}^{\ast}\boldsymbol{\zeta}\right] - \left[\boldsymbol{\varepsilon}\right]^{-1} \left[\boldsymbol{\zeta}^{\ast}\boldsymbol{\zeta}\right], \left[\frac{1}{\varepsilon}\right] \left[\left[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\ast}\right] + \left[\boldsymbol{\varepsilon}\right]^{-1} \left[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\ast}\right] \right] \end{bmatrix}. \end{split} \tag{67}$$

In the case of the hexagonal periodicity we define u and the other periodic quantities inside one hexagon (half the area of the rectangular unit cell) where we can use formally the same equations as above, except for

$$D(\phi) = \frac{\Lambda/2}{\max\_{n=0,\ldots,5} \left[ \cos \left( \phi - \frac{n\pi}{3} \right) \right]},\tag{68}$$

which is now the distance from the hexagon's center to its edge. Here Λ is the hexagon's shortest width (equal to the width Λ*<sup>x</sup>* of the rectangular cell).

Γ

chosen eigenmodes

Sample S2).

norm. frequency

 

(a) Sample S1

Fig. 6. Photonic band structures of two samples

**-**

**-**

> **--**

S1-Γ2 S1-Γ3

H-M1

S2-X3

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 333

Fig. 7. Amplitude distributions of the scaled magnetic field |*Hz*| within the unit cell of four

and the number corresponds to the band. The amplitude distributions of the modes are displayed in Fig. 7, and the corresponding convergence properties of the eigenfrequencies calculated by the above described models are shown in Fig. 8 (Model C' is only compared for

The result of a careful comparison of the numerical efficiencies of all the factorization models can be summarized as follows. Models B, C, and C' always converge considerably faster than Model A. Model C converges faster (having usually one order higher precision) than Model B with two exceptions. The first exception is the case such as in Fig. 7 (S1-Γ2) and Fig. 8(a), where the discontinuities of the polarization transformation matrix **F** coincide with a nearly zero amplitude of the field, so that the discontinuities do not manifest themselves. The second exception is the case such as Fig. 7 (S2-X3) and Fig. 8(d), where the elements are densely arranged which causes rapid variations of the ellipticity between two adjacent elements (which are very close to each other); this requires more Fourier components than the weak discontinuity of the linear polarization u in Model B. The problem is solved in Model C',

which obviously converges fastest among all the four models applied to Sample S2.

Γ **M K** Γ

(b) Sample H

**<sup>Γ</sup>** *kx*

*ky* **M K**



 

ΓΓ

#### **4.2.4 Modified method for densely arranged elements (Model C')**

To analyze a more complicated situation, we consider a photonic crystal with square periodicity where circular elements are densely arranged near each other, i.e., where the radius *R* is almost the half width Λ/2 of the periodic cell. Then the convergence properties of **F** becomes worse, which affects all the derived quantities. For this reason we again redefine the polarization distribution. For the modified Model C' we define u to be still same inside the circle (*r* < *R*), but different outside. Assuming the rotation and ellipticity along the boundary of the square cell

$$\theta\_{\mathsf{b}}(\phi) = \theta\left(D(\phi), \phi\right) = \frac{\pi}{2} \mathbf{round}\left(\phi/\frac{\pi}{2}\right),\tag{69}$$

$$\epsilon\_{\mathbf{b}}(\phi) = \epsilon \left( D(\phi) , \phi \right) = \frac{\pi}{8} (1 - \cos 4\phi) \tag{70}$$

(where "**round**" denotes rounding towards the nearest integer), we define the rotation and ellipticity outside the circle (*r* > *R*) as

$$\theta(r,\phi) = \frac{1}{2}\left\{\theta\_{\mathsf{b}}(\phi) + \phi + [\theta\_{\mathsf{b}}(\phi) - \phi]\cos\frac{\pi[r + D(\phi) - 2\mathsf{R}]}{D(\phi) - \mathsf{R}}\right\},\tag{71}$$

$$\epsilon(r,\phi) = \frac{\epsilon\_{\rm b}(\phi)}{2} \left\{ 1 + \cos \frac{\pi [r + D(\phi) - 2\mathbb{R}]}{D(\phi) - \mathbb{R}} \right\}. \tag{72}$$

Assuming otherwise the same Equations 59, 61, 62, and 65, we obtain for [[η]]C� formally the same matrix as in Equation 67, except that the functions *ξ* and *ζ* are now derived from different azimuth and ellipticity distributions of u. Note that u is again continuous along the cell's boundaries; to evaluate its precise limits [when *x* → ±*D*(0), *y* = const or *y* → ±*D*(*φ*/2), *x* = const] would now be more complicated. The distribution of the basis polarization vector u within the periodic cell is depicted in Fig. 5(d) together with dimensions.

#### **4.3 Numerical examples**

We examine the numerical performances of all factorization models presented in Section 4.2 on three samples of 2D photonic crystals, for which we calculate the eigenfrequencies *ωκ* (where the band number *κ* = 1 stands for the lowest eigenfrequency, *κ* = 2 for the second lowest, etc.) and the corresponding eigenvectors [*Hz*]*<sup>κ</sup>* of Equation 48. All convergences will be presented according to the maximum Fourier harmonics retained inside the periodic medium, which will be kept same for the *x* and *y* directions (*M* = *N*).

First, Sample S1 is a square array of cylindrical rods of the circular cross-section with the diameter 2*R* = 500 nm, square period Λ = 1000 nm, relative permittivity of the rods *ε*<sup>1</sup> = 9, and relative permittivity of the surrounding medium corresponding to vacuum (*ε*<sup>2</sup> = 1). Its dispersion relation is displayed in Fig. 6(a). Similarly, for Sample S2 we assume exactly the same parameters except the diameter of the rods, now being 2*R* = 900 nm. This corresponds to densely arranged elements (the distance between two adjacent rods is only 100 nm). Finally, for Sample H we consider a hexagonal array of cylindrical holes of the circular cross-section with the diameter 2*R* = 600 nm, hexagonal periodicity Λ = 1000 nm (corresponding to the rectangular cell of the dimensions <sup>Λ</sup>*<sup>x</sup>* <sup>=</sup> <sup>1</sup> *<sup>μ</sup>*m, <sup>Λ</sup>*<sup>y</sup>* <sup>=</sup> <sup>√</sup><sup>3</sup> *<sup>μ</sup>*m), relative permittivity of the holes corresponding to vacuum (*ε*<sup>1</sup> = 1), and relative permittivity of the substrate medium (surrounding holes) *ε*<sup>2</sup> = 12. The dispersion relation of Sample H is displayed in Fig. 6(b).

For our analysis we choose the eigenmodes Γ2 and Γ3 of Sample S1, the eigenmode X3 of Sample S2, and the eigenmode M1 of Sample H, where the letter denotes a point of symmetry

Fig. 6. Photonic band structures of two samples

14 Will-be-set-by-IN-TECH

To analyze a more complicated situation, we consider a photonic crystal with square periodicity where circular elements are densely arranged near each other, i.e., where the radius *R* is almost the half width Λ/2 of the periodic cell. Then the convergence properties of **F** becomes worse, which affects all the derived quantities. For this reason we again redefine the polarization distribution. For the modified Model C' we define u to be still same inside the circle (*r* < *R*), but different outside. Assuming the rotation and ellipticity along the boundary

(where "**round**" denotes rounding towards the nearest integer), we define the rotation and

1 + cos *<sup>π</sup>*[*r*+*D*(*φ*)−2*R*] *D*(*φ*)−*R*

Assuming otherwise the same Equations 59, 61, 62, and 65, we obtain for [[η]]C� formally the same matrix as in Equation 67, except that the functions *ξ* and *ζ* are now derived from different azimuth and ellipticity distributions of u. Note that u is again continuous along the cell's boundaries; to evaluate its precise limits [when *x* → ±*D*(0), *y* = const or *y* → ±*D*(*φ*/2), *x* = const] would now be more complicated. The distribution of the basis polarization vector

We examine the numerical performances of all factorization models presented in Section 4.2 on three samples of 2D photonic crystals, for which we calculate the eigenfrequencies *ωκ* (where the band number *κ* = 1 stands for the lowest eigenfrequency, *κ* = 2 for the second lowest, etc.) and the corresponding eigenvectors [*Hz*]*<sup>κ</sup>* of Equation 48. All convergences will be presented according to the maximum Fourier harmonics retained inside the periodic medium, which

First, Sample S1 is a square array of cylindrical rods of the circular cross-section with the diameter 2*R* = 500 nm, square period Λ = 1000 nm, relative permittivity of the rods *ε*<sup>1</sup> = 9, and relative permittivity of the surrounding medium corresponding to vacuum (*ε*<sup>2</sup> = 1). Its dispersion relation is displayed in Fig. 6(a). Similarly, for Sample S2 we assume exactly the same parameters except the diameter of the rods, now being 2*R* = 900 nm. This corresponds to densely arranged elements (the distance between two adjacent rods is only 100 nm). Finally, for Sample H we consider a hexagonal array of cylindrical holes of the circular cross-section with the diameter 2*R* = 600 nm, hexagonal periodicity Λ = 1000 nm (corresponding to the rectangular cell of the dimensions <sup>Λ</sup>*<sup>x</sup>* <sup>=</sup> <sup>1</sup> *<sup>μ</sup>*m, <sup>Λ</sup>*<sup>y</sup>* <sup>=</sup> <sup>√</sup><sup>3</sup> *<sup>μ</sup>*m), relative permittivity of the holes corresponding to vacuum (*ε*<sup>1</sup> = 1), and relative permittivity of the substrate medium (surrounding holes) *ε*<sup>2</sup> = 12. The dispersion relation of Sample H is displayed in Fig. 6(b). For our analysis we choose the eigenmodes Γ2 and Γ3 of Sample S1, the eigenmode X3 of Sample S2, and the eigenmode M1 of Sample H, where the letter denotes a point of symmetry

<sup>2</sup> **round**

*<sup>θ</sup>*b(*φ*) + *<sup>φ</sup>* + [*θ*b(*φ*) <sup>−</sup> *<sup>φ</sup>*] cos *<sup>π</sup>*[*r*+*D*(*φ*)−2*R*]

*φ*/ *<sup>π</sup>* 2 

*D*(*φ*)−*R*

, (69)

, (71)

<sup>8</sup> (1 − cos 4*φ*) (70)

. (72)

**4.2.4 Modified method for densely arranged elements (Model C')**

*θ*b(*φ*) = *θ* (*D*(*φ*), *φ*) = *<sup>π</sup>*

u within the periodic cell is depicted in Fig. 5(d) together with dimensions.

*�*b(*φ*) = *�* (*D*(*φ*), *φ*) = *<sup>π</sup>*

of the square cell

ellipticity outside the circle (*r* > *R*) as

**4.3 Numerical examples**

*θ*(*r*, *φ*) = <sup>1</sup>

*�*(*r*, *φ*) = *�*b(*φ*)

will be kept same for the *x* and *y* directions (*M* = *N*).

2 

> 2

Fig. 7. Amplitude distributions of the scaled magnetic field |*Hz*| within the unit cell of four chosen eigenmodes

and the number corresponds to the band. The amplitude distributions of the modes are displayed in Fig. 7, and the corresponding convergence properties of the eigenfrequencies calculated by the above described models are shown in Fig. 8 (Model C' is only compared for Sample S2).

The result of a careful comparison of the numerical efficiencies of all the factorization models can be summarized as follows. Models B, C, and C' always converge considerably faster than Model A. Model C converges faster (having usually one order higher precision) than Model B with two exceptions. The first exception is the case such as in Fig. 7 (S1-Γ2) and Fig. 8(a), where the discontinuities of the polarization transformation matrix **F** coincide with a nearly zero amplitude of the field, so that the discontinuities do not manifest themselves. The second exception is the case such as Fig. 7 (S2-X3) and Fig. 8(d), where the elements are densely arranged which causes rapid variations of the ellipticity between two adjacent elements (which are very close to each other); this requires more Fourier components than the weak discontinuity of the linear polarization u in Model B. The problem is solved in Model C', which obviously converges fastest among all the four models applied to Sample S2.

(a) Model A (b) Model B (c) Model C

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 335

<sup>4</sup> , *<sup>π</sup>*

This means that u is always linear, perpendicular to the permittivity discontinuities, and its discontinuities (along the lines *x* = *y* and *x* = −*y*) have no concurrent discontinuities with the permittivity function except those at the four corners of the square. Hence, Model B here fulfills nearly the same conditions for the application of the factorization rules as demanded

For Model C we divide the unit cell into four areas in the same manner. This time, however, the basis vector u(*r*, *φ*) depends on both polar coordinates and the corresponding polarization is in general elliptic. In analogy with Section 4 we want u to be perpendicular to the permittivity discontinuities and to remove its discontinuities as much as possible. The most simple way how to do this, although the discontinuities will not be completely removed, is

the remaining angles, and to use Equation 64 for the ellipticity distribution, where *D*(*φ*) now

*N* 

Fig. 10. Convergence of the factorization methods for rods with the square cross-section

<sup>4</sup> ) ∪ � <sup>3</sup>*<sup>π</sup>*

<sup>4</sup> , <sup>5</sup>*<sup>π</sup>*

<sup>4</sup> ), and ˆy for the remaining angles.

<sup>4</sup> , *<sup>π</sup>*

<sup>4</sup> ) ∪ � <sup>3</sup>*<sup>π</sup>*

<sup>4</sup> , <sup>5</sup>*<sup>π</sup>*

<sup>4</sup> ), and *<sup>π</sup>*

<sup>2</sup> for

Fig. 9. Polarization distribution of the factorization methods for periodic squares

to set the azimuth *<sup>θ</sup>*(*φ*) of the polarization to zero for *<sup>φ</sup>* ∈ �−*<sup>π</sup>*

on the azimuthal angle, is set ˆ<sup>x</sup> for *<sup>φ</sup>* ∈ �−*<sup>π</sup>*

in Section 4 for circular elements.

corresponds to the square element.

normalized frequency

#### **5. Two-dimensional photonic crystals with non-circular elements**

In this section we briefly show numerical efficiencies of three factorization models derived above (Models A, B, and C) applied to 2D photonic crystals made as long elements of other shapes, namely rods with the square cross-section and tubes with the ring and split-ring cross-sections, arranged with the square periodicity. For all the three samples we choose the permittivity of the elements *ε* = 9 and permittivity of vacuum *ε* = 1 for the surrounding medium.

#### **5.1 Periodic rods with the square cross-section**

For the photonic crystals made of square-sectioned rods, the period is chosen Λ = 1000 nm, and the width of the square *d* = 600 nm. The distribution of the basis polarization vector u, analogously to Section 4, is displayed in Fig. 9 for all the three compared models.

For Model A, of course, u = xˆ as visible in Fig. 9(a). For Model B, we divide the unit cell into four parts by its diagonals (the lines *x* = *y* and *x* = −*y*). The vector u(*φ*), depending only 16 Will-be-set-by-IN-TECH

Γ2

M1

0.7176

0.442

0.4421

normalized frequency

In this section we briefly show numerical efficiencies of three factorization models derived above (Models A, B, and C) applied to 2D photonic crystals made as long elements of other shapes, namely rods with the square cross-section and tubes with the ring and split-ring cross-sections, arranged with the square periodicity. For all the three samples we choose the permittivity of the elements *ε* = 9 and permittivity of vacuum *ε* = 1 for the surrounding

For the photonic crystals made of square-sectioned rods, the period is chosen Λ = 1000 nm, and the width of the square *d* = 600 nm. The distribution of the basis polarization vector u,

For Model A, of course, u = xˆ as visible in Fig. 9(a). For Model B, we divide the unit cell into four parts by its diagonals (the lines *x* = *y* and *x* = −*y*). The vector u(*φ*), depending only

analogously to Section 4, is displayed in Fig. 9 for all the three compared models.

0.4422

5 10 15 20 25 30

5 10 15 20 25 30

*N*

(d) Sample S2, point X3

*N* (b) Sample S1, point Γ3

model A model B

model C

Γ3

X3

model B model C model C'

0.718

0.7184

normalized frequency

0.7188

model A model B model C

5 10 15 20 25 30

*N*

5 10 15 20 25 30

*N* (c) Sample H, point M1

**5.1 Periodic rods with the square cross-section**

model A model B model C

Fig. 8. Convergences of the Fourier factorization methods for three samples

**5. Two-dimensional photonic crystals with non-circular elements**

(a) Sample S1, point Γ2

0.52576 0.52577 0.52578 0.52578 0.52579 0.5258 0.52581 0.52582

0.1837

0.1838

0.1839

normalized frequency

medium.

0.184

0.1841

0.1842

normalized frequency

Fig. 9. Polarization distribution of the factorization methods for periodic squares

on the azimuthal angle, is set ˆ<sup>x</sup> for *<sup>φ</sup>* ∈ �−*<sup>π</sup>* <sup>4</sup> , *<sup>π</sup>* <sup>4</sup> ) ∪ � <sup>3</sup>*<sup>π</sup>* <sup>4</sup> , <sup>5</sup>*<sup>π</sup>* <sup>4</sup> ), and ˆy for the remaining angles. This means that u is always linear, perpendicular to the permittivity discontinuities, and its discontinuities (along the lines *x* = *y* and *x* = −*y*) have no concurrent discontinuities with the permittivity function except those at the four corners of the square. Hence, Model B here fulfills nearly the same conditions for the application of the factorization rules as demanded in Section 4 for circular elements.

For Model C we divide the unit cell into four areas in the same manner. This time, however, the basis vector u(*r*, *φ*) depends on both polar coordinates and the corresponding polarization is in general elliptic. In analogy with Section 4 we want u to be perpendicular to the permittivity discontinuities and to remove its discontinuities as much as possible. The most simple way how to do this, although the discontinuities will not be completely removed, is to set the azimuth *<sup>θ</sup>*(*φ*) of the polarization to zero for *<sup>φ</sup>* ∈ �−*<sup>π</sup>* <sup>4</sup> , *<sup>π</sup>* <sup>4</sup> ) ∪ � <sup>3</sup>*<sup>π</sup>* <sup>4</sup> , <sup>5</sup>*<sup>π</sup>* <sup>4</sup> ), and *<sup>π</sup>* <sup>2</sup> for the remaining angles, and to use Equation 64 for the ellipticity distribution, where *D*(*φ*) now corresponds to the square element.

Fig. 10. Convergence of the factorization methods for rods with the square cross-section

normalized frequency

approximation would be insufficient.

with some deviations visible in the inset of Fig. 14.

*N* 

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 337

(a) Model A (b) Model B (c) Model C

Since the area of splitting is quite small compared to the full area of the unit cell, we have chosen the polarization distributions of all the three models exactly same as for the symmetric rings in Section 5.2. The condition for normal u and tangential v is not satisfied on the surface proportional to the length 2 × 120 = 240 nm, but is satisfied on the surface proportional to 0.9 × 2*π*(*R*<sup>1</sup> + *R*2) ≈ 7 464 nm, which is 30 times higher area, justifying this negligence. Of course, for modes with most of electromagnetic field resonating in the critical area this

The numerical efficiencies of all the three models are displayed in Fig. 14, again for the 3rd band of the symmetry point Γ. We can describe these performances by exactly the same conclusion as for the symmetric rings in Section 5.2. Models B and C converge similarly and both considerably faster than Model A. The discontinuities of u at the center and along the boundaries of the periodic cell again coincide with the zero amplitude of the mode, though

Fig. 13. Polarization distribution of the factorization methods for split hollow cylinders

Fig. 12. Convergence of the factorization methods for hollow cylinders

The numerical efficiencies of all the three models are displayed in Fig. 10 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, with Model C being a little better. Although the improvement from Model B towards Model C is not so distinct, it could be improved by a more careful choice of a completely continuous distribution of the polarization basis.

#### **5.2 Periodic hollow cylinders**

For the photonic crystals made of periodic hollow cylinders (tubes with the symmetric ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner circular hole) *R*<sup>1</sup> = 400 nm, and the outer diameter *R*<sup>2</sup> = 680 nm. The distribution of the basis polarization vector u for all the three compared models is displayed in Fig. 11.

Fig. 11. Polarization distribution of the factorization methods for hollow cylinders

For Model A again u = xˆ. For both Models B and C we use the same polarization distributions as in Section 4 for the inner area (*r* < *R*1) and for the outer area (*r* > *R*2). For the annulus area (*R*<sup>1</sup> < *r* < *R*2) we simply choose the polarization distribution same as that in Model B. These distributions are obviously the most straightforward analogies of the distributions used in Section 4 for circular elements.

The numerical efficiencies of all the three models are displayed in Fig. 12 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, but now Model C exhibits no improvement against Model B. This is because the discontinuities of u at the center and along the boundaries of the periodic cell quite well coincide with the zero amplitude of the mode, as visible in the inset of Fig. 12, so that the discontinuities do not manifest themselves in the calculations.

#### **5.3 Periodic split hollow cylinders**

For the photonic crystals made of split hollow cylinders (tubes with an asymmetric, split ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner semi-circular hole) *R*<sup>1</sup> = 600 nm, the outer diameter *R*<sup>2</sup> = 720 nm, and the relative azimuthal length of the ring *w<sup>φ</sup>* = 0.9 (where *w<sup>φ</sup>* = 1 means the complete, symmetric ring). The distributions of the basis polarization vector u for all the three compared models are displayed in Fig. 13.

18 Will-be-set-by-IN-TECH

The numerical efficiencies of all the three models are displayed in Fig. 10 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, with Model C being a little better. Although the improvement from Model B towards Model C is not so distinct, it could be improved by a more careful choice of a completely

For the photonic crystals made of periodic hollow cylinders (tubes with the symmetric ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner circular hole) *R*<sup>1</sup> = 400 nm, and the outer diameter *R*<sup>2</sup> = 680 nm. The distribution of the basis polarization vector u for all the three compared models is displayed in Fig. 11.

(a) Model A (b) Model B (c) Model C

For Model A again u = xˆ. For both Models B and C we use the same polarization distributions as in Section 4 for the inner area (*r* < *R*1) and for the outer area (*r* > *R*2). For the annulus area (*R*<sup>1</sup> < *r* < *R*2) we simply choose the polarization distribution same as that in Model B. These distributions are obviously the most straightforward analogies of the distributions used

The numerical efficiencies of all the three models are displayed in Fig. 12 for the 3rd band of the symmetry point Γ. As clearly visible, Models B and C converge considerably faster than Model A, but now Model C exhibits no improvement against Model B. This is because the discontinuities of u at the center and along the boundaries of the periodic cell quite well coincide with the zero amplitude of the mode, as visible in the inset of Fig. 12, so that the

For the photonic crystals made of split hollow cylinders (tubes with an asymmetric, split ring cross-section) we choose the period Λ = 1000 nm, the inner diameter (the diameter of the inner semi-circular hole) *R*<sup>1</sup> = 600 nm, the outer diameter *R*<sup>2</sup> = 720 nm, and the relative azimuthal length of the ring *w<sup>φ</sup>* = 0.9 (where *w<sup>φ</sup>* = 1 means the complete, symmetric ring). The distributions of the basis polarization vector u for all the three compared models are

Fig. 11. Polarization distribution of the factorization methods for hollow cylinders

discontinuities do not manifest themselves in the calculations.

continuous distribution of the polarization basis.

**5.2 Periodic hollow cylinders**

in Section 4 for circular elements.

**5.3 Periodic split hollow cylinders**

displayed in Fig. 13.

Fig. 12. Convergence of the factorization methods for hollow cylinders

Fig. 13. Polarization distribution of the factorization methods for split hollow cylinders

Since the area of splitting is quite small compared to the full area of the unit cell, we have chosen the polarization distributions of all the three models exactly same as for the symmetric rings in Section 5.2. The condition for normal u and tangential v is not satisfied on the surface proportional to the length 2 × 120 = 240 nm, but is satisfied on the surface proportional to 0.9 × 2*π*(*R*<sup>1</sup> + *R*2) ≈ 7 464 nm, which is 30 times higher area, justifying this negligence. Of course, for modes with most of electromagnetic field resonating in the critical area this approximation would be insufficient.

The numerical efficiencies of all the three models are displayed in Fig. 14, again for the 3rd band of the symmetry point Γ. We can describe these performances by exactly the same conclusion as for the symmetric rings in Section 5.2. Models B and C converge similarly and both considerably faster than Model A. The discontinuities of u at the center and along the boundaries of the periodic cell again coincide with the zero amplitude of the mode, though with some deviations visible in the inset of Fig. 14.

To apply the factorization correctly, we must again separate the normal and tangential components of all fields for which we can use the correct rules. Let us define a

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 339

where *Eu* and *Ev* are the normal and tangential components of the vector E to all

Analogously to Section 4, for the case of circular elements we can choose the polarization basis

⎤ ⎥

where the matrix elements have the same meaning as in Section 4. In the new coordinates we

Now let us separate two sets of quantities, those which are continuous (*D*˜ *<sup>u</sup>*, *Ev*, *Ez*) and those which are not continuous (*D*˜ *<sup>v</sup>*, *D*˜ *<sup>z</sup>*, *Eu*) to the discontinuities of the permittivity. Expressing

> ⎡ ⎢ ⎢ ⎢ ⎣

1 *<sup>ε</sup>uu* <sup>−</sup> *<sup>ε</sup>uv*

*εvu*

*εzu*

⎡ ⎢ ⎣

*εuu εuv εuz εvu εvv εvz εzu εzv εzz* ⎤ ⎥ ⎦

*<sup>ε</sup>uu <sup>ε</sup>vv* <sup>−</sup> *<sup>ε</sup>vuεuv*

*<sup>ε</sup>uu <sup>ε</sup>zv* <sup>−</sup> *<sup>ε</sup>zuεuv*

[*D*˜ *<sup>u</sup>*] [*Ev*] [*Ez*]

> ⎡ ⎢ ⎣

[*Eu*] [*Ev*] [*Ez*] ⎤ ⎥

⎤ ⎥

⎡ ⎣ *Eu Ev Ez*

⎤

*<sup>ε</sup>uu* <sup>−</sup> *<sup>ε</sup>uz*

*<sup>ε</sup>uu <sup>ε</sup>vz* <sup>−</sup> *<sup>ε</sup>vuεuz*

*<sup>ε</sup>uu <sup>ε</sup>zz* <sup>−</sup> *<sup>ε</sup>zuεuz*

*εuu*

*εuu*

⎤ ⎥ ⎥ ⎥

<sup>⎦</sup> . (82)

<sup>⎦</sup> , (83)

<sup>⎦</sup> , (81)

*εuu*

⎤ ⎥

cos *φ* − sin *φ* 0 sin *φ* cos *φ* 0 0 01

> *ξ* −*ζ*<sup>∗</sup> 0 *ζ ξ*∗ 0 001

> > ⎡ ⎢ ⎣

⎡ ⎣ *Eu Ev Ez* ⎤

⎦ , (77)

<sup>⎦</sup> , (Model B) (78)

<sup>⎦</sup> , (Model C) (79)

⎦ . (80)

⎡ ⎣ *Ex Ey Ez*

⎤ ⎦ = **F**

space-dependent matrix transform **F**(*x*, *y*) so that

distribution corresponding to Models B and C as

**F** =

⎡ ⎢ ⎣

**F** =

⎡ ⎣ *D*˜ *u D*˜ *v D*˜ *z*

⎡ ⎣ *D*˜ *u Ev Ez*

⎤

⎡ ⎢ ⎣

⎡ ⎢ ⎣ [*D*˜ *<sup>u</sup>*] [*D*˜ *<sup>v</sup>*] [*D*˜ *<sup>z</sup>*] ⎤ ⎥

[*Eu*] [*D*˜ *<sup>v</sup>*] [*D*˜ *<sup>z</sup>*] ⎤ ⎥ <sup>⎦</sup> = [[**G**]]

<sup>⎦</sup> = [[ε{*uvz*}]]B, C

⎦ , **G** =

the second set according to the first one yields

⎤ ⎦ = **G**

for which we can simply use the Laurent rule,

From this we express [*D*˜ *<sup>j</sup>*] according to [*Ej*],

⎡ ⎣ *Eu D*˜ *v D*˜ *z* ⎡ ⎢ ⎣

⎤ ⎦ =

discontinuities of the permittivity.

can write

Fig. 14. Convergence of the factorization methods for split hollow cylinders

#### **6. Anisotropic photonic crystals**

In this section we will briefly demonstrate the application of the factorization rules to 2D photonic crystals made of anisotropic materials, again with the plane of incidence parallel to the *z* axis (with geometry of Fig. 3). Unlike the isotropic crystals, now the TE and TM polarizations are not separable. Instead of the scalar permittivity we define and expand the components of the relative permittivity tensor function

$$\varepsilon\_{jk}(x,y) = \sum\_{m,n=-\infty}^{+\infty} \varepsilon\_{jk,mn} e^{-i(mpx+nqy)}\,,\tag{73}$$

where *εjk*,*mn* are the Fourier coefficients. The wave equation for a generally anisotropic medium, now described by the permittivity tensor

$$
\mathfrak{e} = \begin{bmatrix}
\mathfrak{e}\_{xx} \ \mathfrak{e}\_{xy} \ \mathfrak{e}\_{xz} \\
\mathfrak{e}\_{yx} \ \mathfrak{e}\_{yy} \ \mathfrak{e}\_{yz} \\
\mathfrak{e}\_{zx} \ \mathfrak{e}\_{zy} \ \mathfrak{e}\_{zz}
\end{bmatrix} \tag{74}
$$

is the operator equation

$$
\hat{\mathbf{C}}E = \frac{\omega^2}{c^2} \mathbf{E}, \quad \hat{\mathbf{C}} = \varepsilon^{-1} \begin{bmatrix} -\partial\_y^2 \,\partial\_x \partial\_y & 0 \\ \partial\_x \partial\_y \, -\partial\_x^2 & 0 \\ 0 & 0 & -\partial\_x^2 - \partial\_y^2 \end{bmatrix}. \tag{75}
$$

Similarly as above, Model A assumes all components of the permittivity tensor treated by the Laurent rule, i.e.,

$$\mathbb{E}[\tilde{D}\_j] = \sum\_k \mathbb{[e}\_{jk}] \mathbb{I}[E\_k]. \tag{76}$$

20 Will-be-set-by-IN-TECH

*N* 

In this section we will briefly demonstrate the application of the factorization rules to 2D photonic crystals made of anisotropic materials, again with the plane of incidence parallel to the *z* axis (with geometry of Fig. 3). Unlike the isotropic crystals, now the TE and TM polarizations are not separable. Instead of the scalar permittivity we define and expand the

*εjk*,*mne*

where *εjk*,*mn* are the Fourier coefficients. The wave equation for a generally anisotropic

*εxx εxy εxz εyx εyy εyz εzx εzy εzz*

> ⎡ ⎢ ⎢ ⎣

Similarly as above, Model A assumes all components of the permittivity tensor treated by the

[*D*˜ *<sup>j</sup>*] = ∑ *k*

<sup>−</sup>*∂*<sup>2</sup>

*<sup>∂</sup>x∂<sup>y</sup>* <sup>−</sup>*∂*<sup>2</sup>

⎤ ⎥ ⎥

*<sup>y</sup> ∂x∂<sup>y</sup>* 0

0 0 <sup>−</sup>*∂*<sup>2</sup>

*<sup>x</sup>* 0

*<sup>x</sup>* <sup>−</sup> *<sup>∂</sup>*<sup>2</sup> *y* ⎤ ⎥ ⎥

[[*εjk*]][*Ek*]. (76)

−*i*(*mpx*+*nqy*)

, (73)

<sup>⎦</sup> . (75)

<sup>⎦</sup> , (74)

+∞ <sup>∑</sup>*m*,*n*=−<sup>∞</sup>

ε =

*<sup>c</sup>*<sup>2</sup> <sup>E</sup>, **Cˆ** <sup>=</sup> <sup>ε</sup>−<sup>1</sup>

⎡ ⎢ ⎢ ⎣

Fig. 14. Convergence of the factorization methods for split hollow cylinders

normalized frequency

**6. Anisotropic photonic crystals**

is the operator equation

Laurent rule, i.e.,

components of the relative permittivity tensor function

medium, now described by the permittivity tensor

**Cˆ** <sup>E</sup> <sup>=</sup> *<sup>ω</sup>*<sup>2</sup>

*εjk*(*x*, *y*) =

To apply the factorization correctly, we must again separate the normal and tangential components of all fields for which we can use the correct rules. Let us define a space-dependent matrix transform **F**(*x*, *y*) so that

$$
\begin{bmatrix} E\_x \\ E\_y \\ E\_z \end{bmatrix} = \mathbf{F} \begin{bmatrix} E\_u \\ E\_v \\ E\_z \end{bmatrix}'\tag{77}
$$

where *Eu* and *Ev* are the normal and tangential components of the vector E to all discontinuities of the permittivity.

Analogously to Section 4, for the case of circular elements we can choose the polarization basis distribution corresponding to Models B and C as

$$\mathbf{F} = \begin{bmatrix} \cos \phi - \sin \phi & 0 \\ \sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \end{bmatrix} \text{ (Model B)}\tag{78}$$

$$\mathbf{F} = \begin{bmatrix} \xi & -\xi^\* & 0 \\ \zeta & \xi^\* & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \text{(Model } \mathbb{C}) \tag{79}$$

where the matrix elements have the same meaning as in Section 4. In the new coordinates we can write

$$
\begin{bmatrix} \mathcal{D}\_{\iota} \\ \tilde{D}\_{\upsilon} \\ \tilde{D}\_{z} \end{bmatrix} = \begin{bmatrix} \varepsilon\_{\iota \iota \iota} & \varepsilon\_{\iota \upsilon} & \varepsilon\_{\iota \upsilon} \\ \varepsilon\_{\upsilon \iota} & \varepsilon\_{\upsilon \upsilon} & \varepsilon\_{\upsilon z} \\ \varepsilon\_{z \iota} & \varepsilon\_{z \upsilon} & \varepsilon\_{zz} \end{bmatrix} \begin{bmatrix} E\_{\iota} \\ E\_{\upsilon} \\ E\_{z} \end{bmatrix} . \tag{80}
$$

Now let us separate two sets of quantities, those which are continuous (*D*˜ *<sup>u</sup>*, *Ev*, *Ez*) and those which are not continuous (*D*˜ *<sup>v</sup>*, *D*˜ *<sup>z</sup>*, *Eu*) to the discontinuities of the permittivity. Expressing the second set according to the first one yields

$$\begin{bmatrix} \begin{matrix} E\_{\boldsymbol{u}} \\ \tilde{D}\_{\boldsymbol{\upsilon}} \\ \tilde{D}\_{\boldsymbol{z}} \end{matrix} = \mathbf{G} \begin{bmatrix} \tilde{D}\_{\boldsymbol{u}} \\ E\_{\boldsymbol{\upsilon}} \\ E\_{\boldsymbol{z}} \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} \frac{1}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} & -\frac{\xi\_{\boldsymbol{u}\boldsymbol{u}}}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} & -\frac{\xi\_{\boldsymbol{u}\boldsymbol{u}}}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} \\\\ \frac{\xi\_{\boldsymbol{u}\boldsymbol{u}}}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} \, \varepsilon\_{\boldsymbol{\upsilon}\boldsymbol{u}} - \frac{\xi\_{\boldsymbol{u}\boldsymbol{u}} \varepsilon\_{\boldsymbol{u}\boldsymbol{u}}}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} \, \varepsilon\_{\boldsymbol{\upsilon}\boldsymbol{u}} - \frac{\varepsilon\_{\boldsymbol{\upsilon}\boldsymbol{u}} \varepsilon\_{\boldsymbol{u}\boldsymbol{u}}}{\xi\_{\boldsymbol{u}\boldsymbol{u}}} \end{bmatrix}, \tag{81}$$

for which we can simply use the Laurent rule,

$$
\begin{bmatrix} \left[ E\_{\mathcal{U}} \right] \\ \left[ \tilde{D}\_{\mathcal{U}} \right] \\ \left[ \tilde{D}\_{z} \right] \end{bmatrix} = \begin{bmatrix} \mathbf{G} \\ \mathbf{G} \end{bmatrix} \begin{bmatrix} \left[ \tilde{D}\_{\mathcal{U}} \right] \\ \left[ E\_{\mathcal{U}} \right] \\ \left[ E\_{z} \right] \end{bmatrix} . \tag{82}
$$

From this we express [*D*˜ *<sup>j</sup>*] according to [*Ej*],

$$
\begin{bmatrix}
\left[\tilde{D}\_{\boldsymbol{u}}\right] \\
\left[\tilde{D}\_{\boldsymbol{v}}\right] \\
\left[\tilde{D}\_{\boldsymbol{z}}\right]
\end{bmatrix} = \left[\boldsymbol{\varepsilon}\_{\{\boldsymbol{u}\boldsymbol{v}\boldsymbol{z}\}}\right] \mathbb{I}\_{\mathrm{B},\mathsf{C}} \begin{bmatrix}
\left[E\_{\boldsymbol{u}}\right] \\
\left[E\_{\boldsymbol{v}}\right] \\
\left[E\_{\boldsymbol{z}}\right]
\end{bmatrix} \text{/} \tag{83}
$$

Obviously, u is normal on the whole sphere's surface. The vector v corresponds to the left circular polarization at the sphere's poles (*z* = ±*R*), to the vertical linear polarization on the sphere's equator (*z* = 0), and is elliptical and continuous in the intermediate ranges. The vector w is simply chosen orthogonal to both u and v. Both v and w are tangential to the

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 341

If we extend the radial dependence of the vectors defined by Equations 89–91 to the entire cubic cell (formally replacing *R* by *r*), then we obtain some analogy with Model B. The obtained matrix **F** has no concurrent discontinuities with *ε*, but it is discontinuous at the center and on the boundaries of the cubic cell. To remove the discontinuity at the cell's center, we

√2

√2

<sup>2</sup> (u<sup>0</sup> <sup>−</sup> <sup>u</sup>*R*) cos *<sup>π</sup><sup>r</sup>*

where *D*(*ϑ*, *φ*) is the distance from the cell's center to its boundary along the ray determined

is a scalar function ensuring that u becomes a unit vector everywhere. We could extend v and

but it is not clear whether vext or wext are both perpendicular to u and mutually. To ensure it

where **<sup>P</sup>**<sup>u</sup> <sup>=</sup> uu† is the projector into the space of vectors proportional to <sup>u</sup>, so that **<sup>1</sup>** <sup>−</sup> **<sup>P</sup>**<sup>u</sup> is

where **1** − **P**<sup>u</sup> − **P**<sup>v</sup> is the projector into the space of vectors perpendicular to both u and v. Here *A*<sup>v</sup> and *A*<sup>w</sup> are analogously chosen scalar functions ensuring the unit size of the

<sup>2</sup> (v<sup>0</sup> <sup>−</sup> <sup>v</sup>*R*) cos *<sup>π</sup><sup>r</sup>*

<sup>2</sup> (w<sup>0</sup> <sup>−</sup> <sup>w</sup>*R*) cos *<sup>π</sup><sup>r</sup>*

<sup>2</sup> (v<sup>0</sup> <sup>−</sup> <sup>v</sup>*R*) cos *<sup>π</sup>*(*r*+*D*−2*R*)

<sup>2</sup> (w<sup>0</sup> <sup>−</sup> <sup>w</sup>*R*) cos *<sup>π</sup>*(*r*+*D*−2*R*)

<sup>2</sup> (u<sup>0</sup> <sup>−</sup> <sup>u</sup>*R*) cos *<sup>π</sup>*(*r*+*D*−2*R*)

v<sup>0</sup> = v(0, *ϑ*, *φ*) = zˆ, (93)

(xˆ − *i*yˆ), (92)

(xˆ + *i*yˆ), (94)

*<sup>D</sup>*−*<sup>R</sup>* (*<sup>R</sup>* <sup>&</sup>lt; *<sup>r</sup>* <sup>≤</sup> *<sup>D</sup>*). (95)

*<sup>D</sup>*−*<sup>R</sup>* (*<sup>R</sup>* <sup>&</sup>lt; *<sup>r</sup>* <sup>≤</sup> *<sup>D</sup>*). (97)

*<sup>D</sup>*−*<sup>R</sup>* (*<sup>R</sup>* <sup>&</sup>lt; *<sup>r</sup>* <sup>≤</sup> *<sup>D</sup>*). (98)

*<sup>D</sup>*−*<sup>R</sup>* )]−1/2 for *<sup>r</sup>* <sup>&</sup>gt; *<sup>R</sup>*

*<sup>R</sup>* (*r* ≤ *R*)

u = *A*uuext, (96)

<sup>2</sup> (<sup>1</sup> <sup>+</sup> cos2 *<sup>π</sup>*(*r*+*D*−2*R*)

*<sup>R</sup>* (*r* ≤ *R*)

*<sup>R</sup>* (*r* ≤ *R*)

v = *A*v(**1** − **P**u)vext, (99)

w = *A*w(**1** − **P**<sup>u</sup> − **P**v)wext, (100)

As the second step, we define the basis vectors at the center of the sphere,

<sup>2</sup> (u<sup>0</sup> <sup>+</sup> <sup>u</sup>*R*) + <sup>1</sup>

<sup>2</sup> (u<sup>0</sup> <sup>+</sup> <sup>u</sup>*R*) + <sup>1</sup>

<sup>2</sup> (v<sup>0</sup> <sup>+</sup> <sup>v</sup>*R*) + <sup>1</sup>

<sup>2</sup> (v<sup>0</sup> <sup>+</sup> <sup>v</sup>*R*) + <sup>1</sup>

<sup>2</sup> (w<sup>0</sup> <sup>+</sup> <sup>w</sup>*R*) + <sup>1</sup>

<sup>2</sup> (w<sup>0</sup> <sup>+</sup> <sup>w</sup>*R*) + <sup>1</sup>

1

1

1

the projector to its orthogonal complement. Similarly,

<sup>u</sup><sup>0</sup> <sup>=</sup> <sup>u</sup>(0, *<sup>ϑ</sup>*, *<sup>φ</sup>*) = <sup>1</sup>

<sup>w</sup><sup>0</sup> <sup>=</sup> <sup>w</sup>(0, *<sup>ϑ</sup>*, *<sup>φ</sup>*) = <sup>1</sup>

which are again mutually orthogonal. Then we extend u onto the whole cubic cell,

by the spherical angles. The desired unit vector of the polarization basis is then

*<sup>R</sup>* )]−1/2 for *<sup>r</sup>* <sup>&</sup>lt; *<sup>R</sup>* and *<sup>A</sup>*<sup>u</sup> = [ <sup>1</sup>

sphere's surface.

where *A*<sup>u</sup> = [ <sup>1</sup>

we define

corresponding vectors.

w in a similar way,

will proceed differently.

<sup>u</sup>ext(*r*, *<sup>ϑ</sup>*, *<sup>φ</sup>*) = <sup>1</sup>

<sup>2</sup> (<sup>1</sup> <sup>+</sup> cos<sup>2</sup> *<sup>π</sup><sup>r</sup>*

<sup>v</sup>ext(*r*, *<sup>ϑ</sup>*, *<sup>φ</sup>*) = <sup>1</sup>

<sup>w</sup>ext(*r*, *<sup>ϑ</sup>*, *<sup>φ</sup>*) = <sup>1</sup>

with

[[ε{*uvz*}]]B, C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ �� <sup>1</sup> *<sup>ε</sup>uu* ��−<sup>1</sup> , �� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuv εuu* �� , �� *εvu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−<sup>1</sup> , �� *<sup>ε</sup>vv* <sup>−</sup> *<sup>ε</sup>vuεuv <sup>ε</sup>uu* �� + �� *εvu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuv εuu* �� , �� *εzu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−<sup>1</sup> , �� *<sup>ε</sup>zv* <sup>−</sup> *<sup>ε</sup>zuεuv <sup>ε</sup>uu* �� + �� *εzu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuv εuu* �� , �� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuz εuu* �� �� *<sup>ε</sup>vz* <sup>−</sup> *<sup>ε</sup>vuεuz <sup>ε</sup>uu* �� + �� *εvu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuz εuu* �� �� *<sup>ε</sup>zz* <sup>−</sup> *<sup>ε</sup>zuεuz <sup>ε</sup>uu* �� + �� *εzu <sup>ε</sup>uu* ���� <sup>1</sup> *<sup>ε</sup>uu* ��−1�� *εuz εuu* �� ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . (84)

Finally we obtain the matrix of permittivity in the Cartesian coordinates via the formula

$$\left\| \varepsilon\_{\{xyz\}} \right\|\_{\mathcal{B},\mathcal{C}} = \left\| \mathbf{F} \right\| \left\| \varepsilon\_{\{uvz\}} \right\|\_{\mathcal{B},\mathcal{C}} \left\| \mathbf{F}^{-1} \right\|, \tag{85}$$

where the index B or C corresponds to the chosen model.

#### **7. Distribution of polarization basis for 3D structures**

In this section we will suggest how to create the polarization basis distribution for a 3D photonic crystal. Unlike the previous cases, we now have a fully vectorial wave Equation 1. The corresponding material equation D˜ = *ε*(*x*, *y*, *z*)E must be changed to separate the normal and tangential components of the fields to the *ε* discontinuities, which are now surfaces in the 3D space. For simplicity we assume a 3D photonic crystal made as spheres (with the radius R) arranged in the space with the cubic periodicity. To make a 3D analogy with Model C described in the previous sections, we must find a matrix transform **F**(*x*, *y*, *z*) whose columns, denoted u, v, and w, are complex vectorial functions of space, mutually orthonormal and continuous at all points, where u is the normal vector at each point of the sphere's surface and v and w are tangential.

As the first step we choose the distribution of these vectors on the sphere. Defining the spherical coordinates

$$
\dot{x} = r \sin \theta \cos \phi,\tag{86}
$$

$$y = r \sin \theta \sin \phi\_{\prime} \tag{87}$$

$$z = r \cos \theta \tag{88}$$

and the corresponding unit vectors ˆr, ϑˆ, and φˆ (pointing along the increase of the corresponding coordinate) helps us to define the vectors

$$u\_R(\vartheta, \phi) = u(R, \vartheta, \phi) = e^{i(\vartheta + \phi)} \hat{r},\tag{89}$$

$$\upsilon\_R(\theta,\phi) = \upsilon(R,\theta,\phi) = e^{i(\theta+\phi)} \frac{1}{\sqrt{1+\cos^2\theta}} (\hat{\vartheta} + i\hat{\phi}\cos\theta),\tag{90}$$

$$w\_R(\theta,\phi) = w(R,\theta,\phi) = e^{i(\theta+\phi)} \frac{1}{\sqrt{1+\cos^2\theta}} (\hat{\theta}\cos\theta - i\hat{\phi}).\tag{91}$$

22 Will-be-set-by-IN-TECH

��−1�� *εuv εuu* �� ,

��

��

B, C = [[**F**]][[ε{*uvz*}]]

In this section we will suggest how to create the polarization basis distribution for a 3D photonic crystal. Unlike the previous cases, we now have a fully vectorial wave Equation 1. The corresponding material equation D˜ = *ε*(*x*, *y*, *z*)E must be changed to separate the normal and tangential components of the fields to the *ε* discontinuities, which are now surfaces in the 3D space. For simplicity we assume a 3D photonic crystal made as spheres (with the radius R) arranged in the space with the cubic periodicity. To make a 3D analogy with Model C described in the previous sections, we must find a matrix transform **F**(*x*, *y*, *z*) whose columns, denoted u, v, and w, are complex vectorial functions of space, mutually orthonormal and continuous at all points, where u is the normal vector at each point of the sphere's surface and

As the first step we choose the distribution of these vectors on the sphere. Defining the

and the corresponding unit vectors ˆr, ϑˆ, and φˆ (pointing along the increase of the

*i*(*ϑ*+*φ*)

*<sup>i</sup>*(*ϑ*+*φ*) 1 √

*<sup>i</sup>*(*ϑ*+*φ*) 1 √

1 + cos<sup>2</sup> *ϑ*

1 + cos2 *ϑ*

Finally we obtain the matrix of permittivity in the Cartesian coordinates via the formula

*<sup>ε</sup>vz* <sup>−</sup> *<sup>ε</sup>vuεuz εuu* �� + �� *εvu εuu* ���� <sup>1</sup> *εuu*

*<sup>ε</sup>zz* <sup>−</sup> *<sup>ε</sup>zuεuz εuu* �� + �� *εzu εuu* ���� <sup>1</sup> *εuu*

��−1�� *εuv εuu* �� ,

��−1�� *εuv εuu* �� ,

> �� <sup>1</sup> *εuu*

��−1�� *εuz εuu* ��

*x* = *r* sin *ϑ* cos *φ*, (86) *y* = *r* sin *ϑ* sin *φ*, (87) *z* = *r* cos *ϑ* (88)

rˆ, (89)

(ϑˆ + *i*φˆ cos *ϑ*), (90)

(ϑ<sup>ˆ</sup> cos *<sup>ϑ</sup>* <sup>−</sup> *<sup>i</sup>*φˆ). (91)

��−1�� *εuz εuu* �� ⎤ ⎥ ⎥ ⎥ ⎥ ⎦

. (84)

��−1�� *εuz εuu* ��

B, C[[**F**−1]], (85)

�� <sup>1</sup> *εuu*

*<sup>ε</sup>vv* <sup>−</sup> *<sup>ε</sup>vuεuv εuu* �� + �� *εvu εuu* ���� <sup>1</sup> *εuu*

*<sup>ε</sup>zv* <sup>−</sup> *<sup>ε</sup>zuεuv εuu* �� + �� *εzu εuu* ���� <sup>1</sup> *εuu*

[[ε{*xyz*}]]

where the index B or C corresponds to the chosen model.

**7. Distribution of polarization basis for 3D structures**

corresponding coordinate) helps us to define the vectors

u*R*(*ϑ*, *φ*) = u(*R*, *ϑ*, *φ*) = *e*

v*R*(*ϑ*, *φ*) = v(*R*, *ϑ*, *φ*) = *e*

w*R*(*ϑ*, *φ*) = w(*R*, *ϑ*, *φ*) = *e*

with

[[ε{*uvz*}]]

�� <sup>1</sup> *εuu* ��−<sup>1</sup> ,

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

�� *εvu εuu* ���� <sup>1</sup> *εuu* ��−<sup>1</sup> , ��

�� *εzu εuu* ���� <sup>1</sup> *εuu* ��−<sup>1</sup> , ��

v and w are tangential.

spherical coordinates

B, C =

Obviously, u is normal on the whole sphere's surface. The vector v corresponds to the left circular polarization at the sphere's poles (*z* = ±*R*), to the vertical linear polarization on the sphere's equator (*z* = 0), and is elliptical and continuous in the intermediate ranges. The vector w is simply chosen orthogonal to both u and v. Both v and w are tangential to the sphere's surface.

If we extend the radial dependence of the vectors defined by Equations 89–91 to the entire cubic cell (formally replacing *R* by *r*), then we obtain some analogy with Model B. The obtained matrix **F** has no concurrent discontinuities with *ε*, but it is discontinuous at the center and on the boundaries of the cubic cell. To remove the discontinuity at the cell's center, we will proceed differently.

As the second step, we define the basis vectors at the center of the sphere,

$$u\_0 = u(0, \theta, \phi) = \frac{1}{\sqrt{2}}(\dot{x} - i\dot{y}),\tag{92}$$

$$v\_0 = v(0, \theta, \phi) = \mathbf{\hat{z}},\tag{93}$$

$$w\_0 = w(0, \theta, \phi) = \frac{1}{\sqrt{2}}(\dot{x} + i\dot{y}),\tag{94}$$

which are again mutually orthogonal. Then we extend u onto the whole cubic cell,

$$u\_{\rm ext}(r,\theta,\phi) = \begin{cases} \frac{1}{2}(u\_0 + u\_R) + \frac{1}{2}(u\_0 - u\_R)\cos\frac{\pi r}{R} & (r \le R) \\\frac{1}{2}(u\_0 + u\_R) + \frac{1}{2}(u\_0 - u\_R)\cos\frac{\pi(r+D-2R)}{D-R} & (R < r \le D) \end{cases} \tag{95}$$

where *D*(*ϑ*, *φ*) is the distance from the cell's center to its boundary along the ray determined by the spherical angles. The desired unit vector of the polarization basis is then

$$
\mu = A\_{\mathfrak{u}} u\_{\text{ext}\prime} \tag{96}
$$

where *A*<sup>u</sup> = [ <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>+</sup> cos<sup>2</sup> *<sup>π</sup><sup>r</sup> <sup>R</sup>* )]−1/2 for *<sup>r</sup>* <sup>&</sup>lt; *<sup>R</sup>* and *<sup>A</sup>*<sup>u</sup> = [ <sup>1</sup> <sup>2</sup> (<sup>1</sup> <sup>+</sup> cos2 *<sup>π</sup>*(*r*+*D*−2*R*) *<sup>D</sup>*−*<sup>R</sup>* )]−1/2 for *<sup>r</sup>* <sup>&</sup>gt; *<sup>R</sup>* is a scalar function ensuring that u becomes a unit vector everywhere. We could extend v and w in a similar way,

$$\boldsymbol{v}\_{\text{ext}}(\boldsymbol{r},\boldsymbol{\theta},\boldsymbol{\phi}) = \begin{cases} \frac{1}{2}(\boldsymbol{v}\_{0} + \boldsymbol{v}\_{R}) + \frac{1}{2}(\boldsymbol{v}\_{0} - \boldsymbol{v}\_{R})\cos\frac{\pi\boldsymbol{r}}{R} & (\boldsymbol{r} \le R) \\\ \frac{1}{2}(\boldsymbol{v}\_{0} + \boldsymbol{v}\_{R}) + \frac{1}{2}(\boldsymbol{v}\_{0} - \boldsymbol{v}\_{R})\cos\frac{\pi(\boldsymbol{r} + D - 2R)}{D - R} & (R < \boldsymbol{r} \le D) . \end{cases} \tag{97}$$

$$w\_{\rm ext}(r,\theta,\phi) = \begin{cases} \frac{1}{2}(w\_0 + w\_R) + \frac{1}{2}(w\_0 - w\_R)\cos\frac{\pi r}{R} & (r \le R) \\\frac{1}{2}(w\_0 + w\_R) + \frac{1}{2}(w\_0 - w\_R)\cos\frac{\pi(r + D - 2R)}{D - R} & (R < r \le D) \end{cases} \tag{98}$$

but it is not clear whether vext or wext are both perpendicular to u and mutually. To ensure it we define

$$
\Delta v = A\_v (\mathbf{1} - \mathbf{P}\_u) v\_{\text{ext}\prime} \tag{99}
$$

where **<sup>P</sup>**<sup>u</sup> <sup>=</sup> uu† is the projector into the space of vectors proportional to <sup>u</sup>, so that **<sup>1</sup>** <sup>−</sup> **<sup>P</sup>**<sup>u</sup> is the projector to its orthogonal complement. Similarly,

$$
\Delta w = A\_{w} (\mathbf{1} - \mathbf{P}\_{u} - \mathbf{P}\_{v}) w\_{\text{ext}} \tag{100}
$$

where **1** − **P**<sup>u</sup> − **P**<sup>v</sup> is the projector into the space of vectors perpendicular to both u and v. Here *A*<sup>v</sup> and *A*<sup>w</sup> are analogously chosen scalar functions ensuring the unit size of the corresponding vectors.

Bonod, N., Popov, E. & Neviere, M. (2005b). Light transmission through a subwavelength

Fourier Factorization in the Plane Wave Expansion Method in Modeling Photonic Crystals 343

Boyer, P., Popov, E., Neviere, M. & Tayeb, G. (2004). Diffraction theory in TM polarization:

Chernov, B., Neviere, M. & Popov, E. (2001). Fast Fourier factorization method applied to

Datta, S., Chan, C. T., Ho, K. M. & Soukoulis, C. M. (1993). Effective dielectric constant of

David, A., Benisty, H. & Weisbuch, C. (2006). Fast factorization rule and plane-wave expansion

Dossou, K. B., Poulton, C. G., Botten, L. C., Mahmoodian, S., McPhedran, R. C. & de Sterke,

Genereux, F., Leonard, S. W., van Driel, H. M., Birner, A. & Gosele, U. (2001). Large

Gotz, P., Schuster, T., Frenner, K., Rafler, S. & Osten, W. (2008). Normal vector method for the RCWA with automated vector field generation, *Opt. Express* 16: 17295–17301. Ho, K. M., Chan, C. T. & Soukoulis, C. M. (1990). Existence of a photonic gap in periodic

Hugonin, J. P., Lalanne, P., White, T. P. & Krauss, T. F. (2007). Coupling into slow-mode

Joannopoulos, J. D., Meade, R. D. & Winn, J. N. (1995). *Photonic Crystals: Molding the Flow of*

Joannopoulos, J. D., Villeneuve, P. R. & Fan, S. (1997). Photonic crystals: Putting a new twist

Kinoshita, S. & Yoshioka, S. (2005). Structural colors in nature: The role of regularity and

Kosaka, H., Kawashima, T., Tomita, A., Notomi, M., Tamamura, T., Sato, T. & Kawakami, S.

Krokhin, A. A., Halevi, P. & Arriaga, J. (2002). Long-wavelength limit (homogenization) for

Krokhin, A. A. & Reyes, E. (2004). Homogenization of magnetodielectric photonic crystals,

Li, L. (1996). Use of Fourier series in the analysis of discontinuous periodic structures, *J. Opt.*

Li, L. (1997). New formulation of the Fourier modal method for crossed surface-relief gratings,

Li, L. (1998). Reformulation of the Fourier modal method for surface-relief gratings made with

Li, L. (2003). Fourier modal method for crossed anisotropic gratings with arbitrary

(1998). Superprism phenomena in photonic crystals, *Phys. Rev. B* 58: R10096–R10099.

arbitrary cross section, *J. Opt. Soc. Am. A* 21: 2146–2153.

periodic composite structures, *Phys. Rev. B* 48: 14936–14943.

dielectric structures, *Phys. Rev. Lett.* 65: 3152–3155.

photonic crystal waveguides, *Opt. Lett.* 32: 2638–2640.

irregularity in the structure, *ChemPhysChem* 6: 1442–1459.

two-dimensional photonic crystals, *Phys. Rev. B* 65: 115208.

permittivity and permeability tensors, *J. Opt. A* 5: 345–355.

245: 355–361.

*B* 73: 075107.

*Commun.* 194: 289–297.

*Light*, Princeton Univ.

*Phys. Rev. Lett.* 93: 023904.

*Soc. Am. A* 13: 1870–1876.

*J. Opt. Soc. Am. A* 14: 2758–2767.

anisotropic materials, *J. Mod. Opt.* 45: 1313–1334.

on light, *Nature (London)* 386: 143–149.

crystals, *Phys. Rev. A* 80: 013826.

microstructured aperture: electromagnetic theory and applications, *Opt. Commun.*

application of the fast Fourier factorization method to cylindrical devices with

modal analysis of slanted lamellar diffraction gratings in conical mountings, *Opt.*

method for two-dimensional photonic crystals with arbitrary hole-shape, *Phys. Rev.*

C. M. (2009). Modes of symmetric composite defects in two-dimensional photonic

birefringence in two-dimensional silicon photonic crystals, *Phys. Rev. B* 63: 161101(R).

Finally, the matrix of permittivity in the reciprocal space with correct Fourier factorization, directly applicable to the wave Equation 1, becomes

$$\begin{array}{c} \begin{bmatrix} \mathbf{e} \end{bmatrix} \mathbf{C} = \begin{bmatrix} \mathbf{F} \end{bmatrix} \begin{bmatrix} \begin{bmatrix} \frac{1}{\varepsilon} \end{bmatrix} \end{bmatrix}^{-1} \begin{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \\\begin{bmatrix} \begin{bmatrix} \varepsilon \end{bmatrix} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \end{array} \begin{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \\\begin{bmatrix} \mathbf{0} \end{bmatrix} \end{bmatrix} \begin{bmatrix} \mathbf{F}^{-1} \end{bmatrix} \\\begin{bmatrix} \mathbf{0} \end{bmatrix} \end{array} \begin{bmatrix} \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{0} \end{bmatrix} \end{array} \tag{101}$$

because the first element on the diagonal corresponds to the normal (*u*) components of the fields and the other two correspond to the tangential components (*v*, *w*) of the fields.

#### **8. Conclusion**

We have derived the methodology how to apply the Fourier factorization rules of Li (1996) to various photonic crystals. For the case of 1D crystals there is clear consistency of the correct rules with the classical theory of Yeh et al. (1977). For 2D crystals the convergence properties strongly depend on the chosen distribution of the polarization basis; we have shown that it is desirable to choose a distribution as smooth as possible. The method is also usable for periodic elements of any shape, where complicated shapes require complicated distributions of polarization bases. We can also use it to simulate 2D periodic elements made of anisotropic materials, as well as 3D periodic crystals. Moreover, the method can also be used to photonic devices such as photonic crystal waveguides by applying the demonstrated methodology to the device supercell. It is particularly advantageous for devices where high accuracy is required, e.g., for analyzing defect modes near photonic band edges [Dossou et al. (2009); Mahmoodian et al. (2009)], and for large devices for which the available computer memory enables calculations with only a few Fourier components (photonic crystals fibers with large cladding, or asymmetric 3D crystals).

#### **9. Acknowledgments**

This work is part of the research plan MSM 0021620834 financed by the Ministry of Education of the Czech Republic and was supported by a Marie Curie International Reintegration Grant (no. 224944) within the 7th European Community Framework Programme and by the Grant Agency of the Czech Republic (no. 202/09/P355 and P204/10/P346).

#### **10. References**


24 Will-be-set-by-IN-TECH

Finally, the matrix of permittivity in the reciprocal space with correct Fourier factorization,

because the first element on the diagonal corresponds to the normal (*u*) components of the

We have derived the methodology how to apply the Fourier factorization rules of Li (1996) to various photonic crystals. For the case of 1D crystals there is clear consistency of the correct rules with the classical theory of Yeh et al. (1977). For 2D crystals the convergence properties strongly depend on the chosen distribution of the polarization basis; we have shown that it is desirable to choose a distribution as smooth as possible. The method is also usable for periodic elements of any shape, where complicated shapes require complicated distributions of polarization bases. We can also use it to simulate 2D periodic elements made of anisotropic materials, as well as 3D periodic crystals. Moreover, the method can also be used to photonic devices such as photonic crystal waveguides by applying the demonstrated methodology to the device supercell. It is particularly advantageous for devices where high accuracy is required, e.g., for analyzing defect modes near photonic band edges [Dossou et al. (2009); Mahmoodian et al. (2009)], and for large devices for which the available computer memory enables calculations with only a few Fourier components (photonic crystals fibers with large

This work is part of the research plan MSM 0021620834 financed by the Ministry of Education of the Czech Republic and was supported by a Marie Curie International Reintegration Grant (no. 224944) within the 7th European Community Framework Programme and by the Grant

Antos, R. (2009). Fourier factorization with complex polarization bases in modeling optics of

Antos, R. & Veis, M. (2010). Fourier factorization with complex polarization bases in the

Benisty, H. (2009). Dark modes, slow modes, and coupling in multimode systems, *J. Opt. Soc.*

Block, I. D., Ganesh, N., Lu, M. & Cunningham, B. T. (2008). Sensitivity model for predicting

Bonod, N., Popov, E. & Neviere, M. (2005a). Fourier factorization of nonlinear Maxwell

plane-wave expansion method applied to two-dimensional photonic crystals, *Opt.*

equations in periodic media: application to the optical Kerr effect, *Opt. Commun.*

discontinuous bi-periodic structures, *Opt. Express* 17: 7269–7274.

photonic crystal biosensor performance, *IEEE Sens. J.* 8: 274–280.

Agency of the Czech Republic (no. 202/09/P355 and P204/10/P346).

fields and the other two correspond to the tangential components (*v*, *w*) of the fields.

<sup>−</sup><sup>1</sup> [[0]] [[0]] [[0]] [[*ε*]] [[0]] [[0]] [[0]] [[*ε*]]

⎤

<sup>⎦</sup> [[**F**−1]], (101)

⎡ ⎣ [[ 1 *ε* ]]

directly applicable to the wave Equation 1, becomes

**8. Conclusion**

cladding, or asymmetric 3D crystals).

*Express* 18: 27511–27524.

*Am. B* 26: 718–724.

244: 389–398.

**9. Acknowledgments**

**10. References**

[[ε]]<sup>C</sup> = [[**F**]]


26 Will-be-set-by-IN-TECH

344 Photonic Crystals – Introduction, Applications and Theory

Lin, S.-Y., Chow, E., Hietala, V., Villeneuve, P. R. & Joannopoulos, J. D. (1998). Experimental

Mahmoodian, S., McPhedran, R. C., de Sterke, C. M., Dossou, K. B., Poulton, C. G. &

Meade, R. D., Rappe, A. M., Brommer, K. D., Joannopoulos, J. D. & Alerhand, O. L.

Noda, S. (2006). Recent progresses and future prospects of two- and three-dimensional

Popov, E. & Neviere, M. (2000). Grating theory: new equations in Fourier space leading to fast converging results for TM polarization, *J. Opt. Soc. Am. A* 17: 1773–1784. Reyes, E., Krokhin, A. A. & Roberts, J. (2005). Effective dielectric constants of photonic crystal

Schuster, T., Ruoff, J., Kerwien, N., Rafler, S. & Osten, W. (2007). Normal vector method for

Skivesen, N., Tetu, A., Kristensen, M., Kjems, J., Frandsen, L. H. & Borel, P. I. (2007).

Sozuer, H. S., Haus, J. W. & Inguva, R. (1992). Photonic bands: Convergence problems with

Vukusic, P. & Sambles, J. R. (2003). Photonic structures in biology, *Nature (London)*

Watanabe, K. (2002). Numerical integration schemes used on the differential theory for

Watanabe, K., Petit, R. & Neviere, M. (2002). Differential theory of gratings made of

Yablonovitch, E. (1987). Inhibited spontaneous emission in solid-state physics and electronics,

Yariv, A. & Yeh, P. (1977). Electromagnetic propagation in periodic stratified media. II. Birefringence, phase matching, and x-ray lasers, *J. Opt. Soc. Am.* 67: 438–447. Yeh, P., Yariv, A. & Hong, C.-S. (1977). Electromagnetic propagation in periodic stratified

Photonic-crystal waveguide biosensor, *Opt. Express* 15: 3169–3176.

photonic crystal band gaps, *Phys. Rev. A* 79: 013814.

photonic crystals, *J. Lightwave Technol.* 24: 4554–4567.

the plane-wave method, *Phys. Rev. B* 45: 13962–13972.

anisotropic gratings, *J. Opt. Soc. Am. A* 19: 2245–2252.

anisotropic materials, *J. Opt. Soc. Am. A* 19: 325–334.

media. I. General theory, *J. Opt. Soc. Am.* 67: 423–438.

crystal, *Science* 282: 274–276.

nanotubes, *Phys. Rev. B* 72: 155118.

*Phys. Rev. Lett.* 58: 2059–2062.

48: 8434–8437.

24: 2880–2890.

424: 852–855.

demonstration of guiding and bending of electromagnetic waves in a photonic

Botten, L. C. (2009). Single and coupled degenerate defect modes in two-dimensional

(1993). Accurate theoretical analysis of photonic band-gap materials, *Phys. Rev. B*

of aligned anisotropic cylinders and the optical response of a periodic array of carbon

convergence improvement using the RCWA for crossed gratings, *J. Opt. Soc. Am. A*
