## Abstract

A simple implementation of plane wave method is presented for modeling photonic crystals with arbitrary shaped ‘atoms’. The Fourier transform for a single ‘atom’ is first calculated either by analytical Fourier transform or numerical FFT, then the shift property is used to obtain the Fourier transform for any arbitrary supercell consisting of a finite number of ‘atoms’. To ensure accurate results, generally, two iterating processes including the plane wave iteration and grid resolution iteration must converge. Analysis shows that using analytical Fourier transform when available can improve accuracy and avoid the grid resolution iteration. It converges to the accurate results quickly using a small number of plane waves. Coordinate conversion is used to treat non-orthogonal unit cell with non-regular ‘atom’ and then is treated by standard numerical FFT. MATLAB source code for the implementation requires about less than 150 statements, and is freely available at http://www.lions.odu.edu/~sguox002.

© 2002 Optical Society of America

## 1. Introduction

The plane wave method (PWM) is often used for photonic crystal modeling since it can yield accurate and reliable results. This method requires intensive computations for complicated systems, involving thousands of plane waves, and places a high demand on computer resources and time [1]. There is a tradeoff between the computing time and accuracy. In this paper, we present a simple and fast implementation method using MATLAB, utilizing its abundant functions for numerical analysis and graphics. The whole program may have less than 150 statements for performing the calculations and graphical output. However, the computing time and accurately can be much improved. The key point is to obtain the Fourier transform of the unit cell as accurate as possible. For commonly used ‘atoms’ with regular shapes, such as square, rectangular, circular cylinders, cubes and spheres, we use the analytical Fourier transform; there is no need to do a numerical Fast Fourier Transform (FFT), avoiding the step of dividing the space into small grids. The computation time is effectively reduced and accuracy can be improved. For a finite sized supercell, the Fourier transform is obtained using the shift property.

## 2. Theory of PWM

The PWM is illustrated in several papers [2–5]. Here, we summarize the theory very briefly. Maxwell’s equations in a transparent, time-invariant, source free and non-permeable (µ=µ_{0}) space can be rewritten as Helmholz’s equation:

where ε(r) is the dielectric function, ω is the angular frequency and *c* is the speed of light in vacuum.

In an infinite periodic photonic crystal, using Bloch’s theorem, a mode in a periodic structure can be expanded as a sum of infinite number of plane waves:

where λ=1, 2, k⃗ is the wave vector of the plane wave, G⃗ is the reciprocal lattice vector, e_{̂} represents the two unit axis perpendicular to the propagation direction k⃗ + G⃗. (*e*̂_{1}, *e*̂_{2}, k⃗ + G⃗) are perpendicular to each other. *h*
_{G,λ} is the coefficient of the *H* component along the axes *e*̂_{λ}.

Using the Fourier transform, the dielectric function can be written as:

where Ω is the unit cell and V is the volume of the unit cell.

Finally, Helmholz’s equation can be transformed to an algebraic form [6]:

This is a standard eigenvalue problem and it can be solved using a standard eigen-solver. For 1D and 2D cases, the equation can be simplified.

## 3. Implementation

For simplicity, we assume there is only one kind of ‘atom’ in a photonic crystal and the number of ‘atoms’ in the unit cell or supercell is finite. Once the unit cell or supercell is determined, the set of reciprocal lattice vector G⃗ and unit vectors *e*̂_{1}, *e*̂_{2} can be calculated easily. The key part is to obtain the Fourier coefficient matrix ε(G⃗ - G⃗′) according to Eq.(3). It is obtained by calculating the Fourier coefficient of a single atom first using analytical or numerical Fourier transform, and then calculating the Fourier coefficients for the supercell using shift property, finally re-arranging to get the coefficients ε(G⃗ - G⃗′) and inversing to get ε^{-1}(G⃗ - G⃗′).

#### 3.1 ‘Atoms’ with regular shape

It is advantageous to use analytical Fourier transform when available. We illustrate this using the 2D photonic crystals with the most commonly used circular cylinder. The same procedure can be followed for other shapes and their Fourier transforms can be found in Ref. [6].

Assuming the radius of the cylinder is *R*, the dielectric constant for the cylinder is ε_{a}, the background dielectric constant is ε_{b}, the lattice structure can be represented by the two lattice basis vector *a*⃗_{1} and *a*⃗_{2}. The area of the unit cell is calculated as *A* = |*a*⃗_{1} × *a*⃗_{2}|, the Fourier transform of the unit cell is:

where *J*
_{1} is the 1^{st} order Bessel function, *G* is the modulus of G⃗, ƒ is a fraction parameter:

Here the advantage is that in different lattice geometries, the Fourier transform is the same as in Eq. (5) except that the values of ƒ and *G* are different.

#### 3.2 ‘Atoms’ with arbitrary shape

When atom shape is not regular, only numerical FFT can be used. To use FFT in a non-orthogonal lattice, a non-orthogonal unit cell is first converted to an orthogonal cell using coordinate conversion as illustrated below.

In Cartesian coordinate system, the three basis vectors in real space *a*⃗_{1}, *a*⃗_{2}, *a*⃗_{3} are:

$${\overrightarrow{a}}_{2}={a}_{2x}\hat{x}+{a}_{2y}\hat{y}+{a}_{2z}\hat{z},$$

$${\overrightarrow{a}}_{3}={a}_{3x}\hat{x}+{a}_{3y}\hat{y}+{a}_{3z}\hat{z}$$

The dielectric function in the unit cell is ε(*r*) and the column vector is r⃗ = *m**a*⃗_{1} + *n**a*⃗_{2} + *l**a*⃗_{3}, where *m, n* and *l* are coordinates along the basis vectors. In Cartesian coordinates:

So we have:

where

[*g*] is called the metric for the oblique coordinates [7, 8].

In Fig. 1 we show how these conversions are applied in the case of two examples: A circular cylinder in a triangular lattice is converted to an oblique elliptical cylinder in a square lattice as shown in Fig.1(a), while Fig.1(b) shows the conversion of an elliptical cylinder in a triangular lattice. The resulting band diagram for a photonic crystal of elliptical atoms illustrated in Fig.1(b) is shown in Fig.1(c) using appropriate material and structural parameters for GaAs. The conversion can be conveniently applied to treat photonic crystal fiber with non-regular shaped ‘atoms’ such as the one with elliptical air holes.

#### 3.3 Shift property of Fourier transform

Assuming the Fourier transform of a single atom ε(r) is known as ε_{G}, if ε(*r*) is shifted by an amount *r*
_{0}, then its Fourier transform must be multiplied by
${e}^{i\overrightarrow{k}\xb7{\overrightarrow{r}}_{0}}$
[6] which we call shift property here:

Therefore, for a supercell with several atoms in periodic or random positions, the Fourier transform can be obtained using addition and subtraction:

where r⃗
_{i}
is the location of ‘atom’ *i* in the supercell.

Equations (11–12) are especially suitable for supercell method, such as the photonic crystal with defects. We can obtain the accurate Fourier coefficients of the supercell at the required *G* grid points by doing simple additions and subtractions, requiring only the Fourier coefficients of each single kind of atom. This is especially advantageous for a large supercell with many periodic or random atoms in it.

As an example, we show below how the 3D diamond lattice is worked out:

The diamond lattice is a complex FCC lattice with two spherical atoms in the primitive cell. Assuming the length of the simple cubic side is *a*, the primitive lattice vector basis are defined as *a*⃗_{1} = [0,1,1]*a*/2,*a*⃗_{2} = [1,0,1]*a*/2,*a*⃗_{3} = [1,1,0]*a*/2. The locations of the two atoms in the primitive cell are chosen as r⃗0 = [-1,-1,-1]*a*/8 and r⃗1 = [1,1,1]*a*/8 to keep inversion symmetry. The primitive reciprocal lattice vectors b⃗_{1},b⃗_{2},b⃗_{3} are calculated according to (5) in Appendix B of Ref. [10].

Using shift property and Fourier transform for a sphere, the Fourier coefficient at the reciprocal lattice grid is expressed as:

where *R* is the radius of the sphere,
$f=\frac{2\times \frac{4}{3}\pi {R}^{3}}{V}$
*V* = |*a*⃗_{1} · *a*⃗_{2} × *a*⃗_{3}|.

In Fig. 2 we show the band structure for *R*=√3/8*a*, ε_{a}=13, ε_{b}=1 using our simple program with 343 plane waves, which is in excellent agreement with the result in Ref. [3]. In this paper, the set of G is chosen as G⃗=*n*
_{1}b⃗_{1} + *n*
_{2}b⃗ + *n*
_{3}b⃗_{3},|*n*_{i}
|≤*n* and in this example *n*=3.

## 4. Convergence, accuracy and stability

We performed the TM/TE band calculation for an ideal triangular lattice with air holes in GaAs. The radius of the air hole is 0.28*a*, where *a* is the lattice constant, and the dielectric constant for GaAs is 13.0. All the eigen-frequencies with *k*-point located at **M** (see the inset of Fig. 1(c)) are calculated.

Figures 3 and 4 show the convergence for TM and TE modes as a function of the number of plane waves with different mesh resolutions. Two methods are compared: the analytical Fourier transform and FFT with each grid point averaged by a 10×10 submesh. The iteration error is calculated as *norm*[*X*(*n*)- *X*(*n* - 1)], where *X*(n) is the eigen-frequency vector of the first 10 bands for the *n*th iteration, and the number of plane waves used is *N*
_{PW}=(2n+1)^{2}.

In Fig. 3 and Fig. 4, there exist two converging processes: number of plane waves and mesh resolution. These two processes are almost independent of each other. To achieve accurate results, both convergences must be reached. The iteration errors shown in these two figures are not enough to determine whether the accurate values are achieved, since they represent only one process. When a fine enough mesh and enough number of plane waves are used, the frequencies converge to the accurate values. The mesh grid number does not need to be the same as the number of plane waves as in Ref. [4], since that makes the computation too large. In our case, the number of plane waves is always much smaller than the grid number.

Analytical method has only one converging process; numerical FFT and mesh formation procedures are not required. Therefore, analytical method has higher accuracy for the same number of plane waves than FFT; it can be faster and save a lot of memory and computation time. It may be advantageous especially for 3D cases when the problem size is large and computation is long. The reason why analytical method converges a little bit slower than FFT with a certain mesh resolution (see the iteration error) is probably due to the Gibbs’ phenomenon in Fourier series.

Also, Figs. 5 and 6 respectively show the convergence for TM and TE with mesh resolution using FFT approach; both of them show some oscillations, leading to difficulty in convergence. The effect of averaging at each grid point using a finite mesh is also shown. The convergence is improved by some degree using averaging.

In addition, we compare our results using analytical Fourier coefficients with those by Hermann et al. [9], as shown in Table 1. The computing time for their methods is measured on different computers.

The advantage of using analytical Fourier Transform is evident from Table 1. It requires a small number of plane waves, yet produces accurate results, and converges quickly. The method is even better when some very small features exist in the lattice, where FFT needs larger mesh to reflect all the details.

An example of a heavier computation is shown for a defective photonic crystal: a square lattice with alumina rods in air with the center rod absent. The dielectric constant for alumina is 8.9, and the radius of the rod is 0.2*a*, where *a* is the lattice constant. A 7×7 supercell is used to approximate the crystal structure. The defect frequency in the band gap at k=(0,0,0) was calculated using different number of plane waves. The Fourier transform of this 7×7 supercell is easily obtained using Eq. (5) and Eqs. (11–12). The calculated defect frequencies are listed in Table 2. The iteration error is calculated as *norm*[*X*(*n*)- *X*(*n* - 1)], where *X*(*n*) is the eigenfrequency vector of the first 50 bands for the nth iteration. The band structure of a 7×7 supercell is folded 7^{2} times and the defect band is band 49 in this case.

As illustrated in Table 2, the convergence vs ${N}_{\mathit{\text{PW}}}^{1/2}$ for a large supercell is much slower. Much more plane waves are needed than that for the 1×1 ‘supercell’ to achieve the same accuracy and the computing time increases exponentially with the number of plane waves. Analytical method may save plenty of time by eliminating the converging process of grid resolution.

The defect mode fields *H* or *E* can be obtained conveniently at the same time using the calculated eigen-vector *h*
_{G,λ} and Eq. (2). Defects with one or more cylinders of different sizes [11] or dielectric constants can be treated as different ‘atoms’ and their Fourier transform are obtained the same way as Eq. (5) using different *R* or ε and then added to the supercell according to Eqs. (11–12). Similar mode fields as in Ref. [11] can be obtained, but are not shown here. Other defects such as waveguides can be treated in the same way.

To reduce the interaction of the neighboring defects, a large supercell is generally needed, but the computation will increase accordingly. The convergence curve of the defect frequency using different supercell size is shown in the Fig. 7. The defect frequency for a point defect should be independent of k-vectors for an infinitely large supercell; however for small supercells, coupling between neighboring defects will lead to a finite width of the defect frequency. For large supercells, the interaction between neighboring defects becomes smaller.

## 5. Discussion

Though plane wave method is quite successful in PBG calculations, it has several limitations. The computation grows exponentially when the problem size increases. For complicated problems, such as 3D PBG calculations, the computation is intensive. A less computationally intensive eigen-solver is critical to reduce computation effectively. A few good eigen solvers like Lanzcos and subspace methods may be able to greatly improve the performance from *O*(*n*
^{3}) to *O*(*n*
^{2}) or less, a good example is Ref. [1].

In conclusion, we provide a simple implementation of the plane wave method using MATLAB. Our implementation yields good convergence and accurate results, and the programming effort is minimized. We demonstrated how using analytical Fourier coefficients could improve accuracy and reduce the computation time. The shift property is employed to get the Fourier transform of a supercell, thereby reducing computation and increasing flexibility in treating different problems while maintaining accuracy. Using coordinate conversion, ‘atoms’ with arbitrary shapes can be dealt with easily.

## Acknowledgments

This research is supported by NASA Langley Research Center through NASA-University Photonics Education and Research Consortium (NUPERC).

## References and links

**1. **S. G. Johnson and J. D. Joannopoulos, “Block iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express **8**, 173–190 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173 [CrossRef]

**2. **K. M. Leung, “Plane wave calculation of photonic band structures” in *Photonic band gaps and localizations*, C. M. Soukoulis. ed. (Plenum Press NY1993).

**3. **K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. **65**, 3152–3155 (1990). [CrossRef]

**4. **R. D. Meade and A. M. Rappe et a1., “Accurate theoretical analysis of photonic band gap materials,” Phys. Rev. B **48**, 8434–8437 (1993). [CrossRef]

**5. **K. M. Leung and Y. F. Liu, “Full vector wave calculation of photonic band structures in FCC dielectric media,” Phys. Rev. Lett. **65**, 2646–2649 (1990). [CrossRef]

**6. **D. C. Champeney, *Fourier transforms and their physical applications* (Academic Press, 1973) Chap. 3.

**7. **Geoge Arfken, *Mathematical methods for physicists*, 3^{rd} edition, (Academic Press, 1985), Chap. 2.

**8. **A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. **43**, 773–793 (1996) [CrossRef]

**9. **D. Hermann et al., “Photonic band structure computations,” Opt. Express **8**, 167–172 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-167 [CrossRef]

**10. **J. D. Joannopoulos et al., *Photonic crystals - Molding the flow of light* (Princeton University Press1995).

**11. **P. R. Villeneuve and S. Fan et al., “Microcavities in photonic crystals: mode symmetry, tunability, and coupling efficiency,” Phys. Rev. **B 54**, 7837–7842 (1996).