An introduction to the plane wave expansion method for calculating photonic crystal band diagrams


Aaron J. Danner

Created at the University of Illinois at Urbana-Champaign, Urbana, IL 61801 in 2002

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:







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:






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.)












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: 



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.






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.




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




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 :






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: 




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): 




This is equivalent to the following, where the function : 




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






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










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: 





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: 




The matrix is given by the following, where :





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




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




In the derivation of Equation (32), the following property has been used, where  is the nth-order Bessel function: 






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]





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





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: 





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: 





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





(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: 






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 Crystals:  Molding the Flow of Light. Princeton, NJ:  Princeton University Press, 1995.

[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, University of Michigan, 2001.

[8]   K. Sakoda, Optical Properties of Photonic Crystals. Berlin, Germany:  Springer, 2001.

[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 Crystals:  The Road from Theory to Practice. Boston, MA:  Kluwer Academic Publishers, 2002.

[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, Quebec, 2002.

[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 Jilin University, China kindly provided corrections to Equations 18 and 19 and the text after Equation 15.