An introduction to the plane wave expansion method for calculating photonic crystal band diagrams
Aaron J. Danner
Created
at the
Last edited on January 31, 2011
Abstract: Starting from basic undergraduate-level electromagnetics (Maxwell’s equations) a simple method of finding band diagrams is described. Using the information here, it should be possible to fully understand what band diagrams mean, what they describe, and how they can be calculated. A simple program could easily be written. While not the most efficient method computationally, this is a good introduction to the business of band diagrams and is probably the easiest method to understand.
1. Useful downloads
Right click to download a Mathematica program useful in performing planewave expansion method calculations.
Left click to go to the download site for MathReader, free software which allows Mathematica code to be viewed for study even if Mathematica is not installed.
2. Basic description of the plane wave expansion
method
In order to design photonic
crystals to take advantage of their unique properties, a calculation method is
necessary to determine how light will propagate through a particular crystal
structure. Specifically, given any
periodic dielectric structure, we must find the allowable frequencies
(eigenfrequencies) for light propagation in all crystal directions and be able
to calculate the field distributions in the crystal for any frequency of
light. There are several capable
techniques, but one of the most studied and reliable methods is the plane wave
expansion method. It was used in some of
the earliest studies of photonic crystals [1-4] and is simple enough to be
easily implemented. The method allows
the computation of eigenfrequencies for a photonic crystal to any prescribed
accuracy, commensurate with computing time.
For photonic crystal applications in semiconductors such as GaAs or
dielectrics, Maxwell’s equations, which will govern all field simulations that
follow, take the following forms, where 
, 
, 
,
and 
are the field vectors, 
is the current density, 
is time, and 
is the charge density:
|
|
(1) |
|
|
(2) |
|
|
(3) |
|
|
(4) |
Assuming
perfect dielectric materials ( the relative permeability
)
in a source-free region (
and 
),
Maxwell’s equations can be reduced to four equations, each involving only one type
of field. This decoupling of the fields can be accomplished by taking the curl
of both sides of Equation (2),
and substituting from Equation (4)
to give the two electric field equations.
An equivalent process can be carried out in the opposite order to give
the two magnetic field equations. If it
is assumed that the fields are time-harmonic, then 
and
the decoupled equations can be expressed as the following:
|
|
(5) |
|
|
(6) |
|
|
(7) |
|
|
(8) |
Note that the forms for and
are identical.
This is expected because 
is constant in the equations. However, the relative permittivity 
is not constant in our structures (it is
periodic), so the placement of the 
term is strict. The goal here is to find the energies and
electromagnetic field configurations that are allowed to exist in a periodic
structure. Essentially, what we are
given in the problem is the ,
which will be a function of location, and we need to solve for 
and
the fields.
There
are essentially three different choices of procedure at this point. All four equations, given a dielectric
function, will yield one set of field distributions. (The and
expansions give identical results.) After
that, the other fields can simply be deduced from Maxwell’s equations. The question of which equation to solve
depends on several factors. First, the equations for the magnetic fields
(Equations (7)
and (8))
are in a Hermitian form. Strictly
speaking, the operator 
is Hermitian (see [5] for a detailed description of this property). Hermicity establishes
that the eigenvalues 
are real, and that field distributions with
the same eigenfrequency must be orthogonal.
Usually, Hermitian eigenvalue problems are less complex computationally
to solve [6], but the other forms should not be immediately overlooked as will be clear
in the development which follows.
Each
of the decoupled equations above will yield three component equations if the
vector operations are carried out. In
Cartesian coordinates, they can be expressed as follows for the , ,
and expansions, respectively. Expansions for are
not simplified because in their full forms the extra terms generated by the
inner make the expressions very long. (As seen in Equations (9) - (11), each
equation of consists of four terms on the left side; each
equation of or consists of eight terms; and each equation of would
consist of sixteen terms. The chain rule
applied repeatedly to the inner products creates the extra terms.)
|
|
(9) |
|
|
(10) |
|
|
(11) |
The
fields themselves and the dielectric function can be expanded in Fourier series
along the directions in which they are periodic. This Fourier expansion will be
truncated to a fixed number of terms, limiting the accuracy of the calculation. The truncated problem will yield an
eigenvalue equation for the fields which will allow calculation of the
dispersion curves. It must be pointed
out that regardless of which of the four decoupled equations is solved, the
eigenvalues will be the same. For a
fixed number of terms, the accuracy can be improved by a proper choice. For example, when solving a problem of air
spheres embedded in a dielectric, using the
expansion
would yield much better convergence than the others, while using the or expansions would yield better results for
dielectric spheres in air [6]. The analogy with the two-dimensional structures discussed here
suggests that air cylinders drilled in a dielectric background may be a problem
better suited for calculation using the expansion, in terms of matrix size. The accuracy differences among the three
expansions result from the different resultant spatial orientations and
positions of each field.
The
basic approach for calculating the field distribution and eigenfrequency given
a dielectric function and propagation vector is to first expand and
the three components of the appropriate field vector in Fourier series. These series are then substituted into the
decoupled Maxwell’s equations and the terms are reorganized into an ordinary
eigenvalue problem. When the eigenvalues are calculated employing standard
numerical methods (using a finite-sized matrix formed when the Fourier
expansions are truncated), it is straightforward to use the eigenvalues to find
the allowed propagation frequencies, and the eigenvectors to calculate the
field distributions. The process is best
illustrated by a simple example.
3.
Example: One-dimensional photonic
crystal
The simplest example of a
photonic crystal is a one-dimensional array of air slabs penetrating a
dielectric background. Figure 1 shows
the relevant axes. In this case, we will
consider only waves propagating in the +z direction. In most photonic crystal dispersion curves,
it is usually difficult to distinguish curves as “transverse magnetic"
(TM-like) or “transverse electric" (TE-like), but in this simple case
there are two basic polarizations, viz., 
and
.
Here we consider the . case only, and begin the problem by assuming
that the only field components present are , , and . Although the justification for this may not be immediately apparent,
the symmetry in the problem permits this.
There is also nothing wrong with using all components of and ;
the mode separation is then easily seen. (Indeed, this is the method used in
higher-dimensional photonic crystal problems.) Here, the purpose of the early
simplification is to more clearly illustrate the method.
Figure 1: One-dimensional photonic crystal consisting
of air slabs of width d embedded in a dielectric background with a
periodicity of a.
With only one component,
only a single line in Equation (9) remains, which demonstrates the
justification for using the expansion:
|
|
(12) |
Now, Fourier series
expansions for the field and dielectric can be applied. In this case, a Fourier expansion for the
inverse dielectric function is used. Equivalently, the constant can be
moved to the right side of the equation and could be expanded. This would form a generalized Hermitian
eigenvalue problem, or an ordinary eigenvalue problem if an additional matrix
inversion were carried out in the subsequent step. In the notation that follows, 
will represent all Fourier coefficients. The indices m and n are
integers. The variable 
means
“Fourier coefficients, indexed by the integer n, for the y-component of
the electric field” and the variable 
means “Fourier coefficients, indexed by the
integer m, for .” Ideally, the summations should be infinite, but will be truncated for
computation purposes.
|
|
(13) |
|
|
(14) |
Note that if propagation in
a direction other than z had been included in the formulation, then the
additional terms 
would
have been included in Equation (14).
After
the Fourier expansions are substituted into Equation (12), the initial eigenvalue
equation is obtained.
|
|
(15) |
To simplify, each side of
this equation is multiplied by an orthogonal function 
, where p is an integer, and integrated over a unit
cell, i.e., 
. For a nontrivial solution, p can take
only one value so one summation on each side of the equation can be
dropped. After reorganizing terms and
renaming the sums to use only the letters m and n, the
eigenequation takes its final form of
|
|
(16) |
This forms an ordinary
eigenvalue problem, where the integers m and n are truncated
symmetrically about zero as is appropriate for this type of Fourier
expansion. This corresponds to including
only lower order plane waves in two dimensions.
For example, if m and n were truncated to five terms (-2,
-1, 0, 1, 2), then the full eigenvalue problem would appear as follows, using
the notation 
:
|
|
(17) |
The matrix Q can be diagonalized
using a variety of software packages and numerical methods, and the details
will not be discussed here. After
diagonalization, the eigenvalues and eigenvectors 
will
be known. The eigenvalues give the
dispersion diagram and the eigenvectors can be substituted back into the
Fourier expansion for to find the field distribution at any given
frequency.
The
only remaining problem is to find the dielectric coefficients 
,
which can be obtained using the inverse Fourier transform:
|
|
(18) |
If the integral is split
properly, then 
will be constant in the integration
range. Depending upon where z = 0
is defined, the form of the result may take slightly different forms (but will
make no difference as long as it is defined consistently):
|
|
(19) |
This is equivalent to the
following, where the function 
:
|
|
(20) |
All information is now
available to solve the Q matrix for the eigenvalues. Figure 2 shows the results for nine plane
waves (n and m are integers between -4 and 4 inclusive). In this
case, a structure was chosen with a unit period, 
,
and 
.
Several bandgaps are clearly visible for propagation in this direction. In this example we have examined only the case for propagation in the z
direction. The interested reader is
referred to [7] for a discussion of the case and off-axis propagation in the
one-dimensional structure.
In
studies of photonic crystals, the interest is usually not in the electric or
magnetic field forms themselves. It is
the eigenvalues that carry information on the location of the
modes in momentum space. In the general
case, varying values of , , and allows construction of a complete band diagram. In more complicated structures, the band
diagram is usually constructed at the boundaries of the Brillouin zone.
Figure 2: Dispersion curve for propagation in the z direction for the
one-dimensional photonic crystal structure. Note the presence of several
bandgaps.
4. Fully
vectorial, three-dimensional structures
Essentially, the method
remains the same for more complicated structures. Because of our interest in two-dimensional
structures, we examine here the case of the triangular lattice of air holes
embedded in a dielectric background.
Using fabrication techniques described in the next chapter, arrays of
holes can be created using electron beam lithography and etching methods. The result is a two-dimensional array of air
holes in a semiconductor substrate.
Therefore, for our interests the dielectric function will be periodic
only in the xy plane (uniform in the z direction). This results
in some simplification for the expansion. Although the
structure studied is two-dimensional, propagation in all directions (including
the out-of-plane propagation case) will be considered. Extension of this method to three-dimensional
structures is straightforward and will be explained. In the equations used, a is the
lattice spacing of a unit cell; the lattice itself is triangular within a
medium with dielectric constant 
perforated by infinite air holes (atoms)
of diameter d. The 2D triangular lattice unit cell has been widely
covered in literature [2,3,5,8] and the method described here has been tested and gives equivalent
results for the same problems. In this development,
the triangular lattice supercell is further generalized to the N x N
case and method accuracy is treated.
First, we discuss the unit cell.
The
general formulas for the Fourier expansions in the three-dimensional case,
assuming use of the 
expansion,
are
|
|
(21) |
|
|
(22) |
where the 
vectors are related to the directions of
periodicity. They are actually the
collection of reciprocal lattice vectors; the 
vectors
represent the real lattice vectors, and their relationship is defined by 
[8]. For the structure studied here, the real and reciprocal lattice
vectors are shown in Figures 3 and 4, respectively. Because it is a
two-dimensional structure, there are two reciprocal lattice vectors, 
and 
.
The lattice vectors shown can be expressed as
|
|
(23) |
|
|
(24) |
|
|
(25) |
|
|
(26) |
Figure
3: Real lattice vectors for the 2D
triangular lattice.
Figure 4: Reciprocal lattice vectors for the 2D
triangular lattice.
Equations (21) and (22) now
become the following, specifically for the triangular lattice structure:
|
|
(27) |
|
|
(28) |
At
this point the same process is used as in the one-dimensional case discussed in
Section 2. Equations (27) and (28) are substituted back into the appropriate
decoupled Maxwell equation (Equation (9)). Again, both sides of the resulting
equation are multiplied by an orthogonal function 
and
integrated over a unit cell (see Figure 5). The rectangular area in Figure 5
represents one possible area of integration for the unit cell. After algebraic simplification, the result is
an ordinary eigenvalue equation with form similar to the one encountered in the
one-dimensional case:
|
|
(29) |
The matrix is given by the following, where 
:
|
|
(30) |
Although
Equations (29) and (30) are more complicated to convert into matrix form suitable
for diagonalization, the process is carried out exactly as in the
one-dimensional case and will not be discussed further here. Certainly, more plane waves must be used to
maintain accuracy than in the one-dimensional case, but the method is the same. The calculation of the Fourier dielectric
coefficients 
is also more complicated due to the area
integral, but the method for doing so is equivalent. The general formula for the dielectric
coefficient is given by
|
|
(31) |
where the area of integration
is represented by the rectangular area within the solid lines in Figure 5.
After changing to cylindrical coordinates, the integration is easy to split
into two parts (inside the air holes and in the dielectric background). The
dielectric function in each part then becomes a constant and the integral thus
simplifies for the unit cell to
|
|
(32) |
In
the derivation of Equation (32),
the following property has been used, where 
is the nth-order Bessel function:
|
|
(33) |
Figure 5: Unit cell.
The dotted line represents the actual cell and the solid line represents
the area covered by the integral in the dielectric Fourier coefficient
computation.
Using another property 
results in the simplification [8]
|
|
(34) |
In
previous studies of the triangular lattice unit cell, 140 to 225 plane waves
were used [1,2,3] to calculate accurate dispersion curves. Villeneuve and Piché [2] tested the convergence to 841 plane waves and found that only 225
were necessary for good convergence. The
accuracy of any given set of curves is difficult to predict because the
convergence rates can change between differing structures.
The
plane wave expansion method was carried out to construct dispersion curves for along the Brillouin zone shown in Figure 4. The method described using
the expansion creates spurious modes with zero
frequency, which were removed. Also, in
the in-plane propagation case modes can be separated by polarization into
TE-like (
is in the xy plane) and TM-like (
is in the xy plane) modes. Two examples have been carried out with 
=
13.2 using 441 plane waves. TE-like
modes were separated from the result and are shown in Figures 6 and 7 for
values of and along the Brillouin zone. The
gap between the first and second bands changes with the lattice dimensions,
where d refers to the diameter of the air holes, and a to the
lattice spacing. Figure 8 shows the
variation.
Figure 6: TE-like modes for air holes embedded in a
background of =
13.2 with d/a = 0.5.
Figure 7: TE-like modes for air holes embedded in a
background of =
13.2 with d/a = 0.8.
Figure 8: Variation of TE mode gap with lattice
parameters. The two bands plotted are
the two bands with the lowest eigenfrequencies (=
13.2).
5. Supercell techniques
If a defect is introduced
into the otherwise periodic structure then defect modes can arise in the
photonic band structure. To study defect
modes, the same plane wave expansion method can be used. The basic idea is to replace the unit cell by
a more complicated unit cell and preserve the periodicity. For example, a 4 x 4 supercell with a central
defect can give reasonable accuracy because the missing holes are spaced four
lattice units apart. As long as confined
modes do not couple to one another, then the results of the calculation should
equally apply to the case of an isolated defect (missing hole) in a large array
of perfect photonic crystal. Supercells
are often used to calculate defect states in photonic crystals [9,10], although different authors choose to use different sizes. Although a 4 x 4 supercell is a reasonable
size cell for most calculations, in order to study certain modes with more
accuracy larger supercell structures may be needed. In this thesis, the dielectric coefficients
for the general case of an N x N supercell with a point defect
have been derived.
The
most basic example of a supercell is the unit cell itself (see Figure 5). The
situation for a supercell is shown in Figure 9. The lattice spacing is now
given by Na. For example, a 4 x 4 supercell has an overall periodicity
of 4a with a defect appearing in the lattice once per period (every four
unit cells). The eigenvalue equation taking this into account is identical to
the unit cell case (Equation 29), except
the matrix now includes the effect of the N x N
supercell,
Figure
9: Example of a 4 x 4 supercell. The dotted line represents the supercell
itself, and the solid line represents the area covered by the integral in the
dielectric Fourier coefficient computation.
where .
Analogous to Equation (30) we obtain
|
|
(35) |
The Area factor that
appears outside the integral in Equation (31) is now dependent on the size of
the supercell. For example, in a 4 x 4
supercell this becomes 
.
In addition, it is now more complicated to integrate in cylindrical
coordinates, as many air holes are distributed at points other than the center
of the coordinate system. This can be
accounted for by making the following substitutions into Equation (31):
where 
represents
the position of a hole in the supercell.
This creates constants which can be taken from the main integral and summed
over all the positions of the photonic atoms as follows:
|
|
(36) |
The integral over q depends on the position of
each photonic atom. As illustrated in
Figure 10, depending on the supercell size N, several holes will be cut by
the area of integration, and thus the integral about q will not be from 0 to 2p in each case. Fortunately, the cut atoms can be combined
due to their symmetry such that only one integration needs to be carried
out. To obtain the full equation, it is
helpful to visualize the supercell in three parts: two offset rectangular lattices (each with
period a in x, in y) intermeshed to form the whole photonic atoms around the defect,
and an even number of half atoms at the edges. The numbers and positions of the
atoms change with the supercell size.
One of the rectangular lattices has an even number of atoms on a side
and the other has an odd number of atoms on a side. The odd mesh includes the central whole atoms
which will be removed later. The
following quantities give these numbers, less one, as a function of N:
|
|
(37) |
|
|
(38) |
In Equations (37) and (38),
the Floor function returns the largest integer less than or equal to the argument.
The crucial summation over these photonic atoms and
Figure 10: N x N supercell area of
integration. For a given N, the
four corner atoms that the rectangle cuts will not be present; the whole atoms
inside the rectangle must be included (except the central defect).
the combined half-atoms is
given by
|
|
(39) |
(The latter two large terms
are the combined half-atoms, and the -1 removes the central atom.) When
combined with the main equation, the final coefficients are obtained with
complete generality of supercell size.
When determining the accuracy of the supercell method for a given number
of plane waves, it is sometimes useful to create a supercell with no defects and
compare the results to those of a unit cell.
In this case, a +2 term can be added to Equation (39) to give the proper
summation. (This effectively reinserts
the central atom and the four quarter-atoms at the corners of the region of
integration.) The Bessel function simplifications have been carried out as in
the unit cell case:
|
|
(40) |
An
example of a supercell calculation was carried out for a defect in a 2D
photonic crystal for the out-of-plane propagation case using a 4 x 4
supercell. Figure 11 shows the
fundamental mode in this case. This is a
plot of , or the time-averaged electric field.
Because of the tight confinement of the mode around the defect, a larger
supercell would not give significantly different results in this case. It does become important for higher-order
mode calculations or calculations at lower frequencies, where the confinement
is not so tight. In this case, the confinement
of the energy within a one lattice unit radius is 98.19%. The plot was generated
by substituting the eigenvectors back into the Fourier expansion for the electric field.
Figure 11: Near
field plot of lowest eigenmode (
)
for a structure of d/a = 0.3, 
=
0, 
=
0, 
, =
12.25. 1089 plane waves were used in each direction.
Figure
12 shows the in-plane defect mode calculated using a 4 x 4 supercell. In this case, the defect mode is superimposed
on the TE-like modes of the unit cell.
Folded bands from the supercell itself are not shown for clarity. Note that the defect band appears almost in
the middle of the gap. This demonstrates
that adding a defect introduces a localized confined state that is not present
in the bulk photonic crystal. The
dispersion curve of the mode is independent of frequency because of this
localization.
Figure 12:
In-plane defect mode of a 4 x 4 supercell (with central defect)
superimposed on TE-like modes of a unit cell for air holes embedded in a
background of = 13.2 with d/a = 0.8. (Folded
bands are not shown.)
6. Control of Accuracy
As increasing numbers of
planewaves are used, the eigenvalues approach the correct values
asymptotically. For the 4 x 4 supercell,
Figure 13 shows this behavior. For this
case, i = 12 should give sufficient accuracy for band diagrams around
the point of calculation shown. (This
should correspond to i = 3 for a unit cell structure.) At high
frequencies the accuracy for a given number of planewaves decreases, and the
size of the supercell, if too small, will give incorrect defect mode eigenvalues
because of coupling between adjacent supercells.
Figure 13: Plots
of the four lowest eigenfrequencies (bands 1, 3, 5, 7 in ascending order) as a
function of the (i x i) submatrix size, which is related to the
number of plane waves given by 
.
The structure is a 4 x 4 supercell with d/a = 0.20, = 12.25, and the position of calculation is =
0, =
0, 
.
With increasing numbers of plane waves used, each of the four curves approaches
its actual value asymptotically.
Still,
the question remains which expansion (, ,
or )
will give greater accuracy for a fixed matrix size. The relationship was analyzed in [11] and it was found that for air holes in a square lattice the expansion
gave consistently better convergence results, even when large numbers of plane
waves (1000) were used. It is evident that for structures presented here, the expansion should be used for this method.
Other
methods exist for solving the eigenequation.
Instead of solving for all the eigenfrequencies at once, iterative
techniques can be used [5,9,12] to find eigenvalue and eigenvectors pairs one at a time. These methods seem to rely on ordinary
Hermitian eigenvalue problems, so the is
exclusively used. Computing time can be
saved by use of the fast Fourier transform to carry out the 
operation
in Fourier space, instead of real space as was done here [9]. These methods usually suffer from poor convergence times at high
frequencies [13].
The
plane wave method presented here can also be extended to calculate transmission
spectra [1,8,14], as well as modal characteristics [15,16].
7. References
[1] M.
Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional
systems: The triangular lattice,” Phys. Rev.
B, vol. 44, no. 16, pp.
8565-8571, 1991.
[2] P.
R. Villeneuve and M. Piché, “Photoinc band gaps in two-dimensional square and
hexagonal lattices,” Phys. Rev. B, vol.
46, no. 8, pp. 4969-4972, 1992.
[3] R.
D. Meade, K. D. Brommer, A. M. Rappe, and J.D. Joannopoulos, “Existence of a
photonic band gap in two dimensions,” Appl.
Phys. Lett., vol. 61, no.
4, pp. 495-497, 1992.
[4] K. M.
Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic
dielectric structures,” Phys. Rev.
Lett., vol. 65, no. 25, pp.
3152-3155, 1990.
[5] J.
D. Joannopoulos, R. D. Meade and J.N. Winn, Photonic
[6] H.
S. Sözüer and J. W. Haus, “Photonic bands:
Convergence problems with the plane-wave method,” Phys. Rev.
B, vol. 45, no. 24, pp.
13962-13972, 1992.
[7] J.
D. Shumpert, “Modeling of periodic dielectric structures (electromagnetic
crystals),” Ph.D. dissertation,
[8] K.
Sakoda, Optical Properties of Photonic
[9] T.
Sřndergaard, “Spontaneous emission in two-dimensional photonic crystal
microcavities,” IEEE J. of Quantum Elect., vol. 36, no.
4, pp. 450-457, 2000.
[10] S.
G. Johnson and J. D. Joannopoulos, Photonic
[11] Z. Y.
Yuan, J. W. Haus, and K. Sakoda, “Eigenmode symmetry for simple cubic lattices
and the transmission spectra,” Optics Express, vol. 3, no.
1, pp. 19-27, 1998.
[12] S.
G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods
for Maxwell’s equations in a planewave basis,” Optics Express, vol. 8, no.
3, pp. 173-190, 2001.
[13] S.
G. Johnson (private communication), June 25, 2002.
[14] A.
Barra, D. Cassagne, and C. Jouanin, “Existence of two-dimensional absolute
photonic band gaps in the visible,” Appl.
Phys. Lett., vol. 72, no.
6, pp. 627-629, 1998.
[15] N. Yokouchi, A. J. Danner, and K. D. Choquette, “Effective index
model of 2D photonic crystal confined VCSELs,” presented at LEOS VCSEL Summer
Topical, Mont Tremblant,
[16] J.
C. Knight, T. A. Birks, R. F. Cregan, P. Russell, J.-P. de Sandro, “Photonic
crystals as optical fibres - physics and applications,” Optical Materials.
vol. 11, pp. 143-151, 1999.
8. Acknowledgement
October,
2007: Jack (Zetao) Ma of Shizuoka University, Japan kindly provided a
correction to Equation 39.
January,
2011: Wang Wei of