OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 11, Iss. 2 — Jan. 27, 2003
  • pp: 167–175
« Show journal navigation

Simple plane wave implementation for photonic crystal calculations

Shangping Guo and Sacharia Albin  »View Author Affiliations


Optics Express, Vol. 11, Issue 2, pp. 167-175 (2003)
http://dx.doi.org/10.1364/OE.11.000167


View Full Text Article

Acrobat PDF (151 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

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] [PubMed]

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

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

5

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] [PubMed]

]. 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:

×1ε(r)×H(r)=ω2c2H(r)
(1)

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:

H(r)=G,λhG,λêλei(k+G)·r
(2)

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

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

ε(r)=GεGeiG·rεG=1VΩε(r)eiG·rdΩ
(3)

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

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

]:

Gk+Gk+Gε1(G-G)[ê2·ê2ê2·ê1ê1·ê2ê1·ê1][h1,Gh2,G]=ω2c2[h1,Gh2,G].
(4)

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 ê1, ê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

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

].

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 a1 and a2. The area of the unit cell is calculated as A = |a1 × a2|, the Fourier transform of the unit cell is:

ε(G)=εbδ(G)+(εaεb)2πR2AJ1(GR)GR=εbδ(G)+2(εaεb)fJ1(GR)GR
(5)

where J 1 is the 1st order Bessel function, G is the modulus of G⃗, ƒ is a fraction parameter:

f=VolatomVolcell.
(6)

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 a1, a2, a3 are:

a1=a1xx̂+a1yŷ+a1zẑ
a2=a2xx̂+a2yŷ+a2zẑ,
a3=a3xx̂+a3yŷ+a3zẑ
(7)

The dielectric function in the unit cell is ε(r) and the column vector is r⃗ = ma1 + na2 + la3, where m, n and l are coordinates along the basis vectors. In Cartesian coordinates:

r=(ma1x+na2x+la3x)x̂+(ma1y+na2y+la3y)ŷ+(ma1z+na2z+la3z)ẑ
(8)

So we have:

r2=rT[g]r
(9)

where

[g]=[a1xa1ya1za2xa2ya2za3xa3ya3z][a1xa2xa3xa1ya2ya3ya1za2za3z]
(10)

[g] is called the metric for the oblique coordinates [7

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

, 8

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

].

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.

Fig. 1. (a) Conversion of a triangular lattice with a circular cylinder ‘atom’; (b) Conversion of a triangular lattice with an elliptical cylinder ‘atom’; (c) TE band structure of a 2D triangular lattice with elliptical air holes in GaAs. Data used: Rx = 0.28a, Ry = 0.14a, εa=13, εb=1.0.

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 eik·r0 [6

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

] which we call shift property here:

ε(r+r0)eiG·r0εG.
(11)

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

riε(r+ri)rieiG·riεG,
(12)

where r⃗ i is the location of ‘atom’ i in the supercell.

Equations (1112) 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 a1 = [0,1,1]a/2,a2 = [1,0,1]a/2,a3 = [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

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

].

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

ε(G)=3f(εaεb)(sinGRGRcosGR(GR)3)cos(G·r0)
(13)

where R is the radius of the sphere, f=2×43πR3V V = |a1 · a2 × a3|.

In Fig. 2 we show the band structure for R=√3/8a, εa=13, εb=1 using our simple program with 343 plane waves, which is in excellent agreement with the result in Ref. [3

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] [PubMed]

]. In this paper, the set of G is chosen as G⃗=n 1b⃗1 + n 2b⃗ + n 3b⃗3,|ni |≤n and in this example n=3.

Fig. 2. Band structure of a 3D diamond lattice using 343 plane waves for this calculation, the inset shows the unit cell of the diamond lattice.

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.28a, 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 nth iteration, and the number of plane waves used is N PW=(2n+1)2.

Fig. 3. Convergence of TM mode. (a) Convergence of the first band. (b) The iteration errors for the first 10 bands. A uniform mesh with different resolution is used to represent the unit cell, and each grid is averaged by a 10×10 submesh.
Fig. 4. Convergence of TE mode. (a) Convergence of the first band. (b) The iteration errors for the first 10 bands. A uniform mesh with different resolution is used to represent the unit cell, and each grid is averaged by a 10×10 submesh.

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

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]

], 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.

Fig. 5. Eigen frequency convergence as a function of grid resolution for TM mode in a 2D triangular lattice. 225 plane waves are used for this calculation. Line with ‘o’: grid is averaged by a 10×10 submesh; line with ‘+’:not averaged.
Fig. 6. Eigen frequency convergence as a function of grid resolution for TE mode in a 2D triangular lattice. 225 plane waves are used for this calculation. Line with ‘o’: grid is averaged by a 10×10 submesh; line with ‘+’: not averaged.

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

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] [PubMed]

], as shown in Table 1. The computing time for their methods is measured on different computers.

Table 1. Comparison of several methods

table-icon
View This Table
| View All Tables

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.2a, 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. (1112). 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 72 times and the defect band is band 49 in this case.

Table 2. Defect frequency of TM mode in a 2D square lattice using a 7×7 supercell

table-icon
View This Table
| View All Tables

As illustrated in Table 2, the convergence vs NPW1/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

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

] 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. (1112). 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.

Fig. 7. Convergence of defect frequency for TM mode using different supercell size in a square lattice with the center rod being removed

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

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] [PubMed]

].

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] [PubMed]

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] [PubMed]

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] [PubMed]

6.

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

7.

Geoge Arfken, Mathematical methods for physicists, 3rd 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] [PubMed]

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

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(350.3950) Other areas of optics : Micro-optics

ToC Category:
Research Papers

History
Original Manuscript: December 19, 2002
Revised Manuscript: January 21, 2003
Published: January 27, 2003

Citation
Shangping Guo and Sacharia Albin, "Simple plane wave implementation for photonic crystal calculations," Opt. Express 11, 167-175 (2003)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-11-2-167


Sort:  Journal  |  Reset  

References

  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), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173</a> [CrossRef] [PubMed]
  2. K. M. Leung, "Plane wave calculation of photonic band structures" in Photonic band gaps and localizations, C. M. Soukoulis. ed. (Plenum Press NY 1993).
  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] [PubMed]
  4. R. D. Meade, A. M. Rappe et al., �??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] [PubMed]
  6. D. C. Champeney, Fourier transforms and their physical applications (Academic Press, 1973) Chap. 3.
  7. Geoge Arfken, Mathematical methods for physicists, 3rd edition, (Academic Press, 1985), Chap. 2.
  8. A. J. Ward, J. B. Pendry, �??Refraction and geometry in Maxwell�??s equations,�?? J. Mod. Opt. 43, 773-793 (1996) [CrossRef]
  9. D. Hermann, M. Frank, K. Busch, and P. Wolfle,�??Photonic band structure computations,�?? Opt. Express 8, 167-172 (2001), <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-167">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-167</a> [CrossRef] [PubMed]
  10. J. D.. Joannopoulos et al., Photonic crystals - Molding the flow of light ( Princeton University Press 1995).
  11. P. R. Villeneuve, S. Fan et al., �??Microcavities in photonic crystals: mode symmetry, tunability, and coupling efficiency,�?? Phys. Rev. B 54, 7837-7842 (1996).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited