**1. Introduction**

A photonic crystal (PC) is a periodic structure whose refraction index of the material is periodically modulated on the wavelength scale to affect the electromagnetic wave propagation by creating photonic band gaps. In 1887, Lord Rayleigh is the first to show a band gap in one-dimensional periodic structures *i.e.* a Bragg mirror. In 1987, Eli Yablonovitch and Sajeev John have extended the band gap concept to the two and threedimensional structures and for the first time, they use the term "photonic crystal" (Yablonovitch, 1987; John, 1987).

Progress in computational methods for the photonic crystals is understood through an historical review (Oyhenart, 2005). At the beginning of research in the photonic crystals, the purpose was to find a structure with complete band gap by improving the computational methods. In 1988, John shows theoretically by the scalar method of Korringa-Kohn-Rostoker (KKR) that the face centered cubic lattice (FCC) has a complete band gap between the second and the third band. One year later, Yablonovitch builds this structure and finds a band gap experimentally but the W-point raises a problem. In 1990, Satpathy et al. and Leung et al. confirm the complete band gap by the scalar plane wave method (PWM). A few months later, these two teams improve their methods to obtain vectorial PWM on **D** and **E**  fields. They find that FCC structure does not have complete band gap because W-point and U-point are degenerate. With these results, the editor of the journal "Nature" writes *"Photonic Crystals bite the dust"* (Maddox, 1990). Only two weeks later, Ho et al. created the vectorial PWM on **H** and they do not find the complete band gap in FCC structure but they show a complete band gap in the diamond lattice. In 1992, Sözuer et al. improve convergence of the PWM and they obtain a complete band gap for FCC lattice between 8th and 9th band. This structure that has caused many discussions has a complete band gap but not where it was expected.

To study and understand the propagation of the electromagnetic fields in the photonic crystals, computational methods were improved by using their symmetries and periodicities. We will study the classical methods for microwave devices such as the finite element method and the finite difference time domain. After some modifications of these methods, we obtain the band structure of PC which can be calculated by the methods from the solid state physics. For example the plane wave method, the tight binding method and the multiple-scattering theory will be studied. All these computational

Overview of Computational Methods for Photonic Crystals 269

incident wave. In the first case, we apply a TM incident wave on the lattice with a horizontal symmetry and the tangential magnetic field vanishes on the axis of symmetry. We put perfect magnetic conductor (PMC) as boundary condition on the symmetry axis and we study only half top or bottom of the problem. In the second case, we apply a TE incident wave. Similarly, we reduce the problem with a perfect electric conductor condition (PEC).

On figure 1, the geometry has also a vertical symmetry. However, the electromagnetic field is not symmetrical because the incident wave comes only from the left side. A solution is to divide the incident wave into an even mode and an odd mode (figure 2). For the even mode, we have a vertical symmetry plane of the electric field. Only the left or right part is solved with a perfect magnetic conductor condition on the axis of symmetry. The antisymmetric problem is reduced in the same way with a perfect electric conductor condition. The sum of the symmetric and antisymmetric problem provides the solution of the total problem. The reflection of total problem is the sum *reven* + *rodd* and the transmission is the subtraction *reven rodd*. Two half-problems are solved more quickly than the whole problem because the

Photonic crystals, like the familiar crystals of atoms have discrete translational symmetry **T**. The dielectric function is periodic therefore the electric and magnetic fields can be written as

( ) ( ) ( ) ( ) *<sup>i</sup>*

**Fr ur ur ur T** (5)

**K r**

. with

*n*

() ( )

**r rT**

*e*

computation time increases exponentially with the size of the problem.

Fig. 1. Both kinds of lateral symmetries

Fig. 2. Transverse symmetries of the problem

the product of a plane wave envelope and a periodic function:

 

**2.3 Periodicities of photonic crystals** 

methods will be presented in this article and a PC example will be studied to compare these methods.

#### **2. Equations, symmetries and periodicities in photonic crystals**

#### **2.1 Equations in photonic crystals**

The Maxwell equations without sources control the electromagnetic wave propagation in PC.

$$\begin{aligned} \nabla \times \mathbf{E} \left( \mathbf{r}, t \right) + \frac{\partial \mathbf{B} \left( \mathbf{r}, t \right)}{\partial t} &= \mathbf{0} & \nabla \times \mathbf{H} \left( \mathbf{r}, t \right) - \frac{\partial \mathbf{D} \left( \mathbf{r}, t \right)}{\partial t} &= \mathbf{0} \\ \nabla \cdot \mathbf{D} \left( \mathbf{r}, t \right) &= 0 & \nabla \cdot \mathbf{B} \left( \mathbf{r}, t \right) &= 0 \end{aligned} \tag{1}$$

In these laws, the same physical behavior is observed if we change simultaneously the wavelength and the structure dimensions in the same proportions. Therefore, it is convenient to introduce a normalized wavelength λ0/*a* and a normalized frequency *af*/*c*= *a*/λ0, with *a* the lattice constant of the photonic crystal (Joannopoulos et al., 2008).

Some methods do not solve the Maxwell equations directly but they use the Helmholtz equations, for example the **E**-wave equation or the **H**-wave equation:

$$\begin{aligned} \nabla \times \left[ \nabla \times \mathbf{E} \left( \mathbf{r} \right) \right] &= -\nabla^2 \mathbf{E} \left( \mathbf{r} \right) = k\_0^2 \varepsilon\_r \left( \mathbf{r} \right) \mathbf{E} \left( \mathbf{r} \right) \\ \nabla \times \left[ \varepsilon\_r^{-1} \left( \mathbf{r} \right) \nabla \times \mathbf{H} \left( \mathbf{r} \right) \right] &= k\_0^2 \mathbf{H} \left( \mathbf{r} \right) \end{aligned} \quad \text{with} \quad \begin{cases} \mathbf{E} \left( \mathbf{r}, t \right) = \mathbf{E} \left( \mathbf{r} \right) e^{-i \alpha t} \\ \mathbf{H} \left( \mathbf{r}, t \right) = \mathbf{H} \left( \mathbf{r} \right) e^{-i \alpha t} \end{cases} \tag{2}$$

Equations number to be solved depends on dimension of the photonic crystal. In the twodimensional case, the problem is simplified. It is assumed that the materials are uniform along *z*-axis. It follows that the fields are uniform and the partial derivatives with respect to the variable z vanish. The previous equation is simplified and split-up into TE-polarization and TM-polarization; we have a scalar equation for each polarization:

$$\frac{\partial^2 F\left(\mathbf{r}\right)}{\partial \mathbf{x}^2} + \frac{\partial^2 F\left(\mathbf{r}\right)}{\partial y^2} + k\_o^2 \varepsilon\_r\left(\mathbf{r}\right) F\left(\mathbf{r}\right) = 0 \quad \text{with} \begin{cases} F = E\_z \text{ for T\!i-polarization} \\ F = H\_z \text{ for TM-polarization} \end{cases} \tag{3}$$

The computation time decreases exponentially compared to the three-dimensional case. In the one-dimensional case, the problem is even more simplified. If the materials are uniform along *x*-axis and *y*-axis, we solve analytically one second-order equation:

$$\frac{d^2F\left(\mathbf{r}\right)}{dz^2} + k\_0^2 \varepsilon\_r\left(\mathbf{r}\right)F\left(\mathbf{r}\right) = 0 \quad \text{with} \quad F = E\_{x'}, E\_{y'}, H\_x \text{ or } H\_y \tag{4}$$

#### **2.2 Symmetries of photonic crystals**

Like all sets of differential equations, Maxwell's equations cannot be uniquely solved without a suitable set of boundary conditions. Photonic crystals have symmetries which define boundary conditions. Symmetries of the structure are not a sufficient condition to reduce the computational domain; the electromagnetic field must be also symmetrical. On figure 1, we study symmetries on a lattice of cylinders for different polarization of the

methods will be presented in this article and a PC example will be studied to compare

The Maxwell equations without sources control the electromagnetic wave propagation

 

, , , ,

*t t*

In these laws, the same physical behavior is observed if we change simultaneously the wavelength and the structure dimensions in the same proportions. Therefore, it is convenient to introduce a normalized wavelength λ0/*a* and a normalized frequency *af*/*c*= *a*/λ0, with *a* the lattice constant of the photonic crystal (Joannopoulos et al., 2008).

Some methods do not solve the Maxwell equations directly but they use the Helmholtz

Equations number to be solved depends on dimension of the photonic crystal. In the twodimensional case, the problem is simplified. It is assumed that the materials are uniform along *z*-axis. It follows that the fields are uniform and the partial derivatives with respect to the variable z vanish. The previous equation is simplified and split-up into TE-polarization

The computation time decreases exponentially compared to the three-dimensional case. In the one-dimensional case, the problem is even more simplified. If the materials are uniform

Like all sets of differential equations, Maxwell's equations cannot be uniquely solved without a suitable set of boundary conditions. Photonic crystals have symmetries which define boundary conditions. Symmetries of the structure are not a sufficient condition to reduce the computational domain; the electromagnetic field must be also symmetrical. On figure 1, we study symmetries on a lattice of cylinders for different polarization of the

with

0 *<sup>z</sup>*

 

*z*

**<sup>r</sup> r r** with , , or *xy x y FEEH H* (4)

**r r r r** (3)

**D r B r**

, 0 , 0

**B r D r E r 0 Hr 0**

> *t t t t*

*t t*

with

 

, , *t t*

**Er Er Hr Hr**

for TE-polarization

for TM-polarization

  *i t i t*

*e e*

 

(1)

(2)

**2. Equations, symmetries and periodicities in photonic crystals** 

equations, for example the **E**-wave equation or the **H**-wave equation:

 

**r Hr Hr**

and TM-polarization; we have a scalar equation for each polarization:

2

*k*

*r*

along *x*-axis and *y*-axis, we solve analytically one second-order equation:

*k F*

2 <sup>0</sup> <sup>2</sup> <sup>0</sup> *<sup>r</sup>*

 

*F F F E*

*x y F H F* 

0 2 2

1 2

 

*r*

2 2

2

*d F*

**2.2 Symmetries of photonic crystals** 

*dz*

2 2 0

**Er Er rEr**

0

*r*

*k k*

these methods.

in PC.

**2.1 Equations in photonic crystals** 

incident wave. In the first case, we apply a TM incident wave on the lattice with a horizontal symmetry and the tangential magnetic field vanishes on the axis of symmetry. We put perfect magnetic conductor (PMC) as boundary condition on the symmetry axis and we study only half top or bottom of the problem. In the second case, we apply a TE incident wave. Similarly, we reduce the problem with a perfect electric conductor condition (PEC).

Fig. 1. Both kinds of lateral symmetries

On figure 1, the geometry has also a vertical symmetry. However, the electromagnetic field is not symmetrical because the incident wave comes only from the left side. A solution is to divide the incident wave into an even mode and an odd mode (figure 2). For the even mode, we have a vertical symmetry plane of the electric field. Only the left or right part is solved with a perfect magnetic conductor condition on the axis of symmetry. The antisymmetric problem is reduced in the same way with a perfect electric conductor condition. The sum of the symmetric and antisymmetric problem provides the solution of the total problem. The reflection of total problem is the sum *reven* + *rodd* and the transmission is the subtraction *reven rodd*. Two half-problems are solved more quickly than the whole problem because the computation time increases exponentially with the size of the problem.

Fig. 2. Transverse symmetries of the problem

#### **2.3 Periodicities of photonic crystals**

Photonic crystals, like the familiar crystals of atoms have discrete translational symmetry **T**. The dielectric function is periodic therefore the electric and magnetic fields can be written as the product of a plane wave envelope and a periodic function:

$$\begin{aligned} \boldsymbol{\varepsilon}(\mathbf{r}) &= \boldsymbol{\varepsilon}(\mathbf{r} + n\mathbf{T}) \\ \mathbf{F}(\mathbf{r}) &= e^{i\mathbf{K}\cdot\mathbf{r}} \mathbf{u}(\mathbf{r}) \quad \text{with } \mathbf{u}(\mathbf{r}) = \mathbf{u}(\mathbf{r} + \mathbf{T}) \end{aligned} \tag{5}$$

Overview of Computational Methods for Photonic Crystals 271

Let us consider a partial differential equation with the differential operator of order *n*

1 *N j j j c* 

are the chosen expansion functions and *<sup>j</sup> <sup>c</sup>* are constant coefficients to be determined. The

1

If we use the Galerkin method, the weight functions *wi* are the functions of previous

 

So now, we have to solve a linear system with *cj* unknowns where most of the entries of the matrix *L* are zero. Such matrices are known as sparse matrices, there are efficient solvers for

The basic idea of the finite element method is to divide the computation domain into small subdomains, which are called finite elements, and then use simple functions, such as linear and quadratic functions, to approximate the unknown solution over each element. For plane geometries, the domain is divided into finite triangular sub-domains. For three-dimensional problems, the sub-domains are tetrahedra. These two- and three-dimensional finite elements

> <sup>2</sup> <sup>0</sup> *<sup>r</sup>* **Er rEr** *k*

For most of the finite elements computer programs, the frequency is fixed and the electric field is the unknown (Massaro et al., 2008). We use a commercial software, Ansys HFSS (High Frequency Structure Simulator). This three-dimensional computer program builds an adaptive tetrahedral mesh to model microwave devices, for example microstrips and antennas. These devices have a characteristic length lower than the wavelength of study. It is more difficult to study PC because it has a periodicity close to the wavelength and the matrix of calculation is large. To simplify calculations, we will study a periodic PC according to two directions of space. The directions where PC is infinitely periodic require

with *ij i j L d*

are widely used because it is a variable mesh and adapts to curved structures.

The finite element method usually solves the **E**-wave equation for PC (equation 5):

*N i i ji j j R w rd c w f d*

and *i i <sup>b</sup>*

*f* (6)

(7)

is weakest on all points

(9)

(10)

0

 *fd* 

(8)

applied to a function *φ* and a source function *f*:

*j* 

such problems.

The first step is to expand the function *φ* on a set of functions:

best solution of the equation 6 is obtained when the residue *r f*

interpolations and we write the problem in matrix form:

*j ij i*

*cL b*

1 *N*

**3.1 Transmission calculation of a photonic crystal** 

*j*

of the domain . The weighted residual method requires this condition:

The result above is commonly known as Bloch's theorem in solid state physics (Kittel, 2005).

If we apply these conditions to an infinite lattice of cylinders (figure 3a), the calculation of band structure is reduced to the study of a single cylinder.

Fig. 3. Periodic boundary on photonic crystals

Figure 3b is a lattice of cylinders, finite along *y*-axis and infinite along *x*-axis for studying the transmission of the structure. If we know the transfer matrix of one layer, it is easy to find the transmission of the total structure thanks to the transfer matrices method (TMM). It reduces the computational domain to one layer. It will be detailed in the next sections.

By mixing TMM along *y*-axis and Bloch's conditions along *x*-axis, the structure of figure 3 can be reduced to a single cylinder. Figure 4 summarizes these techniques in the 3D case.

Fig. 4. Reduction of the calculation domain in PC
