**3. Finite element method (FEM)**

The finite element method is very popular in mechanics and civil engineering. It was originally developed in the 40s to solve problems of mechanical structures. A few years later, it has been applied to electromagnetism. Since then, the finite element method extends to all branches of physics and engineering where there exists a partial differential equation (PDE) with boundary conditions. It can be formulated from the variational method or the weighted residual method. We will present in outline the second method which is simpler.

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

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.

The finite element method is very popular in mechanics and civil engineering. It was originally developed in the 40s to solve problems of mechanical structures. A few years later, it has been applied to electromagnetism. Since then, the finite element method extends to all branches of physics and engineering where there exists a partial differential equation (PDE) with boundary conditions. It can be formulated from the variational method or the weighted residual method. We will present in outline the second method

band structure is reduced to the study of a single cylinder.

Fig. 3. Periodic boundary on photonic crystals

Fig. 4. Reduction of the calculation domain in PC

**3. Finite element method (FEM)** 

which is simpler.

Let us consider a partial differential equation with the differential operator of order *n* applied to a function *φ* and a source function *f*:

$$
\mathfrak{T}\varphi = f \tag{6}
$$

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

$$\varphi = \sum\_{j=1}^{N} c\_j \nu\_j \tag{7}$$

*j* are the chosen expansion functions and *<sup>j</sup> <sup>c</sup>* are constant coefficients to be determined. The best solution of the equation 6 is obtained when the residue *r f* is weakest on all points of the domain . The weighted residual method requires this condition:

$$R\_i = \int\_{\Omega} w\_i r d\Omega = \sum\_{j=1}^{N} c\_j \int\_{\Omega} w\_i \left(\mathfrak{T}\,\nu\_j - f\right) d\Omega = 0 \tag{8}$$

If we use the Galerkin method, the weight functions *wi* are the functions of previous interpolations and we write the problem in matrix form:

$$\sum\_{j=1}^{N} c\_j L\_{ij} = b\_i \quad \text{with} \quad L\_{ij} = \underset{\Omega}{\int} \nu\_i \mathfrak{Z} \,\nu\_j d\Omega \quad \text{and} \quad b\_i = \underset{\Omega}{\int} \nu\_i f d\Omega \tag{9}$$

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 such problems.

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 are widely used because it is a variable mesh and adapts to curved structures.

#### **3.1 Transmission calculation of a photonic crystal**

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

$$\nabla \times \left[\nabla \times \mathbf{E}\left(\mathbf{r}\right)\right] = k\_0^2 \varepsilon\_r\left(\mathbf{r}\right) \mathbf{E}\left(\mathbf{r}\right) \tag{10}$$

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

Overview of Computational Methods for Photonic Crystals 273

constant). For the band structure, the number of layer is infinite, and for transmission calculation, we study four and eight layers. The structure is deliberately simple to be studied by all the methods. The band structure, the transmission coefficient in normal incidence and some modes of dielectric PC are plotted on figure 7. For normal incidence, the

The first method with incident waves takes about 7 hours and 960 Mo of working memory on a personal computer to calculate the transmission of 8 layers *i.e.* 8 spheres. Whereas if we use the second method with the wave ports, the memory is reduced to 360 Mo and the CPU time is reduced to 8 min. To obtain this optimization, we reduced the geometry to four quarter-spheres thanks to symmetries and we use the Padé interpolation on frequencies. It is

first two band gaps are found in the band structure and transmission curves.

necessary to make optimizations if you want to use FEM.

Fig. 7. Band structure and transmission coefficient of dielectric PC

Maxwell in a conducting, isotropic and homogeneous medium:

Finite difference method is a numerical method for approximating differential equations using finite difference equations to approximate derivatives. This is very simple to implement but it has some limitations in the mesh geometries. In 1966, Yee proposed a finite difference scheme applied to electromagnetism. The FDTD was born (Taflove & Hagness, 2005). The success of this method is due to the scheme based on the Taylor series of second-order. For example, this method is used to model the effects of cellphones on the human body, antennas and printed circuits. The FDTD uses the two structural equations of

1

*t µ*

*t*

**H**

**E**

 

1

**H E**

**E**

(11)

**4. Finite-difference time-domain method (FDTD)** 

to study only one period according to this direction. Several methods exist to model these conditions.

In the first method, the source is a plane wave and we use periodic boundary conditions (PBC) on lateral faces and absorbing boundary conditions (ABC) on bases. We study only one lateral period (figure 5). If the source has an oblique incidence, the phase shift is easily taken into account by PBC because the E-field is a complex vector in FEM.

Fig. 5. Calculation of E-field in the PC with PBC, ABC and incident wave

The second method requires a very different source, a wave port. This source is a semiinfinite waveguide whose cross section is drawn on bases of structure. Propagative modes of this fictitious guide will be the source of the structure. For a plane wave source, we apply to lateral faces and on the fictitious guide the conditions PEC and PMC (figure 6). The source is transverse and has a normal incidence. It is not possible to change the angle of incidence without change the boundary conditions of the lateral faces. The PC studied on figure 5 can be reduced to quarter-spheres thanks to the symmetry of the fields at normal incidence.

Fig. 6. Calculation of E-field in the PC with PEC, PMC and wave ports

## **3.2 Band structure calculation**

The band structure shows the states which propagate in a PC. These states are differentiated by their frequency and their Bloch wave vector. FEM sets the wave vector and solves the wave equation to find the frequencies. There is no source, only boundary conditions to set because it is an eigenvalue problem. In the case of a cavity resonator, the boundary conditions of the domain are PEC. Whereas, for PC, we choose the Bloch conditions on all faces of the unit cell to set the wave vector. The phase shift of the Bloch conditions is set easily because the fields are complex vectors. The FEM calculates the band structure of dielectric with or without losses, metallic, and metallodielectric PC. Any material can be used by this method, it is the main advantage.

In this chapter, we choose a photonic crystal to study, a cubic lattice with several layers of dielectric spheres which have a permittivity equal to 5.1 and radius equal to 0.4\**a* (*a*: lattice

to study only one period according to this direction. Several methods exist to model these

In the first method, the source is a plane wave and we use periodic boundary conditions (PBC) on lateral faces and absorbing boundary conditions (ABC) on bases. We study only one lateral period (figure 5). If the source has an oblique incidence, the phase shift is easily

The second method requires a very different source, a wave port. This source is a semiinfinite waveguide whose cross section is drawn on bases of structure. Propagative modes of this fictitious guide will be the source of the structure. For a plane wave source, we apply to lateral faces and on the fictitious guide the conditions PEC and PMC (figure 6). The source is transverse and has a normal incidence. It is not possible to change the angle of incidence without change the boundary conditions of the lateral faces. The PC studied on figure 5 can be reduced to quarter-spheres thanks to the symmetry of the fields at normal

PBC on lateral faces

ABC k0

The band structure shows the states which propagate in a PC. These states are differentiated by their frequency and their Bloch wave vector. FEM sets the wave vector and solves the wave equation to find the frequencies. There is no source, only boundary conditions to set because it is an eigenvalue problem. In the case of a cavity resonator, the boundary conditions of the domain are PEC. Whereas, for PC, we choose the Bloch conditions on all faces of the unit cell to set the wave vector. The phase shift of the Bloch conditions is set easily because the fields are complex vectors. The FEM calculates the band structure of dielectric with or without losses, metallic, and metallodielectric PC. Any material can be

**Wave port2** 

In this chapter, we choose a photonic crystal to study, a cubic lattice with several layers of dielectric spheres which have a permittivity equal to 5.1 and radius equal to 0.4\**a* (*a*: lattice

taken into account by PBC because the E-field is a complex vector in FEM.

**Source: incident wave** 

Fig. 5. Calculation of E-field in the PC with PBC, ABC and incident wave

Fig. 6. Calculation of E-field in the PC with PEC, PMC and wave ports

 PEC on up and down faces PMC on left and right faces

Lateral boundary conditions:

conditions.

AB

**E** 

incidence.

**Source Wave port1** 

**3.2 Band structure calculation** 

used by this method, it is the main advantage.

constant). For the band structure, the number of layer is infinite, and for transmission calculation, we study four and eight layers. The structure is deliberately simple to be studied by all the methods. The band structure, the transmission coefficient in normal incidence and some modes of dielectric PC are plotted on figure 7. For normal incidence, the first two band gaps are found in the band structure and transmission curves.

The first method with incident waves takes about 7 hours and 960 Mo of working memory on a personal computer to calculate the transmission of 8 layers *i.e.* 8 spheres. Whereas if we use the second method with the wave ports, the memory is reduced to 360 Mo and the CPU time is reduced to 8 min. To obtain this optimization, we reduced the geometry to four quarter-spheres thanks to symmetries and we use the Padé interpolation on frequencies. It is necessary to make optimizations if you want to use FEM.

Fig. 7. Band structure and transmission coefficient of dielectric PC

#### **4. Finite-difference time-domain method (FDTD)**

Finite difference method is a numerical method for approximating differential equations using finite difference equations to approximate derivatives. This is very simple to implement but it has some limitations in the mesh geometries. In 1966, Yee proposed a finite difference scheme applied to electromagnetism. The FDTD was born (Taflove & Hagness, 2005). The success of this method is due to the scheme based on the Taylor series of second-order. For example, this method is used to model the effects of cellphones on the human body, antennas and printed circuits. The FDTD uses the two structural equations of Maxwell in a conducting, isotropic and homogeneous medium:

$$\begin{aligned} \frac{\partial \mathbf{E}}{\partial t} &= \frac{1}{\mathcal{E}} \nabla \times \mathbf{H} - \frac{\sigma}{\mathcal{E}} \mathbf{E} \\ \frac{\partial \mathbf{H}}{\partial t} &= -\frac{1}{\mu} \nabla \times \mathbf{E} \end{aligned} \tag{11}$$

Overview of Computational Methods for Photonic Crystals 275

The band structure is calculated from the eigenmodes but FDTD is not an eigenvalue problem. We use the unit cell of PC with boundary conditions of Bloch and several Gaussian functions for the source. We only need one calculation to apply the Bloch conditions. After 100000 time-steps of calculation, propagative waves are amplified and evanescent waves vanished. If the number of time steps is too small, the transmission peaks are widened, therefore imprecise. To reduce the number of time-steps without affecting the accuracy, we can use the Padé approximation or ADI-FDTD formulation (Taflove & Hagness, 2005). On Figure 10, we plot the electric field amplitude of the unit

0 0.2 0.4 0.6 0.8 1

af/c

Fig. 10. Electric field amplitude of dielectric PC at the R-point of the Brillouin zone

This curve is similar to the diffraction pattern in solid state physics. It is plotted for all points of the first Brillouin zone to view the band structure. If we compare the curve of figure 10 with the band structure calculated by the FEM, the mode 0.57 is much attenuated and the following mode does not exist on figure 10. In fact, the source was not correctly selected to excite these modes, and so we must ensure that the source excites all

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

af/c

Fig. 9. Transmission coefficient of the dielectric PC

4 layers 8 layers

**4.2 Calculation of band structure** 

cell of PC.

available modes.

5

10

Amplitude

15

20

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Transmission

The FDTD uses an approximation of derivatives by centered finite differences.

$$\frac{d}{d\mathbf{x}\_i}\mu\left(\mathbf{x}\_i\right) = \frac{\mu\left(\mathbf{x}\_i + h\right) - \mu\left(\mathbf{x}\_i - h\right)}{2h} + o\left(h^2\right) \tag{12}$$

Let us apply this approximation to one-dimensional case in order to understand the FDTD principle. The **E** and **H**-field are stepped in time and space. By replacing curls and derivatives, we obtain one-dimensional Yee scheme:

$$\begin{aligned} \frac{E\_k^{n+1/2} - E\_k^{n-1/2}}{\Delta t} &= -\frac{1}{\mathcal{E}} \frac{H\_{k+1/2}^n - H\_{k-1/2}^n}{\Delta \mathbf{x}}\\ \frac{E\_{k+1/2}^{n+1} - E\_{k+1/2}^n}{\Delta t} &= -\frac{1}{\mu} \frac{H\_{k+1}^{n+1/2} - H\_k^{n+1/2}}{\Delta \mathbf{x}} \end{aligned} \quad \text{with} \quad \begin{aligned} E\_k^u &= \text{E}\{n.\Delta t, k\Delta \mathbf{x}\} = \text{E}\{t, \mathbf{x}\} \\ H\_k^u &= \text{H}\{n.\Delta t, k\Delta \mathbf{x}\} = \text{H}\{t, \mathbf{x}\} \end{aligned} \tag{13}$$

The **E** and **H**-field are staggered and updated step by step in time. **E**-field updates are conducted midway during each time-step between successive **H**-field updates, and conversely. This explicit time-stepping scheme avoids the need to solve simultaneous equations, and furthermore it is order *N i.e.* proportional to the size of the system to model. This scheme can be generalized for two-dimensional and three-dimensional problems. The Yee scheme is stable if the wave propagates from one cell to another with a speed less than the light (the Current-Friedrichs-Lewy condition).

#### **4.1 Calculation of transmission coefficient**

FDTD is used to study PC. Calculation domain of finite PC is surrounded by absorbing boundary conditions (ABC). The periodic infinite PC uses periodic boundary conditions (PBC) on the lateral faces (Figure 8). The source is Gaussian with a spectrum which extends on the frequency range to study. The fields are calculated on the time domain and we use the Fourier transform to convert them on frequency domain.

Fig. 8. Source, probe and boundary conditions for infinite PC

For the non-normal incidence, a simple periodic boundary condition cannot be applied but many methods exist to solve this problem. One method is to calculate the computing problem twice, with the sources cos(t) and sin(t). The addition of two calculations is the solution of the source exp(-it). Thanks to this source, we can apply the Bloch conditions on complex vector fields. This method is simplified for the study of band structure.

Figure 9 plot the transmission coefficient for the previous dielectric PC. The transmission calculation of 8 layers is calculated in 32 minutes with 58 Mo of working memory. FDTD is faster than FEM but less accurate.

 <sup>2</sup> 2 *i i*

with

*n k n k*

(12)

( . , ) (, ) ( . , ) (, )

(13)

*E En t k x Et x H Hn tk x Htx* 

PC Probe PBC

*<sup>d</sup> ux h ux h u x <sup>h</sup>*

Let us apply this approximation to one-dimensional case in order to understand the FDTD principle. The **E** and **H**-field are stepped in time and space. By replacing curls and

The **E** and **H**-field are staggered and updated step by step in time. **E**-field updates are conducted midway during each time-step between successive **H**-field updates, and conversely. This explicit time-stepping scheme avoids the need to solve simultaneous equations, and furthermore it is order *N i.e.* proportional to the size of the system to model. This scheme can be generalized for two-dimensional and three-dimensional problems. The Yee scheme is stable if the wave propagates from one cell to another with a speed less than

FDTD is used to study PC. Calculation domain of finite PC is surrounded by absorbing boundary conditions (ABC). The periodic infinite PC uses periodic boundary conditions (PBC) on the lateral faces (Figure 8). The source is Gaussian with a spectrum which extends on the frequency range to study. The fields are calculated on the time domain and we use

For the non-normal incidence, a simple periodic boundary condition cannot be applied but many methods exist to solve this problem. One method is to calculate the computing problem twice, with the sources cos(t) and sin(t). The addition of two calculations is the solution of the source exp(-it). Thanks to this source, we can apply the Bloch conditions on

Figure 9 plot the transmission coefficient for the previous dielectric PC. The transmission calculation of 8 layers is calculated in 32 minutes with 58 Mo of working memory. FDTD is

complex vector fields. This method is simplified for the study of band structure.

ABC ABC

1/2 1/2

The FDTD uses an approximation of derivatives by centered finite differences.

*dx h*

*i i*

1 1/2 1/2

1

1

*n n n n k k k k*

*E E H H t x E E H H t x*

*n n n n k k k k*

 

derivatives, we obtain one-dimensional Yee scheme:

1/2 1/2 1

the light (the Current-Friedrichs-Lewy condition).

the Fourier transform to convert them on frequency domain.

Fig. 8. Source, probe and boundary conditions for infinite PC

faster than FEM but less accurate.

Plane wave

**4.1 Calculation of transmission coefficient** 

1/2 1/2

Fig. 9. Transmission coefficient of the dielectric PC

#### **4.2 Calculation of band structure**

The band structure is calculated from the eigenmodes but FDTD is not an eigenvalue problem. We use the unit cell of PC with boundary conditions of Bloch and several Gaussian functions for the source. We only need one calculation to apply the Bloch conditions. After 100000 time-steps of calculation, propagative waves are amplified and evanescent waves vanished. If the number of time steps is too small, the transmission peaks are widened, therefore imprecise. To reduce the number of time-steps without affecting the accuracy, we can use the Padé approximation or ADI-FDTD formulation (Taflove & Hagness, 2005). On Figure 10, we plot the electric field amplitude of the unit cell of PC.

Fig. 10. Electric field amplitude of dielectric PC at the R-point of the Brillouin zone

This curve is similar to the diffraction pattern in solid state physics. It is plotted for all points of the first Brillouin zone to view the band structure. If we compare the curve of figure 10 with the band structure calculated by the FEM, the mode 0.57 is much attenuated and the following mode does not exist on figure 10. In fact, the source was not correctly selected to excite these modes, and so we must ensure that the source excites all available modes.

Overview of Computational Methods for Photonic Crystals 277

The reflection and transmission calculation is performed using the transfer matrices. It is more interesting to calculate the elements of transmission matrix than the transmitted field in some points. The incident, reflected and transmitted waves are expanded on a plane-

**T0** define a plane-wave

**T** (19)

**T** is not compatible with the new plane-wave

**T0** operator is not self-adjoint. They are

**5.2 Calculation of transmission coefficient** 

wave basis sets thanks to the transfer matrices.

The right and left eigenvectors are different because ˆ

are in a vacuum layer. The transfer matrix ˆ

**6. Finite Integration Technique (FIT)** 

integral form to a set of staggered grids like FDTD (figure 12).

basis set, we convert this matrix:

**T0** is the transfer matrix of a vacuum layer. The eigenvectors of <sup>ˆ</sup>

ˆ

**T**

**0**

with the same propagation direction. This matrix is divided into four blocks:

**T**

ˆ

**0**

*i i ik c i i ik c i i*

*r er l el* 

expanded on the reciprocal space *kx* and *ky* and translate on the direct space with the Fourier transform because the transfer matrix is known in direct space. *kz* is found easily because we

> <sup>ˆ</sup> *ji j i i j ji l r l T r T lr T*

To calculate the transmission coefficient, this matrix is arranged to group the eigenvectors

11 12 21 22 **T T**

**T T**

TMM calculates the transmission coefficient of several layers from one layer (see 1D-MST). A great number of layers can be calculated without increasing the computing

Figure 11 plot the band structure and the transmission coefficient calculated by FDFD, for the previous dielectric PC. If we compares with figure 7, the first bands are correct but inaccuracies on the higher band after 0.6-frequency are due to the weak mesh (7x7x7 cells). If we increase the number of cells, accuracy increases. Calculation is more difficult for high frequencies because the E- and H-field are functions that oscillate more. The transmission calculation of 8 layers is calculated in 22 seconds with 3.3 Mo of working memory for 7x7x7 cells. The computing time is low because only one sphere is actually calculated thanks

In 1977, Weiland proposes a spatial discretization scheme to solve the integral equations of Maxwell (Weiland, 1977). This scheme called Finite Integration Technique (FIT) can be applied to many electromagnetism problems, in time and frequency domain, from static up to high frequency. The basic idea of this approach is to apply the Maxwell equations in

 

   

**<sup>T</sup>** (20)

(21)

ˆ

basis set:

power.

TMM.

#### **5. Finite-difference frequency-domain method (FDFD)**

We will study a finite difference method combined with transfer matrix method and the Fourier transform. The finite difference method substitutes the derivatives in PDE to obtain finite difference schemes. In electromagnetism, FDFD uses the structural equations of Maxwell on space (**k**, ) and applies approximations on the wave vector:

$$\begin{cases} \mathbf{k} \times \mathbf{E} = \alpha \mu\_0 \mathbf{H} \\ \mathbf{k} \times \mathbf{H} = -\alpha \nu \varepsilon \mathbf{E} \end{cases} \text{ with } \begin{cases} k\_x \approx \left( e^{ik\_z a} - 1 \right) / \ i a \\ k\_y \approx \left( e^{ik\_y b} - 1 \right) / \ i b \\ k\_z \approx \left( e^{ik\_z c} - 1 \right) / \ i c \end{cases} \tag{14}$$

*a*,*b*,*c* are lattice constants along *x*, *y* and *z*-axis. We obtain discrete equations which are applied to one layer of the structure on a cubic lattice. In 1992, Pendry and MacKinnon (Pendry & MacKinnon, 1992) used this method with transfer matrix method (TMM) which extended the solution of one layer to the total structure. To apply TMM, the discrete equations must be written on real space. We obtain a set of equations where *z*-component of the **E** and **H**-fields can be removed. The six equations are reduced to four equations which are written in a matrix form:

$$\mathbf{F(r+c)} = \sum\_{i} \hat{\mathbf{T}}(\mathbf{r}, \mathbf{r'}) \mathbf{F(r')} \text{ with } \mathbf{F(r)} = \begin{bmatrix} E\_x(\mathbf{r}) & E\_y(\mathbf{r}) & H\_x(\mathbf{r}) & H\_y(\mathbf{r}) \end{bmatrix}^\mathrm{T} \tag{15}$$

#### **5.1 Calculation of band structure**

The band structure is calculated from unit cell of structure by setting frequency and calculating wave vectors propagating in the unit cell. Bloch's theorem applies to vector **F**:

$$\begin{aligned} \mathbf{F}\left(\mathbf{r} + a'\right) &= e^{ik\_{\circ}a'} \mathbf{F}\left(\mathbf{r}\right) \\ \mathbf{F}\left(\mathbf{r} + b'\right) &= e^{ik\_{\circ}b'} \mathbf{F}\left(\mathbf{r}\right) \\ \mathbf{F}\left(\mathbf{r} + c'\right) &= e^{ik\_{\circ}c'} \mathbf{F}\left(\mathbf{r}\right) \end{aligned} \tag{16}$$

*a ab b c c* ' , ' ' and are the mesh size along *x*, *y* and *z*-axis. The transfer matrix of the unit cell is written from the unit mesh:

$$\mathbf{F}\left(\mathbf{r}+\mathbf{c}'\right) = \sum\_{\mathbf{r}'} \hat{\mathbf{T}}\left(\mathbf{c}',0\right) \mathbf{F}\left(\mathbf{r}'\right) \text{ with } \hat{\mathbf{T}}\left(\mathbf{c}',0\right) = \prod\_{j=1}^{N} \hat{\mathbf{T}}\left(\mathbf{r}\_{j},\mathbf{r}\_{j-1}\right) \tag{17}$$

If the above equations are joined together, we have an eigenvalue problem:

$$\sum\_{\mathbf{r}'} \hat{\mathbf{T}}\left(\mathbf{c}', \mathbf{0}\right) \mathbf{F}\left(\mathbf{r}'\right) = e^{i\mathbf{k}\_1 \mathbf{c}'} \mathbf{F}\left(\mathbf{r}\right) \tag{18}$$

We set the wave number *k* and the wavelength *k0* to solve the eigenvalue problem. The eigenvalues *kz* having an imaginary part are eliminated because they are not propagating waves. The remaining eigenvalues gives us the band structure *kz*(*k0*).

#### **5.2 Calculation of transmission coefficient**

276 Photonic Crystals – Introduction, Applications and Theory

We will study a finite difference method combined with transfer matrix method and the Fourier transform. The finite difference method substitutes the derivatives in PDE to obtain finite difference schemes. In electromagnetism, FDFD uses the structural equations of

0 with

*k e ia*

*x y*

*ik b*

*ik a*

*x*

 

*a*,*b*,*c* are lattice constants along *x*, *y* and *z*-axis. We obtain discrete equations which are applied to one layer of the structure on a cubic lattice. In 1992, Pendry and MacKinnon (Pendry & MacKinnon, 1992) used this method with transfer matrix method (TMM) which extended the solution of one layer to the total structure. To apply TMM, the discrete equations must be written on real space. We obtain a set of equations where *z*-component of the **E** and **H**-fields can be removed. The six equations are reduced to four equations which

<sup>T</sup> <sup>ˆ</sup> , with

The band structure is calculated from unit cell of structure by setting frequency and calculating wave vectors propagating in the unit cell. Bloch's theorem applies to vector **F**:

> 

<sup>1</sup>

ˆ ,0 *<sup>z</sup> ik c*

*c e*

eigenvalues *kz* having an imaginary part are eliminated because they are not propagating

*r j cc c* 

If the above equations are joined together, we have an eigenvalue problem:

*r*

waves. The remaining eigenvalues gives us the band structure *kz*(*k0*).

with ˆ ˆˆ ,0 ,0 ,

*a e b e c e*

**Fr Fr Fr Fr Fr Fr**

 

*x y z*

 

and are the mesh size along *x*, *y* and *z*-axis. The transfer matrix of

1

and the wavelength *k0* to solve the eigenvalue problem. The

**T Fr Fr** (18)

*j j*

*N*

**Fr T Fr T Trr** (17)

(16)

*ik a ik b ik c*

*xy x y c EEHH*

*y*

*z*

1 /

1 /

**kH E** (14)

1 /

 

*k e ib*

*z*

*k e ic*

*ik c*

**F r T rr F r F r rr r r** (15)

**5. Finite-difference frequency-domain method (FDFD)** 

Maxwell on space (**k**, ) and applies approximations on the wave vector:

**kE H**

are written in a matrix form:

*a ab b c c* ' , ' '

We set the wave number *k*

**5.1 Calculation of band structure** 

*r*

 

the unit cell is written from the unit mesh:

 

The reflection and transmission calculation is performed using the transfer matrices. It is more interesting to calculate the elements of transmission matrix than the transmitted field in some points. The incident, reflected and transmitted waves are expanded on a planewave basis sets thanks to the transfer matrices.

ˆ **T0** is the transfer matrix of a vacuum layer. The eigenvectors of <sup>ˆ</sup> **T0** define a plane-wave basis set:

$$\begin{aligned} \hat{\mathbf{T}\_{\mathbf{0}}} \left| r\_{i} \right\rangle &= e^{ik\_{i}c} \left| r\_{i} \right\rangle\\ \left\langle l\_{i} \right| \hat{\mathbf{T}\_{\mathbf{0}}} &= e^{ik\_{i}c} \left\langle l\_{i} \right| \end{aligned} \tag{19}$$

The right and left eigenvectors are different because ˆ **T0** operator is not self-adjoint. They are expanded on the reciprocal space *kx* and *ky* and translate on the direct space with the Fourier transform because the transfer matrix is known in direct space. *kz* is found easily because we are in a vacuum layer. The transfer matrix ˆ **T** is not compatible with the new plane-wave basis set, we convert this matrix:

$$\left\langle l\_{j}\left|\hat{\mathbf{T}}\right|r\_{i}\right\rangle = \left\langle l\_{j}\left|\sum\_{\alpha}\tilde{T}\_{\alpha i}\right|r\_{\alpha}\right\rangle = \sum\_{\alpha} \tilde{T}\_{\alpha i}\left\langle l\_{j}\left|r\_{\alpha}\right\rangle = \tilde{T}\_{ji} \tag{20}$$

To calculate the transmission coefficient, this matrix is arranged to group the eigenvectors with the same propagation direction. This matrix is divided into four blocks:

$$
\tilde{\mathbf{T}} = \begin{pmatrix} \tilde{\mathbf{T}}\_{11} & \tilde{\mathbf{T}}\_{12} \\ \tilde{\mathbf{T}}\_{21} & \tilde{\mathbf{T}}\_{22} \end{pmatrix} \tag{21}
$$

TMM calculates the transmission coefficient of several layers from one layer (see 1D-MST). A great number of layers can be calculated without increasing the computing power.

Figure 11 plot the band structure and the transmission coefficient calculated by FDFD, for the previous dielectric PC. If we compares with figure 7, the first bands are correct but inaccuracies on the higher band after 0.6-frequency are due to the weak mesh (7x7x7 cells). If we increase the number of cells, accuracy increases. Calculation is more difficult for high frequencies because the E- and H-field are functions that oscillate more. The transmission calculation of 8 layers is calculated in 22 seconds with 3.3 Mo of working memory for 7x7x7 cells. The computing time is low because only one sphere is actually calculated thanks TMM.

### **6. Finite Integration Technique (FIT)**

In 1977, Weiland proposes a spatial discretization scheme to solve the integral equations of Maxwell (Weiland, 1977). This scheme called Finite Integration Technique (FIT) can be applied to many electromagnetism problems, in time and frequency domain, from static up to high frequency. The basic idea of this approach is to apply the Maxwell equations in integral form to a set of staggered grids like FDTD (figure 12).

Overview of Computational Methods for Photonic Crystals 279

solid-state physics, **H**-field is expanded on localized wave functions. They are calculated

*BZ*

The Wannier functions are a complete set of orthogonal functions which is defined for each band and each unit cell. The Wannier functions only depend on the quantity (**r**-**R**), **r** is the space position and **R** is any lattice vector. The reverse relation is written as

> 1 2 3 2 ( ,) , 2

*n n e* 

**H***n*(**k**,**r**) functions are the eigenvectors of **H**-wave equation. PWM solve this equation to get the **H***n*(**k**,**r**) functions for the PC without defect. To study a defect in PC, the **H**-field is

( ) , *n n*

Green's functions method solve the **H**-wave equation for the photonic crystal with defect.

Plane wave method is used to study the band structure of PC. It comes from the solid state physics where the electronic wave functions are scalar whereas the electromagnetic fields are vectors. A scalar approximation of fields is not enough to describe band structure. This method is modified to take the vectorial nature of the fields into account. Three computational methods of vectorial PWM were created quasi-simultaneously: the E-field method (Leung & Liu, 1990), the D-field method (Zhang & Satpathy, 1990) and the H-field method (Ho et al., 1990). The **D** and **H**-field are continuous in PC unlike the **E**-field. Moreover, only the differential operators of the wave equation on **E-**field and **H**-field are self-adjoint. We present the H-field method because of these properties**.** PWM expands the field and permittivity on plane-wave basis set. The permittivity is a periodic function with translational symmetry **T**, so Bloch's theorem can be applied to the

*n <sup>c</sup>* **R**

All periodic functions can be expanded on reciprocal space with Fourier series:

*r*

( )

**r**

( )

 **k R G**

*i*

*e d*

*i*

**H kr a Rr** (24)

**Hr R a Rr** (25)

with . ( ) ( ) ( ) ( ) *<sup>i</sup> e* **K r Hr hr hr hr T** (26)

**hr h** (27)

.

**G r G G**

*i*

*e*

 

*i*

*e*

.

**G r G G**

**k R a Rr H kr k** (23)

1 2 3 2 , ( ,) 2

from the Wannier functions **a***n*(**R**,**r**) for a PC without defect:

expanded on the previous Wannier functions:

**8. Plane Wave Method (PWM)** 

follows:

**H**-field:

The **H**-field is written:

*n n*

Fig. 11. Band structure and transmission coefficient of dielectric PC

Fig. 12. Tension and flux component on the mesh

The spatial discretization process is applied to the integral form of the Faraday's law:

$$\begin{aligned} \int\_{\mathcal{S}} \mathbf{E}(\mathbf{r}, t) \cdot d\mathbf{l} &= -\iint\_{\mathcal{S}} \frac{\partial}{\partial t} \mathbf{B}(\mathbf{r}, t) \cdot d\mathbf{S} \\\\ \hat{e}\_x(\mathbf{i}, j, k) + \hat{e}\_y(\mathbf{i} + \mathbf{1}, j, k) - \hat{e}\_x(\mathbf{i}, j + \mathbf{1}, k) - \hat{e}\_y(\mathbf{i}, j, k) &= -\frac{\partial}{\partial t} \hat{\bar{b}}\_z(\mathbf{i}, j, k) \end{aligned} \tag{22}$$

This last equation is an exact form without approximation. This process is applied to the other Maxwell equations in integral form. CST Microwave Studio is the software based on the FIT. Unlike FDTD, we use the local integral form and we can apply the technique of Perfect Boundary Approximation (PBA) which decreases the meshes on the boundaries.

#### **7. The tight binding method (TB)**

The tight binding method (Lidorikis et al., 1998) is less used than the PWM although this method is fast for the calculation of defect states in a PC. By analogy with the TB model in

Fig. 11. Band structure and transmission coefficient of dielectric PC

The spatial discretization process is applied to the integral form of the Faraday's law:

*t*

(22)

, , 1, , , 1, , , , ,

This last equation is an exact form without approximation. This process is applied to the other Maxwell equations in integral form. CST Microwave Studio is the software based on the FIT. Unlike FDTD, we use the local integral form and we can apply the technique of Perfect Boundary Approximation (PBA) which decreases the meshes on the boundaries.

The tight binding method (Lidorikis et al., 1998) is less used than the PWM although this method is fast for the calculation of defect states in a PC. By analogy with the TB model in

*e i j k ei j k ei j k ei j k bi j k*

*xy x y z*

Fig. 12. Tension and flux component on the mesh

*S S*

**7. The tight binding method (TB)** 

(,) (,)

**Er l Br S**

*td td t*

solid-state physics, **H**-field is expanded on localized wave functions. They are calculated from the Wannier functions **a***n*(**R**,**r**) for a PC without defect:

$$\mathbf{a}\_{\boldsymbol{n}}\left(\mathbf{R},\mathbf{r}\right) = \frac{\boldsymbol{\Omega}^{\boldsymbol{\mathbb{X}}}}{\left(2\pi\right)^{\boldsymbol{\mathbb{X}}}} \int\_{\boldsymbol{\mathbb{R}}\mathcal{Z}} \mathbf{H}\_{\boldsymbol{n}}(\mathbf{k},\mathbf{r}) e^{-i\mathbf{k}\cdot\mathbf{R}} d\mathbf{k} \tag{23}$$

The Wannier functions are a complete set of orthogonal functions which is defined for each band and each unit cell. The Wannier functions only depend on the quantity (**r**-**R**), **r** is the space position and **R** is any lattice vector. The reverse relation is written as follows:

$$\mathbf{H}\_n(\mathbf{k}, \mathbf{r}) = \frac{\Omega^{\frac{\mathbf{Z}}{\mathbf{r}}}}{\left(2\pi\right)^{\frac{\mathbf{Z}}{\mathbf{r}}}} \sum\_{\mathbf{G}} \mathbf{a}\_n\left(\mathbf{R}, \mathbf{r}\right) e^{i\mathbf{k} \cdot \mathbf{R}} \tag{24}$$

**H***n*(**k**,**r**) functions are the eigenvectors of **H**-wave equation. PWM solve this equation to get the **H***n*(**k**,**r**) functions for the PC without defect. To study a defect in PC, the **H**-field is expanded on the previous Wannier functions:

$$\mathbf{H(r)} = \sum\_{n} \sum\_{\mathbf{R}} c\_{n} \left( \mathbf{R} \right) \mathbf{a}\_{n} \left( \mathbf{R}, \mathbf{r} \right) \tag{25}$$

Green's functions method solve the **H**-wave equation for the photonic crystal with defect.

#### **8. Plane Wave Method (PWM)**

Plane wave method is used to study the band structure of PC. It comes from the solid state physics where the electronic wave functions are scalar whereas the electromagnetic fields are vectors. A scalar approximation of fields is not enough to describe band structure. This method is modified to take the vectorial nature of the fields into account. Three computational methods of vectorial PWM were created quasi-simultaneously: the E-field method (Leung & Liu, 1990), the D-field method (Zhang & Satpathy, 1990) and the H-field method (Ho et al., 1990). The **D** and **H**-field are continuous in PC unlike the **E**-field. Moreover, only the differential operators of the wave equation on **E-**field and **H**-field are self-adjoint. We present the H-field method because of these properties**.** PWM expands the field and permittivity on plane-wave basis set. The permittivity is a periodic function with translational symmetry **T**, so Bloch's theorem can be applied to the **H**-field:

$$\mathbf{H(r)} = e^{i\mathbf{K} \cdot \mathbf{r}} \mathbf{h(r)} \quad \text{with} \quad \mathbf{h(r)} = \mathbf{h(r+T)} \tag{26}$$

All periodic functions can be expanded on reciprocal space with Fourier series:

$$\begin{aligned} \varepsilon\_r(\mathbf{r}) &= \sum\_{\mathbf{G}} \varepsilon\_{\mathbf{G}} e^{i\mathbf{G} \cdot \mathbf{r}} \\ \mathbf{h}(\mathbf{r}) &= \sum\_{\mathbf{G}} \mathbf{h}\_{\mathbf{G}} e^{i\mathbf{G} \cdot \mathbf{r}} \end{aligned} \tag{27}$$

The **H**-field is written:

Overview of Computational Methods for Photonic Crystals 281

If we study one or two-dimensional PC, the central equation is written in another form and

22 21

**KG KG KG KG**

ˆˆ ˆˆ

**ee ee ee ee**

In the two-dimensional case, we choose the constant permittivity along *z*-axis. The vectors **K** and **G** are in the *xy* plane. The vectors of the central equation are indicated on figure 14.

**G G G**

**KG KG KG KG**

<sup>0</sup> 12 11

<sup>2</sup> ˆ *z*ˆ **e K G**

**K**+**G**

**K GK G** (31)

ˆˆ ˆˆ *<sup>k</sup>*

2

1 1

**G G**

*h h h h*

(32)

(34)

 

**K GK G** (33)

0

0

2 2

**8.1 One and two-dimensional photonic crystals** 

Fig. 14. Vectors definition of the central equation

z

In the central equation, the matrix of the scalar products is simplified:

22 21 12 11

**ee ee ee ee**

**9. One-dimensional multiple-scattering theory (1D-MST)** 

**G**

 

**G**

**G**

dimensional PC *i.e.* multilayer structures.

solve *N* equations.

 

**KG KG KG KG KG KG KG KG**

ˆˆ ˆˆ 1 0 ˆˆ ˆˆ 0 cos

1 ˆ**K G** <sup>y</sup>**e**

The 2*N* equations split into two parts, the TE-polarization and the TM-polarization. On the TM-polarization, the **H**-field vanishes on the *xy*-plane and the central equation is written:

1 1 21

1 1 1 2 22

**G G KG KG G G**

In the one-dimensional case, the equations for TE and TM-polarization become similar. We

The multiple-scattering studies the interaction between objects using the analytical solution for each object taken individually. In the next section, we will apply the MST to the cylinders and the spheres. Before, we will establish the analytical solution of one-

*h kh*

0

ˆ ˆ *h kh*

 *h kh* 

**G G G G**

On the TE-polarization, the **H**-field is parallel to *x*-axis and the equation is written:

**KG KG**

1 2 22

**K GK G e e**

**G G G G**

will be simplified:

1

x

**G G**

$$\mathbf{H(r)} = \sum\_{\mathbf{G}} \mathbf{h}\_{\mathbf{G}} e^{i(\mathbf{K} + \mathbf{G}) \cdot \mathbf{r}} = \sum\_{\mathbf{G}, \lambda = 1, 2} h\_{\mathbf{G}}^{\lambda} e^{i(\mathbf{K} + \mathbf{G}) \cdot \mathbf{r}} \hat{\mathbf{e}}\_{\perp \mathbf{K} + \mathbf{G}}^{\lambda} \tag{28}$$

As the **H**-field is transverse, each plane wave is perpendicular to the propagation vector **K+G**. The transverse plane of the propagation vector is described by the unit vectors <sup>1</sup> <sup>ˆ</sup> **K G <sup>e</sup>** and <sup>2</sup> <sup>ˆ</sup> **K G <sup>e</sup>** . The set of vectors 1 2 ˆ ˆ , , **KG KG e e KG** represents an orthonormal basis. We only need to storage two vectors instead of three, consequently data storage is reduced. The Fourier series expansion is replaced in the **H**-wave equation:

$$\begin{split} \nabla \times \left[ \boldsymbol{\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) \\ \nabla \times \left[ \sum\_{\mathbf{G}'} \boldsymbol{\varepsilon}\_{\mathbf{G}}^{-1} \boldsymbol{\varepsilon}^{\prime \mathbf{G}' \mathbf{r}} \nabla \times \sum\_{\mathbf{G}, \lambda} h\_{\mathbf{G}}^{\lambda} e^{i \left( \mathbf{K} \ast \mathbf{G} \right) \mathbf{r}} \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^{\lambda} \right] &= k\_{0}^{2} \sum\_{\mathbf{G}, \lambda} h\_{\mathbf{G}}^{\lambda} e^{i \left( \mathbf{K} \ast \mathbf{G} \right) \mathbf{r}} \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^{\lambda} \end{split} \tag{29}$$

Some algebraic calculations simplify this equation and we get for every vector **G**, the central equation of the photonic crystals:

$$\sum\_{\mathbf{G}',\mathbf{k}'} \mathbf{s}\_{\mathbf{G}\text{-}\mathbf{G}'}^{-1} \left[ (\mathbf{K} + \mathbf{G}) \times \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^{\lambda} \right] \cdot \left[ (\mathbf{K} + \mathbf{G}') \times \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^{\lambda'} \right] \mathbf{h}\_{\mathbf{G}'}^{\lambda'} = \mathbf{k}\_0^2 \mathbf{h}\_{\mathbf{G}}^{\lambda} \tag{30}$$

This equation is an eigenvalue problem which is solved with classical methods. The Bloch vector **K** is set and we try to find the eigenvalues *k0*. The calculation convergence depends on the *N* number of reciprocal lattice vectors **G**. A minimum number of vectors **G** is necessary to describe correctly the permittivity of the PC. The convergence of the problem is rather slow. As the **H**-field is transverse, the number of equations is decreased from 3*N* to 2*N*. The PWM is difficult to apply to the materials whose permittivity depends on the frequency like metals. On figure 13, we plot the band structure of the dielectric PC. It is calculated in 12 seconds with 10 Mo of working memory for 343 plane-waves. Calculation is fast because the structure is simple. If we compares with figure 7, we obtain the same result.

Fig. 13. Band structure of the previous dielectric PC

( ) ˆ *i i e h e*

**G G**

and <sup>2</sup> <sup>ˆ</sup> **K G <sup>e</sup>** . The set of vectors 1 2

*r*

,

Fig. 13. Band structure of the previous dielectric PC

0 0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

k0a / 2

equation of the photonic crystals:

Fourier series expansion is replaced in the **H**-wave equation:

 

**r Hr Hr**

 

0

**G G G**

*k*

1 2

 . . , 1,2

**G G K G**

 **KGr KGr**

As the **H**-field is transverse, each plane wave is perpendicular to the propagation vector **K+G**. The transverse plane of the propagation vector is described by the unit vectors <sup>1</sup> <sup>ˆ</sup> **K G <sup>e</sup>**

need to storage two vectors instead of three, consequently data storage is reduced. The

1 . . . 2

**G G K G G K G**

*i i i*

 

*e he k he* 

**G r KGr KGr**

, ,

Some algebraic calculations simplify this equation and we get for every vector **G**, the central

1 2

 

**G G K G KG G G <sup>G</sup>**

This equation is an eigenvalue problem which is solved with classical methods. The Bloch vector **K** is set and we try to find the eigenvalues *k0*. The calculation convergence depends on the *N* number of reciprocal lattice vectors **G**. A minimum number of vectors **G** is necessary to describe correctly the permittivity of the PC. The convergence of the problem is rather slow. As the **H**-field is transverse, the number of equations is decreased from 3*N* to 2*N*. The PWM is difficult to apply to the materials whose permittivity depends on the frequency like metals. On figure 13, we plot the band structure of the dielectric PC. It is calculated in 12 seconds with 10 Mo of working memory for 343 plane-waves. Calculation is fast because the structure is simple. If we compares with figure 7, we obtain the same result.

X M R

ka / 2

ˆ ˆ *h kh*

0

ˆ ˆ

**e e**

 

 

**KG e KG e** (30)

0

(29)

**Hr h e** (28)

 

ˆ ˆ , , **KG KG e e KG** represents an orthonormal basis. We only

#### **8.1 One and two-dimensional photonic crystals**

If we study one or two-dimensional PC, the central equation is written in another form and will be simplified:

$$\sum\_{\mathbf{G}'} \mathbf{c}\_{\mathbf{G} - \mathbf{G}'}^{-1} \left| \mathbf{K} + \mathbf{G} \right| \left| \mathbf{K} + \mathbf{G}' \right| \begin{pmatrix} \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^2 \cdot \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^2 & -\hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^2 \cdot \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^1 \\ -\hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^1 \cdot \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^2 & \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^1 \cdot \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^1 \end{pmatrix} \cdot \begin{pmatrix} h\_{\mathbf{G}'}^1 \\ h\_{\mathbf{G}'}^2 \end{pmatrix} = k\_0^2 \begin{pmatrix} h\_{\mathbf{G}}^1 \\ h\_{\mathbf{G}}^2 \end{pmatrix} \tag{31}$$

In the two-dimensional case, we choose the constant permittivity along *z*-axis. The vectors **K** and **G** are in the *xy* plane. The vectors of the central equation are indicated on figure 14.

Fig. 14. Vectors definition of the central equation

In the central equation, the matrix of the scalar products is simplified:

$$
\begin{pmatrix}
\hat{\mathbf{e}}\_{\perp\mathbf{K}+\mathbf{G}}^{2} \cdot \hat{\mathbf{e}}\_{\perp\mathbf{K}+\mathbf{G}'}^{2} & -\hat{\mathbf{e}}\_{\perp\mathbf{K}+\mathbf{G}}^{2} \cdot \hat{\mathbf{e}}\_{\perp\mathbf{K}+\mathbf{G}'}^{1} \\
\end{pmatrix} = \begin{pmatrix}
1 & 0 \\
0 & \cos(\theta - \theta')
\end{pmatrix} \tag{32}
$$

The 2*N* equations split into two parts, the TE-polarization and the TM-polarization. On the TM-polarization, the **H**-field vanishes on the *xy*-plane and the central equation is written:

$$\sum\_{\mathbf{G}'} \varepsilon\_{\mathbf{G}-\mathbf{G}'}^{-1} \left| \mathbf{K} + \mathbf{G} \right| \left| \mathbf{K} + \mathbf{G}' \right| h\_{\mathbf{G}'}^{1} = k\_0^2 h\_{\mathbf{G}}^{1} \tag{33}$$

On the TE-polarization, the **H**-field is parallel to *x*-axis and the equation is written:

$$\begin{aligned} \sum\_{\mathbf{G}'} \boldsymbol{\varepsilon}\_{\mathbf{G} - \mathbf{G}'}^{-1} \left| \mathbf{K} + \mathbf{G} \right| \left| \mathbf{K} + \mathbf{G}' \right| \left( \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}}^{1} \cdot \hat{\mathbf{e}}\_{\perp \mathbf{K} \ast \mathbf{G}'}^{1} \right) h\_{\mathbf{G}'}^{2} = k\_0^2 h\_{\mathbf{G}}^2 \\ \sum\_{\mathbf{G}'} \boldsymbol{\varepsilon}\_{\mathbf{G} - \mathbf{G}'}^{-1} \left( \mathbf{K} + \mathbf{G} \right) \cdot \left( \mathbf{K} + \mathbf{G}' \right) h\_{\mathbf{G}'}^{2} = k\_0^2 h\_{\mathbf{G}}^2 \end{aligned} \tag{34}$$

In the one-dimensional case, the equations for TE and TM-polarization become similar. We solve *N* equations.

#### **9. One-dimensional multiple-scattering theory (1D-MST)**

The multiple-scattering studies the interaction between objects using the analytical solution for each object taken individually. In the next section, we will apply the MST to the cylinders and the spheres. Before, we will establish the analytical solution of onedimensional PC *i.e.* multilayer structures.

Overview of Computational Methods for Photonic Crystals 283

2 1 11 12 2

The transmission of the layer is calculated from inversion of transfer matrix. The inverse will be calculated directly from phase shift matrices and interface matrix to avoid

<sup>1</sup> 1 2 <sup>1</sup> <sup>2</sup> <sup>2</sup>

1 2 2 1 *aa a bb b* 

1 11 12 1 21 11 22 12 21 1 <sup>1</sup> <sup>1</sup> 2 21 22 2 <sup>11</sup> <sup>2</sup> <sup>12</sup> 1

*b S S a T TT TT a a SS b T T b*

1 21 11 21 1 1 11 11 1 *<sup>T</sup>*

*S S*

1

*T T*

The Kronig-Penney model evaluates the electronic levels of a crystal structure in a onedimensional periodic potential (Kittel, 2005). This model has been modified to be used in PC and takes into account the oblique incidences (Mishra & Satpathy, 2003). We use TMM and the field expansion on backward and forward wave. The structure is periodic according to *z*-

a2

b2

The following equation converts the transfer matrix to scattering matrix:

Reflection and transmission coefficients are equal respectively to:

Several layers are calculated similarly.

Fig. 16. Refraction index of one dimensional PC

0 a

**E**<sup>1</sup> **E**<sup>2</sup>

a

n(z)

a1

b1

**9.2 Kronig-Penney model** 

axis (figure 16).

1 1

1 1

**T Φ Λ** (41)

(43)

(42)

z

1 11 11

**T T Λ Φ** (40)

The product of previous matrices provides the transfer matrix from layer 1 to layer 2:

2 1 21 22 with . *a a TT b b TT*

inaccuracies:

#### **9.1 Transfer-matrix method (TMM)**

Transfer-matrix method is known since many years (Born & Wolf, 1999). It is essential for the study of PC. TMM reduces the computational domain. In this section, transfer matrix will be applied to one-dimensional PC. For oblique incidences, we solve the **E**-wave equation for TE-polarization and the **H**-wave equation for TM-polarization. Two polarizations are separated and we use the same steps of calculation for two polarizations. We will study only the **E**-wave equation:

$$
\nabla^2 \mathbf{E}(\mathbf{r}) + k\_0^2 \varepsilon\_r(\mathbf{r}) \mathbf{E}(\mathbf{r}) = 0 \tag{35}
$$

Let us suppose that the layers of 1D-PC are stacked up along ez. The PC is uniform according to e┴ and e║ (figure 15). Because of these symmetries, the **E**-field and **E**-wave equation are simplified:

$$\begin{aligned} \mathbf{E}(\mathbf{r}) &= E(r\_{\parallel})E(r\_{\perp})E(z)\mathbf{e}\_{\perp} = e^{i\mathbf{k}\_{\parallel} \cdot \mathbf{r}\_{\parallel}}E(z)\mathbf{e}\_{\perp} \\ \frac{d^2E(z)}{dz^2} + k\_z^2(z)E(z) &= 0 \quad \text{with } k\_z^2(z) = k\_0^2 \varepsilon\_r(z) - k\_{\parallel}^2 \end{aligned} \tag{36}$$

We use the boundary conditions of each layer to solve the previous equation. The electric field *En* for a layer *n* is written with a forward wave and a backward wave (figure 15):

 *z n*, , *z n ik z ik z E z ae be nn n* (37)

Fig. 15. Field expansion on backward wave and forward wave

We calculate the transfer matrix between layer 1 and 2. The phase shift between the beginning and the end of layer 1 can be written in a transfer matrix form:

$$
\begin{pmatrix} a\_{1,end} \\ b\_{1,end} \end{pmatrix} = \boldsymbol{\Phi}\_n \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} \text{ with } \boldsymbol{\Phi}\_n = \begin{pmatrix} e^{ik\_{z,n}a} & 0 \\ 0 & e^{-ik\_{z,n}a} \end{pmatrix} \tag{38}
$$

Boundary conditions at the interface of layer 1 and layer 2 provide the following expression:

$$
\begin{pmatrix} a\_2 \\ b\_2 \end{pmatrix} = \boldsymbol{\Lambda}\_1^2 \begin{pmatrix} a\_{1,end} \\ b\_{1,end} \end{pmatrix} \text{ with } \boldsymbol{\Lambda}\_n^{n+1} = \frac{1}{2} \begin{pmatrix} 1 + k\_{z,n} \begin{/k\_{z,n+1}} & 1 - k\_{z,n} \begin{/k\_{z,n+1}} \\ k\_{z,n+1} \end{pmatrix} \\ 1 - k\_{z,n} \begin{cases} k\_{z,n+1} & 1 + k\_{z,n} \begin{/k\_{z,n+1}} \end{pmatrix} \end{pmatrix} \tag{39}
$$

Transfer-matrix method is known since many years (Born & Wolf, 1999). It is essential for the study of PC. TMM reduces the computational domain. In this section, transfer matrix will be applied to one-dimensional PC. For oblique incidences, we solve the **E**-wave equation for TE-polarization and the **H**-wave equation for TM-polarization. Two polarizations are separated and we use the same steps of calculation for two polarizations.

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

() ( )( )() ()

**k r E r e e**

*Er Er Ez e Ez*


2

*dz*

Fig. 15. Field expansion on backward wave and forward wave

**E**<sup>1</sup>

ez

1, 1 1, 1

*a a b b* 

*end*

*end*

*end end*

2 2 1, 1 2 1,

*a a <sup>b</sup> <sup>b</sup>*

beginning and the end of layer 1 can be written in a transfer matrix form:

*n*

**Φ** with

*n*

Let us suppose that the layers of 1D-PC are stacked up along ez. The PC is uniform according to e┴ and e║ (figure 15). Because of these symmetries, the **E**-field and **E**-wave

( ) ( ) ( ) 0 ( ) ( )

*dEz k zEz k z k z k*

*z z r*

We use the boundary conditions of each layer to solve the previous equation. The electric field *En* for a layer *n* is written with a forward wave and a backward wave (figure 15):

*z n*, , *z n ik z ik z E z ae be nn n*

**E**<sup>2</sup>

a1

b1

We calculate the transfer matrix between layer 1 and 2. The phase shift between the

Boundary conditions at the interface of layer 1 and layer 2 provide the following expression:

**Λ** with <sup>1</sup> , ,1 , ,1

,

*ik a n ik a e*

a2

b2

0 *z n*

1 1 1 2 1 1 *n zn zn zn zn*

, 0

*e* 

, ,1 , ,1

**Λ** (39)

*zn zn zn zn kk kk kk kk* 

 

*z n*

**Φ** (38)


(37)

2 22 2 <sup>2</sup> with 0 ||

0 1 2 n

*i*

(35)

(36)

an

bn

**9.1 Transfer-matrix method (TMM)** 

We will study only the **E**-wave equation:

equation are simplified:

e║

e┴

The product of previous matrices provides the transfer matrix from layer 1 to layer 2:

$$\mathbf{T}\begin{pmatrix} a\_2\\ b\_2 \end{pmatrix} = \mathbf{T}\begin{pmatrix} a\_1\\ b\_1 \end{pmatrix} \text{ with } \mathbf{T} = \begin{pmatrix} T\_{11} & T\_{12} \\ T\_{21} & T\_{22} \end{pmatrix} = \mathbf{A}\_1^2, \mathbf{O}\_1 \tag{40}$$

The transmission of the layer is calculated from inversion of transfer matrix. The inverse will be calculated directly from phase shift matrices and interface matrix to avoid inaccuracies:

$$
\begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} = \mathbf{T}^{-1} \begin{pmatrix} a\_2 \\ b\_2 \end{pmatrix} = \left[ \boldsymbol{\Phi}\_1 \right]^{-1} \left[ \boldsymbol{\Lambda}\_1^2 \right]^{-1} \begin{pmatrix} a\_2 \\ b\_2 \end{pmatrix} \tag{41}
$$

The following equation converts the transfer matrix to scattering matrix:

$$
\begin{pmatrix} b\_1 \\ a\_2 \end{pmatrix} = \begin{pmatrix} S\_{11} & S\_{12} \\ S\_{21} & S\_{22} \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_2 \end{pmatrix} = \frac{1}{T\_{11}^{-1}} \begin{pmatrix} T\_{21}^{-1} & T\_{11}^{-1}T\_{22}^{-1} - T\_{12}^{-1}T\_{21}^{-1} \\ 1 & -T\_{12}^{-1} \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_2 \end{pmatrix} \tag{42}
$$

Reflection and transmission coefficients are equal respectively to:

$$S\_{11} = \frac{T\_{21}^{-1}}{T\_{11}^{-1}} \qquad S\_{21} = \frac{1}{T\_{11}^{-1}} \tag{43}$$

Several layers are calculated similarly.

#### **9.2 Kronig-Penney model**

The Kronig-Penney model evaluates the electronic levels of a crystal structure in a onedimensional periodic potential (Kittel, 2005). This model has been modified to be used in PC and takes into account the oblique incidences (Mishra & Satpathy, 2003). We use TMM and the field expansion on backward and forward wave. The structure is periodic according to *z*axis (figure 16).

Fig. 16. Refraction index of one dimensional PC

Overview of Computational Methods for Photonic Crystals 285

The multiple-scattering is an analytical theory which calculates the scattering of *N* objects from the scattering of each object independently (Felbacq et al., 1994). In the two-

<sup>222</sup> 1 *E kE k E <sup>r</sup>*

Outside the cylinders, the wave equation can be split into two equations, one for the

*s*

*ss r*

Using Green's theorem, the scattered field outside the cylinder can be written as follows:

The function *H* is a Hankel function. The surface integral can be restricted to the *Cj*

*E H k E dS*

**r rr r**

The cylinders interact to give the scattered field by the structure. To better understand the process of multiple-scattering between the cylinders, we will study a simple example of cylinders with circular section aligned along *x*-axis and excited by an incident field

0 0 *ikx l*

We define the total incident field around the cylinder *i* which takes into account other

*l l*

*l l l*

*E e i J kr a J kr* **<sup>r</sup>** (55)

C2 Cn

…

4 *<sup>j</sup>*

C1

4 *s r ik <sup>E</sup> H k*

<sup>2</sup> 1

**r**

*j*

*<sup>s</sup> <sup>C</sup>*

*E kE k E*

**r r rr**

1

1

*E dS* **<sup>r</sup> rr r r** (53)

0

**r r rr** (51)

(52)

(54)

dimensional case, objects are cylinders and the theory uses the scalar wave equation:

 

**r rr r r**

0 2 2 0 0 222

 

*EE E E kE*

2

**r r**

*s s j*

*E E*

*r j*

*E*

*ik*

*j*

**10. Two-dimensional multiple-scattering theory (2D-MST)** 

incident field *E*0 and the other one for the scattered field *Es*:

cylinders. The scattered field is written as a sum:

propagating along *x*-axis (figure 17).

y

cylinders:

Fig. 17. Incident field on *n* cylinders aligned

x

The incident field is expanded into Bessel functions:

The field *E*(*z*) of the previous section is a Bloch function on *z*. The continuity of tangential fields and the Bloch theorem give the two following relations:

$$\begin{cases} E\_2(a) = e^{i\mathcal{K}a} E\_1(0) \\ \left. \frac{dE\_2}{dz} \right|\_a = e^{i\mathcal{K}a} \left. \frac{dE\_1}{dz} \right|\_0 \end{cases} \tag{44}$$

*K* is the wave vector of Bloch. In the unit cell, the permittivity is not uniform according to *z*axis. It is splitted into sub-cells with a constant permittivity. We calculate the transfer matrix of the unit cell from the expressions of the previous section:

$$
\begin{pmatrix} a\_2 \\ b\_2 \end{pmatrix} = \begin{pmatrix} T\_{11} & T\_{12} \\ T\_{21} & T\_{22} \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} \tag{45}
$$

We expand equation 44 on forward and backward waves:

$$\begin{cases} a\_2 + b\_2 = e^{iKa} \left( a\_1 + b\_1 \right) \\ a\_2 - b\_2 = e^{iKa} \left( a\_1 - b\_1 \right) \end{cases} \tag{46}$$

Equation 46 is replaced in equation 45:

$$
\begin{pmatrix} T\_{11} & T\_{12} \\ T\_{21} & T\_{22} \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} \tag{47}
$$

The Bloch factor *eiKa* and the complex conjugate value *e-iKa* are the two eigenvalues of the transfer matrix because the determinant of this matrix is equal to one. The trace of the transfer matrix is equal to the sum of the eigenvalues:

$$T\_{11} + T\_{22} = e^{i\mathbb{K}a} + e^{-i\mathbb{K}a} \tag{48}$$

If the materials of the photonic crystal do not absorb, we have \* *T T* 11 22 . Similarly to the Kronig-Penney model, the above relation is written from the transmission coefficient:

$$\begin{aligned} \left| S\_{21} = \left| S\_{21} \right| e^{i\phi\_{S\_{21}}} = \left( T\_{11}^\* \right)^{-1} = \left( T\_{22} \right)^{-1} \\ \frac{e^{i\phi\_{S\_{21}}}}{\left| S\_{21} \right|} + \frac{e^{-i\phi\_{S\_{21}}}}{\left| S\_{21} \right|} = e^{iKa} + e^{-iKa} \end{aligned} \tag{49}$$

After simplification, the transcendent equation is written:

$$\frac{\cos\left(\varphi\_{\mathbb{S}\_{21}}\right)}{\left|\mathbb{S}\_{21}\right|} = \cos\left(\mathbb{K}a\right) \tag{50}$$

We get an equation similar to electronic case. The transmission coefficient is different in the TE and TM-polarization. To plot band structure *f( k║*,*K*, *k0*) = 0, we set the wave numbers *k0* and *k║* and calculate the Bloch number *K*.

The field *E*(*z*) of the previous section is a Bloch function on *z*. The continuity of tangential

2 1 2 1

*dE dE <sup>e</sup> dz dz*

*K* is the wave vector of Bloch. In the unit cell, the permittivity is not uniform according to *z*axis. It is splitted into sub-cells with a constant permittivity. We calculate the transfer matrix

> 2 11 12 1 2 21 22 1 *a TT a b TT b*

> > 22 11 22 11

11 12 1 1 21 22 1 1 <sup>=</sup> *iKa TT a a <sup>e</sup> TT b b* 

The Bloch factor *eiKa* and the complex conjugate value *e-iKa* are the two eigenvalues of the transfer matrix because the determinant of this matrix is equal to one. The trace of the

If the materials of the photonic crystal do not absorb, we have \* *T T* 11 22 . Similarly to the

21 21 11 22

*S Se T T*

<sup>21</sup>

*iKa iKa*

<sup>21</sup>

cos cos *<sup>S</sup> Ka*

We get an equation similar to electronic case. The transmission coefficient is different in the TE and TM-polarization. To plot band structure *f( k║*,*K*, *k0*) = 0, we set the wave numbers *k0*

1 1 \*

11 22

21 21

*S S*

*i i*

After simplification, the transcendent equation is written:

and *k║* and calculate the Bloch number *K*.

21 21

Kronig-Penney model, the above relation is written from the transmission coefficient:

*S*

*i*

*e e e e S S*

21

*S* 

 

*iKa iKa ab e ab ab e ab* 

*Ea e E*

 

*a*

( ) (0) *iKa iKa*

0

 

(44)

(46)

(49)

(50)

(45)

(47)

*iKa iKa TT e e* (48)

fields and the Bloch theorem give the two following relations:

of the unit cell from the expressions of the previous section:

We expand equation 44 on forward and backward waves:

transfer matrix is equal to the sum of the eigenvalues:

Equation 46 is replaced in equation 45:

 

#### **10. Two-dimensional multiple-scattering theory (2D-MST)**

The multiple-scattering is an analytical theory which calculates the scattering of *N* objects from the scattering of each object independently (Felbacq et al., 1994). In the twodimensional case, objects are cylinders and the theory uses the scalar wave equation:

$$
\nabla^2 E(\mathbf{r}) + k^2 E(\mathbf{r}) = k^2 \left[ 1 - \varepsilon\_r(\mathbf{r}) \right] E(\mathbf{r}) \tag{51}
$$

Outside the cylinders, the wave equation can be split into two equations, one for the incident field *E*0 and the other one for the scattered field *Es*:

$$\begin{aligned} E(\mathbf{r}) &= E\_0(\mathbf{r}) + E\_s(\mathbf{r})\\ \nabla^2 E\_0(\mathbf{r}) + k^2 E\_0(\mathbf{r}) &= 0\\ \nabla^2 E\_s(\mathbf{r}) + k^2 E\_s(\mathbf{r}) &= k^2 \left[ 1 - \varepsilon\_r(\mathbf{r}) \right] E(\mathbf{r}) \end{aligned} \tag{52}$$

Using Green's theorem, the scattered field outside the cylinder can be written as follows:

$$E\_s\left(\mathbf{r}\right) = -\frac{ik^2}{4} \iint H\left(k\left|\mathbf{r} - \mathbf{r}\prime\right|\right) \left[1 - \varepsilon\_r\left(\mathbf{r}\right)\right] E\left(\mathbf{r}\prime\right) dS\prime \tag{53}$$

The function *H* is a Hankel function. The surface integral can be restricted to the *Cj* cylinders. The scattered field is written as a sum:

$$\begin{aligned} E\_s(\mathbf{r}) &= \sum\_j E\_s^j(\mathbf{r})\\ E\_s^j(\mathbf{r}) &= \frac{ik^2 \left(\varepsilon\_r^j(\mathbf{r}) - 1\right)}{4} \iint\_{C\_j} H\left(k|\mathbf{r} - \mathbf{r}'|\right) E(\mathbf{r}')dS' \end{aligned} \tag{54}$$

The cylinders interact to give the scattered field by the structure. To better understand the process of multiple-scattering between the cylinders, we will study a simple example of cylinders with circular section aligned along *x*-axis and excited by an incident field propagating along *x*-axis (figure 17).

Fig. 17. Incident field on *n* cylinders aligned

The incident field is expanded into Bessel functions:

$$E\_0\left(\mathbf{r}\right) = e^{ikx} = \sum\_l i^l J\_l\left(kr\right) = \sum\_l a\_{0l} J\_l\left(kr\right) \tag{55}$$

We define the total incident field around the cylinder *i* which takes into account other cylinders:

Overview of Computational Methods for Photonic Crystals 287

The three-dimensional multiple-scattering theory use same principles as the twodimensional case but three-dimensional case is more complex because the Helmholtz equation is vectorial. The 3D-MST is separated into two parts, the calculation for one sphere and the generalization to *N* spheres (Oyhenart & Vignéras, 2007). The first part calculates the scattered wave by one sphere with Mie theory. The base of spherical harmonics is

1, ,

functions. The incident field and the scattered field are expanded on spherical wave

0 (1)

 

*s s s lm lm lm l lm*

 

*<sup>i</sup> z kr Y*

1 ˆ

*l m l lm*

 

<sup>4</sup> () ,

(3) 0 ,, ,, , , , ,

with

*E E SE*

The elements of **S**-matrix are Mie's coefficients (Bohren & Huffman, 1998). Figure 19

The second part of the method is an iterative algorithm to calculate scattered field for *N* spheres from one sphere (figure 20). For the first order, we calculate the scattering of the incident field for each sphere. For the second order, the scattered field of first order for one sphere becomes the incident field for the *N*-1 other spheres. With this new incident field, the scattered field is calculated as at first order and so on, for higher orders. This iterative process stops when it is converged. The total scattered field is the contribution of all spheres and all orders. The material and the size of these *N* spheres can be different. Moreover, spheres may be put in a random way, without symmetry conditions on spheres positions.

*E* 

*k* 

 

> 

  (60)

are scalar spherical wave

(61)

1, , 1

*l m M*

**Z Z**

0 ,, ,, , ,

**r Z**

( ) *lm lm l m*

*E*

*Zl* are the Hankel and Bessel spherical functions and , *lm <sup>Y</sup>*

, ,

**E**

*l m*

presents an example of the incident and scattered fields by a sphere.

Fig. 19. Pictures of incident and scattered fields by a metallic sphere

**E Z**

*l*

**Z r**

*i ll i k*

**11. Three-dimensional multiple-scattering theory (3D-MST)** 

written as follows:

functions:

$$E\_i^{\,}\left(\mathbf{r}\right) = \sum\_l a\_{il} I\_l\left(kr\right) \tag{56}$$

In the case of the cylinders with circular section, the integral 54 is calculated analytically:

$$E\_s^j\left(\mathbf{r}\right) = \sum\_l b\_{jl} H\_l^{(1)}\left(kr\right) \text{ with } b\_{jl} = S\_{jl} a\_{jl} \tag{57}$$

*Sjl* is the scattering coefficients (Bohren & Huffman, 1998). The field *<sup>j</sup> Es* **<sup>r</sup>** scattered by a cylinder *j* is shifting on another cylinder *i* to become an incident field *Gjil.bjl* on this cylinder with *Gjil* the translation coefficients (Felbacq et al., 1994). If we shift all the scattered fields on the cylinder *i*, and if we add the initial incident field to it, we get the total incident field on the cylinder *i*:

$$a\_{il} = \sum\_{j \neq i} G\_{jil} b\_{jl} + a\_{0l} \tag{58}$$

By using the equation 57, we get the multiple-scattering equation:

$$a\_{il} = \sum\_{j \neq i} G\_{jil} S\_{jl} a\_{jl} + a\_{0l} \tag{59}$$

*Gjil* and *Sjl* are the translation and scattering matrix coefficients. To get the total incident field and the field scattered by the structure, it is necessary to calculate *Gjil* and *Sjl* coefficients and to solve the multiple-scattering equation. This method is suitable for the study of defects in PC because it does not impose any condition on the position and the material of cylinders. On figure 18, we study a 2D-PC with 80 cylinders doped by a microcavity.

Fig. 18. Transmission coefficient of triangular 2D-PC with a defect

 *<sup>i</sup> i il l l*

cylinder *j* is shifting on another cylinder *i* to become an incident field *Gjil.bjl* on this cylinder with *Gjil* the translation coefficients (Felbacq et al., 1994). If we shift all the scattered fields on the cylinder *i*, and if we add the initial incident field to it, we get the total incident field

> *il jil jl l* 0 *j i a Gb a*

*il jil jl jl l* 0

*Gjil* and *Sjl* are the translation and scattering matrix coefficients. To get the total incident field and the field scattered by the structure, it is necessary to calculate *Gjil* and *Sjl* coefficients and to solve the multiple-scattering equation. This method is suitable for the study of defects in PC because it does not impose any condition on the position and the material of cylinders.

*j i a G Sa a* 

On figure 18, we study a 2D-PC with 80 cylinders doped by a microcavity.

Fig. 18. Transmission coefficient of triangular 2D-PC with a defect

In the case of the cylinders with circular section, the integral 54 is calculated analytically:

 *j* (1) *s jl l l*

*Sjl* is the scattering coefficients (Bohren & Huffman, 1998). The field

By using the equation 57, we get the multiple-scattering equation:

on the cylinder *i*:

*E a J kr* **<sup>r</sup>** (56)

(58)

(59)

*<sup>j</sup> Es* **<sup>r</sup>** scattered by a

*E b H kr* **<sup>r</sup>** with *jl jl jl b Sa* (57)

### **11. Three-dimensional multiple-scattering theory (3D-MST)**

The three-dimensional multiple-scattering theory use same principles as the twodimensional case but three-dimensional case is more complex because the Helmholtz equation is vectorial. The 3D-MST is separated into two parts, the calculation for one sphere and the generalization to *N* spheres (Oyhenart & Vignéras, 2007). The first part calculates the scattered wave by one sphere with Mie theory. The base of spherical harmonics is written as follows:

$$\begin{aligned} \left| \mathbf{Z}\_{-1,l,m} \right\rangle &= \frac{i^l \sqrt{4\pi}}{i\sqrt{l(l+1)}} z\_l(kr) \mathbf{r} \times \nabla Y\_{lm} \left( \theta, \phi \right) \\\\ \left| \mathbf{Z}\_{1,l,m} \right\rangle &= -\frac{i}{k} \nabla \times \hat{\mathbf{Z}}\_{M=-1} \end{aligned} \tag{60}$$

*Zl* are the Hankel and Bessel spherical functions and , *lm <sup>Y</sup>* are scalar spherical wave functions. The incident field and the scattered field are expanded on spherical wave functions:

$$\begin{aligned} \mathbf{E}\_0(\mathbf{r}) &= \sum\_{\sigma,l,m} \mathbf{E}^0\_{\sigma,l,m} \Big| \mathbf{Z}^{(1)}\_{\sigma,l,m} \\ \mathbf{E}\_s &= \sum\_{\sigma,l,m} \mathbf{E}^s\_{\sigma,l,m} \Big| \mathbf{Z}^{(3)}\_{\sigma,l,m} \Big> \text{with } \mathbf{E}^s\_{\sigma,l,m} = \mathbf{S}\_{\sigma l} \mathbf{E}^0\_{\sigma,l,m} \end{aligned} \tag{61}$$

The elements of **S**-matrix are Mie's coefficients (Bohren & Huffman, 1998). Figure 19 presents an example of the incident and scattered fields by a sphere.

Fig. 19. Pictures of incident and scattered fields by a metallic sphere

The second part of the method is an iterative algorithm to calculate scattered field for *N* spheres from one sphere (figure 20). For the first order, we calculate the scattering of the incident field for each sphere. For the second order, the scattered field of first order for one sphere becomes the incident field for the *N*-1 other spheres. With this new incident field, the scattered field is calculated as at first order and so on, for higher orders. This iterative process stops when it is converged. The total scattered field is the contribution of all spheres and all orders. The material and the size of these *N* spheres can be different. Moreover, spheres may be put in a random way, without symmetry conditions on spheres positions.

Overview of Computational Methods for Photonic Crystals 289

The computational methods are studied and compared on table 1. The first four methods are three-dimensional numerical methods coming from electromagnetism. The FEM gives very precise results but it requires many resources systems. Another method, the FDTD is very used and converges without an excessive mesh thanks to its formulation. If you want more precise results, the mesh becomes heavy and you need to use FEM. The FDFD studies photonic crystals with a high number of layers by modeling only one layer. The FIT is a method which studies numerically PC with a high number of objects, without holding excessive resources system but the result is approximate. These methods calculate dielectric and metallic PC to obtain reflection and transmission coefficients and the band structure.

Then, we have others methods resulting from the solid state physics which require a very small computing time. They were adapted from the scalar methods of the solid state physics. The PWM calculate only the band structure. The tight binding method is less used than the PWM although this method is fast for the calculation of defect states in a PC. The multiple-scattering theory which is also used in optics studies analytically large finite PC and it easily takes the defects into account in the PC. For three-dimensional PC, there is no simple method, fast and accurate. We must remove a feature. For example, the 3D-

Origin E.M. E.M. E.M. E.M. Q.M. Q.M. Q.M. E.M. both Maxwell's equations Freq. Time Freq. Time Freq. Freq. Freq. Freq. Freq. Calculation Num. Num. Num. Num. Num. Num. Analyt. Analyt. Analyt. Geometries 3D 3D 3D 3D 2D/3D 3D 1D Cyl. Sphere

Expansion, series X X X X X X Free space, finite PC X X X X X X Infinite periodic PC X X X X X X X X Band structure X X X X X X X X Transmission/reflection X X X X X X X Metallic PC X X X X X Single defect in finite-PC X X X X X X Periodic defect X X X X X X X X X computing speed slow med. fast fast fast med. fast fast fast data storage large low low low low low low low low

Free software X X X X X X

FEM FDTD FDFD FIT TB PWM 1D-MST 2D-MST 3D-MST

Analyt. : analytical method Num. : numerical method

Med. : medium Cyl. : cylinder

**12. Conclusion: Comparison of computational method** 

MST is fast and accurate but the computer program is complex to write.

commercial software X X X

Popular method in PC X X X

Q.M. : quantum mechanics and solid state physics

Table 1. Specifications of the computational methods

Abbreviations:

E.M. : electromagnetism Freq. : frequency domain Time : time domain

Discrete equations X X X X

These two last remarks show all the interest of this method in calculation of PC with defects and random structures.

Fig. 20. Block diagram of multiple-scattering method, the total scattered field is the sum of the scattered fields for all orders.

For periodic structures, calculation is simplified. Figure 21 plot the transmission coefficient for the previous infinite dielectric PC. The transmission calculation of 8 layers is calculated in 11 minutes with 63 Mo of working memory. By using the principle of KKR-method of solid state physics, MST can also calculate the band structure (Wang and Al, 1993). MST is the fastest method for the finite structures and random structures.

Fig. 21. Transmission coefficient of the infinite periodic dielectric PC

These two last remarks show all the interest of this method in calculation of PC with defects

Fig. 20. Block diagram of multiple-scattering method, the total scattered field is the sum of

For periodic structures, calculation is simplified. Figure 21 plot the transmission coefficient for the previous infinite dielectric PC. The transmission calculation of 8 layers is calculated in 11 minutes with 63 Mo of working memory. By using the principle of KKR-method of solid state physics, MST can also calculate the band structure (Wang and Al, 1993). MST

is the fastest method for the finite structures and random structures.

Fig. 21. Transmission coefficient of the infinite periodic dielectric PC

and random structures.

the scattered fields for all orders.
