Numerical Study of Three-Dimensional Photonic Crystals with Large Band Gaps

全文

(1)Journal of the Physical Society of Japan Vol. 73, No. 3, March, 2004, pp. 727–737 #2004 The Physical Society of Japan. Numerical Study of Three-Dimensional Photonic Crystals with Large Band Gaps R. L. C HERN1;2 , C. Chung C HANG2 , Chien C. C HANG2 y and R. R. H WANG1;3 1 Institute of Physics, Academia Sinica, Taipei 115, Taiwan, Republic of China Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan, Republic of China 3 Department of System Engineering and Naval Architecture, National Taiwan Ocean University, Keelung 202, Taiwan, Republic of China (Received November 5, 2003) 2. The study presents a finite difference formulation for efficiently computing band structures of threedimensional photonic crystals. First of all, we will show how to correctly discretize the double-curl equation for the magnetic field so that the transversality condition is exactly satisfied in the discrete sense. The first few branches of nontrivial eigenfrequencies that determine the major full band gaps of photonic crystals are computed by interlacing an inverse method with conjugate gradient projection and full multigrid acceleration. The presently developed method is applied to compute band structures of photonic crystals with modified simple cubic lattice, tetragonal square spiral structure (direct and inverse structures), and diamond structure with sp3 -like configuration. The computed results for the modified simple cubic and square spiral structures are in close agreement with those obtained by previous authors. Moreover, the sp3 -like configuration made of silicon and air is reported to have a large band gap which is larger than the largest reported elsewhere in the literature. KEYWORDS: photonic crystals, large band gaps, diamond structure, multigrid acceleration DOI: 10.1143/JPSJ.73.727. 1.. is compact. For a comparison, the finite difference method (FDM) leads to a sparse matrix and ordinary matrix eigenvalue problem, whilst PWE leads to ordinary matrix problem but the matrix is always full, and FEM leads to sparse matrices, yielding a generalized eigenvalue problem. The latter statement is true of other formulations of finite element methods.15) There are several other methods using different analytic/discrete formulations of Maxwell’s equations, such as transfer matrix method,16,17) multiple scattering method18) and finite-difference time-domain (FDTD) method.19,20) These different methods lead to different complexities to computing photonic band structures. However, finite difference analysis of photonic band structures in three dimensions based on direct discretization of the vector eigenvalue problem (1) has not been established in the literature. In the present study, three important issues are addressed for the finite difference method. (i) First of all, the FDM does not use series-expansion (like PWE) nor uses any basis functions (like FEM). It is therefore essentially important that how the formulation fulfills the transversality condition. Actually, from eq. (1), r H ¼ 0 implies that the left hand side of eq. (1) must be correctly discretized so that the resultant matrix is divergence-free in the discrete sense. (ii) If FDM is correctly formulated, the next issue of importance is the appearance of many zero eigenfrequencies, as already noted in the literature of finite element analysis.21) The number of zero eigenfrequencies was found to be equal to the number of internal nodes of finite elements.22) In the present finite difference formulation, it is shown to be exactly equal to the number of grid points. If the size of the matrix for the eigenvalue problem is large, it is impractical to apply a direct eigenvalue solver (such as QR algorithm23)) or to deflate the many zeros in an iterative method (such as inverse iteration). In the study of photonic band gaps, the first few branches of non-zero eigenfrequency are of primary interest. It is therefore important that the first few non-zero frequency bands can be computed without the necessity of deflating the zero frequencies. This can be achieved by. Introduction. In recent years, photonic crystals have received extensive studies because of their academic and practical interests. The most attractive/distinguished feature of photonic crystals is their full band gap structures.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 important applications such as defect cavities,7) optical waveguides,8) defect-mode photonic crystal lasers,9) and feedback mirror in laser diodes.10) Photonic crystals, also named photonic band-gap materials, are electromagnetic analogy of electronic structures in quantum mechanics. But there are different difficulties with eigenvalue problems of the two analogous fields. It has been a challenging work to compute efficiently band structures of three-dimensional photonic crystals. In computing band structures of three-dimensional photonic crystals, we are faced with a vector (double-curl) eigenvalue problem, say, in H field formulation, 2 1 ! r rH ¼ H: ð1Þ " c It is a nontrivial problem to formulate a numerical method so that the transversality condition r H ¼ 0 is correctly enforced.11) The transversality condition can be easily done in the widely used method of plane wave expansions (PWE) by choosing, for each component hG eiðkþGÞr two orthogonal polarizations hG so that ðk þ GÞ hG ¼ 0.12) In the finite element method (FEM), the so-called vector edge elements with each individual basis function satisfying the transversality condition have been proposed for different shapes of elements.13) Nonfulfillment of the transversality condition would lead to nonphysical spurious solutions.14) On the other hand, finite difference is the most straightforward method to discretize the double-curl equation and the resultant matrix . E-mail: [email protected] E-mail: [email protected]. y. 727.

(2) 728. J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. Fig. 1. Modified simple cubic lattice comprising dielectric spheres and connecting thin circular cylinders.25). employing projections onto the range space of the discrete matrix. (However, it is noted that the problem of many zeros does not occur in the PWE as in each plane wave component the number of polarizations is chosen to be the independent two rather three.) (iii) As the third issue of importance, it is of great interest to develop a fast algorithm for the eigenvalue problem. An inverse method, accelerated by multigrid technique with use of projection is proposed for this purpose. The method exploits the sparsity of the matrix for the eigenvalue problem in the finite difference formulation. Because of the above mentioned difficulties, the present method is a nontrivial extension of a similar method, recently developed by the authors24) for computing photonic band structures in two dimensions. In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.25) Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.26) The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,27,28) as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.29) The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp3 -like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.12,28) In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp3 structure of the electrons of diamond atoms. Recently, submicron diamond-lattice photonic crystals have been successfully produced by two-photon laser nanofabri-. R. L. C HERN et al.. Fig. 2. Tetragonal square spiral structure comprising circular cylinders.27,28). Fig. 3. diamond structure with sp3 -like configuration comprising dielectric spheres and connecting spheroids.. cation (photopolymerization).32) The order of presentation of the paper is organized as follows. In §2, we show how to correctly formulate the finite difference method for the double curl operator of the photonic eigenvalue problem. In §3, we develop the numerical method (inverse iteration with the full multigrid acceleration) and present the fast algorithm, in which two alternative methods of projection are proposed to avoid the necessity of deflating zeros). In §4, we first present numerical results that illustrate the efficiency of the presently developed method. Then, the band structures are computed for the modified simple cubic lattice, the tetragonal square spiral structure (direct and inverse structure) and the diamond structure with sp3 -like configuration. Finally, concluding remarks with a summary of results are drawn in §5. 2.. Basic Equations and Finite Difference Formulation. The electromagnetic waves propagating in the photonic crystals are well described by Maxwell’s equations. For linear isotropic and frequency-independent dielectric materials with permeability close to one, the time-harmonic.

(3) J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. CHERN et al.. modes of electromagnetic waves are described by 2 ! LH ¼ H c. ð2Þ. with the operator L¼r. 1 r ". ð3Þ. along with r H ¼ 0. In the above equations, " ¼ "ðrÞ is the dielectric function, ! is the angular frequency of the timeharmonic modes, and c is the speed of light in vacuum. In component form, eq. (2) is written as 2 32 3 2 3 Lxx Lxy Lxz Hx 2 H x ! 6 7 6L 76 7 ð4Þ 4 yx Lyy Lyz 54 Hy 5 ¼ 4 Hy 5; c Lzx Lzy Lzz Hz Hz where. . Lxx Lyy Lzz Lxy Lyx Lzx. @ 1 @ @ 1 @ þ ; ¼ @y " @y @z " @z @ 1 @ @ 1 @ þ ; ¼ @x " @x @z " @z @ 1 @ @ 1 @ þ ; ¼ @x " @x @y " @y @ 1 @ @ 1 @ ; Lxz ¼ ; ¼ @y " @x @z " @x @ 1 @ @ 1 @ ; Lyz ¼ ; ¼ @x " @y @z " @y @ 1 @ @ 1 @ ; Lzy ¼ ; ¼ @x " @z @y " @z. l¼1. ð5Þ. Sl. where Sl , l ¼ 1 to 6 are the six surfaces enclosing Vcell . Next,. Hz(i,j,k+1/2) Z Hx(i-1/2,j,k) Hy(i,j-1/2,k). Hy(i,j+1/2,k) Y. Hx(i+1/2,j,k) Hz(i,j,k-1/2). we approximate the right hand side of eq. (6) by 6 Z 6 X X 1 ðLHÞ nda r r H nA; 0¼ " Cl l¼1 Sl l¼1. ð7Þ. where Cl denotes the center of Sl and A is the area of Sl . Now we approximate the double-curl operator on the right hand side of eq. (7) at the centers Cl of six surfaces Sl by a secondorder central difference scheme. For example, consider Lxx Hx and Lxy Hy at the center of front surface of Vcell . First, the first-order derivatives at the middle points ði þ 1=2; j þ 1=2; kÞ of the H field are approximated by @Hx 1 ðHx Þiþ1=2; jþ1;k ðHx Þiþ1=2; j;k ; ð8Þ @y iþ1=2; jþ1=2;k h where h is the grid spacing. Then the second-order derivatives are approximated using the above derivatives, @ 1 @Hx @y " @y iþ1=2; j;k ! ð9Þ 1 1 @Hx 1 @Hx : h " @y iþ1=2; jþ1=2;k " @y iþ1=2; j1=2;k Substituting (8) for the derivatives in (9) eliminates the middle fields at ði þ 1=2; j 1=2; kÞ. The complete finite difference formulation for the diagonal term Lxx Hx is. and Hx , Hy and Hz are components of the H field. This section is aimed to discretize the operator L on the left hand side of eq. (2) so that the transversality condition is correctly imposed in discrete sense. For this purpose, we first take the divergence of eq. (2) and integrate over a unit grid cell Vcell , which is shown in Fig. 4. In this figure, Hx , Hy and Hz are interleaved to occupy the sites of front and rear, left and right, top and bottom sites, respectively. The reason of choosing this becomes clear shortly. Applying the divergence theorem, we have Z 6 Z X r ðLH Þd ¼ ðLH Þ nda; ð6Þ Vcell. 729. X. Fig. 4. The unit grid cell for finite difference formulation. Hx , Hy and Hz are interleaved to occupy the front and rear, left and right, top and bottom sites, respectively.. ðLxx Hx Þiþ1=2; j;k iþ1=2; j;k 1 Hx Hxiþ1=2; j;k 2 þ h "iþ1=2; j1=2;k "iþ1=2; jþ1=2;k iþ1=2; j;k 1 Hx Hxiþ1=2; j;k þ þ 2 h "iþ1=2; j;k1=2 "iþ1=2; j;kþ1=2 1 Hxiþ1=2; j1;k H iþ1=2; jþ1;k þ x 2 h "iþ1=2; j1=2;k "iþ1=2; jþ1=2;k 1 Hxiþ1=2; j;k1 Hxiþ1=2; j;kþ1 þ 2 h "iþ1=2; j;k1=2 "iþ1=2; j;kþ1=2 and the off-diagonal term Lxy Hy is 1 Lxy Hy iþ1=2; j;k 2 h 1 2 h. Hyi; j1=2;k "iþ1=2; j1=2;k. . Hyi; jþ1=2;k "iþ1=2; jþ1=2;k. Hyiþ1;j1=2;k. ð10Þ !. "iþ1=2; j1=2;k . Hyiþ1;jþ1=2;k "iþ1=2; jþ1=2;k. ! : ð11Þ. Likewise, we discretize the other terms in eq. (5) in a similar way, and collect the results into eq. (7). Altogether, there are 96 terms in these approximations. The 96 terms cancel each other exactly, showing that the transversality condition is satisfied in a discrete form of eq. (6). The discrete transversality condition is satisfied even when the grid cell cross the interface of different dielectric materials. The detailed proof of the discrete transversality condition is given in Appendix A. Indeed, our numerical experiments further confirm this condition by giving 64 zero eigenvalues for matrix of size 192 192 (grid size 4 4 4), and 512 zeros for matrix of size 1536 1536 (grid size 8 8 8). It is now cleared that the grid arrangement in Fig. 4 is sufficient to fulfill the transversality condition in the present finite difference formulation. Namely, we use Hx only at half-integer points i þ 1=2 and integer points j and k, Hy at half-integer points j þ 1=2 and integer points i and k, and Hz.

(4) 730. J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. C HERN et al.. at half-integer points k þ 1=2 and integer points i and j. For simplicity, we denote ði þ 1=2; j; kÞ by a=2, ði; j þ 1=2; kÞ by b=2 and ði; j; k þ 1=2Þ by c=2. To be consistent with the above finite difference formulation, the double-curl operator in eq. (3) is discretized to yield a discretization matrix A as follows: 2 3 2 a=2 3 Hx Lxx Lxy Lxz 6 L 7 6 7 ½AHi; j;k 4 yx Lyy Lyz 5 4 Hyb=2 5; ð12Þ Lzx Lzy Lzz C.D. Hzc=2 where ‘‘C.D.’’ denotes that the central difference is applied. The matrix A has 39 finite difference terms. Appendix B contains all the details. It is also noticed that the space arrangement for the H field is equivalent to that in Yee’s cell,33) but the formulation, which contains second-order derivatives, is different from Yee’s scheme. In the finite element method, the vector edge elements have been proposed to satisfy the transversality condition.13) In their formulations, the transversality condition is satisfied on each individual basis function while in the present method, the transversality condition is satisfied in a less obvious way. If the matrix size is small, the eigenvalues can be obtained by direct methods. However, the operation count is of order N 3 , which becomes prohibitively large for large N. The difficulty of large operation counts can be alleviated by an iterative solver, but the large number of zero eigenvalues causes another difficulty. Zero modes appear before any nontrivial eigenmodes when we solve the eigenfrequencies from the smallest one. Therefore, it is impractical and inaccurate to deflate a large number of zero eigenmodes. In the present study, deflation of zero modes is avoided by introducing a projection operator onto the range space of the matrix A, which will be discussed in the next section. Next, it is very important to choose the domain of computation that is advantageous to the finite difference formulation. For the three types of configuration, the domains of computation are chosen to be; (i) a cubic coincident with the primitive cell for the modified simple cubic lattice, as shown in Fig. 5, (ii) a tetragon coincident with the primitive cell for the tetragonal square spiral. s. a3. c. a1. a2. a. a. Fig. 6. Domain of computation for the tetragonal square spiral structure in Fig. 2 is a tetragon with square of side length a and height c. The square spiral structure is composed of circular cylinder with radius r, length L and pitch c.. r a3 b. √ 23 a′ a2. a1. √3 2. a′. 3 Fig. 7. Domain of computation for the diamond structure with pffiffisp ffi -like 0 0 3 =2 and configuration in Fig. 3 is a tetragon with length a , width a pffiffiffiffiffiffiffiffi pffiffiffi height a0 2=3, where a0 ¼ a= 2 with a the lattice constant. The radius of the dielectric sphere is r. The connecting spheroid has minor axis length b, and the foci located at the centers of the spheres.. structure, as shown in Fig. 6, and (iii) a tetragon with one edge aligned with one of the lattice translation vectors for the diamond structure with sp3 -like configuration. This domain has and must have the same volume of the primitive cell as shown in Fig. 7. Finally, Bloch’s theorem is applied at the boundary of domain of computation:. r. a3. H k ðr þ ai Þ ¼ eikai H k ðrÞ;. a. a1. a2. a. a′. a. Fig. 5. Domain of computation for the modified simple cubic lattice in Fig. 1 is a cube with side length a. The radius of the dielectric sphere is r and the radius of the connecting cylinder is s.. ð13Þ. where Hk is the Bloch function for magnetic field associated with the wave vector k in the first Brillouin zone. The letters ai (i ¼ 1; 2; 3) denote the lattice translation vectors; (i) for the modified simple cubic lattice, a1 ¼ að1; 0; 0Þ, a2 ¼ að0; 1; 0Þ and a3 ¼ að0; 0; 1Þ, (ii) for the tetragonal square spiral structure, a1 ¼ að1; 0; 0Þ, a2 ¼ að0; 1; 0Þ and a3 ¼ cð0; 0; 1Þ, andp(iii) a0 ð1; pffiffiffiffiffiffiffi ffi 0; 0Þ, ffiffiffi for the diamond0 structure,paffiffi1ffi ¼ 0 a2 ¼ a ð1=2; 3=2; 0Þ and a3 ¼ a ð1=2; 1=2 3; 2=3Þ with p ffiffi ffi a0 ¼ a= 2. In the last case, since a2 and a3 have the component a0 =2 in the x-direction, application of Bloch’s theorem in the y- and z-directions should be additionally.

(5) J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. CHERN et al.. InverseEigen { for n ¼ 1 to S Initial guess b do Solve ðA IÞx ¼ b by LUD or PCG Deflate x by q1 to qn1 Set b ¼ x=kxk Rayleigh Quotient n ¼ hb; Abi until kðA n IÞbk2 < Set qn ¼ b end }. 0. 20. 40. 60. 80. 100. 120. where is the parameter for inverse iteration, 1 2 S is the sequence of S smallest eigenvalues, and q1 to qS are the corresponding eigenvectors. The inner product h; i is defined as X xi yi ð16Þ hx; yi ¼. 140. 160. 180 0. 20. 40. 60. 80. 100 nz = 2496. 120. 140. 160. 180. i. Fig. 8. 13-band structure of the discretization matrix A of size 192 192.. accompanied by the shift of one-half cell p inffiffiffithe x-direction. pffiffiffi Furthermore, a3 has the component a0 =2 3 ¼ ða0 3=2Þ=3 in the y-direction, application of Bloch’s condition in the zdirection should be further accompanied by the shift of onethird cell in the y-direction. 3.. 731. Inverse Method With Multigrid Acceleration and Methods of Projection. In eq. (12) we obtain the discretization matrix A of the vector differential operator L. The matrix structure is shown in Fig. 8. Matrix A has a sparsity pattern very different from that for the two-dimensional operators (either in transverseelectric or magnetic modes). For the two-dimensional case, there are only 5 nonzero entries in each row of A, but there are 13 nonzero entries in each row for the three-dimensional case. For example, the number of nonzero entries of matrix A of size 192 192 (grid size 4 4 4) is 192 13 ¼ 2496. 3.1 Inverse iteration From a practical point of view, the first few branches of eigenvalues are of primary interest. A natural choice for this purpose is the method of inverse iteration, applied to obtain the smallest eigenvalues as well as the corresponding eigenvectors. For the eigenvalue problem Ax ¼ x;. ð14Þ. the method of inverse iteration finds the smallest eigenvalue, which is closest to , by iteratively solving ðA IÞx ¼ b:. ð15Þ. The second eigenvalue is obtained by repeating the same method with deflation of the first eigenvector. The third eigenvalue is obtained by repeating the same method with deflating the first and second eigenvectors, and so on. In particular, the Gram–Schmidt orthogonalization is used to deflate, or to eliminate the components of previously obtained eigenvectors. The inverse algorithm23,34) is presented as follows:. in view of the orthonormalization property of the H field: Z 1 H ðrÞ H n ðrÞdr ¼ mn ; ð17Þ Vcell Vcell m where H n is the nth eigenmode of the magnetic field, Vcell is the volume of the primitive cell and mn denotes the Kronecker delta. LUD denotes the lower–upper decomposition solver for matrix inversion, and PCG the preconditioned conjugate gradient solver. 3.2 Methods of projection In §2, we mentioned that matrix A has N=3 zero modes, where N ¼ 3n3 for the grid size n n n. Since only the first few nontrivial eigenmodes are of practical interests, it is important to remove from the iterated solutions the components of the null (zero-mode) space NðAÞ ¼ fxjAx ¼ 0g. But it would be impractical and inaccurate to delate all the N=3 zero modes. Instead, we propose two methods of projection onto the range space RðAÞ ¼ fAxg that dispenses with the necessity of deflating the many zero modes. 3.2.1 Projection by multiplying A In the original inverse algorithm, we iteratively solve x¼. 1 b ðA IÞ. ð18Þ. to obtain the eigenvector x. Let z be any solution corresponding to zero frequency, that is, Az ¼ 0. Taking inner product with x, and using the Hermitian property of A, we obtain hx; Azi ¼ hAx; zi ¼ 0. The result indicates that Ax is orthogonal to z. We can therefore multiply the solution x by A to eliminate the components of zero modes z. This idea is implemented by modifying the inverse algorithm as x¼. 1 Ab: ðA IÞ2. ð19Þ. This is equivalent to applying a forward iteration without shift followed by two inverse iterations with shift . The overall effect is one inverse iteration, as is the original inverse algorithm. The new algorithm, however, removes the components of zero modes at each iteration, keeping the iterative solutions falling into the range space of A..

(6) 732. J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. C HERN et al.. 3.2.2 Projection by conjugate gradient method The above method exhibits a somewhat oscillatory behavior, though the projection is successful and stable. A better idea of projecting the solution x onto the range space of RðAÞ is to find w which minimizes the object function f ¼. 1 kAw xk2 : 2. ð20Þ. The solution which does not contain zero-mode components is then obtained by replacing x by Aw. The minimization of f is performed by equating the derivatives of f with respect to the real and imaginary parts of w to zero, respectively. This is equivalent to solving the normal equation ðA AÞw ¼ A x:. eration is to embed the multigrid solver in the inverse iteration loop, namely, multigrid takes place of the expensive LUD or less expensive PCG. A smart and novel idea for acceleration is to embed inverse iteration in the multigrid acceleration structure. In this idea, the Rayleigh quotient is updated at each level of multigrid acceleration, which results in a more efficient algorithm. Let there be L levels of grids, with m ¼ 1 the coarsest level and m ¼ L the finest level. On each grid level m, the differential operator is discretized to form the matrix AðmÞ . The inverse method with the full multigrid acceleration implemented with the good idea and conjugate gradient projection is designed as follows:. ð21Þ. FMGInverse { Call InverseEigen at m ¼ 1 for m ¼ 2 to L for n ¼ 1 to S Interpolate qnðm1Þ to bðmÞ do Solve ðAðmÞ IÞxðmÞ ¼ bðmÞ by MGV Project xðmÞ onto Range(AðmÞ ) by CG ðmÞ Deflate xðmÞ by qðmÞ 1 to qn1 ðmÞ Set b ¼ xðmÞ =kxðmÞ k Rayleigh Quot. nðmÞ ¼ hbðmÞ ; AðmÞ bðmÞ i until kðAðmÞ nðmÞ IÞbðmÞ k2 < ðmÞ Set qðmÞ n ¼ b end end }. In practice, it is not necessary to solve w explicitly and multiply it by A, but to solve A y ¼ A x;. ð22Þ. keeping in mind that y ¼ Aw falls into the range space of A. It is noticed that A is singular and the above equation may not be solvable. However, if y is kept falling into the range space RðAÞ, a unique solution can be obtained. This is achieved by applying the conjugate gradient method. For completeness, the conjugate gradient algorithm35,36) for solving Ax ¼ b is given as follows: CG { n ¼ 0; x0 ¼ 0; r0 ¼ b; p0 ¼ r0 do n¼nþ1 ¼ rn1 rn1 =pn1 Apn1 xn ¼ xn1 þ pn1 rn ¼ rn1 Apn1 ¼ rn rn =rn1 rn1 pn ¼ rn1 þ pn1 until krn k2 < } If A lies in the range space of A, then by construction each rn (pn ) also belongs to the range space of A. Starting from x0 ¼ 0, the conjugate gradient algorithm uses consecutive expansion of orthogonal spaces (Krylov spaces) to ensure that each iterative solution xn would fall into the range space of A. Projection is achieved by applying the algorithm to eq. (22) with b replaced by A x. Numerical results indicates that this algorithm is a more stable and efficient method for projecting the iterative solutions. 3.3 Multigrid acceleration For the problem of arbitrary shape and/or with large dielectric contrast, high resolution grid is necessary. Therefore, the dimension of A is usually very large. The matrix inverse problem may be solved in two ways: direct method and iterative method. For direct solvers such as LU decomposition or QR factorization, the operation count is of order N 3 , which is computationally expensive. In this study, we choose the multigrid method to accelerate convergence of matrix inversion in the algorithm. This extends our previous algorithm for two-dimensional problems.24) Multigrid methods use multilevels of grids to attenuate the errors of different frequencies.36,37) A plain idea for accel-. where 1ðmÞ 2ðmÞ SðmÞ. ð23Þ. is the sequence of S smallest eigenvalues on the grid level m, ðmÞ and qðmÞ are the corresponding eigenvectors. 1 qS ðmÞ RangeðA Þ denotes the range space of AðmÞ and MGV is the multigrid V-cycle solver.36,37) In particular, SSOR method with Chebyshev acceleration38) is used for relaxing the errors in the smooth step of the multigrid V-cycle. 4.. Results and Discussion. As an example of test, we compute the first six frequency bands at the point Xð; 0; 0Þ in the first Brillouin zone for a simple cubic lattice of dielectric sphere with radius r=a ¼ 0:345 for dielectric contrast "="0 ¼ 13. Table I shows a comparison of computing times in central processing unit (CPU) seconds between the present method (FMGI), conjugate gradient method and plane wave expansion method.. Table I. Computing times in CPU seconds for the first six frequency bands at the point X for a simple cubic lattice of dielectric sphere with radius r=a ¼ 0:345 for dielectric contrast "="0 ¼ 13. Ngrid =Nwave. FMGI. PCG. 888. 3.7. 3.0. PWE 5.1. 10 10 10. —. —. 32.3. 12 12 12. —. —. 155.8. 14 14 14. —. —. 628.6. 16 16 16. 25.3. 32.4. —. 32 32 32. 280.5. 636.1. —. 64 64 64. 3416.5. —. —.

(7) J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. CHERN et al.. 733. 0.6. 104. Gap-Midgap Ratio = 0.1363. Frequency (ωa/2πc). Computing Time (CPU sec). 0.5 103. 102. 10. 0.3. 0.2. R. 1. X. FMGI CG PWE 10. 0.4. 0.1. 0. 102. 103. 104. 105. M. 0. 106. Number of Grid Points or Waves. The plot of the computing times is shown in Fig. 9. All the computations are performed on a Pentium 4-2.8 GHz PC. The computing time for the present method is of order N 1:18 , which is slightly larger than order N for the two-dimensional problems.24) This is because additional computations are due to the conjugate gradient projection, which is in general slower than the pure multigrid method. If the multigrid solver is replaced by the conjugate gradient method, the computing time is of order N 1:43 . On the other hand, the direct plane wave expansion has the computational cost of about the order N 2:87 . However, it must be noticed that iterative solvers with suitable preconditioners may significantly reduce the operation counts for PWE.39) The accuracy of the present method is illustrated by the convergence test against the grid size for the band edges and gap–midgap ratios for the optimized diamond structure with sp3 -like configuration. Computations are performed on three different grids: N ¼ 24 24 24, 48 48 48 and 96 96 96. Table II lists the numerics of the location of the band edge for the diamond structure with sp3 -like configuration in Fig. 2. Basically, the computed results agree with each other to three significant digits. In other words, the relative differences of the computed results with the different grids are less than 0.1%. Figure 10 shows the band structure for the modified simple cubic lattice in Fig. 1 at r=a ¼ 0:345 and s=a ¼ 0:11 for "="0 ¼ 13. The complete band gap lies between the 5th and 6th branches with the lower edge located at point X and the upper edge at point M. The gap–midgap ratio is 0:1363 with center frequency 0:4501. The modified simple cubic lattice was previously investigated by Biswas et al., who obtained a gap–midgap ratio 0:12 for "="0 ¼ 12:96.25) At the same geometric parameters, the overall band structure of the. Table II. Convergence test against the grid size for the band edges and gap–midgap ratios of the diamond structure with sp3 -like configuration in Fig. 3 for r=a ¼ 0:12 and b=a ¼ 0:11. Ngrid. 24 24 24. 48 48 48. 96 96 96. !up. 0.6890. 0.6879. 0.6874. !low. 0.5016. 0.5013. 0.5010. !mid. 0.5953. 0.5946. 0.5942. !. 0.1874. 0.1866. 0.1864. !=!mid. 0.3147. 0.3139. 0.3137. X. M. Γ. R. Fig. 10. Band structure computed with 48 48 48 grid for the modified simple cubic lattice in Fig. 1 at r=a ¼ 0:345 and s=a ¼ 0:11 for "="0 ¼ 13.. 0.5. 0.4. Frequency (ωa/2πc). Fig. 9. Comparison of the present method, conjugate gradient method and plane wave expansions in computing time.. Γ. Gap-Midgap Ratio = 0.162. 0.3. 0.2 R. Z A. 0.1. X. M. 0. Γ. Z. A. M. Γ. X. R. Z. Fig. 11. Band structure computed with 48 48 48 grid for the tetragonal square spiral structure in Fig. 2 at r=a ¼ 0:13, c=a ¼ 1:3 and L=a ¼ 1:65 for "="0 ¼ 13.. present study is in reasonably good agreement with their results. However, our gap–midgap ratio is somewhat larger. One possible reason for this could be that there are relatively smaller number of plane waves used in ref. 25, compared to the presently larger number of grid points N ¼ 48 48 48. Figure 11 shows the band structure for the tetragonal square spiral structure in Fig. 2 at r=a ¼ 0:13, c=a ¼ 1:3 and L=a ¼ 1:65 for "="0 ¼ 13. A photonic band gap with gap– midgap ratio 0:162 opens between the 4th and 5th bands. The upper edge closes at point R and the lower edge at point A. The center frequency of the band gap is 0:3974. Figure 12 shows the band structure for the inverse structure of the tetragonal square spiral structure with r=a ¼ 0:333, c=a ¼ 1:69 and L=a ¼ 1:5 for "="0 ¼ 13. The inverse structure has a similar band structure with the band edges occurring at the same symmetric points: R (upper) and A (lower), but with an even larger band gap. The gap–midgap ratio is 0:2593 with center frequency 0:3826. These configurations were studied by Toader and John, who obtained the gap–midgap ratios 0:162 and 0:236 for direct and inverse structures, respectively, for "="0 ¼ 11:9.28) For the direct structure, their gap–midgap ratio is identical to that obtained in the present study, while for the inverse structure, our result is a little bit larger. Figure 13 shows the band structure for the diamond structure with sp3 -like configuration in Fig. 3 at r=a ¼ 0:12.

(8) 734. J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. R. L. C HERN et al.. 0.5 100. Gap-Midgap Ratio = 0.2593. 0.3. Residue. Frequency (ωa/2πc). 0.4. 10. -1. 10. -2. 10. -3. 10. -4. 0.2 R. Z A. 0.1. X. M. 0. Γ. Z. A. Γ. M. X. Z. R. Fig. 12. Band structure computed with 48 48 48 grid for the inverse structure of tetragonal square spiral at r=a ¼ 0:333, c=a ¼ 1:69 and L=a ¼ 1:5 for "="0 ¼ 13.. No Precond Jacobi SOR SSOR. 100. 101. 102. Iteration Number. Fig. 14. The effect of preconditioners on the convergence of the conjugate gradient method as an eigenvalue solver of two-dimensional photonic crystals.. 0.8. 100 0.7. Gap-Midgap Ratio = 0.3139 10. -1. 10. -2. 10. -3. 10. -4. 0.5. Residue. Frequency (ωa/2πc). 0.6. 0.4 0.3. L U. 0.2. X. W K. 0.1. X. U. L. Γ. X. W. K. Fig. 13. Band structure computed with 48 48 48 grid for the diamond structure with sp3 -like configuration in Fig. 3 at r=a ¼ 0:12 and b=a ¼ 0:11 for "="0 ¼ 13.. and b=a ¼ 0:11 for "="0 ¼ 13. A large photonic band gap occurs between the 2nd and 3rd branches with gap–midgap ratio 0:3139 and center frequency 0:5946. This band gap has the upper edge at point L and the lower edge at point W. The lower edge from W to K is a flat curve, which almost has the edge value at point U. It is rather interesting to notice that the gap is found to be larger than ever reported in the literature for any other three-dimensional photonic crystals. Finally, we illustrate the difficulty in accelerating convergence of numerical methods for computing three-dimensional photonic band structures. For this, we compare convergence of the conjugate gradient solver for computing two- and three-dimensional photonic band structures as well as three-dimensional Laplace equation. The conjugate gradient method is preconditioned with Jacobi, SOR or SSOR preconditioners.38,40) For two-dimensional problem, these preconditioners do help accelerate convergence of the conjugate gradient method for the eigenvalue solver. Figure 14 shows the effect of these preconditioners on the convergence of the eigenvalue solver for the first band at point X for a square lattice. The relaxation parameter has been optimized for best convergence in each case. Figure 15 shows the effect of these preconditioners on the convergence of the three-dimensional Laplace equation. Likewise, the relaxation parameters have been optimized for best convergence. The results show that the SSOR preconditioner accelerates convergence while the other three do not. 5. 10. 15. 20. Iteration Number. Fig. 15. The effect of preconditioners on the convergence of the conjugate gradient method for the three-dimensional Laplace equation.. 100. Residue. 0. No Precond Jacobi SOR SSOR. 10. -1. 10. -2. 10. -3. 10. -4. 100. No Precond Jacobi SOR SSOR. 101. 102. Iteration Number. Fig. 16. The effect of preconditioners on the convergence of the conjugate gradient method as an eigenvalue solver of three-dimensional photonic crystals.. accelerate convergence. In spite of the two successful examples, none of the well-known preconditioners are effective in accelerating convergence of the conjugate gradient method as an eigenvalue solver for three-dimensional photonic crystals. Figure 16 shows the effect of these preconditioners on the first band at point X for the modified simple cubic lattice. The results show that these preconditioners do not accelerate and even slow down the convergence. There might be some other preconditioners which are suitable for the three-dimensional photonic crystals, but.

(9) J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. they would be very much problem-dependent. On the other hand, the presently developed method are applicable to compute band structures of three-dimensional photonic crystals of arbitrary configuration with strong dielectric interfaces. 5.. Conclusions. In the present study, an efficient method in finite difference formulation was developed for computing band structures of three-dimensional photonic crystals. The method is a substantial extension of our two-dimensional version. (i) First, we have shown how the transversality condition can be satisfied for the vector eigenvalue problem in the finite difference formulation. (ii) Second, fulfillment of the transversality condition results in a large number of unwanted zero-frequency modes. The number of zero frequencies is exactly equal to the number of the internal nodes. The method of inverse iteration has been applied to obtain the first few nontrivial eigenmodes. Projection by the matrix itself, or conjugate gradient method has been employed to avoid the necessity of deflating a large number of zero eigenmodes. (iii) Third, fast computation of photonic band structures is achieved by embedding inverse iteration in the full multigrid algorithm with use of the method of projection. The present method was shown to have a computational cost of order N 1:18 , which is slightly larger than the linear order N for the two-dimensional problems.24) This is due to the more expensive part of the conjugate gradient projection. Nevertheless, the method is much faster than the most commonly used method of plane wave expansions, which has been demonstrated to be of order N 2:87 . The accuracy of the present method has been tested against various grid sizes from N ¼ 24 24 24 to 96 96 96 96, and basically the results agree to each other to three significant digits. Examples of applications include (i) the modified simple cubic lattice, (ii) the tetragonal square spiral structure (direct and inverse structure), and (iii) the diamond structure with sp3 -like configuration. The first two examples were previously studied by other authors.25,27,28) Our numerical results are in close agreement with those obtained by previous authors. However, with the fast algorithm the band structures presented in this study can be easily computed with the grid of size N ¼ 48 48 48, compared to the relatively small N ¼ 140028) used in plane wave expansions. Both the second and third examples may be considered as modifications of the diamond structure. The second example is amenable to the successful fabrication technique: GLAD (GLancing Angle Deposition). The third example with sp3 -like configuration made of silicon and air proposed by the present authors has a large band gap with gap–midgap ratio 0:3139, which is larger than ever reported in the literature. The latter result is of particular interest in view of recently successful production of diamond-lattice photonic crystals by two-photon laser nanofabrication (photopolymerization).32) In general, our numerical experience indicates that the diamond-like structures and their modifications are the ideal three-dimensional configurations to produce large photonic band gaps. As a matter of fact, if square spirals are replaced by triangular or hexagonal spirals, we do not find any significant band gaps for these. R. L. CHERN et al.. 735. configurations. These examples have demonstrated that the presently developed method is capable of resolving the multieigenvalue band structures of three-dimensional photonic crystals of arbitrary configuration with strong dielectric interfaces. The method is currently being extended to be applicable to metallic and metallodielectric photonic crystals; the results will be reported elsewhere. Acknowledgments This work was supported in part by the National Science Council of the Republic of China under Contract No. NSC 91-2212-E-002-072. Appendix A:. Proof of the Discrete Transversality Condition. Recall that we denote ði þ 1=2; j; kÞ by a=2, ði; j þ 1=2; kÞ by b=2 and ði; j; k þ 1=2Þ by c=2. The complete discretized finite difference formulation of eq. (7) can be written as 6 X 1 r r H nA " Cl l¼1 b=2 c=2 a=2 Gb=2 Gc=2 ; ¼ Ga=2 x þ Gy þ Gz Gx y z. Ga=2 x . . Hxa=2. . H a=2 Hxa=2 H a=2 þ x þ þ x "a=2;b=2 "a=2;c=2 "a=2;c=2. "a=2;b=2 a=2;b Hx H a=2;b H a=2;c H a=2;c þ x þ x þ x "a=2;b=2 "a=2;b=2 "a=2;c=2 "a=2;c=2 ! Hyb=2 Hya;b=2 Hyb=2 Hya;b=2 þ þ "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2 c=2 Hz Hza;c=2 Hzc=2 Hza;c=2 þ ; þ "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2. Gb=2 y Hyb=2. . "a=2;b=2. þ. Hya;b=2. . "a=2;b=2. Hyb=2 "a=2;b=2 þ. þ. Hya;b=2 "a=2;b=2. Hyb=2 "b=2;c=2 þ. þ. Hyb=2;c "b=2;c=2. !. Hyb=2 "b=2;c=2 þ. Hyb=2;c. !. "b=2;c=2. . Hxa=2 H a=2;b H a=2 H a=2;b x x þ x "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2 c=2 Hz Hzb;c=2 Hzc=2 Hzb;c=2 þ þ ; "b=2;c=2 "b=2;c=2 "b=2;c=2 "b=2;c=2 þ. Gc=2 z . . Hzc=2. þ. Hzc=2 Hzc=2 H c=2 þ þ z "a=2;c=2 "b=2;c=2 "b=2;c=2. . "a=2;c=2 a;c=2 Hz Hza;c=2 Hzb;c=2 Hzb;c=2 þ þ þ "a=2;c=2 "a=2;c=2 "b=2;c=2 "b=2;c=2 a=2 Hx Hxa=2;c Hxa=2 Hxa=2;c þ þ "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2 ! Hyb=2 Hyb=2;c Hyb=2 Hyb=2;c þ þ ; "b=2;c=2 "b=2;c=2 "b=2;c=2 "b=2;c=2.

(10) 736. J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004. Gxa=2 . . Hxa=2. þ. R. L. C HERN et al.. . Hxa=2 Hxa=2 H a=2 þ þ x "a=2;b=2 "a=2;c=2 "a=2;c=2. "a=2;b=2 a=2;b Hx Hxa=2;b Hxa=2;c Hxa=2;c þ þ þ "a=2;b=2 "a=2;b=2 "a=2;c=2 "a=2;c=2 ! Hya;b=2 Hyb=2 Hya;b=2 Hyb=2 þ þ "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2 a;c=2 Hz Hzc=2 H a;c=2 Hzc=2 z þ þ ; "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2. Gyb=2 Hyb=2. . "a=2;b=2. þ. Hya;b=2. . "a=2;b=2. Hyb=2 "a=2;b=2 þ. þ. Hya;b=2 "a=2;b=2. Hyb=2 "b=2;c=2 þ. þ. Hyb=2;c "b=2;c=2. Lyy Hy b=2. "b=2;c=2 Hyb=2;c. . Hzc=2. þ. !. Summing these terms together shows that all the 96 elements exactly cancel out to zero. The discrete transversality condition of the present formulation is thus proved. Appendix B:. . Discretized Equations for the Eigenvalue Problem. The complete finite difference formulation of the discretization matrix A in eq. (12) is presented as follows: ðLxx Hx Þa=2 1 Hxa=2 H a=2 Hxa=2 H a=2 2 þ x þ þ x h "a=2;b=2 "a=2;b=2 "a=2;c=2 "a=2;c=2 1 Hxa=2;b Hxa=2;b Hxa=2;c Hxa=2;c 2 þ þ þ ; h "a=2;b=2 "a=2;b=2 "a=2;c=2 "a=2;c=2 Lxy Hy a=2 ! Hyb=2 Hya;b=2 Hyb=2 Hya;b=2 1 2 þ ; h "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2 Lxz Hz a=2 c=2 1 Hz Hza;c=2 Hzc=2 Hza;c=2 2 þ ; h "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2. Lzy Hy. "a=2;b=2. þ. þ. Hya;b=2 "a=2;b=2. "b=2;c=2 þ. þ. Hyb=2;c "b=2;c=2. Hyb=2. !. "b=2;c=2 þ. Hyb=2;c. ! ;. "b=2;c=2. c=2. 1 2 h. . "a=2;c=2 a;c=2 Hz Hza;c=2 Hzb;c=2 Hzb;c=2 þ þ þ "a=2;c=2 "a=2;c=2 "b=2;c=2 "b=2;c=2 a=2;c Hx Hxa=2 Hxa=2;c Hxa=2 þ þ "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2 ! Hyb=2;c Hyb=2 Hyb=2;c Hyb=2 þ þ : "b=2;c=2 "b=2;c=2 "b=2;c=2 "b=2;c=2. Hya;b=2. "a=2;b=2. Hyb=2. Lzx Hx c=2 a=2 1 Hx Hxa=2;c Hxa=2 Hxa=2;c 2 þ ; h "a=2;c=2 "a=2;c=2 "a=2;c=2 "a=2;c=2. "b=2;c=2. Hzc=2 Hzc=2 H c=2 þ þ z "a=2;c=2 "b=2;c=2 "b=2;c=2. þ. Hyb=2. Lyz Hz b=2 c=2 1 Hz Hzb;c=2 Hzc=2 Hzb;c=2 2 þ ; h "b=2;c=2 "b=2;c=2 "b=2;c=2 "b=2;c=2. Hxa=2;b Hxa=2 Hxa=2;b Hxa=2 þ þ "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2 b;c=2 Hz Hzc=2 Hzb;c=2 Hzc=2 þ þ ; "b=2;c=2 "b=2;c=2 "b=2;c=2 "b=2;c=2. . "a=2;b=2. 1 2 h. . Gzc=2. Hyb=2. 1 2 h. !. Hyb=2. þ. Lyx Hx b=2 a=2 1 Hx Hxa=2;b Hxa=2 Hxa=2;b 2 þ ; h "a=2;b=2 "a=2;b=2 "a=2;b=2 "a=2;b=2. Hyb=2 "b=2;c=2. . Hyb=2;c "b=2;c=2. . Hyb=2 "b=2;c=2. þ. Hyb=2;c "b=2;c=2. ! ;. Lzz Hz c=2 1 Hzc=2 H c=2 Hzc=2 H c=2 2 þ z þ þ z h "a=2;c=2 "a=2;c=2 "b=2;c=2 "b=2;c=2 1 Hza;c=2 Hza;c=2 Hzb;c=2 Hzb;c=2 2 þ þ þ : h "a=2;c=2 "a=2;c=2 "b=2;c=2 "b=2;c=2. 1) 2) 3) 4) 5) 6) 7). 8) 9) 10) 11). 12) 13). 14). 15) 16). E. Yablonovitch: Phys. Rev. Lett. 58 (1987) 2059. S. John: Phys. Rev. Lett. 58 (1987) 2486. N. Garcia and A. Z. Genack: Phys. Rev. Lett. 66 (1991) 1850. S. John: Phys. Today 44 (1991) No. 5, 32. S. Lin and G. Arjavalingam: J. Opt. Soc. Am. B 11 (1994) 2124. T. F. Krauss, R. M. De La Rue and S. Brand: Nature 383 (1996) 699. J. S. Foresi, P. R. Villeneuve, J. Ferrera, E. R. Thoen, G. Steinmeyer, S. Fan, J. D. Joannopoulos, L. C. Kimerling, H. I. Smith and E. P. Ippen: Nature 390 (1997) 143. A. Mekis, J. C. Chen, I. Kurland, S. Fan, P. R. Villeneuve and J. D. Joannopoulos: Phys. Rev. Lett. 77 (1996) 3787. O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus and I. Kim: Science 284 (1999) 1819. D. L. Bullock, C. Shih and R. S. Margulies: J. Opt. Soc. Am. B 10 (1993) 399. S. G. Johnson and J. D. Joannopoulos: The MIT Photonic Crystals Tutorial home page http://ab-initio.mit.edu/photons/ tutorial/ (2003). K. M. Ho, C. T. Chan and C. M. Soukoulis: Phys. Rev. Lett. 65 (1990) 3152. A. Bossavit and J. C. Verite: IEEE Trans. Magn. 18 (1982) 431; M. Hano: IEEE Trans. Microwave Theory Tech. 32 (1984) 1275; G. Mur and A.T. de Hoop: IEEE Trans. Magn. 21 (1985) 2188. C. W. Crowley, P. P. Silvester and H. Hurwitz Jr.: IEEE Trans. Magn. 24 (1988) 397; K. Ise, K. Inoue and M. Koshiba: IEEE Trans. Microwave Theory Tech. 38 (1990) 1352; K. D. Paulsen and D. R. Lynch: ibid. 39 (1991) 395. D. C. Dobson, J. Gopalakrishnan and J. E. Pasciak: J. Comput. Phys. 161 (2000) 668. J. B. Pendry and A. MacKinnon: Phys. Rev. Lett. 69 (1992) 2772..

(11) J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 17) Y. C. Tsai, J. B. Pendry and K. W.-K. Shung: Phys. Rev. B 59 (1999) R10401. 18) X. Wang, X.-G. Zhang, Q. Yu and B. N. Harmon: Phys. Rev. B 47 (1993) 4161. 19) C. T. Chan, Q. L. Yu and K. M. Ho: Phys. Rev. B 51 (1995) 16635. 20) A. J. Ward and J. B. Pendry: Phys. Rev. B 58 (1998) 7252. 21) J. M. Jin: The Finite Element Method in Electromagnetics (John Wiley & Sons, New York, 2002) 2nd ed. 22) A. Chatterjee, J. M. Jin, J. L. Volakis: IEEE Trans. Microwave Theory Tech. 40 (1992) 2106. 23) B. N. Parlett: The Symmetric Eigenvalue Problem (SIAM, Philadelphia, 1998). 24) R. L. Chern, C. C. Chang, C. C. Chang and R. R. Hwang: Phys. Rev. E 68 (2003) 26704. 25) R. Biswas, M. M. Sigalas, K. M. Ho and S. Y. Lin: Phys. Rev. B 65 (2002) 205121. 26) S. Lin, J. G. Fleming, R. Lin, M. M. Sigalas, R. Biswas and K. M. Ho: J. Opt. Soc. Am. B 18 (2001) 32. 27) O. Toader and S. John: Science 292 (2001) 1133. 28) O. Toader and S. John: Phys. Rev. E 66 (2002) 16610.. R. L. CHERN et al. 29) 30) 31) 32) 33) 34) 35) 36) 37) 38) 39) 40). 737. A. Chutinan and S. Noda: Phys. Rev. B 57 (1998) 2006. K. Robbie, M. J. Brett and A. Lakhtakia: Nature 384 (1996) 616. S. Kennedy, M. Brett, O. Toader and S. John: Nano Lett. 2 (2002) 59. K. Kaneko, H. B. Sun, X. M. Duan and S. Kawata: Appl. Phys. Lett. 83 (2003) 2091. K. S. Yee: IEEE Trans. Antennas Propagat. 14 (1966) 302. G. H. Golub and C. F. Van Loan: Matrix Computations (Johns Hopkins University Press, London, 1996) 3rd ed. L. N. Trefethen and D. Bau: Numerical Linear Algebra (SIAM, Philadelphia, 1997). J. W. Demmel: Applied Numerical Linear Algebra (SIAM, Philadelphia, 1997). U. Trottenberg, C. Oosterlee and A. Schu¨ller: Multigrid (Academic Press, London, 2001). D. Kincaid and W. Cheney: Numerical Analysis (Brooks/Cole Publishing, California, 1996) 2nd ed. S. G. Johnson and J. D. Joannopoulos: Photonic Crystals: The Road from Theory to Practice (Kluwer Academic, Boston, 2002). L. A. Hageman and D. M. Young: Applied Iterative Methods (Academic Press, San Diego, 1981)..

(12)