Large full band gaps for photonic crystals in two dimensions computed by an inverse method
with multigrid acceleration
R. L. Chern,1,2C. Chung Chang,2Chien C. Chang,2and R. R. Hwang1 1Institute of Physics, Academia Sinica, Taipei 115, Taiwan, Republic of China
2Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan, Republic of China 共Received 21 April 2003; revised manuscript received 12 June 2003; published 29 August 2003兲 In this study, two fast and accurate methods of inverse iteration with multigrid acceleration are developed to compute band structures of photonic crystals of general shape. In particular, we report two-dimensional pho-tonic crystals of silicon air with an optimal full band gap of gap-midgap ratio⌬/mid⫽0.2421, which is 30% larger than ever reported in the literature. The crystals consist of a hexagonal array of circular columns, each connected to its nearest neighbors by slender rectangular rods. A systematic study with respect to the geometric parameters of the photonic crystals was made possible with the present method in drawing a three-dimensional band-gap diagram with reasonable computing time.
DOI: 10.1103/PhysRevE.68.026704 PACS number共s兲: 02.70.Bf, 42.70.Qs, 78.20.Bh
I. INTRODUCTION
Photonic crystals with large full band gaps are of aca-demic and practical interest关1,2兴. Large full band gaps allow strong photon localization with the gap 关3,4兴, and a detailed manipulation of photonic defect states 关5,6兴. They have im-portant applications such as defect cavities 关7兴, optical waveguides关8兴, defect-mode photonic crystal lasers 关9兴, and feedback mirror in laser diodes关10兴.
Because of the diversity of photonic crystals, there is a great demand for a general method for fast and accurate pre-diction of their band-gap structures. The most widely used method for this purpose is plane wave expansion 关11–13兴, which has indeed provided useful information to band-gap structures. But there are several known disadvantages with this method: slow convergence关14,15兴, the need to interpo-late the dielectric property关16,17兴, and sophisticated skills to diagonalize a dense matrix关18兴. The purpose of this study is therefore twofold: one is to develop fast algorithms for com-puting band structures of photonic crystals, the other is to propose possible large full band-gap structures, the analysis of which requires major computational efforts.
Among various other possible techniques—transfer ma-trix method 关19,20兴, multiple scattering method 关21,22兴, fi-nite difference time domain method 关23,24兴, and finite ele-ment method 关25,26兴—finite difference multigrid method
关27,28兴 has the potential to meet the above standards. On the
other hand, the studies 关29,30兴 provided an interesting algo-rithm toward optimization of band structures, but only in E polarization or H polarization. In the present study, a highly fast and accurate inverse method with multigrid acceleration is developed to be applicable to photonic crystals comprising of general-shape dielectric-dielectric materials. The opera-tion count is a good order N 共the matrix size兲, where the proportional constant is dependent upon the shape and prop-erty of the materials used. Moreover, the developed algo-rithm is capable of resolving multieigenvalue band struc-tures, and applicable to photonic crystals with interfaces of strong contrast. These salient features of the present method make it a highly efficient computational technique to analyze
photonic crystals with two or more geometric parameters. In this study, we further propose to consider a hexagonal lattice shown in Fig. 1. The reason of choosing this geometry is: previous studies indicate that band gaps for E polarization are favored in a lattice of isolated high- region, and band gaps for H polarization are favored in a connected lattice
关31兴. A compromise must therefore be reached between the
sizes of the dielectric column and the connecting rod in order to have a large full band gap for a given column radius. This is a problem with two geometrical parameters, i.e., the radius of the circular columns and the width of the connecting rods. Searching for their optimal values is often computationally demanding, but with the presently developed methods, this can be done with great efficiency. The rest of the paper pro-vides a detailed description of the numerical algorithms and results of the computed band structures.
II. BASIC EQUATIONS AND NUMERICAL METHODS
For linear isotropic and frequency-independent dielectric materials with permeability close to one, the time-harmonic modes in two dimensions for E polarization 共TM兲 can be written as ⫺1
冉
2 x2⫹ 2 y2冊
E⫽冉
c冊
2 E 共1兲and H polarization共TE兲 as ⫺
冋
x冉
1 x冊
⫹ y冉
1 y冊册
H⫽冉
c冊
2 H, 共2兲where E and H are the electric and magnetic field intensities, respectively, and ⫽(r) is the dielectric function. To dis-cretize Eqs. 共1兲 and 共2兲, a second-order central difference scheme is used. At the grid point near the interface, the di-electric function for E polarization is interpolated as关16,18兴
⫽ faa⫹ fbb 共3兲
and for H polarization as 1 ⫽fa 1 a⫹ fb 1 b , 共4兲
whereaandbare dielectric constants of materials a and b,
respectively, and fa and fb are fractions of the grid cell
which contain a andb, respectively. The domain of
com-putation is chosen as a rectangle with the same area of one primitive cell of the hexagonal lattice as shown in Fig. 2. Bloch’s theorem is applied at the domain boundary:
Ek共r⫹ai兲⫽eik•aiEk共r兲, 共5兲
Hk共r⫹ai兲⫽eik•aiHk共r兲, 共6兲
where Ek and Hk are the Bloch functions for electric and magnetic fields, respectively, associated with the wave vector
k in the first Brillouin zone, and ai(i⫽1,2) is lattice
transla-tion vector.
From a practical point of view, the first few branches of eigenvalues are of primary interest. As a first step, it is natu-ral for us to propose the method of inverse iteration to com-pute the eigenvalues as well as eigenvectors.
Let A be the discretization matrix of the differential op-erator in Eqs. 共1兲 and 共2兲. The basic idea is to solve (A
⫺I)x⫽b by inverse iteration 关32,33兴, whereis chosen to be close to the eigenvalue one wishes to compute. The in-verse algorithm is as follows:
InverseEigenˆ for n⫽1 to S
Initial guess b
do
Solve (A⫺I)x⫽b by LUD or PCG
Deflate x by q1 to qn⫺1
Set b⫽x/储x储
Rayleigh Quotient n⫽
具
b,Ab典
until储(A⫺nI)b储2⬍⑀ Set qn⫽b
end ‰
where 1⭐2⭐•••⭐S is the sequence of smallest
eigen-values, and q1, . . . ,qS are the corresponding eigenvectors.
The inner product
具
•,•典
is defined as具
x,y典
⫽兺
i
wi¯xiyi,wi⫽
再
共ri兲 for E polarization
1 for H polarization in view of the different orthogonal properties of E field and
H field: 1 V
冕
V Em*共r兲共r兲En共r兲dr⫽␦mn, 共7兲 1 V冕
V Hm*共r兲Hn共r兲dr⫽␦mn, 共8兲 where En and Hn are the nth eigenmodes of the electric and magnetic fields, respectively. V denotes the volume of the primitive cell and ␦mnis the Kronecker delta. LUD denotesthe lower-upper decomposition solver for matrix inversion and PCG the preconditioned conjugate gradient method.
TABLE I. Computing times in CPU seconds for six bands.
Ngrid/Nwave FMGI PWE
32⫻32 0.03 5.80 48⫻48 57.06 64⫻64 0.17 339.53 128⫻128 0.74 12836共est.兲 256⫻256 2.74 512⫻512 9.39 1024⫻1024 33.59
FIG. 3. Comparison of PWE, a plain idea and a good idea. The plain idea embeds the multigrid solver in the inverse iteration loop, while the good idea embeds the inverse iteration in the multigrid acceleration structure.
Next, we propose two ideas to accelerate the convergence of matrix inversion in the algorithm. A plain idea is to embed a multigrid solver in the inverse iteration loop, namely, mul-tigrid takes place of the expensive LUD. A good idea is to embed inverse iteration in the multigrid acceleration struc-ture. In this idea, the Rayleigh quotient is updated at each grid level, which results in a more efficient algorithm. The basic idea of multigrid method关34,35兴 is to solve the matrix problem by approximating the solution on a fine grid, solv-ing the residue on a coarse grid, and then improvsolv-ing the solution on the fine grid. Successively applying this idea on each level of grids causes relaxation of errors on different resolution, and hence accelerates the convergence.
Let there be L levels of grids. On each grid level m, the differential operator is discretized to form the matrix A(m). The multilevel algorithm based on the good idea is as fol-lows: FMGInverse ˆ Call InverseEigen at m⫽1 for m⫽2 to L for n⫽1 to S Interpolate qn(m⫺1) to b(m) do Solve (A(m)⫺I)x(m)⫽b(m) by MGV Deflate x(m) by q1(m) to qn(m)⫺1 Set b(m)⫽x(m)/储x(m)储 Rayleigh Quot.n(m)⫽
具
b(m),A(m)b(m)典
until储(A(m)⫺n(m)I)b(m)储2⬍⑀Set qn(m)⫽b(m)
end end ‰
where MGV denotes the well-known multigrid V-cycle solver in numerical linear algebra 关34,35兴.
The whole idea of proposing an inverse method with mul-tigrid acceleration for computing photonic band structures has now become clear. First of all, the interface interpola-tions 共3兲 and 共4兲 are used to take care of strong dielectric contrasts. Second, the deflation in the algorithm of inverse iteration enables singling out the eigenvalues one by one from the smallest even when they may be degenerated. Third, interlacing the inverse iteration and multigrid accel-eration makes the whole method an amazingly fast algorithm for computing photonic band structures, as demonstrated be-low.
III. RESULTS AND DISCUSSION
As an example of test, we compute the first six frequency bands at the point X of a square array of dielectric columns with radius r/a⫽0.45 and dielectric contrast /0⫽13. Table I shows a comparison of computing times in central processing unit 共CPU兲 seconds between the present method
共FMGI兲 and plain plane wave expansion 共PWE兲. The
com-parison is made by averaging the computing times over 68 wave numbers. Both computations are performed on a PC FIG. 4. Band structures by 256⫻256 grid for the hexagonal
lattice in Fig. 1 (/0⫽13, r/a⫽0.155, and d/a⫽0.035).
TABLE II. Convergence test for maximum band gap against the grid size. Ngrid 642 1282 2562 5122 10242 au p/2c 0.4917 0.4940 0.4940 0.4941 0.4940 alow/2c 0.3946 0.3876 0.3873 0.3873 0.3873 amid/2c 0.4431 0.4408 0.4407 0.4407 0.4407 a⌬/2c 0.0971 0.1064 0.1067 0.1068 0.1067 ⌬/mid 0.2191 0.2413 0.2421 0.2423 0.2421
TABLE III. Comparison of the gap-midgap ratios for two-dimensional photonic crystals.
/0 amid/2c a⌬/2c ⌬/mid Present 13 0.4407 0.1067 0.2421 11.4 0.4503 0.0888 0.1972 关31兴 13 0.186 关36兴 11.4 共0.71兲 0.1158 共0.1631兲 关37兴 12.25 0.0872 0.182 关38兴 11.4 (⬃0.45) 0.0762 (⬃0.1693) 关39兴 12.96 0.090 0.171 关40兴 11.4 ⬎0.9 0.0967 ⬍0.1074 ⬎0.8 0.0849 ⬍0.1061
FIG. 5. A map of band gaps by 128⫻128 grid for the hexagonal lattice in Fig. 1 by varying radius r/a of circular columns. The dielectric contrast/0is 13 with fixed d/a⫽0.035.
Pentium 4. The computing time for the present method is of a good order N⫽Ngrid, and tends to be superlinear at large
N. The plane wave expansion, on the contrary, has the
com-putational cost of about the order N2.873. Figure 3 shows a comparison of PWE, the plain idea and the good idea of implementing inverse iteration with multigrid acceleration. The good-idea algorithm is roughly four times faster than the plain-idea algorithm, and even better at larger N. Neverthe-less, both of them are faster than PWE by almost two orders in Ngrid.
The fast algorithm developed above enables efficient search for the optimal geometry of the photonic crystals in Fig. 1. The largest gap-midgap ratio for silicon air (/0
⫽13) corresponds to r/a⫽0.155 and d/a⫽0.035. Figure 4
shows the detailed band structure at these ratios. It is inter-esting to notice that the full band gap is almost the
simulta-neous maximum band gap of both E and H polarizations. In
order to ensure the accuracy of the band gap, computations are performed on five different grids. Table II lists the
nu-merics of the location of the maximum full band gap for different grids for N⭓128⫻128. In this paper, most of the computations are performed with the grid N⫽128⫻128. All the computed results coincide with each other to three sig-nificant digits. Table III has a list of references that reported numerics of large band gaps. It is found that for GaAs air (/0⫽11.4) the present ⌬/mid⫽0.1972 is the largest, while for silicon air the present ⌬/mid⫽0.2421 is not
only the largest but 30% larger than reported in Ref.关31兴.
Figure 5 is a gap map versus r/a 共with fixed d/a
⫽0.035). This map is obtained by computing the band
struc-tures at every 0.0025 of r/a. Figure 6 is a gap map versus
d/a 共with fixed r/a⫽0.155). This gap map is obtained by
computing the band structures at every 0.0025 of d/a. The two maps show clearly the location of the full band gap with maximum gap-midgap ratio. In particular, by increasing d/a from 0.035, the full band gap is gradually embedded into a large band gap of H polarization. By varying r/a near 0.155, the band gaps of E polarization and of H polarization have a good overlap with each other. It is noted that addition of a smaller column in the center of the unit cell of a square/ honeycomb lattice for symmetry reduction can increase the full band gap关40兴. However, the present study does not find further increase of the full band gap by inserting a circular column, no matter what its radius is, in the center of the photonic structure of Fig. 1. Figure 7 shows the gap map. As
ri/a (ri: the radius of the inserted column兲 increase from 0, the full band gap shrinks gradually in size due to the de-crease of the band gap in E polarization. Figure 8 is the result of all elaboration plotted in a three-dimensional band-gap diagram. It gives an overall outlook of the major full band gap, and shows the most sensitive direction to geometric parameters. This plot is composed of as many band gaps as computed for 61 values of r/a and 41 values of d/a, each for
E and H polarizations. There are totally 60 024 eigenvalues
computed.
IV. CONCLUDING REMARKS
In conclusion, we have developed two fast inverse meth-ods with multigrid acceleration for computing photonic band FIG. 6. A map of band gaps by 128⫻128 grid for the hexagonal
lattice in Fig. 1 by varying width d/a of connecting rods. The dielectric contrast/0is 13 with fixed r/a⫽0.155.
FIG. 7. A map of band gaps by 128⫻128 grid for the hexagonal lattice in Fig. 1 by varying radius ri/a of the inserted center col-umns. The dielectric contrast/0is 13 with fixed r/a⫽0.155 and d/a⫽0.035.
FIG. 8. A three-dimensional terrain of the major full band gaps for 0⭐r/a⭐0.3 and 0⭐d/a⭐0.2.
structures in two dimensions. The fast algorithms are capable of resolving multieigenvalue band structures, and applicable to photonic crystals with interfaces of strong contrast. In par-ticular, the good idea algorithm is even faster than the plane idea algorithm. The methods enable us to explore photonic band structures with two or more geometric parameters sys-tematically. In addition, we propose a photonic structure that has significantly larger full band gaps compared to those ever reported in the literature. The study has illustrated that opti-mal design of photonic structures can benefit greatly from physical principle, aided by fast computation. Presently, the methods are applied to two-dimensional photonic crystals. In
principle, they can be extended to problems in three dimen-sions. Three-dimensional problems are, however, more chal-lenging by presenting themselves as eigenproblems of vector operators. One such extension is currently being under inten-sive investigation, and the results will be reported elsewhere.
ACKNOWLEDGMENT
This work was supported in part by the National Science Council of the Republic of China under Contract No. NSC 90-2212-E-002-238.
关1兴 E. Yablonovitch, Phys. Rev. Lett. 58, 2059 共1987兲. 关2兴 S. John, Phys. Rev. Lett. 58, 2486 共1987兲.
关3兴 N. Garcia and A.Z. Genack, Phys. Rev. Lett. 66, 1850 共1991兲. 关4兴 S. John, Phys. Today 44共5兲, 32 共1991兲.
关5兴 S. Lin and G. Arjavalingam, J. Opt. Soc. Am. B 11, 2124 共1994兲.
关6兴 T. Krauss et al., Nature 共London兲 383, 699 共1996兲.
关7兴 J.S. Foresi, P.R. Villeneuve, J. Ferrera, E.R. Thoen, G. Stein-meyer, S. Fan, J.D. Joannopoulos, L.C. Kimerling, H.I. Smith, and E.P. Ippen, Nature共London兲 390, 143 共1997兲.
关8兴 A. Mekis, J.C. Chen, I. Kurland, S. Fan, P.R. Villeneuve, and J.D. Joannopoulos, Phys. Rev. Lett. 77, 3787共1996兲. 关9兴 O. Painter, R.K. Lee, A. Scherer, A. Yariv, J.D. OBrien, P.D.
Dapkus, and I. Kim, Science 284, 1819共1999兲.
关10兴 D.L. Bullock, C. Shih, and R.S. Margulies, J. Opt. Soc. Am. B
10, 399共1993兲.
关11兴 K.M. Leung and Y.F. Liu, Phys. Rev. Lett. 65, 2646 共1990兲. 关12兴 Z. Zhang and S. Satpathy, Phys. Rev. Lett. 65, 2650 共1990兲. 关13兴 K.M. Ho, C.T. Chan, and C.M. Soukoulis, Phys. Rev. Lett. 65,
3152共1990兲.
关14兴 H.S. So¨zu¨er, J.W. Haus, and R. Inguva, Phys. Rev. B 45, 13 962共1992兲.
关15兴 L. Shen and A. He, J. Opt. Soc. Am. A 19, 1021 共2002兲. 关16兴 R.D. Meade, A.M. Rappe, K.D. Brommer, J.D. Joannopoulos,
and O.L. Alherhand, Phys. Rev. B 48, 8434共1993兲. 关17兴 P. Lalanne, Phys. Rev. B 58, 9801 共1998兲.
关18兴 S.G. Johnson and J.D. Joannopoulos, Photonic Crystals 共Klu-wer Academic, Boston, 2002兲.
关19兴 J.B. Pendry and A. MacKinnon, Phys. Rev. Lett. 69, 2772 共1992兲.
关20兴 M.M. Sigalas, C.T. Chan, K.M. Ho, and C.M. Soukoulis, Phys. Rev. B 52, 11 744共1995兲.
关21兴 X. Wang, X.-G. Zhang, Q. Yu, and B.N. Harmon, Phys. Rev. B
47, 4161共1993兲.
关22兴 K.M. Leung and Y. Qiu, Phys. Rev. B 48, 7767 共1993兲. 关23兴 C.T. Chan, Q.L. Yu, and K.M. Ho, Phys. Rev. B 51, 16 635
共1995兲.
关24兴 A.J. Ward and J.B. Pendry, Phys. Rev. B 58, 7252 共1998兲. 关25兴 D.C. Dobson, J. Comput. Phys. 149, 363 共1999兲.
关26兴 W. Axmann and P. Kuchment, J. Comput. Phys. 150, 468 共1999兲.
关27兴 D. Hermann, M. Frank, K. Busch, and P. Wolfle, Opt. Express
8, 167共2001兲.
关28兴 A. Brandt, S. McCormick, and J. Ruge, SIAM 共Soc. Ind. Appl. Math.兲 J. Sci. Stat. Comput. 4, 244 共1983兲.
关29兴 S.J. Cox and D.C. Dobson, SIAM 共Soc. Ind. Appl. Math.兲 J. Appl. Math. 59, 2108共1999兲.
关30兴 S.J. Cox and D.C. Dobson, J. Comput. Phys. 158, 214 共2000兲. 关31兴 J.D. Joannoupoulos, R.D. Meade, and J.N. Winn, Photonic
Crystals共Princeton University Press, Princeton, 1995兲. 关32兴 B.N. Parlett, The Symmetric Eigenvalue Problem 共SIAM,
Philadelphia, 1998兲.
关33兴 G.H. Golub and C.F. Van Loan, Matrix Computations, 3rd ed. 共Johns Hopkins University Press, London, 1996兲.
关34兴 J.W. Demmel, Applied Numerical Linear Algebra 共SIAM, Philadelphia, 1997兲.
关35兴 U. Trottenberg, C. Oosterlee, and A. Schu¨ller, Multigrid 共Aca-demic Press, London, 2001兲.
关36兴 L. Shen, S. He, and S. Xiao, Phys. Rev. B 66, 165315 共2002兲. 关37兴 N. Susa, J. Appl. Phys. 91, 3501 共2002兲.
关38兴 M. Qiu and S. He, J. Opt. Soc. Am. B 17, 1027 共2000兲. 关39兴 X.H. Wang, B.Y. Gu, Z.Y. Li, and G.Z. Yang, Phys. Rev. B 60,
11 417共1999兲.
关40兴 C.M. Anderson and K.P. Giapis, Phys. Rev. Lett. 77, 2949 共1996兲.
